∎

11institutetext: Eugene Ndiaye 22institutetext: Riken AIP
22email<EMAIL_ADDRESS>33institutetext: Ichiro Takeuchi
44institutetext: Riken AIP and Nagoya Institute of Technology
44email<EMAIL_ADDRESS>
# Root-finding Approaches for Computing Conformal Prediction Set

Eugene Ndiaye  Ichiro Takeuchi

###### Abstract

Conformal prediction constructs a confidence set for an unobserved response of
a feature vector based on previous identically distributed and exchangeable
observations of responses and features. It has a coverage guarantee at any
nominal level without additional assumptions on their distribution. Its
computation deplorably requires a refitting procedure for all replacement
candidates of the target response. In regression settings, this corresponds to
an infinite number of model fit. Apart from relatively simple estimators that
can be written as pieces of linear function of the response, efficiently
computing such sets is difficult and is still considered as an open problem.
We exploit the fact that, _often_ , conformal prediction sets are intervals
whose boundaries can be efficiently approximated by classical root-finding
algorithm. We investigate how this approach can overcome many limitations of
formerly used strategies and we discuss its complexity and drawbacks.

###### Keywords:

Prediction Set Uncertainty Quantification Distribution-Free Inference
Conformal Prediction Reliability

##  1 Introduction

(Vovk et al. ,, 2005; Shafer & Vovk,, 2008) introduced Conformal Prediction
(CP) as a general method for predicting a confidence set of a random variable
from its point prediction. Given an observed data set
$\mathcal{D}_{n}=\\{(x_{1},y_{1}),\cdots,(x_{n},y_{n})\\}$ sampled from a
distribution $\mathbb{P}$, it constructs a $100(1-\alpha)\%$ confidence set
that contains the unobserved response $y_{n+1}$ of a new instance $x_{n+1}$.
In this way, it equips traditional statistical learning algorithms with a
confidence value when predicting the response of a new test example. The
general idea is to learn a predictive model on the augmented database
$\mathcal{D}_{n+1}(z)=\mathcal{D}_{n}\cup(x_{n+1},z)$ where $z$ replaces the
unknown response $y_{n+1}$. We can therefore define a prediction loss for each
observation and rank them. A candidate $z$ will be considered as conformal or
typical if the rank of its loss is sufficiently small. The conformal
prediction set will merely collect the most typical $z$ as a confidence set
for $y_{n+1}$. As long as the sequence $\\{(x_{i},y_{i})\\}_{i=1}^{n+1}$ is
exchangeable 111Their joint probability distribution is invariant w.r.t.
permutation of the data. and the predictive model is invariant with respect to
permutation of the data, this method benefits from a strong coverage guarantee
without any assumption on the distribution. This holds for any finite sample
size $n$.

Several extensions and applications of conformal prediction have been
developed for designing uncertainty sets in active learning (Ho & Wechsler,,
2008), anomaly detection (Laxhammar & Falkman,, 2015; Bates et al. ,, 2021b),
image classification (Angelopoulos et al. ,, 2020), few shot learning (Fisch
et al. ,, 2021), time series (Xu & Xie,, 2021) or to infer the performance
guarantee for statistical learning algorithms (Holland,, 2020; Cella & Ryan,,
2020; Bates et al. ,, 2021a). We refer to the extensive reviews in
(Balasubramanian et al. ,, 2014) for other applications to artificial
intelligence.

Despite these attractive properties, the computation of conformal prediction
sets is challenging for regression problems since an infinite number of models
must be fitted with an augmented training set $\mathcal{D}_{n+1}(z)$, for all
possible $z\in\mathbb{R}$. This is not only expensive, it is simply impossible
in most cases. In general, efficiently computing conformal sets with the full
data remains an open problem. The current successful approaches for
calculating the set of conformal predictions are twofold.

  * •

Exhaustive search with a _homotopy continuation_. The fundamental idea is to
rely on the fact that the typicalness function that map each candidate with
the rank of its prediction loss is piecewise constant. As follows, if we
carefully manage to list all its transition points, we can find exactly where
it is above the prescribed confidence level. For estimators that have a
closed-form formula (e.g. Ridge or Lasso), it is possible to draw the
solutions curve w.r.t. the input candidate $z$. They are often pieces of
linear function which enable the exhaustive listing of the change points of
the rank function; see (Nouretdinov et al. ,, 2001) and (Lei,, 2019).

  * •

_Inductive confidence machine_ also called _Splitting_ (Papadopoulos et al. ,,
2002; Lei et al. ,, 2018). The observed dataset is divided into two parts. A
proper training set to fit the regression model and an independent calibration
set to calculate prediction losses and ranks. This method is the most
computationally efficient because it requires only a single model fit on a
sub-part of the data. Separating the roles of the data to build the model and
to evaluate its performance avoids refitting without loss of coverage
guarantees. The use of splitting techniques in statistics can be dated at
least to (Cox,, 1975)

These strategies have some noticeable limitations. The homotopy methods rely
on strong assumptions on the model fit and are numerically unstable due to
multiple matrix inversions that are potentially poorly conditioned. They can
suffer from exponential complexity in the worst cases, and must frequently be
abandoned because of extremely small step sizes (Gärtner et al. ,, 2012;
Mairal & Yu,, 2012). The data splitting approach does not use all the data in
the training phase. It generally results in a less accurate confidence region
i.e. of wider size. As an alternative, a common heuristic unduly restricts the
function evaluations to an arbitrary discrete grid of trial values $z$ and
select the most typical one among them. These strategies might lose the
coverage guarantee and are still computationally inefficient. As a viable
alternative, one relaxes the exact computation of the regression model at
every step and then approximately follows the homotopy continuation path by
tightly controlling the optimization error as in (Giesen et al. ,, 2010;
Ndiaye et al. ,, 2019). (Ndiaye & Takeuchi,, 2019) has shown this is a safer
discretization strategy and that it can cope with more general nonlinear
regressions. Still, it is so far limited to convex problems with strong
regularity assumptions on the model fit and fails to be applicable to most
machine learning prediction methods.

### Summary of the contributions

We build on the striking remark that for common practical situations, the
conformal prediction set is a bounded interval of the real line. Its
boundaries are the roots of coverage level $\alpha$ minus the typicalness
function and these can be efficiently computed by a root-finding algorithm
such as bisection search, with high precision and without suffering from the
limitations mentioned above. Despite its simplicity, it overcomes the
limitations of the aforementioned strategies and significantly improves and
extends the applicability of full conformal prediction to problems where it
was considered intractable so far.

We highlight some advantages of our approach.

  * •

_Efficiency_ : we demonstrate that computing a full conformal prediction set
is tractable under mild assumptions. Relying on a bisection search,
approximating the boundaries of the full exact conformal set at a prescribed
accuracy $\epsilon>0$, requires about $O(\log_{2}({1}/{\epsilon}))$ number of
model fit. The latter, trained on the whole data, allows to obtain a more
accurate confidence set than splitting methods. Accordingly, we maintain both
statistical and computational efficiency.

  * •

_Flexibility_ : our strategy offers considerable freedom on the choice of the
regression estimator. For example, it can be defined as an output of a
gradient descent process that maximizes a likelihood. It can be terminated
when the norm of the gradient is smaller than a tolerance $\epsilon_{0}$ or
after $100$ iterations of the algorithm. Consequently, the estimator can be
parameterized by the number of iterations or the optimization error resulting
from an iterative process as long as the symmetry of the data is preserved.
The proposed root-finding approaches are easily applicable to more
sophisticated recent machine learning techniques, such as deep neural networks
or models involving a non-convex regularization.

  * •

