Universal Prediction Band via Semi-Definite Programming

Universal Prediction Band via Semi-Definite Programming

\affil[1]{University of Chicago, Booth School of Business}

Keywords--- Uncertainty quantification, heteroscedasticity, nonparametric prediction band, semi-definite programming, sum-of-squares.

Introduction

A plausible criticism from the statistics community to modern machine learning is the lack of rigorous uncertainty quantification, with perhaps the exception in conformal prediction vovk2005algorithmic, lei2018DistributionFreePredictive, romano2019ConformalizedQuantile. Instead, the machine learning community would argue that conventional uncertainty quantification based on idealized distributional assumptions may be too restrictive for real data. Nevertheless, without a doubt, uncertainty quantification for predictive modeling is essential to statistics, learning theory, and econometrics. This paper will resolve the above inference dilemma by introducing a new method with provable uncertainty quantification via semi-definite programming. This paper provides an alternative approach to the now-standard conformal prediction for uncertainty quantification, with novel theoretical insights and computational advantages. The proposed method learns a data-adaptive, heteroscedastic prediction band that is: (a) universally applicable without strong distributional assumptions, (b) with desirable theoretical coverage with or without any user-specified predictive model, and (c) easy to implement via standard convex programs (when used in conjunction with a wide range of positive-definite kernels).

Let $(x, y) \in \mathcal{X} \times \mathbb{R}$ be the covariates and response data pair drawn from an unknown probability distribution $\mathcal{P}$. There are plenty of regression or predictive models---denoted by $\mathsf{m}_0(x)$---that estimate $\mathsf{m}(x) := \mathop{\mathbb{E}}[\by|\bx = x]$ sufficiently well with finite data. However, to make downstream decisions reliable, a good prediction band quantifying the uncertainty in $|\by - \mathsf{m}_0(\bx)|$ with provable coverage is urgently needed. The prediction band is of particular relevance to complex machine learning models that construct $\mathsf{m}_0(x)$ in a less transparent way, such as deep neural networks and boosting machines. This paper makes progress in filling in such a gap: we estimate a nonparametric, heteroscedastic prediction band $\widehat{\mathsf{PI}}(x)$ that enjoys provable coverage with minimal data assumptions for any predictive model. Our approach can be viewed as a novel variance interpolation with confidence and leverages techniques from sum-of-squares relaxations for nonparametric variance estimation. On a non-technical level, this paper enriches the toolbox of applied researchers with a theoretically justified new methodology for uncertainty quantification and visualization, as in conformal prediction.

Semi-definite Programs and Prediction Bands

We introduce our procedure for constructing the predictive band in this section. Let $K(\cdot, \cdot): \mathcal{X} \times \mathcal{X} \rightarrow \mathbb{R}$ be a continuous symmetric and positive-definite kernel function. Given $n$ data pairs $\{ (x_i, y_i) \}_{i=1}^n$ and the corresponding kernel matrix $\bK \in \mathbb{S}^{n \times n}$ with $\bK_{ij}=K(x_i, x_j)$, our prediction band is constructed based on the following semi-definite program (SDP)

where the optimization variable $\bB \in \mathbb{S}^{n\times n}$ is a symmetric positive semi-definite (PSD) matrix, $\mathsf{m}_0(\cdot): \mathcal{X} \rightarrow \mathbb{R}$ is a given predictive model (user-specified), and $\sK_i \in \mathbb{R}^n$ denotes the $i$-th column of the kernel matrix $\bK$. Given the estimated $\widehat{\bB}$, the prediction band, $\widehat{\mathsf{PI}}(x)$ that maps each $x$ to an interval, can be constructed accordingly

with $\delta \in \mathbb{R}$ being a scalar quantifying confidence. $\delta$ can be later calibrated for exact coverage control, and may be set as $0$ if $n$ is large. Here $\widehat{\mathsf{v}}(x)$ estimates the variability in the “deviations” $e_i := y_i - \mathsf{m}_0(x_i)$. A few remarks on such deviations are in place.

Let $K^{\mathsf{m}}$ and $K^{\mathsf{v}}$ specify two kernel functions, corresponding to the conditional mean and variance functions respectively. $\bK^{\mathsf{m}}, \bK^{\mathsf{v}} \in \mathbb{S}^{n \times n}$ denote empirical kernel matrices on finite data with size $n$. For any $\gamma\geq 0$, the following convex SDP program constructs the prediction band and the conditional mean function simultaneously

where the optimization variables are $\bB \in \mathbb{S}^{n\times n}$ and $\alpha \in \mathbb{R}^n$. Given the solution $\widehat{\bB}$ and $\widehat{\alpha}$, the $\widehat{\mathsf{PI}}(x)$ is constructed as

A Numerical Illustration

Before diving into the motivations behind the above SDPs (Sec. (ref)) and corresponding theory (Sec. (ref)), let us first visually illustrate the empirical performance of the constructed prediction bands on a toy numerical example. A complete simulation study comparing our methods and conformal predictions will be deferred to Section (ref). Details of the data generating processes will be elaborated therein as well. The quick exercise is to showcase that convex programs (ref) and ((ref)) are easy to implement using standard optimization toolkits (say, CVX cvx), and construct flexible prediction bands with desired coverage properties. As a motivating example, we try out the SDP ((ref)), which simultaneously estimates the conditional mean and variance functions. A minimal 10-line Python implementation is provided in Listing. (ref).

