Universal Prediction Band via Semi-Definite Programming

We propose a computationally efficient method to construct nonparametric, heteroscedastic prediction bands for uncertainty quantification, with or without any user-specified predictive model. Our approach provides an alternative to the now-standard conformal prediction for uncertainty quantification, with novel theoretical insights and computational advantages. The data-adaptive prediction band is universally applicable with minimal distributional assumptions, has strong non-asymptotic coverage properties, and is easy to implement using standard convex programs. Our approach can be viewed as a novel variance interpolation with confidence and further leverages techniques from semi-definite programming and sum-of-squares optimization. Theoretical and numerical performances for the proposed approach for uncertainty quantification are analyzed.


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 [Vovk et al., 2005, Lei et al., 2018, Romano et al., 2019. 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) ∈ X × R be the covariates and response data pair drawn from an unknown probability distribution P. There are plenty of regression or predictive models-denoted by m 0 (x)-that estimate m(x) := E[y|x = x] sufficiently well with finite data. However, to make downstream decisions reliable, a good prediction band quantifying the uncertainty in |y − m 0 (x)| with provable coverage is urgently needed. The prediction band is of particular relevance to complex machine learning models that construct 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 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(·, ·) : X × X → R be a continuous symmetric and positive-definite kernel function. Given n data pairs {(x i , y i )} n i=1 and the corresponding kernel matrix K ∈ S n×n with K ij = K(x i , x j ), our prediction band is constructed based on the following semi-definite program (SDP) min B Tr(KB) B 0 where the optimization variable B ∈ S n×n is a symmetric positive semi-definite (PSD) matrix, m 0 (·) : X → R is a given predictive model (user-specified), and K i ∈ R n denotes the i-th column of the kernel matrix K.
Given the estimated B, the prediction band, PI(x) that maps each x to an interval, can be constructed accordingly and K x := [K(x, x 1 ), . . . , K(x, x n )] ∈ R n , with δ ∈ R being a scalar quantifying confidence. δ can be later calibrated for exact coverage control, and may be set as 0 if n is large.
Here v(x) estimates the variability in the "deviations" e i := y i − m 0 (x i ). A few remarks on such deviations are in place.
• First, e i 's can be computed based on any user-specified predictive model m 0 (x) that estimates the conditional mean m 0 (x) ≈ E[y|x = x], be it accurate or not.
• Second, in the absence of such a predictive model for the conditional mean, one can set m 0 (x) ≡ 0 and learn a conditional second-moment function to assess uncertainty.
• Last, as shown next in (3), in practice, one can simultaneously learn the conditional mean and variance functions using a variant of the above SDP. Therefore, a pre-specified model m 0 (x) is not required.
Let K m and K v specify two kernel functions, corresponding to the conditional mean and variance functions respectively. K m , K v ∈ S n×n denote empirical kernel matrices on finite data with size n. For any γ ≥ 0, the following convex SDP program constructs the prediction band and the conditional mean function simultaneously B 0 where the optimization variables are B ∈ S n×n and α ∈ R n . Given the solution B and α, the PI(x) is constructed as