_Simplicity_ : the proposed methods are straightforward to implement. One
substantially benefits from freely available scientific computing software
packages like scikit-learn (Pedregosa et al. ,, 2011) or scipy (Virtanen et
al. ,, 2020) to adjust models and find the endpoints of the conformal set.

We also introduce an interpolation point of view of grid based approaches that
properly justifies how the coverage guarantee can be maintained along with
reduced computational time. In the case where a piecewise linear (or constant)
interpolation scheme and a simple conformity score (e.g. the absolute value)
is used, the assumption that the conformal set is an interval is not required.
The computations can be easily carried following a homotopy strategy. To
further reduce the number of model evaluations, we additionally provide a
differentiable approximation of the rank function which effectively improve
the computational efficiency of the root-finding solvers. We carefully analyze
its coverage guarantee and point out the trade-off between calibration and
number of model evaluations when such smoothing techniques are used. Such
smoothing is mainly beneficial when a high precision is required.

#### Code.

The source of our implementation is available at

https://github.com/EugeneNdiaye/rootCP

#### Notation.

For a non zero integer $n$, we denote $[n]$ the set $\\{1,\cdots,n\\}$. We
denote by $Q_{1-\alpha}$, the $(1-\alpha)$-quantile of a real valued sequence
$(U_{i})_{i\in[n+1]}$, defined as the variable
$Q_{1-\alpha}=U_{(\lceil(n+1)(1-\alpha)\rceil)}$, where $U_{(i)}$ are the
$i$-th order statistics. The interval $[a-\tau,a+\tau]$ will be denoted
$[a\pm\tau]$. For an index $j$ in $[n+1]$, the rank of $U_{j}$ among
$U_{1},\cdots,U_{n+1}$ is defined as
$\mathrm{Rank}(U_{j})=\sum_{i=1}^{n+1}\mathbb{1}_{U_{i}\leq U_{j}}$.

##  2 Conformal Prediction

We recall the arguments presented in (Vovk et al. ,, 2005; Shafer & Vovk,,
2008; Lei et al. ,, 2018) while underlining in details the intuitions and
principles that sustain the construction and validity of conformal prediction.
Let us consider an input random variable $X$ and output $Y$. The goal is to
construct a confidence set for the variable $Y$ i.e. find a set
$\mathcal{C}(X)$ such that

$\mathbb{P}(Y\in\mathcal{C}(X))\geq
1-\alpha,\quad\forall\alpha\in(0,1)\enspace.$ (1)

Given a prediction function $\mu(\cdot)$ that maps the input to the output
space and a loss measure $S$, one can assess the prediction error as
$E=S(Y,\mu(X))$. It is a random variable with cumulative distribution function
$F$ and quantile $Q$ defined as $F(z)=\mathbb{P}(E\leq z)$ and
$Q(\alpha)=\inf\\{z\in\mathbb{R}:\,F(z)\geq\alpha\\}$. The main tool for
building a set $\mathcal{C}(X)$ that satisfies the probabilistic guarantee in
Equation 1, is the following classical result222By definition, we have
$Q(\alpha)$ is the smallest real value $z$ such that $F(z)\geq\alpha$. Thus
$\alpha\leq F(Q(\alpha))=\mathbb{P}(E\leq
Q(\alpha))=\mathbb{P}(F(E)\leq\alpha)$.:

$\forall\alpha\in(0,1),\qquad\mathbb{P}(F(E)\leq\alpha)\geq\alpha\enspace.$
(2)