The first example is a linear model with heteroscedastic error: the conditional mean $\mathsf{m}(x)$ being a linear function and variance $\mathsf{v}(x)$ being a quadratic (with the conditional variance generated from a uniform distribution). We generate a training dataset of size $n = 40$ and compare the coverage among three methods: (a) SLR, simple linear regression, (b) SDP1, a SDP ((ref)) with linear kernels $K^{\mathsf{m}}$ and $K^{\mathsf{v}}$ for both mean and variance functions, and (c) SDP2, a SDP ((ref)) with a linear $K^{\mathsf{m}}$ and a (degree-$3$) polynomial kernel $K^{\mathsf{v}}$. The coverage is compared on the same test dataset of size $N=800$. Here $\gamma = 0.1$. See Fig. (ref) for details and Table. (ref) for coverage statistics.

The second example is a non-linear, heteroscedastic error model: mean $\mathsf{m}(x)$ and variance $\mathsf{v}(x)$ functions lying in a reproducing kernel Hilbert space (RKHS) with a radial basis function (rbf) kernel. Here $n = 60$ training samples and $N = 800$ test samples are generated. Three methods being compared are: (a) SLR, (b) SDP1, rbf kernel for $K^{\mathsf{m}}$ and linear kernel for $K^{\mathsf{v}}$, and (c) SDP2, rbf kernels for both $K^{\mathsf{m}}$ and $K^{\mathsf{v}}$, summarized in Fig. (ref) and Table. (ref). Here $\gamma = 1$.

These two numerical examples are minimal yet informative. In Fig. (ref), SLR misspends a wide prediction bandwidth on data where the conditional variances are small yet fails to capture the large conditional variance cases, resulting in the overall coverage of $86\%$ and a median bandwidth of $8.21$. SDP1/SDP2 re-distributes the bandwidth budget leveraging the heteroscedastic nature and achieves an improved coverage $91\%$/$94\%$, with a smaller median bandwidth of $7.47$/$7.30$. Such an effect is even more pronounced in Fig. (ref). Observe first that in SDP2, the prediction band constructed by (ref) almost perfectly contours the heteroscedastic variances, thus achieving a $>99\%$ prediction coverage with a merely $3.35$ median bandwidth, in contrast to SLR with a $96\%$ coverage and a $4.80$ bandwidth. Second, a better conditional variance estimate also improves performance in learning the conditional mean, as seen in the differences between \textcolor{blue}{Blue lines} and \textcolor{red}{Red lines}. The errors are also numerically summarized in the column “MSE” of Table (ref). Leveraging the heteroscedasticity in data, our prediction band distributes the bandwidth in a data-adaptive way, thus improving the overall coverage.

Sum-of-Squares, Interpolation and Connections to Literature

The SDPs proposed in (ref) and ((ref)) are inspired by recent advancements in optimization and learning theory. We will elaborate on the connections to related works and explain the innovations in our approach. We start with some basic observations about the SDPs. First, when $\alpha$ is not an optimization variable, ((ref)) recovers ((ref)). Second, the constraints set of ((ref)) is always non-empty since $\alpha = 0, ~\bB = \max_i \|\sK_i^{\mathsf{v}} \|^{-2} y_i^2 \cdot \bI$ is feasible.

\paragraph{Sum-of-Squares and Phase Retrieval} As shown in Prop. (ref), the infinite-dimensional SDP with a nuclear-norm minimization is equivalent to (ref),

Here $\mathcal{H}^{\mathsf{m}}, \mathcal{H}^{\mathsf{v}}$ denote two RKHSs where the conditional mean and variance functions reside. $\phi_{x_i} \in \mathcal{H}$ is the feature map w.r.t. the Hilbert space $\mathcal{H}$ and $\langle \cdot,\cdot \rangle_{\mathcal{H}}$ is the Hilbert space inner-product. We call it the infinite-dimensional SDP since the optimization variables $(\beta, \bA)$ are (function, operator) rather than finite-dimensional (vector, matrix). A few remarks are in place. First, if the kernel $K^{\mathsf{m}}$ is universal, $\langle \phi_{x}^{\mathsf{m}}, \beta \rangle_{\mathcal{H}^{\mathsf{m}}}$ is dense in $L^2$ and hence can universally approximate all conditional mean function $\mathsf{m}(x)$. Second, as for the conditional variance which has positivity constraints over a continuum $x \in \mathcal{X}$ with $0 \leq \mathsf{v}(x) = \big( y - \mathsf{m}(x) \big)^2$, we relax the positivity constraints using a sum-of-squares form

It turns out that when $K^{\mathsf{v}}$ is universal, the above sum-of-squares function can approximate all smooth, positive functions fefferman1978positivity, bagnell2015learning, marteau-ferey2020NonparametricModels, thus explaining the name “universal” in the title. Remark that sum-of-squares optimization lasserre2001global for nonparametric estimation has recently been considered; see bagnell2015learning, marteau-ferey2020NonparametricModels, curmei2020shape. The further relaxation changing from equality in ((ref)) to inequality will be discussed in the next paragraph. Last, the minimum nuclear norm objective translates to a particular form of “minimum bandwidth” in the prediction band as $\mathsf{v}(x) = \langle \phi_x^{\mathsf{v}}, \bA \phi_x^{\mathsf{v}} \rangle_{\mathcal{H}^{\mathsf{v}}}$. In language, for all prediction bands that shelter the data, (ref) aims to find the one with minimum bandwidth.