A Numerical Illustration
Before diving into the motivations behind the above SDPs (Sec. 1.3) and corresponding theory (Sec. 2), 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 3.1. Details of the data generating processes will be elaborated therein as well. The quick exercise is to showcase that convex programs (1) and (3) are easy to implement using standard optimization toolkits (say, CVX [Grant and Boyd, 2014]), and construct flexible prediction bands with desired coverage properties. As a motivating example, we try out the SDP (3), which simultaneously estimates the conditional mean and variance functions. A minimal 10-line Python implementation is provided in Listing. 1. Here the training and test data share the same conditional distribution y|x = x and thus m(x). The training and test data are shared in three plots. A good coverage corresponds to when Blue band covers essentially all Red dots. Statistics are summarized in Table 1.
The first example is a linear model with heteroscedastic error: the conditional mean m(x) being a linear function and variance 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 (3) with linear kernels K m and K v for both mean and variance functions, and (c) SDP2, a SDP (3) with a linear K m and a (degree-3) polynomial kernel K v . The coverage is compared on the same test dataset of size N = 800. Here γ = 0.1. See Fig. 1a for details and Table. 1 for coverage statistics.
The second example is a non-linear, heteroscedastic error model: mean m(x) and variance 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 m and linear kernel for K v , and (c) SDP2, rbf kernels for both K m and K v , summarized in Fig. 1b and Table. 1. Here γ = 1.
These two numerical examples are minimal yet informative. In Fig. 1a, 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. 1b. Observe first that in SDP2, the prediction band constructed by (3) 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 Blue lines and Red lines. The errors are also numerically summarized in the column "MSE" of Table 1. 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 (1) and (3) 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 α is not an optimization variable, (3) recovers (1). Second, the constraints set of (3) is always non-empty since α = 0, Sum-of-Squares and Phase Retrieval As shown in Prop. 3, the infinite-dimensional SDP with a nuclear-norm minimization is equivalent to (3), A 0 Here H m , H v denote two RKHSs where the conditional mean and variance functions reside. φ xi ∈ H is the feature map w.r.t. the Hilbert space H and ·, · H is the Hilbert space inner-product. We call it the infinite-dimensional SDP since the optimization variables (β, A) are (function, operator) rather than finitedimensional (vector, matrix). A few remarks are in place. First, if the kernel K m is universal, φ m x , β H m is dense in L 2 and hence can universally approximate all conditional mean function m(x). Second, as for the conditional variance which has positivity constraints over a continuum we relax the positivity constraints using a sum-of-squares form It turns out that when K v is universal, the above sum-of-squares function can approximate all smooth, positive functions [Fefferman and Phong, 1978, Bagnell and Farahmand, 2015, Marteau-Ferey et al., 2020, thus explaining the name "universal" in the title. Remark that sum-of-squares optimization [Lasserre, 2001] for nonparametric estimation has recently been considered; see [Bagnell and Farahmand, 2015, Marteau-Ferey et al., 2020, Curmei and Hall, 2020. The further relaxation changing from equality in (4) 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 v( In language, for all prediction bands that shelter the data, (3) aims to find the one with minimum bandwidth.
The curious reader may wonder where the nuclear norm A 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 m (x, x ) ≡ 0 and K v (x, x ) = x, x (the linear kernel with φ v x = x), and force the inequality constraints to be equal, our SDP in (3) is equivalent to phase retrievel Conceptually, the minimum nuclear norm procedure estimates the smallest number of factors that could generate the variance.
Min-norm Variance Interpolation with Confidence Now we discuss the tuning parameter γ ∈ [0, ∞] and reveal the connection to the recent min-norm interpolation literature , Bartlett et al., 2020, Ghorbani et al., 2020, Liang and Recht, 2021. In the limit of γ → 0, (3) reduces to the familiar min-norm interpolation with kernel K m (whenever it has full rank, since optimal B = 0) In the limit of γ → ∞, (3) reduces to (since optimal α = 0) Now it is clear what the role of the tuning parameter γ is: it trades off the conditional mean m(x) and variance v(x) to explain the variability in y's witnessed on the data. A small γ aims to use a complex mean m(x) and a parsimonious variance v(x) to explain the overall variability, and vice versa.
From the above discussion, it is also clear that the SDP (3) can be viewed as a min-norm variance interpolation with confidence. Instead of having the typical equality constraints in interpolation 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) As we shall see, the notion of confidence in this variance interpolation is closely related to the notion of margin in classification [Bartlett et al., 1998, Liang andSur, 2020].

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 v (x, x ) = 1(x = x ), then the decision variable B only matters in its diagonal component, and our SDP (3) reduces to the kernel ridge regression Moreover, a slight modification of (3) 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 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.

Conformal Prediction
Based on the exchangeability of data and a user-specified nonconformity measure, Vovk et al. [2005], Shafer and Vovk [2008] 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 Shafer and Vovk [2008] Section 4.3) usually requires to enumerate over all the possibilities of z = (x, y) ∈ X × Y, and for each possibility, calculate n nonconformity measures via the leave-one-out method. Therefore, the total budget is n × |Y| × |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 ∈ 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 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 1: the prediction band has coverage probability > 95% on a new data point (x, y), for 99.9999% dataset {(x i , y i )} n i=1 of size n drawn from the same distribution. Such a distinction on "confidence" vs. "probability" is discussed extensively in Section 2.2 of Shafer and Vovk [2008].
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. Lei et al. [2018] alleviates the computational burden of the conformal prediction by introducing the sample-splitting technique. Remarkably, theory on the bandwidth is also studied in Lei et al. [2018], thus providing an angle to probe the statistical efficiency. Romano et al. [2019] studies the problem that existing conformal methods can form nearly constant or weakly varying bandwidth and provide conservative results. Romano et al. [2019] 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.
Residual Subsampling and Quantile Regression An alternative approach for uncertainty quantification that leverages the metric structure in x ∈ X is to resample the residuals locally. Typically, this is done by first fitting a predictive model m 0 (x), and defining a local neighborhood around a new data x, 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 [Koenker and Bassett Jr, 1978, Koenker and Hallock, 2001, Belloni and Chernozhukov, 2011, Belloni et al., 2019, ξ τ (·) := arg min ξ 1 n n i=1 ρ τ y i − ξ(x i ) where τ ∈ (0, 1) is a quantile parameter, ρ τ (·) : R → R + is the tilted absolute value function, and ξ τ (·) : X → R is the estimated conditional quantile function. However, it is not guaranteed that over all x ∈ X , the estimated conditional quantile function satisfies ξ τ1 (x) < ξ τ2 (x) for two quantiles τ 1 < τ 2 . In other words, it is entirely possible that for several x's, the conditional prediction intervals are empty [Chernozhukov et al., 2010].

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 (x, y) ∼ P. To highlight the main arguments in a simple form, let us consider the setting m 0 (x) ≡ 0. Otherwise, the same proof follows by replacing y with y − m 0 (x). Define the corresponding prediction band, with a confidence parameter δ ∈ (0, 1] We need the following assumptions before stating the theorem, where (x, y) ∼ P and C > 0 denotes a universal constant.
[S1] (Kernel and RKHS) The continuous symmetric kernel K is positive definite and satisfies sup x∈X K(x, x) ≤ C. In addition, eigenvalues of the associated integral operator T : [S3] (Non-wild uncertainty) There exists a constant ω > 0 such that P y 2 > t · K(x, x) < exp(−Ct ω ) for all t ≥ 1.
Discussion of Assumptions All the above assumptions are mild. The eigenvalue decay in [S1] is almost identical to Tr(T ) < ∞ (bounded trace integral operator).
[S2] is also minimal, since it is only not true when there is no variability in y|x = x.
[S3] is the most stringent one, which requires the variability of y to exhibit a certain tail-decay. For small ω ∈ (0, 1), [S3] can be much milder than exponential tail-decay. Bounded or Gaussian y|x = x satisfies [S3] with arbitrarily large ω or ω = 2, respectively. With some extra work, [S3] can be relaxed to the case of a sufficiently rapid polynomial tail-decay.
[S2] can be relaxed to restricting only to x with P y 2 > 0 | x = x = 1.

What does the Theorem Entail
A few remarks are in order before we sketch the proof of Theorem 1.
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 (x, y) ∼ P drawn from the same distribution, with a non-asymptotic coverage probability (on the new data x, y) With δ = 0.5, the bandwidth Length PI(x) = 2.45 v(x) is at a heteroscedastic level adaptive to x with corresponding coverage probability at least 1 − O ( 1 n ). Here O hides polylog factors. The coverage can be arbitrary close to 1 with large n without the need of increasing δ, which is in clear distinction to the conventional wisdom that coverage 1 can only be possible with an increasing δ regardless of n. Again, we emphasize that the above coverage guarantee holds essentially on 99.9999% 1 − n −10 of the datasets In Section 2.2, we propose a rigorous calibration algorithm to choose δ (α) to achieve a constant coverage level 1 − α ∈ (0, 1).
Optimality If one wishes to obtain the classic 95% coverage probability, then choosing δ = O ( 1 n ) suffices, which translates to Recall that in classic simple linear regression, the prediction interval is of length withŝ = Σiê 2 i n−2 being the estimated residual standard error. The fact that (6) and (7) share the 1 n fluctuation seems to indicate the optimality of our Theorem (in terms of the dependence on n).
Data Adaptivity Curiously, the objective value of the convex optimization program quantifies the uncertainty of the prediction band: a smaller Opt n implies (a) a better confidence/coverage guarantee and (b) a narrower prediction band overall. More importantly, the 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 Opt n is a function of the dataset, our Theorem reveals which dataset allows for a better prediction band. Remark that Opt n = v(·) 2 is also a particular norm of the heteroscedastic variance function, quantified by the nuclear norm of the associated PSD operator A 0 with v(x) = φ x , Aφ x H . Curiously, a simpler variance function 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 1 is in an inequality form with a mild dependence on the non-explicit universal constant C γ,ξ,η,ω ; however, the coverage is in a stronger conditional sense conditioned on {(x i , y i )} n i=1 . 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 δ. This section provides a theoretically justified calibration procedure in choosing δ (in the SDP band) to control coverage at a user-specified level 1 − α ∈ (0, 1), as in conformal prediction. Such fine calibration on δ can be helpful numerically in the moderate n and constant α 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 )} n i=1 and the calibration set {x j , y j } m j=1 , with in total n + m data points drawn i.i.d. from P. The training set will be used to construct prediction band PI(·, δ). The calibration data set will be used to choose δ to calibrate coverage.
Calibration Procedure Suppose that there exsits large enough constants ∆, L > 0 such that The calibration uses a dyadic search to select δ (α) ∈ [−1, ∆] with the set {x j , y j } m j=1 . The goal of the calibration is to ensure P(x,y)∼P y ∈ PI(x, δ (α)) ≥ 1 − α.