It implies $F(E)=F(S(Y,\mu(X))\leq 1-\alpha$ with probability larger than
$1-\alpha$. One then defines a confidence set for $Y$ as the collection of
candidate $z$ that satisfy the same inequality i.e.
$\mathcal{C}(X)=\\{z:F(S(z,\mu(X))\leq 1-\alpha\\}$. It turns out that the
same principle can be applied to compute a confidence set for sequential
observations. To do so, the coverage bound in Equation 2 can be extended to
empirical cumulative distribution and empirical quantile functions defined as:

$\displaystyle F_{n+1}(z)$
$\displaystyle=\frac{1}{n+1}\sum_{i=1}^{n+1}\mathbb{1}_{E_{i}\leq z}\enspace,$
$\displaystyle Q_{n+1}(\alpha)$
$\displaystyle=\inf\\{z\in\mathbb{R}:F_{n+1}(z)\geq\alpha\\}\enspace.$

###### Lemma 1

For a sequence of exchangeable random variables $E_{1},\cdots,E_{n+1}$, it
holds $\mathbb{P}(F_{n+1}(E_{n+1})\leq\alpha)\geq\alpha$, for any
$\alpha\in(0,1)$.

###### Proof

We follow the proof in (Romano et al. ,, 2019). By definition of the empirical
quantile, we have

$\alpha\leq
F_{n+1}(Q_{n+1}(\alpha))=\frac{1}{n+1}\sum_{i=1}^{n+1}\mathbb{1}_{E_{i}\leq
Q_{n+1}(\alpha)}\enspace.$

Taking the expectation on both side, we have

$\displaystyle\alpha$
$\displaystyle\leq\frac{1}{n+1}\sum_{i=1}^{n+1}\mathbb{P}(E_{i}\leq
Q_{n+1}(\alpha))$
$\displaystyle=\frac{1}{n+1}\sum_{i=1}^{n+1}\mathbb{P}(F_{n+1}(E_{i})\leq\alpha)\enspace.$

Moreover, we have for any $i$ in $[n]$,
$\mathbb{P}(F_{n+1}(E_{i})\leq\alpha)=\mathbb{P}(F_{n+1}(E_{n+1})\leq\alpha)$
by excheangeability. Hence the result.

Using Lemma 1, we have $F_{n+1}(E_{n+1})\leq 1-\alpha$ with probability larger
than $1-\alpha$. Then, given the $n$ previous observations, one can define a
confidence set for an unobserved variable $E_{n+1}$ as the random set

$\\{z:F_{n+1}(z)\leq 1-\alpha\\}\enspace.$

In supervised statistical learning problems, where we observe both the
responses and the features, one can apply this principle while taking benefits
of an underlying model trained on the observed data. For the augmented dataset
$\mathcal{D}_{n+1}(z)=\mathcal{D}_{n}\cup\\{(x_{n+1},z)\\}$ for
$z\in\mathbb{R}$, an example of predictive model is given as
$\hat{\mu}_{z}(x)=\Phi(x,\hat{\beta}(z))$, where $\Phi$ is a regression model
e.g. a kernel machine or a Deep Neural Network with parameter $\hat{\beta}(z)$
adjusted on the data. For example, by using empirical risk minimization
principle, one defines

$\hat{\beta}(z)\in\mathop{\mathrm{arg\,min}}_{\beta\in\mathbb{R}^{p}}L(\beta|\mathcal{D}_{n+1}(z))+\lambda\Omega(\beta)\enspace,$
(3)

where $\lambda>0$ and

$L(\beta|\mathcal{D}_{n+1}(z))=\sum_{i=1}^{n}\ell(y_{i},\Phi(x_{i},\beta))+\ell(z,\Phi(x_{n+1},\beta))$

is the data fitting term and the regularization function $\Omega$ enforces a
structured solutions e.g. sparsity.

#### Examples.

A popular example of a instance-wise loss function found in the literature is
the power norm, where $\ell(a,b)=|a-b|^{q}$. When $q=2$, this corresponds to
classical linear regression. Cases where $q\in(0,2)$ are common in robust
statistics. In particular, $q=1$ is known as least absolute deviation. The
logcosh loss $\ell(a,b)=\gamma\log(\cosh(a-b)/\gamma)$ is a differentiable
alternative to the $\ell_{\infty}$-norm. One can also have the Linex loss
function (Gruber,, 2010; Chang & Hung,, 2007) which provides an asymmetric
loss $\ell(a,b)=\exp(\gamma(a-b))-\gamma(a-b)-1$, for $\gamma\neq 0$. The
regularization functions $\Omega$ e.g. Ridge (Hoerl & Kennard,, 1970) or
sparsity inducing norms (Bach et al. ,, 2012; Obozinski & Bach,, 2016) can be
considered as well as non convex penalties (Xie & Huang,, 2009).

Given the fitted model $\mu_{z}(\cdot)$ and a loss measure $S$, let us define
the sequence of instance-wise prediction errors as:

$\displaystyle\forall i\in[n],\qquad E_{i}(z)$
$\displaystyle=S(y_{i},\,\mu_{z}(x_{i}))\enspace,$ $\displaystyle E_{n+1}(z)$
$\displaystyle=S(z,\,\mu_{z}(x_{n+1}))\enspace.$

The sequence $\\{E_{1}(y_{n+1}),\cdots,E_{n}(y_{n+1}),E_{n+1}(y_{n+1})\\}$ is
exchangeable as long as the data $\\{(x_{i},y_{i})\\}_{i=1}^{n+1}$ is
exchangeable and the model fit $\mu_{z}(\cdot)$ is invariant w.r.t.
permutation of the data. In this way, we can apply Lemma 1 and then define the
conformal set as

$\Gamma^{(\alpha)}(x_{n+1})=\\{z:F_{n+1}(E_{n+1}(z))\leq 1-\alpha\\}\enspace,$
(4)

where

$F_{n+1}(E_{n+1}(z))=\frac{1}{n+1}\sum_{i=1}^{n+1}\mathbb{1}_{E_{i}(z)\leq
E_{n+1}(z)}\enspace.$ (5)

Lemma 1 implies that the set $\Gamma^{(\alpha)}(x_{n+1})$ is a valid
confidence set for $y_{n+1}$ in the sense of Equation 1 i.e.
$\mathbb{P}(y_{n+1}\in\Gamma^{(\alpha)}(x_{n+1}))\geq 1-\alpha$ for any
$\alpha$ in $(0,1)$. Somehow, the refitting procedure with the extended
dataset $\mathcal{D}_{n+1}(z)$, puts all the variables on equal feet and
preserves the exchangeability of the sequence of prediction errors. One can
rewrite the CP set as

$\Gamma^{(\alpha)}(x_{n+1})=\\{z:\pi(z)\geq\alpha\\}\enspace,$

where $z\mapsto\pi(z)$ is the typicalness function that measures how conformal
a candidate is. It is defined as

$\pi(z)=1-\frac{1}{n+1}\mathrm{Rank}(E_{n+1}(z))\enspace.$ (6)

Lemma 1 reads $\mathbb{P}(\pi(y_{n+1})\leq\alpha)\leq\alpha$ i.e. the random
variable $\pi(y_{n+1})$ takes small values with small probability. Thus, it is
unlikely that $y_{n+1}$ will take the value $z$ when $\pi(z)$ is small. More
precisely, $\pi(y_{n+1})$ is (sub) uniformly distributed as usual for
classical statistics for hypothesis testing. For example $p$-value function
satisfies such a property under the null hypothesis; see (Lehmann & Romano,,
2006, Lemma 3.3.1). One can then interpret the typicalness $\pi(\cdot)$ as a
$p$-value function for testing the null hypothesis $H_{0}:y_{n+1}=z$ against
the alternative $H_{1}:y_{n+1}\neq z$, for $z$ in $\mathbb{R}$. The conformal
prediction set merely corresponds to the collection of candidate $z$ for which
the null hypothesis $H_{0}$ is not rejected.

##  3 Computing Conformal Prediction Set

For regression problem where $y_{n+1}$ lies in a subset of $\mathbb{R}$, one
need to evaluate $\pi(z)$ in Equation 6 and so refitting the model
$\mu_{z}(\cdot)$ for infinitely many candidate $z$. This merely renders the
overall computation challenging and leaves the problem open in general.
Nevertheless, some peculiar regularity structure of the typicalness function
$\pi(\cdot)$ can be exploited. For example, by utilizing the fact that it is
piecewise constant, it is sufficient to enumerate the transition points (when
they are finite) to compute the conformal set. This is possible for a limited
number of cases e.g. Ridge or Lasso where the map $z\mapsto\mu_{z}(\cdot)$ can
be explicitly described. Unfortunately, only a very small class of statistical
learning problems has such nice regularity structure. For other estimators,
the computation of CP set when $y_{n+1}$ can take countless number of values,
is unclear.

#### Splitting.

To overcome this issue, the split conformal prediction set introduced in
(Papadopoulos et al. ,, 2002), separates the model fitting and the score
ranking step. Let us define

  * •

the training set
$\mathcal{D}_{\rm{tr}}=\\{(x_{1},y_{1}),\cdots,(x_{m},y_{m})\\}$ with $m<n$ ,

  * •

the calibration set
$\mathcal{D}_{\rm{cal}}=\\{(x_{m+1},y_{m+1}),\cdots,(x_{n},y_{n})\\}$.

Then the model is fitted on the training set $\mathcal{D}_{\rm{tr}}$ to get
$\mu_{\rm{tr}}(\cdot)$ and define the score function on the calibration set
$\mathcal{D}_{\rm{cal}}$:

$\displaystyle\forall i\in[m+1,n],\,E_{i}^{\rm{cal}}$
$\displaystyle=S(y_{i},\,\mu_{\rm{tr}}(x_{i})),\text{ and }$ $\displaystyle
E_{n+1}^{\rm{cal}}(z)$ $\displaystyle=S(z,\,\mu_{\rm{tr}}(x_{n+1}))\enspace.$

Thus, we obtain the split typicalness function as

$\displaystyle\pi_{\rm{split}}(z)$
$\displaystyle=1-F_{\rm{split}}(E_{n+1}^{\rm{cal}}(z)),\text{ where}$
$\displaystyle F_{\rm{split}}(E_{n+1}^{\rm{cal}}(z))$
$\displaystyle=\frac{1}{n-m+1}\sum_{i=m+1}^{n+1}\mathbb{1}_{E_{i}^{\rm{cal}}\leq
E_{n+1}^{\rm{cal}}(z)}\enspace.$

The latter is proportional to the rank of the $(n+1)th$ score on the
calibration set. Finally, we define

$\displaystyle\Gamma_{\rm{split}}^{(\alpha)}(x_{n+1})$
$\displaystyle=\\{z:\pi_{\rm{split}}(z)\geq\alpha\\}$
$\displaystyle=\\{z:E_{n+1}^{\rm{cal}}(z)\leq
Q_{1-\alpha}^{\rm{cal}}\\}\enspace,$

where $Q_{1-\alpha}^{\rm{cal}}$ is the $(1-\alpha)$ quantile of the
calibration scores $\\{E_{m+1}^{\rm{cal}},\cdots,E_{n+1}^{\rm{cal}}\\}$. When
the score function is the absolute value $S(a,b)=|a-b|$, the split CP set is
the interval
$\Gamma_{\rm{split}}^{(\alpha)}(x_{n+1})=[\mu_{\rm{tr}}(x_{n+1})\pm
Q_{1-\alpha}^{\rm{cal}}]$. While this approach avoids the computational
bottleneck, the statistical efficiency of the model can be reduced due to a
significantly smaller sample available during the training and calibration
phase. Moreover, the length of the split conformal set tends to have a higher
variance.

###  3.1 Approximation to a Prescribed Accuracy

In this paper, we promptly take advantages of the remarkable fact that the
conformal regions are _often_ intervals. We subsequently take an alternative
direction which carefully avoids tracking the integral path of all model fit.
When the $(1-\alpha)$-level set of the repartition function $z\mapsto
F_{n+1}(E_{n+1}(z))$ is convex e.g. Figure 1, we propose to employ a numerical
root-finding solver to approximate the endpoints of the interval.

### Outline of the Algorithm: rootCP

Assuming that the conformal set is a non empty interval of finite length, we
denote

$\Gamma^{(\alpha)}(x_{n+1})=[\ell_{\alpha}(x_{n+1}),u_{\alpha}(x_{n+1})]\enspace.$

Given a tolerance $\epsilon>0$, we proceed as follows:

  1. 1.

find $z_{\min}<z_{0}<z_{\max}$ such that

$\pi(z_{\min})<\alpha<\pi(z_{0})\text{ and }\alpha>\pi(z_{\max})\enspace.$ (7)

  2. 2.

Perform a bisection search in $[z_{\min},z_{0}]$. It will output a point
$\hat{\ell}$ such that $\ell_{\alpha}(x_{n+1})$ belongs to
$[\hat{\ell}\pm\epsilon]$ after at most
$\log_{2}(\frac{z_{0}-z_{\min}}{\epsilon})$ iterations.

  3. 3.

Perform a bisection search in $[z_{0},z_{\max}]$. It will output a point
$\hat{u}$ such that $u_{\alpha}(x_{n+1})$ belongs to $[\hat{u}\pm\epsilon]$
after at most $\log_{2}(\frac{z_{\max}-z_{0}}{\epsilon})$ iterations.

(a) Failed initialization at first trial (b) Successful initialization at
first trial

Figure 1: Illustration of the initialization steps when both the initial
prediction $\mu_{\mathcal{D}_{n}}(x_{n+1})$ and midpoint of split conformal
interval $\mu_{\rm{tr}}(x_{n+1})$ fails to be in the full and exact conformal
prediction set. The synthetic data are generated with sklearn as $X,y=$
make_regression$(n=300,p=50)$. We choose an optimization accuracy of
$\epsilon_{0}=\left\lVert(y_{1},\cdots,y_{n})\right\rVert_{2}^{2}/10^{4}$ for
approximating the ridge estimator. The number of trial points is $d=10$ and
$\displaystyle z^{(\epsilon)}\in\mathop{\mathrm{arg\,max}}_{z\in
C_{d}}\pi^{(\epsilon)}(z)$ is the most conformal trial candidate at precision
$\epsilon\geq 0$.

(a) Boundaries of conformal sets (b) Computational times

Figure 2: Benchmark on ridge regression. Conformal prediction set computed
with various regularization parameter on synthetic dataset generated from
sklearn as $X,y=$ make_regression$(n=1000,p=100)$ with $90$ informative
features. For the splitting method, we average the results of $100$
independent run. For the proposed root finding method, we approximate the
boundaries of the exact set at precision $10^{-12}$.

### Initialization

For the initial lower and upper bounds, we suggest

$z_{\min}=\min_{i\in[n]}y_{i}\text{ and
}z_{\max}=\max_{i\in[n]}y_{i}\enspace.$

For most of the situations encountered in our numerical experiments, we
consistently get $\pi(z_{\min})$ and $\pi(z_{\max})$ both smaller than the
threshold level $\alpha$. Otherwise, we can always take values even farther
apart without affecting the complexity thanks to the logarithmic dependence in
the length of the initialization brackets. This is especially necessary when
the total number of samples $n$ is small333Indeed, we have
$\mathbb{P}(y_{n+1}\in[z_{\min},z_{\max}])\geq 1-\frac{1}{n+1}$. Hence when
$n$ is sufficiently large i.e. $n\geq 1+1/\alpha$, then $[z_{\min},z_{\max}]$
is a $(1-\alpha)$ confidence set.. The most crucial part is to choose $z_{0}$
so that $\pi(z_{0})>\alpha$. It is equivalent to get a point in the interior
of the conformal set itself. In the ideal case where the length of the
conformal set is extremely small, finding an initialization point might be
notoriously hard. Indeed, it corresponds to a _rare event_ equivalent to
sampling a point in a low probability region. We adopted a simple strategy
which consists in estimating $y_{n+1}$ with the observed data
$\mathcal{D}_{n}$. We subsequently denote it

$z_{0}=\mu_{\mathcal{D}_{n}}(x_{n+1})\enspace.$

In our sequence of repetitive numerical experiments, this choice rarely fails.
Naturally, its success depends on the prediction capabilities of the model
fit. In the rare cases where it fails, we propose to test the initialization
condition on some query points selected on an initial estimation
$[z_{\alpha}^{-},z_{\alpha}^{+}]$ of the CP set. This localization step aims
to exploit additional problem structure and can be interpreted as an iterative
importance sampling to maintain a reasonably low computational cost.

  1. 1.

_Localization._ Given an easy to compute estimate set
$[z_{\alpha}^{-},z_{\alpha}^{+}]$ that is potentially larger444One can
slightly enlarge the set by taking
$[z_{\alpha}^{-}-0.5(z_{\alpha}^{+}-z_{\alpha}^{-}),z_{\alpha}^{+}+0.5(z_{\alpha}^{+}-z_{\alpha}^{-})].$
than the targeted conformal set, we select its mid point

$z_{0}=\frac{z_{\alpha}^{+}+z_{\alpha}^{-}}{2}\enspace.$

If $z_{0}$ satisfies $\pi(z_{0})>\alpha$, then we have a valid initialization
by paying only a single model fit. Otherwise, we run the next step on the
bracket $[z_{\alpha}^{-},z_{\alpha}^{+}]$.

For instance, one can use
$[z_{\alpha}^{-},z_{\alpha}^{+}]=\Gamma_{\rm{split}}^{(\alpha)}(x_{n+1})$ or
$[z_{\alpha}^{-},z_{\alpha}^{+}]=\\{z:\pi_{\mathcal{D}_{n}}(z)>\alpha\\}$
where $\pi_{\mathcal{D}_{n}}(\cdot)$ is an unsafe estimation of $\pi(\cdot)$
with $\mu_{z}(x)$ replaced by $\mu_{\mathcal{D}_{n}}(x)$ for any candidate $z$
and any input feature $x$.

  2. 2.

_Sampling._ For a small number $d$ e.g. $d=5$, and given a bracket search
$[z^{-},z^{+}]$, select candidates $C_{d}=\\{z_{1},\cdots,z_{d}\\}$ uniformly.
If there is $\pi(z_{0})>\alpha$ for a $z_{0}$ in $C_{d}$, we have a valid
initialization. Otherwise, we use these query points to interpolate the model
fit as in eq. 8. Thus, by selecting additional points that have a higher
typicalness according to the interpolated model, one can refine the sampling
set $C_{d}$ and repeat the process.

Additionally, we explain below how this model fit interpolation can be used to
obtain an alternative CP set. Note that its midpoint can also be used as a
candidate for initialization. For computational efficiency, one can rely
predominantly on the fact that for the usual prediction problems in machine
learning, it is unneeded to optimize below the inevitable statistical error
(Bousquet & Bottou,, 2008). This means that a high optimization accuracy in
the model fit might be unnecessary to achieve better generalization
performances. Therefore, with a coarse optimization tolerance, we can preview
the final shape of the conformity function. Whence, one can replace $\pi(z)$
by $\pi^{(\epsilon_{0})}(z)$ which is computed with a rough optimization error
$\epsilon_{0}$ in order to guess the shape of the function $\pi(\cdot)$.
Similarly, if $\pi(z_{0})$ fail to be a valid initialization, then decrease
$\epsilon_{0}$ and repeat the process. In all our experiments, it works fine
after very few number of iterations. Nonetheless, we do not have strategies to
avoid worst-case situations or any mathematical guaranty in the total number
of iterations needed to find a valid initialization. We illustrate this
strategy in Figure 1 in both situations of failure and success.

#### Further complexity reduction.

In cases where the regression map $z\mapsto\mu_{z}(x)$ for any feature $x$,
can be traced with homotopy as in Ridge (Nouretdinov et al. ,, 2001) and Lasso
(Lei,, 2019), it takes $O(n^{2})$ to compute the exact conformal set. This can
be reduced to $O(n\log n)$ by sorting the roots of the instance-wise scores
$E_{i}(z)-E_{n+1}(z)$ for $i$ in $[n]$ and cleverly flattening the double loop
when evaluating the ranks of the score functions (Vovk et al. ,, 2005, Chapter
2.3). By relaxing the exactness, none of these two steps is needed in our
approach. We obtain an asymptotic improvement to $O(n\log_{2}(1/\epsilon))$
and an easier to implement algorithm.

When the model fit is parameterized by the solution of optimization problem in
Equation 3, the regularity of the loss function and penalty terms play a major
role in the computational tractability of the full conformal prediction set.
Leveraging smoothness and convexity assumptions on the loss or penalty
functions, it has been shown in (Ndiaye & Takeuchi,, 2019) that approximate
solutions can be used without refitting the model for close candidates $z$.
The resulting conformal set is
$\overline{\Gamma}^{(\alpha,\epsilon)}=\\{z\in\mathbb{R}:\,\overline{\pi}(z,\epsilon)>\alpha\\}$
where the corresponding typicalness function incorporates the optimization
error i.e.

$\overline{\pi}(z,\epsilon)=\frac{1}{n+1}\sum_{i=1}^{n+1}\mathbb{1}_{E_{i}(z)\geq
E_{n+1}(z)-2\sqrt{2\nu\epsilon}}\enspace.$

One can further show that the typicalness function based on exact solution
$\hat{\pi}(\cdot)$ is uniformly upper bounded by
$\overline{\pi}(\cdot,\epsilon)$ and then
$\hat{\Gamma}^{(\alpha)}\subset\overline{\Gamma}^{(\alpha,\epsilon)}$. This
can be equally used to reduce the number of model fit and also to wrap the CP
set based on exact solution by applying rootCP directly to
$\overline{\pi}(\cdot,\epsilon)$ instead of computing a whole approximation
path.

Compared with the homotopy approach, rootCP will always make a smaller number
of model fit. By merely storing each model evaluation, it benefits from warm-
start boosting by employing the solutions of the previous function call. Alas,
the homotopy approaches require either an exact solution or an approximate one
with a strict control of the optimization error. This control is not always
available if one does not provide a computable upper bound of it, for example
by precisely evaluating the duality gap. Such bounds are hardly available in
non-convex setting which greatly reduce its applicability in modern machine
learning techniques. Meanwhile, the complexity of the approximate homotopy
algorithm is $O(\frac{z_{\max}-z_{\min}}{\sqrt{\epsilon_{0}}})$ where
$\epsilon_{0}$ is the optimization error of the model fit. Additionally to the
linear dependence in the initial interval length, it cannot be launched for
small $\epsilon_{0}$ whereas the number of model fit in the root-finding
approach is not degraded. In a nutshell, the proposed method avoids the
computation of the whole path. Hence, it enjoys an exponential improvement
over the homotopy approach w.r.t. to the initial interval length
$z_{\max}-z_{\min}$ and an overall complexity that is independent of
$\epsilon_{0}$. It can then be used with highly optimized model fit where the
homotopy method cannot even be launched.

#### Drawbacks.

Full conformal prediction set is _not always_ an interval. When it is a union
of few well separated intervals, our proposed method cannot be applied without
finely bracketing these intervals. One can include a human in the loop. The
discrete function $\pi^{(\epsilon_{0})}$ offers a cheap pre-visualization of
the landscape of the conformity function that allows to detect these
situations and infer a proper bracketing. At this point, efficiently
enumerating all the roots remain a challenging task that we leave as an open
problem. In the following proposition, we provide a sufficient condition so
that the conformal set is an interval.

###### Proposition 1

If for any $i$ in $[n]$, the difference of instance-wise error function
$z\mapsto E_{i}(z)-E_{n+1}(z)$ is quasi-concave and has two zeros $a_{i}\leq
b_{i}$ such that

$\max_{i\in[n]}a_{i}\leq\min_{i\in[n]}b_{i}\enspace,$

then the typicalness function $\pi(\cdot)$ is quasi-concave and the conformal
prediction set at a level $\alpha\in(0,1)$ is either empty or an interval.

Proposition 1 generalizes (Lei,, 2019, Theorem 3.3) which provides a
sufficient condition so that the Lasso conformal set is an interval.
Unfortunately, such sufficient conditions is not testable for most problems.
Indeed, it requires knowing all the zero crossing points of the function
$z\mapsto E_{i}(z)-E_{n+1}(z)$ for all index $i\in[n]$ which is as hard as
computing the whole function $z\mapsto\pi(z)$. We remind that even when the
typicalness function $\pi(\cdot)$ is not quasi-concave, our algorithm is still
valid as long as the conformal set is an interval; which is a much weaker
assumption than quasi-concavity.

Figure 3: Illustration of the smoothed conformal set with data generated from sklearn as $X,y=$ make_regression$(n=300,p=50)$. The smoothed typicalness function $\pi(\cdot,\gamma)$ is evaluated with several hyperparameter $\gamma$. The underlying estimator is the ridge regressor with parameter $\lambda=p/\lVert\beta_{\mathrm{LS}}\rVert^{2}$ where $\beta_{\mathrm{LS}}$ is the Least-squares estimator on the observed dataset $\mathcal{D}_{n}$.  |  |  | Approximate homotopy / rootCP
---|---|---|---
| oracleCP | splitCP | 1e-2 | 1e-4 | 1e-6
Coverage | 0.94 | 0.90 | 0.95 / 0.93 | 0.94 / 0.94 | 0.94 / 0.94
Length | 0.693 | 1.054 | 0.747 / 0.699 | 0.711 / 0.708 | 0.707 / 0.707
Time (s) | 3.088 | 1.409 | 4.289 / 2.478 | 29.987 / 9.173 | 316.67 / 20.92
Table 1: Computing a conformal set for a Lasso regression problem on a climate
data set with $n=814$ observations and $p=73570$ features. We display the
coverage, length and execution times for different methods averaged over $100$
randomly held-out validation data sets.

###  3.2 Interpolated Conformal Prediction

The full conformal prediction set is computationally expensive since it
requires knowing exactly the map $z\mapsto\mu_{z}(\cdot)$. The splitting
approach does not use all the data in the learning phase but is
computationally efficient since it requires a single model fit. Alternatively,
it was proposed in (Lei et al. ,, 2018) to use an arbitrary discretization and
its theoretical analysis in (Chen et al. ,, 2018) unfortunately failed to
preserve the coverage guarantee. In this section, we argue that grid based
strategy with an interpolation point of view, stands as an ”in-between”
strategy that exploits full data with a restricted computational time while
preserving the coverage guarantee. We propose to compute a conformal
prediction set based on an interpolation of the model fit map given a finite
number of query points. The main insight is that the underlying model fit
plays a minor role in the coverage guarantee; the only requirement is to be
symmetric with respect to permutation of the data. As such, the model path
$z\mapsto\hat{\mu}_{z}(\cdot)$ can be replaced by an interpolated map
$z\mapsto\tilde{\mu}_{z}(\cdot)$ based on query points $z_{1},\cdots,z_{d}$.
It reads to a valid prediction set as long as the interpolation preserves the
symmetry 555Otherwise one can always perform a symmetrization e.g. using a
model parameter
$\tilde{\beta}(z)=\frac{1}{(n+1)!}\sum_{\sigma\in\Sigma_{n+1}}\beta(w_{\sigma(1)},\cdots,w_{\sigma(n)},w_{\sigma(n+1)})$,
where $w_{i}=(x_{i},y_{i})$ if $i$ in $[n]$, $w_{n+1}=(x_{n+1},z)$ and
$\Sigma_{n+1}$ is the group of permutation of $[n+1]$.. For instance, one can
rely on a piecewise linear interpolation

$\tilde{\mu}_{z}=\begin{cases}\frac{z_{1}-z}{z_{1}-z_{\min}}\hat{\mu}_{z_{\min}}+\frac{z_{\min}-z}{z_{1}-z_{\min}}\hat{\mu}_{z_{1}}&\text{
if }z\leq z_{\min}\enspace,\\\
\frac{z-z_{t+1}}{z_{t}-z_{t+1}}\hat{\mu}_{z_{t}}+\frac{z-z_{t}}{z_{t+1}-z_{t}}\hat{\mu}_{z_{t+1}}&\text{
if }z\in[z_{t},z_{t+1}]\enspace,\\\
\frac{z-z_{d}}{z_{\max}-z_{d}}\hat{\mu}_{z_{\max}}+\frac{z_{\max}-z}{z_{\max}-z_{d}}\hat{\mu}_{z_{d}}&\text{
if }z\geq z_{\max}\enspace,\end{cases}$ (8)

where $\hat{\mu}_{z}(x)$ is a prediction map trained on the augmented dataset
$\mathcal{D}_{n+1}(z)$ as in Equation 3. As before, one defines the instance-
wise score functions

$\displaystyle\forall i\in[n],\qquad\tilde{E}_{i}(z)$
$\displaystyle=S(y_{i},\,\tilde{\mu}_{z}(x_{i}))\enspace,$
$\displaystyle\tilde{E}_{n+1}(z)$
$\displaystyle=S(z,\,\tilde{\mu}_{z}(x_{n+1}))\enspace.$

The conformal set based on interpolated model fit is then defined as

$\displaystyle\tilde{\Gamma}^{(\alpha)}(x_{n+1})$
$\displaystyle=\\{z:\tilde{\pi}(z)\geq\alpha\\},\text{ where }$
$\displaystyle\tilde{\pi}(z)$
$\displaystyle=1-\frac{1}{n+1}\sum_{i=1}^{n+1}\mathbf{1}_{\tilde{E}_{i}(z)\leq\tilde{E}_{n+1}(z)}\enspace.$

Since the map $z\mapsto\tilde{\mu}_{z}(\cdot)$ is (or can be made) symmetric,
it is immediate to see that $\tilde{\Gamma}^{(\alpha)}(x_{n+1})$ is a valid
conformal set following the same proof technique.

We remind that the conformal set can be highly concentrated around its mid
point and the typicalness of most candidates is close to zero. Whence, we
suggest to restrict the query points around an estimate of the conformal set
provided by a localization step.

######  Remark 1 (Interpolation of the typicalness map)

Given the query points and their corresponding typicalness
$(z_{1},\pi(z_{1})),\cdots,(z_{d},\pi(z_{d}))$, one can also directly learn a
function that approximate the typicalness $z\mapsto\pi(z)$. However, in this
case, we could not establish the theoretical coverage guarantee of this
method. Moreover, when the conformal set is highly localized, most of the
$\pi(z_{i})$ might be close to zero leading to a flat and poorly interpolated
typicalness map.

Previous discretization approaches did not preserve the coverage guarantee or
did it at expensive cost by approximating a model fit path on a wide range,
with a high precision at every step, and was restricted to convex problems.
The interpolation point of view that we provided allows to compute a valid
conformal prediction set with arbitrary discretization without loss in the
coverage guarantee and without restriction to convex problems. Also note that
depending on the interpolation used, there is no need to assume that the
conforming prediction set is an interval. Indeed, in the case below of
piecewise linear interpolation, one can easily enumerate all the change points
of the conformal function as in homotopy methods. However, in general we also
recommend the use of rootCP. Finally, note that in the case of the ridge
estimator (which is linear in $z$), the exact conformal prediction coincides
with that of the interpolation. In all our experiments, we have chosen $d=8$.

###  3.3 Smoothed Conformal Prediction

Conformal prediction set relies on rank computations. The latter function is
piecewise constant and has no useful first order information in the sense that
it is either null or undefined. We propose a smooth approximation of the
typicalness function to reduce the number of query points. In addition to
exchangeability, we merely use the fact that $F_{n+1}$ is increasing and the
linearity of the sum to obtain the coverage guarantee. Likewise, one should be
able to replace $\mathbb{1}_{E_{i}-z\leq 0}$ with a continuously
differentiable and increasing function $\phi_{\gamma}(E_{i}-z)$. Hence,
replacing the function $z\mapsto F_{n+1}(E_{n+1}(z))$ by a smoother one allows
the use of more efficient gradient or quasi-Newton-based root finding methods.
We further investigate the influence of such smoothing on the coverage
guarantee. In practice, we simply choose the sigmoid function
$\phi_{\gamma}(x)=\frac{\mathrm{e}^{-\gamma x}}{1+\mathrm{e}^{-\gamma x}}$ as
in (Qin et al. ,, 2010). We have

$\displaystyle\mathrm{Rank}(u_{j})$
$\displaystyle:=\sum_{i=1}^{n+1}\mathbb{1}_{u_{i}-u_{j}\leq 0}$
$\displaystyle\approx\sum_{i=1}^{n+1}\phi_{\gamma}(u_{i}-u_{j})=:\mathrm{sRank}(u_{j},\gamma)\enspace.$

The main advantage is that the map $\mathrm{sRank}(\cdot,\gamma)$ improve the
regularity of $\mathrm{Rank}(\cdot)$ and allows faster convergence. The smooth
approximation of the typicalness function is then defined as

$\displaystyle\pi(z)$
$\displaystyle\approx\pi(z,\gamma):=1-\frac{1}{n+1}\mathrm{sRank}(E_{n+1}(z),\gamma)\enspace,$

and the smoothed conformal prediction set (illustrated in Figure 3) as

$\Gamma^{(\alpha,\gamma)}(x_{n+1})=\\{z:\pi(z,\gamma)>\alpha\\}\enspace.$

Now computing an approximation of the conformal prediction set is equivalent
to finding the smallest and largest solution of the equation
$\pi(z,\gamma)=\alpha$ which is often easier to solve than $\pi(z)=\alpha$.
Using different root-finding solvers, we illustrate the computational
advantages by displaying the reduction of the number of model fit in Figure 3.

######  Remark 2 (Gradient Based Solvers)

When, for any feature $x$, the regression map $z\mapsto\mu_{z}(x)$ is
differentiable, the solutions of equation $\pi(z,\gamma)=\alpha$ could be
approximated with more efficient gradient based root-finding algorithm.
However, the function $z\mapsto\pi(z,\gamma)$ is mostly flat except at a tiny
vicinity of the conformal set which makes the convergence difficult unless a
good initialization is found. One could also rely on a regularized version by
minimizing $(\pi(z,\gamma)-\alpha)^{2}+\tau z^{2}$ which requires a proper
tuning of the hyper parameter $\tau$. Both of these strategies turn out to be
less stable and need further investigations.

We analyze the statistical consequences of using a continuous version of the
indicator function. We recall the definition of the smoothed version of the
empirical cumulative distribution and empirical quantile:

$\displaystyle\tilde{F}_{n+1}(z)=\frac{1}{n+1}\sum_{i=1}^{n+1}\phi_{\gamma}(E_{i}-z)\enspace,$
$\displaystyle\tilde{Q}_{n+1}(\alpha)=\inf\\{z\in\mathbb{R}:\tilde{F}_{n+1}(z)\geq\alpha\\}\enspace.$

######  Proposition 2 (Coverage guarantee of the smooth relaxation)

For a sequence of exchangeable random variable $E_{1},\ldots,E_{n+1}$, it
holds for any $\tilde{\alpha}$ in $(0,1)$,

$\mathbb{P}(\tilde{F}_{n+1}(E_{n+1})\leq\tilde{\alpha})\geq\tilde{\alpha}-\Delta(\gamma)\enspace,$

where $\Delta(\gamma)=\sup_{x}(\phi_{\gamma}-\mathbb{1}_{\cdot\leq 0})(x)$.

###### Proof

By definition of $\tilde{F}_{n+1}$ and $\tilde{Q}_{n+1}$, we have

$\displaystyle\tilde{\alpha}$
$\displaystyle\leq\tilde{F}_{n+1}(\tilde{Q}_{n+1}(\tilde{\alpha}))$
$\displaystyle=\frac{1}{n+1}\sum_{i=1}^{n+1}\phi_{\gamma}(E_{i}-\tilde{Q}_{n+1}(\tilde{\alpha}))$
$\displaystyle=\frac{1}{n+1}\sum_{i=1}^{n+1}\mathbb{1}_{E_{i}\leq\tilde{Q}_{n+1}(\tilde{\alpha})}+(\phi_{\gamma}-\mathbb{1}_{\cdot\leq
0})(E_{i}-\tilde{Q}_{n+1}(\tilde{\alpha}))$
$\displaystyle\leq\frac{1}{n+1}\sum_{i=1}^{n+1}\mathbb{1}_{\tilde{F}_{n+1}(E_{i})\leq\tilde{\alpha}}+\Delta(\gamma)\enspace.$

We conclude by taking the expectation on both side along with exchangeability.

In order to display a probabilistic statement, one needs to maintain the
indicator function when defining the typicalness function. Replacing it with a
continuous version will distort the coverage guarantee as described in
Proposition 2. To obtain a well $\alpha$-calibrated confidence set, one must
take into account such approximation error by choosing
$(\tilde{\alpha},\gamma)$ such that

$\tilde{\alpha}-\Delta(\gamma)\geq\alpha\enspace.$

If $\tilde{\alpha}$ is fixed, one needs to be careful when choosing $\gamma$.
Otherwise, we obtain a vacuous upper bound and all the coverage guarantee is
lost. Meanwhile, if $\gamma$ is chosen such that $\phi_{\gamma}$ is a lower
approximation of the indicator function, then $\tilde{\alpha}$ can be taken as
$\alpha$ and there is no calibration loss. However, when $\Delta(\gamma)$ is
close to zero, $\phi_{\gamma}$ will be flat almost everywhere and we will not
get useful first order information. This brings a _trade-off_ between number
of model fitting (which influences the computational time) and efficiency i.e.
length of the interval (wider $\tilde{\alpha}$-level set).

#### Building a gap.

To finely assess how the vanilla conformal and smoothed conformal set can be
related in practice, one can simply design both a lower and upper
approximation of the indicator function i.e. $\phi_{\gamma}^{+}$ and
$\phi_{\gamma}^{-}$. In that case, it is easy to see that

$\ell_{\gamma}^{+}\leq\ell_{\alpha}(x_{n+1})\leq\ell_{\gamma}^{-}\text{ and
}u_{\gamma}^{-}\leq u_{\alpha}(x_{n+1})\leq u_{\gamma}^{+}\enspace,$

which is equivalent to

$\Gamma_{\gamma}^{(\alpha,-)}(x_{n+1})\subset\Gamma^{(\alpha)}(x_{n+1})\subset\Gamma_{\gamma}^{(\alpha,+)}(x_{n+1})\enspace.$

The overall complexity is moderately expanded (we now need to compute two
different conformal prediction sets) and not too time consuming as long as the
underlying model fit is reasonably computable.

##  4 Experiments

Figure 4: Benchmarking conformal sets for different regression models on real
datasets. We display the lengths of the confidence sets over $100$ random
permutation of the data. We denoted $\overline{cov}$ the average coverage and
$\overline{T}$ the average computational time normalized with the average time
for computing oracleCP which requires a single model fit on the whole data.

We numerically examine the performance of the root-finding methods to compute
various conformal prediction sets for regression problems on both synthetic
and real databases. The experiments were conducted with a coverage level of
$0.9$ i.e. $\alpha=0.1$. For comparisons, we run the evaluations on $100$
repetitions of examples and display the average of the following performance
statistics for different methods: 1) the empirical coverage i.e. the
percentage of times the prediction set contains the held-out target $y_{n+1}$;
2) the length of the confidence intervals; 3) the execution time.

