Conformal inference for regression on Riemannian Manifolds

Conformal prediction is a powerful set of statistical tools that operates under minimal assumptions about the underlying model. It is primarily used to construct confidence sets that are applicable to a diverse range of problems in fields such as machine learning and statistics. Is unique in its ability to construct prediction regions that guarantee a specified level of coverage for finite samples, irrespective of the distribution of the data. This is particularly crucial in high-stakes decision-making scenarios where maintaining a certain level of coverage is critical.

Unlike other methods, conformal prediction does not require strong assumptions about the sample distribution, such as normality. The aim is to construct prediction regions as small as possible to yield informative and precise predictions, enhancing the tool’s utility in various applications.

The approach was first proposed by Vovk, Gammerman, and Shafer in the late 1990s, as referenced in [57]. Since its inception, it has been the subject of intense research activity. Originally formulated for binary classification problems, the method has since been expanded to accommodate regression, multi-class classification, functional data, functional time series, and anomaly detection, among others. Several applications of this method can be found in the book [5].

In the context of regression, conformal prediction has proven to be efficient in constructing prediction sets, as evidenced by works such as [61], [62], and [34]. To enhance the performance of these prediction sets, particularly to decrease the length of prediction intervals (when the output is one-dimensional), a combination of conformal inference and quantile regression was proposed in [53]. In [20], computationally efficient conformal inference methods for Bayesian models were proposed.

The method has also been extended to functional regression, where the predictors and responses are functions, rather than vectors, see for instance [41], [21] and [18]. In this case, the prediction regions take the form of functional envelopes that have a high probability of containing the true function.

In the field of classification, conformal prediction has been employed to tackle a broad spectrum of problems. These include image classification, as in [40], and text classification, as in [58]. For multi-class classification, a prevalent approach is to create prediction sets that have a high likelihood of encompassing the correct class. This can be realized via the use of ‘Venn prediction sets’, as outlined in [43]. These sets partition the label space into overlapping regions, each one corresponding to a distinct class.

In summary, conformal prediction provides powerful statistical tools and has been successfully applied to a wide range of problems in machine learning, statistics, and related fields. Its main advantage is its ability to provide distribution-free prediction regions that can be used in the presence of any underlying distribution of the data. As this field of research continues to evolve, it is expected to find even more applications in the future.

Although, as mentioned, it has been extended to various scenarios, there are no proposals or extensions to the case where the output is in a Riemannian manifold. This case is of cumbersome importance because there are certain types of data that, due to their inherent characteristics, must be treated as data on manifolds. A prominent example of this are covariance matrices (for instance, the volatility of a portfolio, see [8] and [28]), which are data in the manifold of positive definite matrices (see [14]). Another example is given by Vectorcardiograms that are condensed into data on the Stiefel manifold (see [13]). The outcome of performing a principal component analysis results in data on the Grassmannian manifold (see [27]), and the wind measurements can be represented as data on the cylinder (see [12]).

Other applications in image analysis are developed in [49] and more generally in machine learning in [24].

The extension to this context is not trivial since when the output $Y$ belongs to a Riemannian manifold $\mathcal{M}$, several arguments used in Euclidean conformal inference are not straightforwardly extendable. For instance, the Frechet mean on the manifold can be defined, but it does not always exist nor is it necessarily unique. Moreover, formulating a regression model on the manifold is not easy because of the absence of any additive structure, although there have been recent advances in this area (see for instance [50]).

We extend the methodologies outlined by [61] to the broader context of pairs $(X,Y)$ where $X\in\mathbb{R}^{d}$ and the response $Y$ has support included on a, smooth enough, $\ell$-dimensional manifold $\mathcal{M}$. [61] is focused on the scenario where $Y$ is a real-valued variable, constructing confidence sets—referred to as confidence bands—via density estimators. One of the key tools to get the consistency of the empirical conformal region to its population counterpart, as in [61], is Theorem 1. It states that the kernel-density estimator is uniformly consistent on manifolds. Once this is established, to adapt the ideas in [61], it must also be proved that the conditional density is uniformly bounded from above, this is done in Lemma 4. Finally, a further distinction from [61] is that—because the manifold may possess a nonempty boundary—points located well inside the manifold are handled differently from those situated near its edge. Here, we obtain an upper bound for the discrepancy between the empirical confidence set $\hat{\mathbf{C}}_{n}(X)$ and its theoretical counterpart $\mathbf{C}_{P}(X)$, in terms of $\nu$, the volume measure on $\mathcal{M}$. More precisely, as detailed in Theorem 2 (see also Remark 1) the probability that this discrepancy being larger than $n^{1/((\ell+3)(d+2))}$ is bounded from above by $A_{\lambda}n^{-\lambda}$, for any $\lambda>0$, $A_{\lambda}$ being a constant that depends only on $\lambda$, where $\ell$ denotes the dimensionality of $\mathcal{M}$ and $X\in\mathbb{R}^{d}$. Our main and stronger theorem is for compact manifolds, however, we considered in section 5 the case of non-compact manifolds.