The curious reader may wonder where the nuclear norm $\| \bA \|_{\star} $ arises from. The first reason is conceptual: the nuclear norm is a relaxation for rank, and the procedure is to minimize the number of factors (rank) that explain the variance. The second reason is a connection to phase retrieval: specify $K^{\mathsf{m}}(x, x') \equiv 0$ and $K^{\mathsf{v}}(x, x') = \langle x, x' \rangle$ (the linear kernel with $\phi_{x}^{\mathsf{v}} = x$), and force the inequality constraints to be equal, our SDP in (ref) is equivalent to phase retrievel

Conceptually, the minimum nuclear norm procedure estimates the smallest number of factors that could generate the variance.

\paragraph{Min-norm Variance Interpolation with Confidence} Now we discuss the tuning parameter $\gamma \in [0, \infty]$ and reveal the connection to the recent min-norm interpolation literature liang2020JustInterpolate, bartlett2020BenignOverfitting, bartlett2021DeepLearning, ghorbani2020LinearizedTwolayers, montanari2020GeneralizationError, liang2021InterpolatingClassifiers. In the limit of $\gamma \rightarrow 0$, (ref) reduces to the familiar min-norm interpolation with kernel $\bK^{\mathsf{m}}$ (whenever it has full rank, since optimal $\bB=0$)

In the limit of $\gamma \rightarrow \infty$, (ref) reduces to (since optimal $\alpha=0$)

Now it is clear what the role of the tuning parameter $\gamma$ is: it trades off the conditional mean $\mathsf{m}(x)$ and variance $\mathsf{v}(x)$ to explain the variability in $y$'s witnessed on the data. A small $\gamma$ aims to use a complex mean $\mathsf{m}(x)$ and a parsimonious variance $\mathsf{v}(x)$ to explain the overall variability, and vice versa.

From the above discussion, it is also clear that the SDP ((ref)) can be viewed as a min-norm variance interpolation with confidence. Instead of having the typical equality constraints in interpolation $$ \big\langle \sK_i^{\mathsf{v}} , \bB \sK_i^{\mathsf{v}} \big\rangle = \big( y_i - \big\langle \sK_i^{\mathsf{m}}, \alpha \big\rangle \big)^2 \;, $$ which violates the disciplined convex programming ruleset (due to the quadratic form on the RHS), we further relax to inequality constraints to incorporate additional “confidence” (and to make the problem convex at the same time) $$ \big\langle \sK_i^{\mathsf{v}} , \bB \sK_i^{\mathsf{v}} \big\rangle \geq \big( y_i - \big\langle \sK_i^{\mathsf{m}}, \alpha \big\rangle \big)^2 \;. $$ As we shall see, the notion of confidence in this variance interpolation is closely related to the notion of margin in classification bartlett1998BoostingMargin, liang2020PreciseHighdimensional.

\paragraph{Support Vector Regression} We illustrate that minor modifications to our SDP formulation lead to other problems, including support vector regression and kernel ridge regression. Specify the variance kernel as the trivial one $K^{\mathsf{v}}(x,x') = \mathbbm{1}(x=x')$, then the decision variable $\mathbf{B}$ only matters in its diagonal component, and our SDP (ref) reduces to the kernel ridge regression

Moreover, a slight modification of (ref) is to use the absolute deviation rather than the squared deviation in the constraints, namely

In this case support vector regression is exactly our procedure with the specification $K^{\mathsf{v}}(x,x') = \mathbbm{1}(x=x')$,

In summary, our SDP generalize beyond support vector machines, with the new non-trivial variance component for rigorous uncertainty quantification for heteroscedastic data.

Literature Review

There are increasingly many approaches proposed to address the uncertainty quantification dilemma in machine learning due to its significance and centrality. However, very few methods are theoretically grounded and universally applicable to the best of our knowledge. Many approaches are merely heuristics or data visualization tools. This section divides related theoretical studies in the literature into two categories and discusses how our method significantly differs from them and could potentially lead to a stronger theory.

\paragraph{Conformal Prediction} Based on the exchangeability of data and a user-specified nonconformity measure, vovk2005algorithmic, shafer2008TutorialConformal pioneered the field of conformal prediction, which uses past data to determine precise levels of confidence in new predictions. To some extent, the elegant theory of conformal prediction, motivated by online learning and sequential prediction, resolved the uncertainty quantification dilemma. The conformal prediction algorithm (see, for instance shafer2008TutorialConformal Section 4.3) usually requires to enumerate over all the possibilities of $z = (x, y) \in \mathcal{X} \times \mathcal{Y}$, and for each possibility, calculate $n$ nonconformity measures via the leave-one-out method. Therefore, the total budget is $n \times |\mathcal{Y}| \times |\mathcal{X}|$, which can be expensive for continuous $y$ and multi-dimensional $x$. Much of the above computation can be saved if additional information about the metric structure in $x \in \mathcal{X}$ can be leveraged. In contrast, our SDP approach constructs the prediction band over all the $x$'s at once, leverages the metric structure in $\mathcal{X}$, and suffers at most a computational budget of $n^2$. An additional key feature in our approach is in the coverage theory established in Theorem (ref): the prediction band has coverage probability $>95\%$ on a new data point $(\bx, \by)$, for $99.9999\%$ dataset $\{(x_i, y_i)\}_{i=1}^n$ of size $n$ drawn from the same distribution. Such a distinction on “confidence” vs. “probability” is discussed extensively in Section 2.2 of shafer2008TutorialConformal.