From Lemma 1, we have $\pi(y_{n+1})\geq\alpha$ with probability larger than
$1-\alpha$. Whence one can define the OracleCP as $\pi^{-1}([\alpha,+\infty))$
where $\pi$ is obtained with a model fit optimized on the oracle data
$\mathcal{D}_{n+1}(y_{n+1})$. In the case where the conformity function is the
absolute value, we obtain the reference prediction set as in (Ndiaye &
Takeuchi,, 2019)

$\texttt{oracleCP:
}[\mu_{y_{n+1}}(x_{n+1})\;\pm\;Q_{1-\alpha}(y_{n+1})]\enspace.$

We remind that the target variable $y_{n+1}$ is not available in practice. In
the case of Ridge regression, _exact_ conformal prediction sets can be
computed by homotopy without data splitting and without additional assumptions
(Nouretdinov et al. ,, 2001). This allows us to finely assess the precisions
of the proposed approaches and illustrate the speed up benefit in Figure 2.

We also illustrate the performance of our approach compared to the approximate
homotopy method for the Lasso problem. The latter uses all the data and does
not loose in statistical efficiency. However, as it has already been noticed
in (Ndiaye & Takeuchi,, 2019), this method becomes unusable when the
optimization accuracy becomes small. We have shown that it does not affect our
method because its number of model fits does not increase when the
optimization error decreases see Table 1.