The kernel-based density estimator proposed in [12] uses the Euclidean distance instead of the geodesic distance on the manifold, which may be unknown in some cases, see Equation (4). Due to the curse of dimensionality, kernel-based density estimators are generally unsuitable for high-dimensional problems. In Section 6, we discuss an approach to overcome this limitation that builds on the ideas introduced in [30]. In essence, instead of estimating the conditional density $p(y|x)$ locally (i.e., using only sample points close to $x$), the proposed method uses all sample points $(X_{i},Y_{i})$ whose estimated conditional densities $\hat{p}(Y_{i}|X_{i})$ are similar (in a sense to be defined).

Conformal inference is tackled in metric spaces in [44]; however, given the generality, convergence rates are not provided, as in the present work, but only consistency in probability, for the case where the confidence bands are constructed through regression.

In the following, $\aleph_{n}=\{(X_{i},Y_{i}):i=1,\dots,n\}$ denotes an i.i.d. sample of $(X,Y)$ with distribution $P$, where $X\in\mathbb{R}^{d}$ and $Y\in\mathcal{M}$. Here $(\mathcal{M},\rho)$ is a compact $\ell$-dimensional submanifold of $\mathbb{R}^{D}$. The case of non compact manifolds is discussed in Section 5. We further denote by $P_{X}$ the marginal distribution of $X$ and by $\text{supp}(P_{X})$ its support. We denote the volume measure on $\mathcal{M}$ by $\nu$ and use $\|\cdot\|$ to denote the Euclidean norm on $\mathbb{R}^{D}$. We denote by $\mu$ the $d$-dimensional Lebesgue measure in $\mathbb{R}^{d}$.

We will now give the set of assumptions that we will require. H1 to H4 are also imposed in [61]. Hypotheses H0 and H5 are imposed to guarantee the uniform convergence of the kernel-based density estimator of the conditional density.

• H0

  $\mathcal{M}$ is a $\mathcal{C}^{2}$ submanifold, and if $\partial\mathcal{M}\neq\emptyset$, then $\partial\mathcal{M}$ is also a $\mathcal{C}^{2}$ submanifold.

• H1

  The joint distribution $P$ has a density $p_{X,Y}(x,y)$ w.r.t. $\mu\times\nu$, and the marginal distribution $P_{X}$ has a density $p_{X}$ (w.r.t. $\mu$). We denote by $p(y|x)=p_{X,Y}(x,y)/p_{X}(x)$ the conditional density, where $p(y|x)=0$ if $p_{X}(x)=0$.

• H2

  $p_{X}$ is such that that there exist $b_{1},b_{2}$ such that $0<b_{1}\leq p_{X}(x)\leq b_{2}<\infty$ for all $x$ in the support of $X$.

• H3

  $p(y|x)$ is Lipschitz continuous as a function of $x$, i.e., there exists a constant $L>0$ such that $\|p(\cdot|x)-p(\cdot|x^{\prime})\|_{\infty}\leq L|x-x^{\prime}|$.

• H4

  There are positive constants $\epsilon_{0}$, $\gamma$, $c_{1}$, and $c_{2}$ such that for all $x$ in the support of $X$,

  $$c_{1}\epsilon^{\gamma}\leq P(\{y:|p(y|x)-t_{x}^{\alpha}|<\epsilon\}|X=x)\leq c_{2}\epsilon^{\gamma},$$

  for all $\epsilon\leq\epsilon_{0}$, where $t_{x}^{\alpha}$ is given by (2). Moreover $\inf_{x}t_{x}^{(\alpha)}\geq t_{0}>0$.

• H5

  $p(y|A_{k})$ has continuous partial second derivatives, for short $p(y|A_{k})\in\mathcal{C}^{2}$.

The hypothesis H0 is satisfied by a very wide range of manifolds, among them, the Stieffel and Grassmannian Manifold. The cone of positive definite matrices, among many others. It is not a restrictive hypothesis in applications.

Assumptions H1 to H5 impose regularity on the joint and conditional densities of the data. In particular, the assumption of Lipschitz continuity H3 ensures that small changes in the input $x$ lead to small changes in the output $y$. Assumption H2 implies that $X$ es compactly supported. Regarding Assumption H4, quoting [61], “is related to the notion of the ‘$\gamma$-exponent’ condition that was introduced by [51], and widely used in the density level set literature ([56, 52]). It ensures that the conditional density $p(\cdot|x)$ is neither too flat nor too steep near the contour at level $t_{x}^{\alpha}$, so the cut-off value $t_{x}^{(\alpha)}$ and the conditional density level set $\mathbf{C}_{P}(x)=L_{x}(t_{x}^{\alpha})$ can be approximated from a finite sample […]. Assumption H4 also requires that the optimal cut-off values $t_{x}^{(\alpha)}$ be bounded away from zero.”. Assumption H5 is required to get the uniform convergence of the Kernel based estimator of $p(y|A_{k})$ on manifolds, see [7]. These assumptions play a crucial role in the development and analysis of conformal prediction methods.