There has been a vast line of recent work on extending the conformal prediction idea further to address the bottlenecks above in the regression setting. The body of work proliferates, and we certainly cannot do justice here. lei2018DistributionFreePredictive alleviates the computational burden of the conformal prediction by introducing the sample-splitting technique. Remarkably, theory on the bandwidth is also studied in lei2018DistributionFreePredictive, thus providing an angle to probe the statistical efficiency. romano2019ConformalizedQuantile studies the problem that existing conformal methods can form nearly constant or weakly varying bandwidth and provide conservative results. romano2019ConformalizedQuantile proposes conformalized quantile regression to address this issue. One shared feature of our SDP approach is that the prediction band is fully adaptive to heteroscedasticity. Finally, we would like to emphasize that conformal prediction constructs prediction bands in a numerical, black-box fashion without a structural understanding of the variance function. In contrast, our SDP approach provides a transparent and efficient way of learning the variance function, a complementary contribution to the conformal literature.

\paragraph{Residual Subsampling and Quantile Regression} An alternative approach for uncertainty quantification that leverages the metric structure in $x \in \mathcal{X}$ is to resample the residuals locally. Typically, this is done by first fitting a predictive model $\mathsf{m}_0(x)$, and defining a local neighborhood around a new data $\bx$, then subsampling the residuals for uncertainty quantification via (conditional) quantiles. The validity of the above approach crucially depends on how many “similar residuals” to pool information from. However, the curse of dimensionality comes in since data points are far from each other in high dimensions, posing challenges in pooling the residuals. One can also use either the obtained residuals or the original responses $y$ to fit a conditional quantile regression model koenker1978regression, koenker2001quantile, belloni2011, belloni2019conditional, $\widehat{\xi}^{\tau}(\cdot) := \arg\min_{\xi} \tfrac{1}{n} \sum_{i=1}^n \rho_{\tau}\big(y_i - \xi(x_i) \big)$ where $\tau \in (0,1)$ is a quantile parameter, $\rho_{\tau}(\cdot):\mathbb{R} \rightarrow \mathbb{R}_+$ is the tilted absolute value function, and $\widehat{\xi}^{\tau}(\cdot): \mathcal{X} \rightarrow \mathbb{R}$ is the estimated conditional quantile function. However, it is not guaranteed that over all $x \in \mathcal{X}$, the estimated conditional quantile function satisfies $\widehat{\xi}^{\tau_1}(x) < \widehat{\xi}^{\tau_2}(x)$ for two quantiles $\tau_1 < \tau_2$. In other words, it is entirely possible that for several $x$'s, the conditional prediction intervals are empty chernozhukov2010quantile.

Theory for Uncertainty Quantification

In this section, we develop a theory for the coverage property of the prediction band constructed above, under the mild assumption that the data are i.i.d. drawn with $(\bx, \by)\sim\mathcal{P}$. To highlight the main arguments in a simple form, let us consider the setting $\mathsf{m}_0(x) \equiv 0$. Otherwise, the same proof follows by replacing $\by$ with $\by - \mathsf{m}_0(\bx)$. Define the corresponding prediction band, with a confidence parameter $\delta \in (0,1]$

We need the following assumptions before stating the theorem, where $(\bx, \by)\sim\mathcal{P}$ and $\text{C}>0$ denotes a universal constant.

\paragraph{Discussion of Assumptions} All the above assumptions are mild. The eigenvalue decay in [S1] is almost identical to $\operatorname{Tr}(\mathcal{T}) < \infty$ (bounded trace integral operator). [S2] is also minimal, since it is only not true when there is no variability in $\by|\bx = x$. [S3] is the most stringent one, which requires the variability of $\by$ to exhibit a certain tail-decay. For small $\omega \in (0,1)$, [S3] can be much milder than exponential tail-decay. Bounded or Gaussian $\by|\bx=x$ satisfies [S3] with arbitrarily large $\omega$ or $\omega = 2$, respectively. With some extra work, [S3] can be relaxed to the case of a sufficiently rapid polynomial tail-decay. \textsf{[S2]} can be relaxed to restricting only to $x$ with $\mathop{\mathbb{P}}\big[ \by^2 > 0 ~ |~ \bx = x \big] = 1$.

What does the Theorem Entail

A few remarks are in order before we sketch the proof of Theorem (ref).