We run experiments on more complex regression model such as orthogonal
matching pursuit (OMP), feedforward neural network also called Multi-Layer
Perceptron (MLP), Random Forest and Gradient Boosting (with warm start). In
most of these settings, the estimator is obtained by approximating a solution
of a non-convex optimization problem where none of the homotopy methods are
available. We can observe in Figure 4 that the root-finding approach computes
a full conformal prediction set while maintaining a reasonable computational
time. In the worst cases observed, it cost about $30$ times a single model fit
which roughly correspond to the $15$ model fit for each root as predicted by
the complexity (e.g. $\epsilon=10^{-4}$). Moreover, this computational time is
significantly reduced by the interpolationCP while achieving a statistical
performance almost identical to the vanilla one in all our simulations.

##  5 Conclusion

Since its introduction, the computation of a confidence region with full
conformal prediction methods has been a major weakness to its adoption by a
broader audience. The algorithms available until now were based on too strong
assumptions that limited them to estimators whose map $z\mapsto\mu_{z}(\cdot)$
can be traced e.g. by using homotopy methods. We have shown that the
limitations of the previous methods can be overcome by directly estimating the
endpoints of the $\alpha$-level set of the typicalness function with a root
finding algorithm. Therefore, it is unnecessary to train the regression
estimator an infinite number of times nor to make strong additional
assumptions on the prediction model. As long as the conformal set represents
an interval containing the point prediction obtained from the observed data,
it can easily be estimated with only about ten number of model fits. The
proposed approach can be readily applied to recent generalizations of the
conformal prediction set beyond the exchangeability assumption as in
(Chernozhukov et al. ,, 2018, 2021).