Let us introduce some important sequences of positive real numbers that will play a key rol all along the manuscript.

1. 1.

   $w_{n}=(\log(n)/n)^{1/(d+2)}$,

2. 2.

   $\gamma_{n}=\lfloor b_{1}nw_{n}^{d}/2\rfloor$, being $b_{1}$ the positive constant introduced in H2. Observe that $\gamma_{n}\to\infty$.

3. 3.

   Given $h_{n},c_{n}\to 0$ as $n\to\infty$, we will consider the subsquences $c_{\gamma_{n}}$, and $h_{\gamma_{n}}$ being $\gamma_{n}=\lfloor b_{1}nw_{n}^{d}/2\rfloor$ as before.

In this section we introduce a slightly modified version of the estimated local marginal density $\hat{p}^{(x,y)}(v|A_{k})$ and of the local conformity rank $\pi_{n,k}(x,y)$ originally introduced in [61]. We make the assumption that $\textrm{supp}(P_{X})=[0,1]^{d}$, and that $\mathcal{A}=\{A_{k},k=1,\dots,T\}$ is a finite partition of $[0,1]^{d}$ consisting of equilateral cubes with sides of length $w_{n}$. This is a quite common and technical assumption in conformal inference, but we can assume that, for instance, $X$ is such that $\textrm{supp}(P_{X})\subset[-R,R]^{d}$, for some $R>0$. It only changes the constants appearing in Theorem 2, which remains true.

Let $n_{k}=\sum_{i=1}^{n}\mathbbm{1}_{\{X_{i}\in A_{k}\}}$. Given a sequence $h_{n}\to 0$, a kernel function $K(\cdot):\mathbb{R}\to\mathbb{R}$, and $h_{n_{k}}$, we define

We aim to prove that $\hat{p}(v|A_{k})$ provides a uniform estimate of $p(y|A_{k})$ across $v$ and $k$. To achieve this, we need to ensure that there are sufficiently many sample points in each $A_{k}$. This is guaranteed by lemma 9 of [61], which states that if we choose $w_{n}=(\log(n)/n)^{1/(d+2)}$, then, with probability one, for all $n$ large enough,

Recall that $b_{1}$ and $b_{2}$ were defined in Assumption H2. In what follows, we assume that $n$ is sufficiently large to ensure (5).

The corresponding augmented estimate, based on $\aleph_{n}\cup\{(x,y)\}$ is, for any $(x,y)\in A_{k}\times\mathcal{M}$,

For any $(X_{n+1},Y_{n+1})=(x,y)\in A_{k}\times\mathcal{M}$, consider the following local conformity rank

As proved in Proposition 2 of [61], the set

has finite sample local validity, i.e., it satisfies (3). The finite sample local validity, as established by (7), does not require any assumptions and it holds under very general conditions.

The following result is the key theorem. Its proof is deferred to the Appendix. The proof follows the ideas used to prove Theorem 1 of [12]. It states that $p(y|A_{k})$ can be estimated uniformly by (4) for all $y$ that are far enough from the boundary of $\mathcal{M}$. Additionally, this estimation can be made uniformly across all $k$. We assume, as in [12], that $K$ is a Gaussian kernel. This restriction is, as in [12], purely technical. More recent references on kernel-density estimation on manifolds are [9, 10], [15, 16, 6, 63].

This result is of fundamental importance in conformal prediction, as it provides a way to estimate the conditional density of $Y$ given $A_{k}$ for any $k$ in a non-parametric way. It implies that the estimation error is uniform across all partitions $A_{k}$, which is a key requirement for conformal prediction methods. In particular, it allows us to construct conformal prediction regions that are valid with a given level of confidence.

The following theorem is the main result of our paper. It states that $\hat{\mathbf{C}}_{n}(X_{i})$ consistently estimates $\mathbf{C}_{P}(X_{i})$ uniformly for all $X_{i}\in\aleph_{n}$. The proof, which is deferred to the Appendix, see Subsection 9.1, is based on some ideas from the proof of Theorem 1 of [61]. Specifically, Lemma 1 and Lemma 3 are adapted from [61], while Lemma 2 is new and is required to adapt the proof of Lemma 3. Lemma 4 is the same as Lemma 8 of [61]. Finally, the proof of Theorem 2 uses the adapted lemmas and considers separately the sample points that are close to $\mathcal{M}$ and those that are far away from this boundary. The first set of points is shown to be negligible with respect to $\nu$ using techniques from geometric measure theory.