Algorithm 1: Calibration for Coverage Control
The following result derives the theoretical ground for the calibration procedure. We defer the proof to Section A.2.
Lemma 2 (Calibration). Consider the calibration procedure in Algorithm 1 and L, ∆, α specified therein. If the size of the calibration set m is large enough such that then the calibrated δ (α) satisfies the coverage control with probability at least 1 − 2m −10 on the calibration set {x j , y j } m j=1 .

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 δ (α) for a fixed coverage level 1−α ∈ (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 1. Observe that by definition There exists a k ∈ N such that f A ∈ F k with 2 k−1 ≤ Opt n < 2 k , and thus we continue to bound .
For term (ii), we will use the high probability symmetrization in Prop. 4. Introduce i.i.d. Rademacher variables { i } n 1 independent of the data. Note that we only need to consider k ≤ k 0 such that 2 k0 = [log(n)] cω , where we use the upper estimate on Opt n obtained in Prop. 5, which is implied by Assumption [S3]. With probability at leat 1 − 2 exp(−t) on the data {z i } n 1 , uniformly over all k ≤ 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. 4, formally upper bounded in (15)), first recall 2 k ≤ 2( Opt n ∨ 1), we know Similarly, by Prop. 4, the deviation term (iii) can be bounded by 6( Opt n ∨ 1) · sup x,y φx y 2 H · k0+t n with k 0 = c ω 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 E φx y ⊗ φx y = 0, and φx y ⊗ φx y op ≤ sup x,y φx y 2 H a.s. , Naively applying the matrix Bernstein inequality, one would expect the term (iv) to behave like sup x,y φx y 2 H · t n ∨ t n with probability 1−dim(φ x )·exp(−t) on { i } n 1 . This is educative yet wrong, since dim(φ) is infinity. To make things rigorous, we rely on a truncation technique to look at a finite-dimensional version φ ≤m x truncated at a level m = poly(n) to apply matrix Bernstein, and then estimate the remaining contribution from φ >m x by the eigenvalue decay in Assumption [S1]. With details given in Prop. 6, we derive that The final piece of the puzzle lies in the term sup (x,y)∈dom(P) φx y 2 H , 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 y 2 |x=x (·) : [0, 1] → R + for the conditional random variable y 2 |x = x. Let's only look at data (x i , y i )'s lying in the region Ω := {(x, y) | y 2 > Q y 2 |x=x (1 − η)}, and denote P| Ω as the conditional distribution of data (x, y) conditioning on the region Ω. Claim that for any PI(x, δ) This is based on two facts. First, regardless of the ordering of Q y 2 |x=x (1 − η) and (1 + δ) v(x). Second, conditioning on Ω does not change the marginal distribution of x due to the quantile construction, namely P x | Ω ≡ P x . Marginalizing (12) over x proves the above claim.

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 [Lei et al., 2018, Algorithms 1 and 2 respectively]. 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 x ∼ Unif([− √ 3, √ 3]), and the conditional distribution y|x = x ∼ · v(x) where the independent error is either drawn from a standard normal N(0, 1) (Fig. 2a) or a uniform distribution Unif([− √ 3, √ 3]) (Fig. 2b). The conditional mean function m(x) = 0. The heteroscedastic variance function scales as v(x) = 1 + x + 4x 2 , depending on x. For each simulated data set, we generate {(x i , y i )} n i=1 i.i.d. from the above data generating process (DGP), with n = 50 as the training data, marked by Blue dots. The test data set are drawn from the same DGP, with N = 500 marked by 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 δ in SDPs as in Algorithm 1, and the homoscedestic conformal bandwidth as in SplitCF. We compare five methods that construct the prediction bands, illustrated by the Blue band, on the Gaussian error data set in Fig. 2a, and Uniform error data set in Fig. 2b. For all methods, the desired coverage is set at 1 − α = 95%. For the SDPs in (3) γ = 10. Table 2 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. 2a-2b and numerically in Table 2, the two conformal methods, SplitCF and CF, estimates a conservative, wide prediction band that is almost homoscedastic. In contrast,

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 Fama and French [1993]. 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 γ in trading off mean and variance.
The first experiment aims to access the prediction coverage in the SDP (3), 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 Blue dots) to construct the prediction bands, each illustrated in Fig. 3a-3c. As for the test data, we use the standardized monthly data (N = 1134, normalized to zero mean and unit standard deviation, shown as 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 3. 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 MKT is high vs. low, so does the conditional expectation. The conditional variability in SMB is roughly minimized when the market is significantly below average. While for the value factor 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 γ, discussed in Sec. 1.3. 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 γ = 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 γ (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 v(x) for a simple/small conditional mean m(x) in explaining y|x = x as γ increases, illustrated by Fig. 4.

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 y−m(x), 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 highdimensional 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 [Caponnetto and De Vito, 2007), or modeling the conditional variance function only (Bayesian justification of kriging/Gaussian processes regression, the covariance function is specified by a kernel, see [Handcock and Stein, 1993, Stein, 2005, 2012) 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 1 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 δ as long as the sample size n is large enough. Second, the data-adaptive quantity Opt n controls both the average bandwidth and the coverage guarantee of the prediction band 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. 1 for interested readers. On the left is SMB as the response and on the right is HML. For each response variable, we run two sub-experiments: the top row corresponds to a linear m(x) and a quadratic v(x), and the bottom row corresponds to a degree-3 polynomial m(x) and a quadratic v(x). Each sub-figure corresponds to a specific γ, noted in its title. The cv γ denotes the cross-validated optimal γ using the validation dataset. Here the (train, valid, test) dataset has size proportional to 1:3:9, denoted by Blue dots, Red pluses, and Red dots respectively.

A.1 Remaining Propositions
In this section, we collect the remaining propositions.
Proposition 3 (Representation). The kernelized version of the SDP as in (3) is equivalent to the following infinite-dimensional SDP Proof. Noticing that the solution to the infinite-dimensional problem must lie in the span of data, namely A = i,j B ij φ v xi ⊗ φ v xj with some PSD B ∈ S n×n , and β = i α i φ m xi with some α ∈ R n . With the above representation, plug in the infinite-dimensional SDP and recall Tr(A) = A , we can derive (3). When β is not a decision variable, this representation theorem applies to (1).
Proposition 4 (Symmetrization). Let F be a class of functions f : Z → R, with sup x∈Z |f (z)| ≤ M . Then with probability at least Proof. First, with McDiarmid's inequality, we know w.h.p.
Apply Giné-Zinn symmetrization to the first term on the RHS, then apply McDiarmid's inequality again, we can establish the claim. See  for details.
Proposition 5 (Objective value estimate). Under [S3], the following holds with probability at least 1 − n −10 , Proof. Apply union bound on the tails given by [S3], with the choice t 0 = [log(n)] cω we know In view of Prop. 3, the above certifies that A := t 0 ·I lies in the feasibility set φ xi , Aφ xi H = t 0 ·K(x i , x i ) ≥ y 2 i , which implies t 0 being an upper bound on Opt n .
Proposition 6 (Operator-norm estimate). Under [S1], for any {x i } n i=1 , the following holds Proof. Recall [S1], due to the Mercer's theorem, one can represent φ x as an infinite-dimensional vector, with each coordinate of φ j x corresponding to the eigenfunction of the integral operator, with j = 1, . . . , ∞ and λ j ≤ Cj −τ . To bound the operator norm, recall the Rayleigh quotient form, for any h ∈ H with h 2 H = 1 h, 1 For the second term, recall the eigenvalue decay λ j ≤ Cj −τ with τ > 1, we know it is upper bounded by Cm −(τ −1) . Choosing log(m) = C τ log(n) with a constant large enough, we know the second term is dominated by the first term. By integrating the tail bound to obtain a bound on the expected value, we complete the proof.

A.3 Remaining Experimental Details
All experiments are conducted using the Python language. The minimal implementation is provided below