Nevertheless, we insist that full conformal prediction set is not always an
interval. In this case, our approach fails without properly bracketing these
intervals. Both, testing the interval assumption or finding a proper
bracketing are still difficult. Another severe and silent disadvantage is that
the conformal set itself might be ill-defined when applied to regression
estimators that depends on solving a non-convex optimization problem or a
stochastic scheme (e.g. stochastic gradient descent). Indeed, in these cases,
given a fixed candidate z, $\pi(z)$ can takes multiples values depending on
the initialization or the random seed. Moreover, the symmetry assumptions
might be violated if the instances are not used evenly (e.g. when using
stochastic gradient descent with importance sampling). As a future work, it
would therefore be interesting to understand how these points can negatively
affect the coverage guarantee and the computational complexity.

###### Acknowledgements.

This work was partially supported by MEXT KAKENHI (20H00601, 16H06538), JST
Moonshot R&D (JPMJMS2033-05), NEDO (JPNP18002, JPNP20006) and RIKEN Center for
Advanced Intelligence Project.

## References

  * Angelopoulos et al. , (2020) Angelopoulos, A., Bates, S., Malik, J., & Jordan, M. 2020.  Uncertainty Sets for Image Classifiers using Conformal Prediction.  ICLR.
  * Bach et al. , (2012) Bach, F., Jenatton, R., Mairal, J., & Obozinski, G. 2012.  Optimization with sparsity-inducing penalties.  Foundations and Trends in Machine Learning.
  * Balasubramanian et al. , (2014) Balasubramanian, V., Ho, S-S., & Vovk, V. 2014.  Conformal prediction for reliable machine learning: theory, adaptations and applications.  Elsevier.
  * Bates et al. , (2021a) Bates, S., Angelopoulos, A., Lei, L., Malik, J., & Jordan, M. 2021a.  Distribution-Free, Risk-Controlling Prediction Sets.  arXiv preprint arXiv:2101.02703.
  * Bates et al. , (2021b) Bates, S., Candès, E., Lei, L., Romano, Y., & Sesia, M. 2021b.  Testing for Outliers with Conformal p-values.  arXiv preprint arXiv:2104.08279.
  * Bousquet & Bottou, (2008) Bousquet, O., & Bottou, L. 2008.  The tradeoffs of large scale learning.  NeurIPS.
  * Cella & Ryan, (2020) Cella, L., & Ryan, R. 2020.  Valid distribution-free inferential models for prediction.  arXiv preprint arXiv:2001.09225.
  * Chang & Hung, (2007) Chang, Y-C., & Hung, W-L. 2007.  LINEX loss functions with applications to determining the optimum process parameters.  Quality & Quantity.
  * Chen et al. , (2018) Chen, W., Chun, K-J., & Barber, R. F. 2018.  Discretized conformal prediction for efficient distribution-free inference.  Stat.
  * Chernozhukov et al. , (2018) Chernozhukov, V., Wüthrich, K., & Zhu, Y. 2018.  Exact and robust conformal inference methods for predictive machine learning with dependent data.  Conference On Learning Theory.
  * Chernozhukov et al. , (2021) Chernozhukov, V., Wüthrich, K., & Zhu, Y. 2021.  An exact and robust conformal inference method for counterfactual and synthetic controls.  Journal of the American Statistical Association.
  * Cox, (1975) Cox, D. R. 1975.  A note on data-splitting for the evaluation of significance levels.  Biometrika.
  * Fisch et al. , (2021) Fisch, A., Schuster, T., Jaakkola, T., & Barzilay, R. 2021.  Few-shot Conformal Prediction with Auxiliary Tasks.  ICML.
  * Gärtner et al. , (2012) Gärtner, B., Jaggi, M., & Maria, C. 2012.  An exponential lower bound on the complexity of regularization paths.  Journal of Computational Geometry.
  * Giesen et al. , (2010) Giesen, J., Jaggi, M., & Laue, S. 2010.  Approximating parameterized convex optimization problems.  European Symposium on Algorithms.
  * Gruber, (2010) Gruber, M. 2010.  Regression estimators: A comparative study.  JHU Press.
  * Ho & Wechsler, (2008) Ho, S-S., & Wechsler, H. 2008.  Query by transduction.  IEEE transactions on pattern analysis and machine intelligence.
  * Hoerl & Kennard, (1970) Hoerl, A. E., & Kennard, R. W. 1970.  Ridge regression: Biased estimation for nonorthogonal problems.  Technometrics.
  * Holland, (2020) Holland, M. J. 2020.  Making learning more transparent using conformalized performance prediction.  arXiv preprint arXiv:2007.04486.
  * Laxhammar & Falkman, (2015) Laxhammar, R., & Falkman, G. 2015.  Inductive conformal anomaly detection for sequential detection of anomalous sub-trajectories.  Annals of Mathematics and Artificial Intelligence.
  * Lehmann & Romano, (2006) Lehmann, E. L., & Romano, J. P. 2006.  Testing statistical hypotheses.  Springer Science & Business Media.
  * Lei, (2019) Lei, J. 2019.  Fast Exact Conformalization of Lasso using Piecewise Linear Homotopy.  Biometrika.
  * Lei et al. , (2018) Lei, J., G’Sell, M., Rinaldo, A., Tibshirani, R. J., & Wasserman, L. 2018.  Distribution-free predictive inference for regression.  Journal of the American Statistical Association.
  * Mairal & Yu, (2012) Mairal, J., & Yu, B. 2012.  Complexity Analysis of the Lasso Regularization Path.  ICML.
  * Ndiaye & Takeuchi, (2019) Ndiaye, E., & Takeuchi, I. 2019.  Computing Full Conformal Prediction Set with Approximate Homotopy.  NeurIPS.
  * Ndiaye et al. , (2019) Ndiaye, E., Le, T., Fercoq, O., Salmon, J., & Takeuchi, I. 2019.  Safe Grid Search with Optimal Complexity.  ICML.
  * Nouretdinov et al. , (2001) Nouretdinov, I., Melluish, T., & Vovk, V. 2001.  Ridge regression confidence machine.  ICML.
  * Obozinski & Bach, (2016) Obozinski, G., & Bach, F. 2016.  A unified perspective on convex structured sparsity: Hierarchical, symmetric, submodular norms and beyond.  HAL Id : hal-01412385, version 1.
  * Papadopoulos et al. , (2002) Papadopoulos, H., Proedrou, K., Vovk, V., & Gammerman, A. 2002.  Inductive confidence machines for regression.  European Conference on Machine Learning.
  * Pedregosa et al. , (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cournapeau, D., Brucher, M., Perrot, M., & Duchesnay, E. 2011.  Scikit-learn: Machine Learning in Python.  J. Mach. Learn. Res.
  * Qin et al. , (2010) Qin, T., Liu, T-Y., & Li, H. 2010.  A general approximation framework for direct optimization of information retrieval measures.  Information retrieval.
  * Romano et al. , (2019) Romano, Y., Patterson, E., & Candes, E. J. 2019.  Conformalized quantile regression.  NeurIPS.
  * Shafer & Vovk, (2008) Shafer, G., & Vovk, V. 2008.  A tutorial on conformal prediction.  Journal of Machine Learning Research.
  * Virtanen et al. , (2020) Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., Evgeni, Burovski, Pearu, Peterson, Warren, Weckesser, Jonathan, Bright, et al. . 2020.  SciPy 1.0: fundamental algorithms for scientific computing in Python.  Nature methods.
  * Vovk et al. , (2005) Vovk, V., Gammerman, A., & Shafer, G. 2005.  Algorithmic learning in a random world.  Springer.
  * Xie & Huang, (2009) Xie, H., & Huang, J. 2009.  SCAD-penalized regression in high-dimensional partially linear models.  The Annals of Statistics.
  * Xu & Xie, (2021) Xu, C., & Xie, Y. 2021.  Conformal prediction interval for dynamic time-series.  ICML.