In this section, we discuss how to remove the compactness assumption imposed in Theorem 2, while retaining all the other hypotheses. We will prove that the set $\hat{\mathbf{C}}_{n}(X_{n+1})$ satisfies asymptotic conditional validity (see [30]). In other words, there exist random sets $\Lambda_{n}$ such that

and

According to Theorem 6 in [30], to obtain asymptotic conditional validity it suffices to prove that

It worth to be mentioned that asymptotic conditional validity is considerably weaker than (9) when $p(y|x)$ is bounded from above.

A different method for estimating the conformal region is introduced in [30]. Instead of using a partition of $\mathbb{R}^{d}$ on cubes to get a partition of $X_{1},\dots,X_{n}$, to build kernel-based estimators, this approach employs a new partition constructed as follows. Consider, as in [30], the functions

where $\hat{p}(y|x)$ is any estimator of $p(y|x)$, not necessarily kernel-based. Define the conditional $\alpha$-quantile of $Y|X$ by $q_{\alpha}(x)=H^{-1}(\alpha|x).$ Similarly, let $\hat{q}_{\alpha}(x)=\hat{H}^{-1}(\alpha|x)$, be the estimate of the conditional $\alpha$-quantile of $Y|X$.

Let $\mathcal{F}$ be a partition of $\mathbb{R}^{+}$. Then we create a partition $\mathcal{A}=\{A_{1},\dots,A_{T}\}$ of $X_{1},\dots,X_{n}$ by assigning $X_{i}$ and $X_{j}$ to the same set $A_{r}\in\mathcal{A}$ if and only if $\hat{q}_{\alpha}(X_{i})$ and $\hat{q}_{\alpha}(X_{j})$ lie in the same element of $\mathcal{F}$. We then build the kernel-based estimator (4) according to this new partition. This approach is particularly suitable for high-dimensional feature spaces.

We propose either the set $\hat{\mathbf{C}}_{n}(x)$ (as defined in (7)) or one of the algorithms presented in [30]. According to Theorem 25 in [30], to establish the asymptotic conditional validity of this proposal (in the sense of (10)), it is necessary to assume some smoothness of $H(z,x)$ (see Assumption 23 in [30]), as well as the compactness of $M$ and the existence of sequences $\eta_{n}=o(1)$ and $\rho_{n}=o(1)$ satisfying

This condition indeed holds in our case; however, proving it in detail would require rewriting the proof of Theorem 1 in [12]. We leave this verification to the reader. As in the non-compact case, the consistency result obtained is considerably weaker than (9).

Since the computational burden of approximating $\hat{\mathbf{C}}_{n}(x)$ is high, [61] proposes replacing that set with $\mathbf{C}_{n}^{+}(x)$, see (11), which, by including $\hat{\mathbf{C}}_{n}(x)$, has coverage of at least $1-\alpha$. Let $A_{x}\in\mathcal{A}$ the element that contains $x$.

1. 1.

   Let $\hat{p}(x,y)$ be the joint density estimator based on the subsample (of cardinality $n_{x}$) of points for which $X_{i}\in A_{x}$

2. 2.

   Let $Z_{i}=(X_{i},Y_{i})$ for $i=1,\dots,n_{x}$, and let $Z_{(1)},Z_{(2)},\dots,Z_{(n_{x})}$ denote the sample ordered increasingly by $\hat{p}(X_{i},Y_{i})$.

3. 3.

   Let $j=\lfloor n_{x}\alpha\rfloor$ and define

   $$\mathbf{C}_{n}^{+}(x)\;=\;\Bigl{\{}\,y:\;\hat{p}(x,y)\;\geq\;\hat{p}\bigl{(}X_{(j)},Y_{(j)}\bigr{)}\;-\;\frac{K(0)^{2}}{n_{x}\,h^{d+1}}\Bigr{\}}.$$

   (11)

We work within the conformal‐inference framework, where the response variable takes values on a Riemannian manifold while the covariates lie in $\mathbb{R}^{d}$. Under assumptions analogous to those in [61], we extend their results to show that the conditional oracle set $\mathbf{C}_{P}(x)$ can be estimated consistently. A central tool in our analysis is the almost‐sure, uniform consistency of the classical kernel density estimator on manifolds, as proved in [12]. We also cover the cases of manifolds with boundary and non‐compact manifolds, and we discuss computational strategies for handling high‐dimensional input data.

Our theoretical contributions are demonstrated through two simulation studies and two real-world data examples.

Finally, one could explore alternative nonconformity scores beyond the KDE-based criterion we employ—for example, quantile-based or regression-based scores—each of which demands a nontrivial adaptation of the Euclidean theory to the manifold setting.