\paragraph{Coverage} First, the above theorem says that the prediction band constructed using the SDP based on a dataset of size $n$, will correctly cover a fresh data point $(\bx, \by) \sim \mathcal{P}$ drawn from the same distribution, with a non-asymptotic coverage probability (on the new data $\bx, \by$) $$ 1 - \delta^{-1} \sqrt{\tfrac{\log^3(n)}{n}} \;. $$ With $\delta = 0.5$, the bandwidth $\text{Length}\big[ \widehat{\mathsf{PI}}(x) \big] = 2.45\sqrt{\widehat{\mathsf{v}}(x)}$ is at a heteroscedastic level adaptive to $x$ with corresponding coverage probability at least $1 - O^\star(\sqrt{\tfrac{1}{n}})$. Here $O^\star$ hides polylog factors. The coverage can be arbitrary close to $1$ with large $n$ without the need of increasing $\delta$, which is in clear distinction to the conventional wisdom that coverage $1$ can only be possible with an increasing $\delta$ regardless of $n$. Again, we emphasize that the above coverage guarantee holds essentially on $99.9999\% \ll 1-n^{-10}$ of the datasets $\{(x_i, y_i\}_{i=1}^n$. In Section (ref), we propose a rigorous calibration algorithm to choose $\delta^\star(\alpha)$ to achieve a constant coverage level $1-\alpha \in (0, 1)$.

\paragraph{Optimality} If one wishes to obtain the classic $95\%$ coverage probability, then choosing $\delta = O^\star(\sqrt{\tfrac{1}{n}})$ suffices, which translates to

Recall that in classic simple linear regression, the prediction interval is of length

with $\hat{s}=\sqrt{\tfrac{\Sigma_i \hat{e}_i^2}{n-2}}$ being the estimated residual standard error. The fact that (ref) and ((ref)) share the $\sqrt{\tfrac{1}{n}}$ fluctuation seems to indicate the optimality of our Theorem (in terms of the dependence on $n$).

\paragraph{Data Adaptivity} Curiously, the objective value of the convex optimization program quantifies the uncertainty of the prediction band: a smaller $\widehat{\mathsf{Opt}}_n$ implies (a) a better confidence/coverage guarantee and (b) a narrower prediction band overall. More importantly, the $\widehat{\mathsf{Opt}}_n$ can be calculated directly from data! We find such an optimization/inference interface exciting: the data-adaptive bound lets us know the coverage guarantee specific to the current dataset. Put differently, the convex program constructs the prediction band via its solution and at the same time, reveals the confidence via its objective value. Since $\widehat{\mathsf{Opt}}_n$ is a function of the dataset, our Theorem reveals which dataset allows for a better prediction band. Remark that $\widehat{\mathsf{Opt}}_n = \| \widehat{\mathsf{v}} (\cdot) \|_{\star}^2$ is also a particular norm of the heteroscedastic variance function, quantified by the nuclear norm of the associated PSD operator $\widehat{\bA} \succeq 0$ with $\widehat{\mathsf{v}}(x) = \langle \phi_x, \bA \phi_x \rangle_{\mathcal{H}}$. Curiously, a simpler variance function $\widehat{\mathsf{v}} (x)$ (with a small norm) will simultaneously result in a narrower band and better coverage. We emphasize that the above discussion is in sharp contrast to the conventional wisdom that a narrow band usually leads to poor coverage guarantees.

Calibration and Coverage Control

One nice feature about conformal prediction is that it directly operates on a user-specified coverage level (e.g., $95\%$), albeit the resulting procedure only achieves some form of coverage guarantee in the marginal sense. In contrast, the coverage level guarantee of the current SDP approach in Theorem (ref) is in an inequality form with a mild dependence on the non-explicit universal constant $C_{\gamma, \xi, \eta, \omega}$; however, the coverage is in a stronger conditional sense conditioned on $\{ (x_i, y_i) \}_{i=1}^n$. The constant won't significantly affect the coverage guarantee in the large $n$ setting; nevertheless, curious readers may wonder if more transparent control could be achieved by tuning $\delta$. This section provides a theoretically justified calibration procedure in choosing $\delta$ (in the SDP band) to control coverage at a user-specified level $1-\alpha \in (0, 1)$, as in conformal prediction. Such fine calibration on $\delta$ can be helpful numerically in the moderate $n$ and constant $\alpha$ setting (e.g., $5\%$).

The calibration idea is based on sample-splitting. Split the samples into two parts, the training set $\{ (x_i, y_i) \}_{i=1}^{n} $ and the calibration set $\{ x'_j, y'_j\}_{j=1}^m$, with in total $n+m$ data points drawn i.i.d. from $\mathcal{P}$. The training set will be used to construct prediction band $\widehat{\mathsf{PI}}(\cdot, \delta)$. The calibration data set will be used to choose $\delta$ to calibrate coverage.

\paragraph{Calibration Procedure} Suppose that there exsits large enough constants $\Delta, L>0$ such that

The calibration uses a dyadic search to select $\delta^\star(\alpha) \in [-1, \Delta]$ with the set $\{ x'_j, y'_j\}_{j=1}^m$. The goal of the calibration is to ensure $\mathop{\mathbb{P}}_{(\bx, \by)\sim \mathcal{P}}\big[ \by \in \widehat{\mathsf{PI}}(\bx, \delta^\star(\alpha)) \big] \geq 1- \alpha$.

The following result derives the theoretical ground for the calibration procedure. We defer the proof to Section (ref).

Intuition and Proof Sketch

We first explain the key intuition before presenting the details of the proof sketch. The proof first leverages a representation theorem that relates the finite-dimensional (kernelized) SDP to an infinite-dimensional SDP to decouple the dependencies among $x_i$'s. Next, we propose to use empirical process theory to analyze the prediction coverage, inspired by the margin-based analysis originally done in analyzing classification. Finally, in controlling the uniform deviations between empirical and population coverage, we use properties of the PSD operators and conditional quantile functions. The additional calibration procedure in fine-tuning the confidence parameter $\delta^\star(\alpha)$ for a fixed coverage level $1-\alpha \in (0, 1)$ also hinges on uniform deviation arguments. Our analysis fundamentally relies on empirical process theory and is crucially different from conformal prediction analysis (based on exchangeability). Since the SDP provides a rigorous coverage guarantee as conformal prediction, we hope the new proof idea opens new doors to study uncertainty quantification.

Now we sketch the proof of Theorem (ref). Observe that by definition

Define a hinge function $h_{\delta}(t): t \mapsto \max\{ \tfrac{1+\delta}{\delta}(1-t) , 0 \}$, we have

and thus

Define a real positive function (indexed by $\bA$) on the data $z=(x, y)$, $ f_{\bA}(z): z \mapsto \big\langle \tfrac{\phi_x}{y}, \bA \tfrac{\phi_x}{y} \big\rangle_{\mathcal{H}} \;. $ Here $\tfrac{\phi_x}{y} \in \mathcal{H}$ lies in the RKHS, and $\bA: \mathcal{H} \rightarrow \mathcal{H}$ is a PSD operator. Define a sequence of function spaces according to its nuclear-norm radius $ \mathcal{F}_{k} := \{ f_{\bA} ~:~ 2^{k-1}< \| \bA \|_\star \leq 2^k \} $ for all $k \in \mathbb{N}$ and $\mathcal{F}_{0} := \{ f_{\bA} ~:~ \| \bA \|_\star \leq 1 \}$.

With the Prop. (ref) establishing the equivalence between the kernelized SDP and the infinite-dimensional SDP, $\widehat{\bA} = \sum_{i,j} \widehat{\bB}_{ij} \phi_{x_i} \otimes \phi_{x_j}$, we know that $y^{-2} \widehat{\mathsf{v}}(x) = f_{\widehat{\bA}}(z)$. There exists a $k \in \mathbb{N}$ such that $f_{\widehat{\bA}} \in \mathcal{F}_k$ with $ 2^{k-1} \leq \widehat{\mathsf{Opt}}_n < 2^k$, and thus we continue to bound

For term (i), recall the optimality condition of (ref),

which further implies $h_{\delta} \circ f_{\widehat{\bA}}(z_i) = 0$ for all $i=1, \cdots, n$. Therefore term (i) is zero.

For term (ii), we will use the high probability symmetrization in Prop. (ref). Introduce i.i.d. Rademacher variables $\{ \epsilon_i\}_{1}^n$ independent of the data. Note that we only need to consider $k \leq k_0$ such that $2^{k_0} = [\log(n)]^{ \text{c}_{\omega}}$, where we use the upper estimate on $\widehat{\mathsf{Opt}}_n$ obtained in Prop. (ref), which is implied by Assumption [S3]. With probability at leat $1 - 2\exp(-t)$ on the data $\{ z_i\}_{1}^n$, uniformly over all $k \leq k_0$

where the last step follows from the duality between the nuclear norm and operator norm. Before getting into the deviation term (iii) (originated by Prop. (ref), formally upper bounded in (ref)), first recall $2^k \leq 2 (\widehat{\mathsf{Opt}}_n \vee 1)$, we know

Similarly, by Prop. (ref), the deviation term (iii) can be bounded by $6 (\widehat{\mathsf{Opt}}_n \vee 1) \cdot \sup_{x,y}\| \tfrac{\phi_x}{y} \|_{\mathcal{H}}^2 \cdot \sqrt{\tfrac{k_0+t}{n}} $ with $k_0 = \text{c}_{\omega} \log\log(n)$.

To bound the expected operator norm for the above random matrix, namely term (iv), we rely on matrix Bernstein's inequality plus a truncation technique. Observe that

and that

Naively applying the matrix Bernstein inequality, one would expect the term (iv) to behave like $ \sup_{x,y}\| \tfrac{\phi_x}{y} \|_{\mathcal{H}}^2 \cdot \big( \sqrt{\tfrac{t}{n}} \vee \tfrac{t}{n} \big) $ with probability $1 - \text{dim}(\phi_x) \cdot \exp(-t)$ on $\{ \epsilon_i \}_{1}^n$. This is educative yet wrong, since $\text{dim}(\phi)$ is infinity. To make things rigorous, we rely on a truncation technique to look at a finite-dimensional version $\phi_x^{\leq m}$ truncated at a level $m = \text{poly}(n)$ to apply matrix Bernstein, and then estimate the remaining contribution from $\phi_x^{> m}$ by the eigenvalue decay in Assumption [S1]. With details given in Prop. (ref), we derive that

The final piece of the puzzle lies in the term $\sup_{(x,y) \in \mathrm{dom}(\mathcal{P})}\| \tfrac{\phi_x}{y} \|_{\mathcal{H}}^2$, which appears in both the main term (iv) and deviation term (iii). It is not true that a.s. for all $x, y$, the above term is bounded. To resolve this issue, we rely on a conditional quantile technique. Introduce the conditional quantile function $Q_{\by^2|\bx = x}(\cdot): [0,1] \rightarrow \mathbb{R}_+$ for the conditional random variable $\by^2|\bx = x$. Let's only look at data $(x_i, y_i)$'s lying in the region

and denote $\mathcal{P}|_{\Omega}$ as the conditional distribution of data $(\bx, \by)$ conditioning on the region $\Omega$. Claim that for any $\widehat{\mathsf{PI}}(\bx, \delta)$

This is based on two facts. First,

regardless of the ordering of $Q_{\by^2|\bx = x}(1-\eta)$ and $(1+\delta) \widehat{\mathsf{v}}(x)$. Second, conditioning on $\Omega$ does not change the marginal distribution of $\bx$ due to the quantile construction, namely $\mathcal{P}_{\bx} |_{\Omega} \equiv \mathcal{P}_{\bx}$. Marginalizing (ref) over $x$ proves the above claim.

The inequality ((ref)) makes the analyses upper bounding $\text{\small (LHS)}$ from (ref)-((ref)) applicable, with the changes: (a) $\mathcal{P}|_{\Omega}$ replacing $\mathcal{P}$, and (b) $\widehat{\mathop{\mathbb{E}}}$ denoting average over data points inside $\Omega$ rather than the whole dataset. With the conditioning on $\Omega$, Assumption [S2] implies $ Q_{\by^2|\bx = x}(1-\eta) \geq \xi \cdot K(x, x), $ and thus

Now, we only need to estimate the effective sample size inside $\Omega$ to complete the analyses. By the quantile construction, $\mathop{\mathbb{P}}_{(\bx, \by) \sim \mathcal{P}}\big[ (\bx, \by) \in \Omega \big] = \eta$, a simple Bernstein's inequality asserts that

with probability at least $1 - \exp(- \text{c}_{\eta} \cdot n)$ on $\{ z_i\}_{1}^n$.

Finally, plug (ref) into upper bounds on terms (iii) and (iv), with $\tfrac{\eta}{2} \cdot n$ replacing $n$ in (ref) and (ref), we have proved that

with probability at least $1 - \exp(- \text{c}_{\eta} \cdot n) - 2\exp(-t)$ on $\{ z_i\}_{1}^n $.

Numerical Studies

We now study the numerical performance of our procedure.

Empirical Example: Comparison to Conformal Prediction

This section compares our SDP methods to the conformal prediction methods, measuring the coverage and statistical efficiency of the prediction bands. We compare five methods on two simulated data sets, including (a) standard prediction band using simple linear regression (SLR), without accounting for heteroscedasticity, as a baseline; (b) SDP prediction bands proposed in this paper, with two specifications of the kernels, denoted as SDP1 and SDP2; (c) conformal prediction bands, including full conformal prediction (CF) and split conformal prediction (SplitCF), see lei2018DistributionFreePredictive. For each method, we report the coverage probability, efficiency statistics (including median length and average length), and finally, the mean squared error (MSE) of the estimated conditional mean function.

We first explain the two simulated data sets. Here the $\bx \sim \text{Unif}([-\sqrt{3}, \sqrt{3}])$, and the conditional distribution $\by| \bx = x \sim \epsilon \cdot \sqrt{\mathsf{v}(x)}$ where the independent error $\epsilon$ is either drawn from a standard normal $\text{N}(0,1)$ (Fig. (ref)) or a uniform distribution $\text{Unif}([-\sqrt{3}, \sqrt{3}])$ (Fig. (ref)). The conditional mean function $\mathsf{m}(x)=0$. The heteroscedastic variance function scales as $\mathsf{v}(x) = 1 + x + 4x^2$, depending on $x$. For each simulated data set, we generate $\{(x_i, y_i)\}_{i=1}^n$ i.i.d. from the above data generating process (DGP), with $n = 50$ as the training data, marked by \textcolor{blue}{Blue dots}. The test data set are drawn from the same DGP, with $N=500$ marked by \textcolor{red}{Red dots}. For the SDP1/SDP2, and SplitCF, an independent calibration data set with $n=50$ is used. The calibration set is used to choose $\delta$ in SDPs as in Algorithm (ref), and the homoscedestic conformal bandwidth as in SplitCF. We compare five methods that construct the prediction bands, illustrated by the \textcolor{blue}{Blue band}, on the Gaussian error data set in Fig. (ref), and Uniform error data set in Fig. (ref). For all methods, the desired coverage is set at $1-\alpha = 95\%$. For the SDPs in (ref) $\gamma = 10$. Table (ref) summarizes the coverage, efficiency, and estimation error.

Nearly all methods achieve $95\%$ desired coverage, with the only exception of SLR. The focus will be on comparing efficiencies, namely, which method estimates a smaller, truly heteroscedastic band in achieving the desired coverage. As seen visually in Fig. (ref)-(ref) and numerically in Table (ref), the two conformal methods, SplitCF and CF, estimates a conservative, wide prediction band that is almost homoscedastic. In contrast, both SDP approaches estimate desirable heteroscedastic bands that are on average much shorter, with the closest ($94.20\%$ and $95.40\%$) to the desired $95\%$ coverage. Finally, we would like to remark that all four methods SLR, SDP1, SDP2 and SplitCF are efficient to compute. Yet, the full conformal method CF involves discretizing input space, which is computationally intensive. Our empirical results show that the two conformal methods can be unnecessarily conservative and form bands of nearly constant width across $x$ romano2019ConformalizedQuantile.

Real Data Example: Fama-French Factors

In this section, we apply our method of constructing the prediction band to the celebrated three-factor dataset created by fama1993CommonRisk. We choose this dataset for three reasons: (a) financial data are known to suffer severe heteroscedasticity, (b) the factors are believed to be different sources explaining returns of diversified portfolios, thus when conditioned on one factor, the other factors should have large, heteroscedastic conditional variability, and (c) the factors---Market, Size, and Value--- correspond nicely to our common sense about the financial market for exploratory data analysis.

Let us first explain the data in plain language. The dataset consists of yearly and monthly observations of four variables from July 1926 to December 2020. The four variables are (a) Risk-free return rate ( RF), the one-month Treasury bill rate (i.e., interest rate), (b) Market factor ( MKT), the excess return on the market (i.e., market return minus interest rate), (c) Size factor ( SMB, Small Minus Big), the average difference in returns between small and big portfolios according to the market capitalization, and (d) Value factor ( HML, High Minus Low), the difference in returns between value and growth portfolios. We design two experiments, one focusing on the prediction coverage and bandwidth, and the other on exploring the role of the tuning parameter $\gamma$ in trading off mean and variance.

The first experiment aims to access the prediction coverage in the SDP (ref), using MKT (as $x$) to predict other variables (as $y$): RF and two other factors SMB and HML. Here we use yearly data ($n = 94$ from 1927-2020, shown as \textcolor{blue}{Blue dots}) to construct the prediction bands, each illustrated in Fig. (ref)-(ref). As for the test data, we use the standardized monthly data ($N = 1134$, normalized to zero mean and unit standard deviation, shown as \textcolor{red}{Red dots}) as a surrogate for test $(x, y)$ pairs. Namely, we match $12$ test data to each training data. We verified that after standardization, the histograms of yearly and monthly data match nicely for all four variables. For each type of response variable, we run two SDPs with different kernels. The summary statistics about the coverage probability, median, and mean bandwidth are given in Table (ref).

We note a few observations regarding the empirical results. First, all models achieve desirable coverage (all $>95\%$). Second, controlling for MKT, all other factors have significant heteroscedastic error left unexplained. For the RF, a high MKT return implies a low expected RF interest, and more importantly, a small variability, compared to the low MKT return case. For the size factor SMB, the conditional variability is much larger when the \textsf{ MKT} is high vs. low, so does the conditional expectation. The conditional variability in \textsf{ SMB} is roughly minimized when the market is significantly below average. While for the value factor \textsf{ HML}, conditional variability is minimized when the market is slightly below its average.

The second experiment aims to verify the mean and variance trade-offs by tuning the parameter $\gamma$, discussed in Sec. (ref). Here we use the monthly return data, and for each sub-experiment, we split the data into (train, valid, test) parts. We train models with different $\gamma = 0.1, 1, 10$ on the training data, then valid their performances on the validation data. We finally evaluate the performance using the test data with the cross-validated optimal $\gamma$ (based on the validation data). A nice feature about this experiment is that, one can visualize how the SDPs trade a complex/large conditional variance $\mathsf{v}(x)$ for a simple/small conditional mean $\mathsf{m}(x)$ in explaining $\by|\bx=x$ as $\gamma$ increases, illustrated by Fig. (ref).

Summary

The current paper progresses to resolve the uncertainty quantification dilemma faced by modern machine learning models. There are two innovative viewpoints we are taking. First, rather than relying on idealized parametric distributional assumptions on error $\by-\mathsf{m}(\bx)$, we make minimal assumptions. Both the conditional mean and variance functions are modeled nonparametrically and can universally approximate all functions. It is worth noting that such flexibility does not hinder computational feasibility due to the sum-of-squares and convex relaxations. The computational complexity and statistical guarantee scale favorably with high-dimensional covariates $x$. Second, rather than modeling the conditional mean only and giving up the variance (Frequentist justification, the conditional mean is assumed inside an RKHS, see caponnetto2007optimal, liang2020JustInterpolate, liang2020MultipleDescent), or modeling the conditional variance function only (Bayesian justification of kriging/Gaussian processes regression, the covariance function is specified by a kernel, see handcock1993bayesian,stein2005space, stein2012interpolation) for the variability in data, we model both the mean and the variance and prove strong, non-asymptotic Frequentist coverage guarantees. Such a modeling advantage enables the uncertainty quantification with or without any black-box predictive model, whether accurate or not.

To conclude, our Theorem (ref) established a strong, non-asymptotic coverage guarantee in the language of Neyman, yet with two distinct new features. First, the coverage probability can go to $1$ with a fixed confidence parameter $\delta$ as long as the sample size $n$ is large enough. Second, the data-adaptive quantity $\widehat{\mathsf{Opt}}_n$ controls both the average bandwidth and the coverage guarantee of the prediction band $\widehat{\mathsf{PI}}(x)$. A small objective value of the SDP makes the prediction band accurate and narrow simultaneously. Finally, our procedure for constructing prediction bands can be viewed as a novel variance interpolation with confidence and further leverages techniques from semi-definite programming and sum-of-squares optimization. We conducted simulated and real data experiments to validate the prediction interval's numerical performance for uncertainty quantification. A minimal 10-line Python implementation is provided in Listing. (ref) for interested readers.

Acknowledgement

Liang acknowledges the generous support from the NSF Career award (DMS-2042473) and the William S. Fishman fellowship. Liang thanks Michael Stein, Ruey Tsay, and Vladimir Vovk for constructive comments. Liang thanks the editor, associate editor, and three anonymous referees for the valuable suggestions that improve the paper.