extfieldsofstudy:
Mathematics
Computer Science
Economics

text:
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


Look Closer to Segment Better: Boundary Patch Refinement for Instance Segmentation

Tremendous efforts have been made on instance segmentation but the mask quality is still not satisfactory. The boundaries of predicted instance masks are usually imprecise due to the low spatial resolution of feature maps and the imbalance problem caused by the extremely low proportion of boundary pixels. To address these issues, we propose a conceptually simple yet effective post-processing refinement framework to improve the boundary quality based on the results of any instance segmentation model, termed BPR. Following the idea of looking closer to segment boundaries better, we extract and refine a series of small boundary patches along the predicted instance boundaries. The refinement is accomplished by a boundary patch refinement network at higher resolution. The proposed BPR framework yields significant improvements over the Mask R-CNN baseline on Cityscapes benchmark, especially on the boundary-aware metrics. Moreover, by applying the BPR framework to the PolyTransform + SegFix baseline, we reached 1st place on the Cityscapes leaderboard.

Introduction
Instance segmentation, which aims to assign a pixelwise instance mask with a category label to each object in an image, has great potential in various computer vision applications, such as autonomous driving and robotics. Mask R-CNN [15] is a prevailing two-stage instance segmentation framework, which first employs a Faster R-CNN [34] detector to detect objects in an image and further per-* Equal contribution. † Corresponding author. forms binary segmentation within each detected bounding box. Other methods [16,27] built upon Mask R-CNN consistently achieve superior performance. Driven by the recent development of one-stage detectors [24,39,55], a number of one-stage instance segmentation frameworks [3,4,8,21,38,42,43,44,48,53] have been proposed. However, the quality of the predicted instance mask is still not satisfactory. One of the most important problems is the imprecise segmentation around instance boundaries. As shown in Figure 1(left), the predicted instance masks of Mask R-CNN are coarse and not well-aligned with the real object boundaries. Empirically, correcting the error pixels near object boundaries can improve the mask quality a lot. We conducted an upper bound analysis in Table 1. A large gain (9.4/14.2/17.8 in AP) can be obtained by simply replacing the predictions with ground-truth labels for pixels within a certain Euclidean distance (1px/2px/3px) to the predicted boundaries, especially for small objects.

Mask R-CNN Ours
We argue that there are two critical issues leading to lowquality boundary segmentation. (1) The low spatial resolution of the output, e.g. 28×28 in Mask R-CNN or at most 1/4 input resolution in some one-stage frameworks [38,43], makes finer details around object boundaries disappear. The predicted boundaries are always coarse and imprecise (see Figure 1 and 4). (2) Pixels around object boundaries only make up a small fraction of the whole image (less than 1% [18]), and are inherently hard to classify. Treating all pixels equally may leads to an optimization bias towards smooth interior areas, while underestimating the boundary pixels.
As a long-standing challenge in dense prediction tasks, many studies have attempted to improve the boundary quality, while the above issues are still not well solved. For example, BMask R-CNN [9] and Gated-SCNN [37] employ an extra branch to enhance the boundary awareness of mask features, which can fix the optimization bias to some extent, while the low resolution issue remains unsolved. Poly-Transform [23] and SegFix [50] act as a post-processing scheme to improve the boundary quality. PolyTransform [23] employs a deforming network with the cropped instance patch to predict the offsets of polygon vertices, while suffering from a large computational overhead. SegFix [50] replaces the coarse predictions of boundary pixels with interior predictions, but it relies on precise boundary predictions. We argue that the instance boundary prediction task shares a similar complexity with instance segmentation.
Considering the human annotation behavior for instance segmentation, the annotators usually first localize and categorize each object in the given image, and then explicitly or implicitly segment some coarse instance masks at a low resolution. Afterwards, to obtain a high-quality mask, the annotators need to repeatedly zoom into the local boundary regions and explore the sharper boundary segmentation at higher resolution. Intuitively, high-level semantics are required to localize and roughly segment objects, while lowlevel details (e.g. colour consistency and contrast) are more critical for segmenting the local boundary regions.
In this paper, motivated by the human segmentation behavior, we propose a conceptually simple yet effective post-processing framework to improve the boundary quality through a crop-then-refine strategy. Specifically, given a coarse instance mask produced by any instance segmentation model, we first extract a series of small image patches along the predicted instance boundaries. After concatenated with mask patches, the boundary patches are fed into a refinement network, which performs binary segmentation to refine the coarse boundaries. The refined mask patches are then reassembled into a compact and high-quality instance mask, shown in Figure 1(right). We termed the proposed framework as BPR (Boundary Patch Refinement). The proposed framework can alleviate the aforementioned issues and improve the mask quality without any modifi- cation or fine-tuning to the segmentation models. Since we only crop around object boundaries, the patches are allowed to be processed with the much higher resolution than previous methods, so that low-level details can be retained better. Concurrently, the fraction of boundary pixels in the small patch is naturally increased, which can alleviate the optimization bias. The proposed BPR framework significantly improves the results of Mask R-CNN baseline (+4.3% AP on Cityscapes dataset), and produces substantially better masks with finer boundaries. We found that the model trained on the results of Mask R-CNN can be easily transferred to refine the results of other instance segmentation models as well, without the need for re-training. We outperform some boundary refinement methods [19,50] and show that these methods are complementary by successfully transferring our model to improve their results. Furthermore, by applying our BPR framework to the "Poly-Transform + SegFix" baseline [50], we established a new state-of-the-art on the Cityscapes test set with AP of 42.7%, and ranked 1 st place on the Cityscapes leaderboard by the CVPR 2021 submission deadline.

Related Work
Instance Segmentation. Recent studies on instance segmentation can be divided into two categories: two-stage and one-stage methods, as briefly reviewed below.
Two-stage methods usually follow the classical detectthen-segment strategy. The dominant method is still Mask R-CNN [15], which inherits from the two-stage detector Faster R-CNN [34] to first detect objects in an image and further performs binary segmentation within each detected bounding box. Following Mask R-CNN, PANet [27] enhances feature representation through bottom-up path augmentation. Mask Scoring R-CNN [16] adds an additional mask-IoU head to re-score the mask predictions. These methods consistently achieve superior performance.
One-stage methods recently attracts more attention due to the rapid development of one-stage detectors [24,39,55]. Some methods [3,4,21,48,53] continue to adapt the detect-then-segment strategy but replace the detectors with the one-stage alternatives. YOLACT [3] achieves real-time speed by learning a set of prototypes and the prototypes are assembled with the learned linear coefficients. BlendMask [4] further improves this idea by assembling with attention maps. Some recent proposed methods [8,38,42,43] eliminate the need for detection by directly segmenting objects in a location-wise manner. CondInst [38] and SOLOv2 [43] achieve remarkable performance with high efficiency. In addition, there are some approaches [11,13,32] built upon the semantic segmentation models, which usually learn the pixel-wise embeddings and then cluster them into instances. Several works [2,33,44,47] replace the pixel-wise instance representation into the contour-based representation.
Our proposed framework is agnostic to the instance segmentation methods, thus it can be applied to refine the results of any instance segmentation model, both one-stage and two-stage methods.
Semantic Segmentation. Modern semantic segmentation approaches are pioneered by fully convolutional networks (FCNs) [29]. Many studies have been proposed on this foundation to improve the segmentation results, such as increasing the resolution of feature maps with dilated/atrous convolutions [6,7], enriching context information [12,49,52,54], using an encoder-decoder architecture [7,17,30,35], or some refinement schemes [20,22,50]. Minaee et al. [31] provided a comprehensive review of these approaches. In this paper, we adopt the prevailing HRNet [41] in our framework, which can maintain high-resolution representation throughout the whole network.
Boundary Refinement for Segmentation. Most recent studies focused on boundary refinement aim at designing a boundary-aware segmentation model by integrating an extra and specialized module to process boundaries. For example, BMask R-CNN [9] and Gated-SCNN [37] employ an extra branch to enhance the boundary awareness of mask features by estimating boundaries directly. PointRend [19] iteratively samples the feature points with unreliable predictions and refines them with a shared MLP. Another line of work attempts to refine the boundaries based on the results of existing segmentation models with a post-processing scheme. SegFix [50] is a general refinement mechanism, which replaces the unreliable predictions of boundary pixels with the predictions of interior pixels. The effectiveness of SegFix highly depends on the accuracy of boundary prediction. However, it is very challenging to directly estimate precise instance boundaries. Intuitively, the instance segmentation task could easily be settled if the precise boundaries are already given. Our method shares more similarities with PolyTransform [23], which transforms the contour of instance into a set of polygon vertices. A Transformer [40] based network is applied to predict the offsets of vertices towards object boundaries. It achieves superior performance while suffering from a large computational overhead due to the large instance patch and the heavy Transformer architecture. Our proposed method is also a post-processing scheme. Different from these methods, we focus on refining the boundary patches to improve the mask quality.

Framework
An overview of the proposed framework is illustrated in Figure 2. As a post-processing mechanism, the proposed framework can be applied to refine the results of any prevailing instance segmentation model, without any modification or fine-tuning to the segmentation models themselves.

Boundary Patch Extraction
Given an instance mask produced by an instance segmentation model, we first need to determine which part of the mask should be refined. Based on the findings of previous works [9,50] and our verification experiments in Table 1, we propose an effective sliding-window style algorithm to extract a series of patches along the predicted instance boundaries. Specifically, we densely assign a group of squared bounding boxes where the central areas of the box should cover the boundary pixels, shown in Figure 2(b). The obtained boxes still contain large overlaps and redundancies, thus we further apply a Non-Maximum Suppression (NMS) algorithm to filter out a subset of patches (Figure 2c). Empirically, with the larger overlaps, the segmentation performance can be boosted, while simultaneously suffering from the larger computational cost. We can adjust the NMS threshold to control the amount of overlap to achieve a better speed/accuracy trade-off. In addition to image patches, we also extract the corresponding binary mask patches from the given instance mask. The concatenated image and mask patches (Figures 2d and 2e) are resized and fed into the following boundary patch refinement network.

Boundary Patch Refinement
Mask Patch. The benefit of the binary mask patch is that it accelerates training convergence and provides location guidance for the instance to be segmented. As discussed in the previous works on semantic segmentation [41,49], context information plays a vital role for pixel-wise classification. Therefore, the cropped image patches are hard to be classified independently due to the limited context information. With the help of location and semantic information provided by the mask patches, the refinement network can eliminate the need for learning instance-level semantics from scratch. Instead, the refinement network only needs to learn how to locate the hard pixels around the decision boundary and push them to the correct side. We believe this goal can be achieved by exploring low-level image properties (e.g. colour consistency and contrast) provided in the local and high-resolution image patches. More importantly, the adjacent instances are likely to share an identical boundary patch, while the learning goals are totally different and ambiguous. Together with different mask patches for each instance, these issues can be avoided. As compared in Table 2, the model has trouble to converge without using the mask patches, examples of which are shown in Figure 3.
Boundary Patch Refinement Network. The role of this refinement network is to perform binary segmentation for each extracted boundary patch individually. Any semantic segmentation model can be employed for this task by simply modifying the input channels to 4 (3 for the RGB image patch and 1 for the binary mask patch) and output classes to 2. For the sake of convenience, we adopt the state-of-theart HRNetV2 [41] as the refinement network in our implementation, which can maintain high-resolution representation throughout the whole network. By increasing the input size appropriately, the boundary patches can be processed with much higher resolution than in previous methods.
Reassembling. The refined boundary patches are reassembled into a compact instance-level mask by replacing their previous predictions. Predictions are unchanged for those pixels without refinement. For the overlapping areas of adjacent patches, the results are aggregated by simply averaging the output logits and applying a threshold of 0.5 to distinguish the foreground and background.

Learning and Inference
The refinement network is trained based on the boundary patches extracted from training images and tested on validation or testing images. We do not directly train or fine-tune the instance segmentation models. During training, we only extract boundary patches from instances whose predicted masks have an Intersection over Union (IoU) overlap larger than 0.5 with the ground-truth masks, while all predicted instances are retained during inference. The model outputs are supervised with the corresponding ground-truth mask patches using the pixel-wise binary cross-entropy loss. We simply fix the NMS eliminating threshold to 0.25 during training, while adopting different thresholds during inference based on the speed requirements. See Appendix A for more implementation details.

Datasets and Metrics
Datasets. We mainly report the results on Cityscapes [10], a real-world dataset with high-quality instance segmentation annotations. We only used the fine data, containing 2, 975/500/1, 525 images for train/val/test, which are collected from 27 cities, with a high resolution of 1024×2048. Eight instance categories are involved, including bicycle, bus, person, train, truck, motorcycle, car, and rider.
Metrics. The COCO-style [25] mask AP (averaged over 10 IoU thresholds ranging from 0.5 to 0.95 in the step of 0.05), AP 50 (AP at an IoU of 0.5) and AP S /AP M /AP L (for small/medium/large instances) were reported in most of our experiments. The official Cityscapes-style AP [10] was only used to report the final results for a fair comparison, which is slightly higher than the COCO-style AP. Similar to [23,37,50], we also used a boundary F-score to evaluate the quality of the predicted boundaries. A mask was considered correct if the boundary is within a certain distance threshold from the ground-truth. We used a threshold of 1px and only compute for true positives that are determined on the same 10 IoU thresholds ranging from 0.5 to 0.95. The boundary F-score was computed in an instance-wise manner and then averaged over them, termed AF.

Ablation Study
We investigated the effectiveness of the proposed framework through extensive ablation experiments on the configurable design choices. We started the refinement with the results of Mask R-CNN ResNet-FPN-50 baseline trained on Cityscapes fine data (with COCO pre-training). We adopted the lightweight HRNetV2-W18-Small as the refinement network in the default setting, with input size equal to 128×128. The boundary patches were extracted with patch size equal to 64×64 without padding, and the inference NMS threshold was set to 0.25 by default.
Effects of Mask Patch. To validate the effect of mask patch for boundary refinement, we made a comparison by eliminating the mask patches while keeping other settings unchanged. As indicated in Table 2, the model trained with image patches solely yielded a terrible result, even worse than the segmentation results before refinement. However, together with mask patches, we achieved a significant improvement (+3.4% in AP, +11.9% in AF) by refining the Mask R-CNN segmentation results. We also show some patch-wise examples in Figure 3. For a simple case with one dominant instance in the image patch (first row), both of the models (w/ or w/o mask patch) produced reasonable results. However, as for cases with multiple instances crowded in the image patch, the model without mask patch (last column) failed to distinguish which object should be segmented, leading to coarse (4th row) or completely wrong (2nd and 3rd rows) predictions. In contrast, with the help of mask patches, we produced high-quality predictions with accurate and distinct boundaries (3rd column).
Patch Size. We increased the boundary patch size by cropping with a larger box and/or with padding. Note that the padded areas were only used to enrich the context and not used for reassembling. As the patch size gets larger, the w/ mask AP AP 50 AF AP S AP M AP L -36.4 60. 8      instance boundaries to achieve better segmentation results.
There are multiple choices about how to extract the boundary patches for refinement. We compared three extraction schemes and listed the results in Table 4. The most straightforward scheme is dividing the input image into a group of patches according to the pre-defined grids, and then picking only the patches that covering the predicted boundaries. We varied the patch size and found the results were consistently worse than our proposed "dense sampling + NMS filtering" scheme. One of the most important reasons is the imbalanced foreground/background ratio. We observed that some extracted patches are almost entirely filled with either foreground or background pixels. These patches are hard to refine due to the lack of context, thus leading to sub-optimal results. In contrast, by restricting the center of patches to cover the boundary pixels, the imbalance problem can be alleviated. Another scheme, similar to some previous works [23,28], is cropping the whole instance based on the detected bounding box and further re-segmenting the instance patch. As shown, even though the input size was increased to 512×512, the results are still sub-optimal, which demonstrated the effectiveness of our local boundary patches. See Appendix B for detailed descriptions. Input Size of the Refinement Network. The extracted patches are resized into a larger scale before refinement. Table 5 shows the impact of input size. We also report the approximate inference speed of the refinement network, with a fixed batch size of 135 (on average 135 patches per image). As the input size increases, the AP and AF scores increase accordingly, and slightly drop after 256. The boundaries can be processed with the higher resolution with the larger input size, thus more details can be retained.  Alternatives of refinement network. We compared different backbones for our refinement network in Table 6. As shown, a stronger backbone usually lead to higher performance, but at the expense of lower speed. Since the model essentially performs binary segmentation for patches, it can further benefit from the advances in semantic segmentation, such as increasing the resolution of feature maps [6,7,41] and more effective backbones [45,51].
NMS Eliminating Threshold. We studied the impact of different NMS eliminating thresholds during inference, shown in Table 7. As the threshold gets larger, the number of boundary patches increases rapidly. The overlap of adjacent patches provides a chance to correct unreliable predictions of the inferior patches. As shown, the resulting boundary quality was consistently improved with a larger threshold, and reached saturation around 0.55. We believe a better speed/accuracy trade-off can be achieved by setting a proper threshold.

Transferability
What the BPR model learned is a general ability to correct error pixels around instance boundaries. We can easily transfer this ability of boundary refinement to refine the results of any instance segmentation model. Specifically, once we get a model trained on the boundary patches extracted from the train-set predictions of Mask R-CNN on Cityscapes, we can make inference to refine any predictions (on Cityscapes train/val/test sets) produced by any models (not only Mask R-CNN), without the need for training from scratch. After training, the BPR model becomes modelagnostic, similar to SegFix [50]. We validated the trans-    Table 10: Results on COCO. AP is measured on the higher-quality LVIS [14] annotations. We improved based on the results of Mask R-CNN ResNeXt-FPN-101 baseline.
ferability by applying the model trained on Mask R-CNN results to refine the predictions of PointRend [19] and Seg-Fix [50]. Note that these two methods are also designed to improve boundary quality in segmentation. As shown in Table 9, the transferred model still improved the results of PointRend and SegFix by a large margin, suggesting that our method is compatible with them.

Overall Results
Comparison with State-of-the-art Methods. We integrated the optimal design choices and hyperparameters found in above ablation experiments into a stronger BPR model. Specifically, we adopted the HRNetV2-W48 as our refinement network, with 256×256 input patches resized from 64×64, and a NMS threshold of 0.55 during inference. We evaluated the framework on Cityscapes val and test sets and compared the performance against some state-ofthe-art methods in Table 8. (1) Compared with the Mask R-CNN baseline, we achieved a significant improvement (+4.3% AP in both val and test). We outperformed Seg-Fix [50] by a large margin, which is also a boundary refinement module applied to the same baseline with ours. Furthermore, by applying our BPR model to the results already refined by SegFix, we can still improved a lot (slightly lower than applying BPR only). (2) We transferred the above BPR model to refine the results of the stronger PolyTransform [23] baseline (1 st place at CVPR 2020). Our "PolyTransform + BPR" consistently improved 2.3% AP on Cityscapes test set and also outperformed "PolyTransform + Seg-Fix" (2 nd place at ECCV 2020) by a large margin (+1.2%). By applying BPR to "PolyTransform + SegFix", we established a new state-of-the-art on Cityscapes test with AP of 42.7%, reaching 1 st place on the Cityscapes leaderboard by the CVPR 2021 submission deadline.
Qualitative Results. We show some qualitative results on Cityscapes val in Figure 4. Compared with the coarse predictions of Mask R-CNN, our BPR framework generated substantially better instance segmentation results with precise and distinct boundaries. It largely alleviated the oversmoothing issues [19] in previous methods caused by the low resolution feature maps. More results are included in Appendix E. In addition, we also provided a detailed limitation analysis in Appendix F.
Speed. Only the speed of refinement network was con-  Table 5 and 6, excluding the patch extraction and reassembling time. As a whole pipeline, it takes about 211ms to process a single Cityscapes image (1024×2048) on a single RTX 2080Ti GPU under the default setting of ablation experiments, which is still much faster than Poly-Transform [23]. The detailed speed calculation and more speed analysis are included in Appendix C.
Results on COCO Dataset. To demonstrate the generality of our framework, we also report the results on the more challenging COCO dataset [25], which contains 80 categories and more images (118k/5k for train/val). It is important to note that the coarse annotations in COCO may not fully reflect the improvements in mask quality [14]. Following PointRend [19], we further report the AP measured using the higher quality LVIS [14] annotations. We randomly sampled about 8% of instances for fast training. As shown in Table 10, we improved the powerful Mask R-CNN ResNeXt-FPN-101 baseline by 0.8% AP and 1.7% AP on val2017. The coarse annotations on COCO train2017 may provide ambiguous optimization objectives, especially for our local boundary patches. It may mislead the learning of our BPR model, leading to suboptimal results. This issue was also observed in some contour-based instance segmen-tation methods [33,44,47]. We believe that training with more instances on higher quality annotations (e.g. LVIS) can further improve the results. More analysis on COCO dataset is included in Appendix D.

Conclusion
In this paper, we propose a conceptually simple yet effective boundary refinement framework to improve the boundary quality for any instance segmentation model. Starting from a coarse instance mask, we extract and refine a series of boundary patches along the predicted instance boundaries through an effective refinement network. The proposed framework achieved consistent and impressive improvements based on different baselines. Qualitative results show that our approach produced high-quality masks with precise and distinct boundaries.
In Section 4.2, we compared the proposed "dense sampling + NMS" scheme with another two patch extraction schemes: pre-defined grid and instance-level patch. Here we provide the implementation details and further analysis of these two schemes. As illustrated in Figure S1b, the predefined grid scheme simply divides the input image into a group of patch candidates according to a pre-defined grid. Candidates that covering both foreground and background pixels are choosen as boundary patches for refinement. This straightforward scheme yields plenty of inferior patches, as indicated by yellow dashed boxes in Figure S1b, which have the imbalanced foreground/background ratio and may lack of real boundary cues, thus leading to sub-optimal results. Another scheme is extracting the instance-level patch (Figure S1c) based on the detected bounding box, which is similar to previous studies [23,28]. This scheme can be viewed as an improved Mask R-CNN equipped with a stand-alone mask head, while still fails to solve the optimization bias issue and the learning process is dominated by interior pixels. Different from these methods, by adaptively extracting patches along the predicted boundaries in a sliding-window style ( Figure S1a) and refining the local boundary regions separately, the above issues can be alleviated.

C. More Speed Analysis
The inference time of our proposed framework is independent of the original instance segmentation models, which consists of three parts: patch extraction, refinement, and reassembling. Note that only the refinement part was considered when calculating the FPS in Table 5 and 6. Besides, the FPS was measured in an imprecise manner by fixing the batch size to 135 (average number of patches per image), while the exact number of patches varies from image to image. Here we report the total inference time, which measured by calculating the exact inference time for each image individually and then averaging them. Taking the default setting (HRNet-W18s with input size of 128×128) in our ablation experiments as example, it takes about 211ms (52ms,81ms,78ms for the above three parts respectively) to process an image (1024×2048) of Cityscapes on a single RTX 2080Ti GPU, which is still much faster than PolyTransform (575ms * per image [23]). Undoubtedly, the network speed can be further improved with more efficient backbones (e.g. MobileNets), smaller input size (e.g. 32×32 or 64×64), and less inference patches (e.g. with lower NMS thresholds or adaptively selecting the most * Measured on a single GTX 1080Ti GPU, which is about 35% slower than our RTX 2080Ti GPU with FP32 training (ref. lambdalabs.com). unreliable patches). Note that the BPR models can still achieve a remarkable performance under these lightweight settings (Tables 5,6,7). The patch extraction and reassembling steps can also be accelerated with more CPU cores.

D. More Analysis on COCO Dataset
In theory, the proposed framework, as a general boundary refinement mechanism, can be applied to any instance segmentation dataset. We achieved impressive performance on Cityscapes, while the AP improvement on COCO dataset was not as high as we got on Cityscapes (see Table 10). The most critical problem is that the coarse polygon-based annotations on COCO dataset yield significantly lower boundary quality [14]. Several examples (which are ubiquitous on COCO) are shown in Figure S2. The misalignment between annotations and real instance boundaries may greatly increase the optimization difficulty of our refinement model. Especially, the coarse annotations may provide ambiguous optimization objectives for our local boundary patches, thus hampering the model convergence. We observed that some contour-based instance segmentation methods [33,44,47], which are sensitive to the quality of boundary annotations, also suffered from this misalignment issue. It seems that the coarse COCO annotations may not friendly to these methods and it is hard to achieve very high AP scores based on these approaches. In spite of this, we still significantly improved the Mask R-CNN results in some cases, shown in Figure S3. Some results are even better than their annotations (the first three examples in Figures S2, S3).

E. More Qualitative Results
We provid more qualitative results on Cityscapes val, including image-level ( Figure S4) and patch-level (Figure S5) results. As shown, our proposed framework consistently improves the instance segmentation results of Mask R-CNN and produces substantially better instance masks with more precise boundaries.

F. Limitation Analysis
The performance of our proposed framework relies on the boundary quality of initial masks. Some failure cases are illustrated in Figure S6. For example, our model failed to produce an optimal mask if the initially predicted boundaries are far from the real object boundaries (1st row), but note that we still refined this case to some extent (IoU was improved). In addition, if the initial mask largely oversegments the neighboring instance, our model may regard the two instances as a whole and further enlarge this error (2nd and 3rd rows) since we only process the local boundary regions without a global view. We analyzed the IoU improvements for all predicted instances on Cityscapes val set, shown in Figure S7. In most cases, our refinement model can effectively improve the mask IoU (red dots above the dash line). However, we found that it's hard to refine instance masks with extremely lower IoU (e.g. < 0.1) due to the poor quality of initial boundaries. In addition, we observed that the improvement for smaller instances (about 2% in AP S ) is not as high as we got for larger instances (about 5% in AP L ). Compared to the upper-bound results (see Table 1), there is still a large step to take for boundary refinement, especially for small instances.

G. More Transferring Results
In Table 9, we verified that the BPR model trained on Mask R-CNN results can be effectively transferred to refine the results of PointRend and SegFix. As an opposite directions with Table 9, we instead trained the BPR model on PointRend or SegFix results and transferred them to refine the Mask R-CNN predictions. As shown in Table S1, the transferring is also workable.   Table S1: BPR models were trained on PointRend and Seg-Fix results respectively, and transferred to refine the Mask R-CNN predictions.


extfieldsofstudy:
Medicine
Physics
Computer Science

text:
Hybrid quantum investment optimization with minimal holding period

In this paper we propose a hybrid quantum-classical algorithm for dynamic portfolio optimization with minimal holding period. Our algorithm is based on sampling the near-optimal portfolios at each trading step using a quantum processor, and efficiently post-selecting to meet the minimal holding constraint. We found the optimal investment trajectory in a dataset of 50 assets spanning a 1 year trading period using the D-Wave 2000Q processor. Our method is remarkably efficient, and produces results much closer to the efficient frontier than typical portfolios. Moreover, we also show how our approach can easily produce trajectories adapted to different risk profiles, as typically offered in financial products. Our results are a clear example of how the combination of quantum and classical techniques can offer novel valuable tools to deal with real-life problems, beyond simple toy models, in current NISQ quantum processors.

The field of quantum computing is living an unprecedented expansion due to recent experimental advances. Quantum technologies as a whole [1][2][3][4][5][6][7][8] are being benefited from recent results both on quantum computer implementations [9][10][11][12][13][14][15] as well as on quantum communications [16][17][18][19][20][21][22] . Due to this, people have started to think seriously about industrial applications of quantum computers 4 . Among the different verticals, finance is among the most promising ones, given the ubiquitousness of intractable mathematical problems. For a detailed description of applications of quantum computing in finance, see Ref. 23 . Among these applications, one of the most prominent is quantum optimization. There are many important optimization problems in finance which can be solved more efficiently using quantum computing. See Refs. [24][25][26][27][28] for some examples. In this setting, the most paradigmatic optimization problem in finance is that of portfolio optimization, both in its static and dynamic versions.
Our aim in this work is to solve the dynamic portfolio optimization problem. An issue that investors often face is that short term investments tend to be taxed much higher than long term investments. It is common for investors to impose a minimal holding period, preventing any purchased asset from being sold before a predetermined period of time. We build upon the work from Ref. 29 , and demonstrate an efficient post-selection protocol to impose the minimal holding constraint. Similarly to Ref. 29 , we further impose that investors must invest in integer bundles, as is typically the case for exchange-traded funds (ETF) shares.
Financial model. Our aim is to find the best investment trajectory for an investor, given a level of risk that they want to take. According to Modern Portfolio Theory, the optimal investment at a defined level of risk is the one which maximizes profit 30 . The portfolio's risk-or volatility-is computed from the assets' covariance matrix. The ratio of returns to risk is the Sharpe ratio, which is our metric for comparing investments. These are defined mathematically in what follows.
Let us define ω tn , the fraction of the total budget invested in asset n at time t. The optimal holding trajectory ω tn minimizes the Modern Portfolio Theory cost function: www.nature.com/scientificreports/ Here, µ t is the vector of the logarithmic returns at time t and t is the covariance matrix. The risk aversion γ controls the portfolios penalty for risk. This determines the amount of risk an investor is willing to take. The forecasted returns and covariance matrices can be deduced from the stocks' prices. This is detailed, for instance, in Ref. 29 . Following this notation, the Sharpe ratio is given by which, as we said, is nothing but the ration between the total return (numerator) and volatility (denominator). At any point in time, we would like the entire available budget to be invested. This constrains the holdings ω t to be normalized at any time t. We enforced this by penalizing portfolios which did not respect this constraint. The cost function we optimized is therefore where the Lagrange multiplier ρ is an hyperparameter of the model. Note also that Eq. (3) can be written as: for some h t , t ∈ {0, N t } . This means that the optimal investment at time t is independent of our investment history, as long as no more constraints are included that correlate trading times among each other. In such a situation it is therefore sufficient to minimize h t at every time t to compute the optimal holdings ω t . This is a good approximation, for instance, for the trading investment funds of an everyday investor, for which transaction costs and market impact may be negligible.
We assume shares can only be sold in large bundles. These constraints imply that our objective variables-the instantaneous asset holdings ω tn -are integer variables. At any trading time, we can encode the qubit values to our qubits using a binary encoding: where x n,q ∈ {0, 1} is the readout value of the qth qubit assigned to institution n, and K is the total investment. This encoding has several consequences. First, the holdings are always bound by: ω tn ∈ [0, 1] . Since ω tn ≥ 0 , investors are not given the option to sell short securities. Allowing short-selling could be an interesting extension for this problem. Second, this encoding allows investors to split their total investment in a maximum of K bundles. This discretization is part of the reason why this problem is so hard to solve on classical computers 29 . From a practical perspective, using the bit variables the cost function can then be written as with x a vector of bit variables, and Q a matrix of real numbers. In the language of combinatorial optimization, this is a quadratic unconstrained binary optimization (QUBO) problem, and is the natural input for the D-Wave quantum annealer.
In Eq. (5), we have introduced the bit depth N q , the number of qubits which encode the nth asset holdings. Typically, we would choose N q such that 2 N q − 1 ≥ K . This choice allows the investor to invest their entire budget into a single asset. We may, however, want to impose a diversification constraint. Typically, large financial institutions are not allowed to invest more than 5% of their total budget in any single asset. Choosing N q such that 2 N q − 1 = 0.05K naturally implements this diversification constraint.

Minimum holding period.
The holding period is the amount of time which elapses between an investment's purchase and its sale (or sale of a security). Because long-term gains are taxed more favourably than short-term gains, it is common to demand that investments meet a minimum holding period. In this study, we imposed a minimal holding period of 7 days, which is a natural timeframe for some realistic banking products.
Here we imposed this constraint by post-selecting investment trajectories which respected this condition. The number of possible trajectories grows exponentially O N K×N t a . It is therefore computationally prohibitive to successively consider each trajectory and verify if they meet the minimal 7 day holding period.
Instead, we build the investment trajectory iteratively, following the flow chart presented Fig. 1. At every time t, we choose holdings ω guess for that instant. We compute them by sampling at each time t the low-energy subspace of h t in Eq. (4) using, in our case, the D-Wave machine. the fact that the problem can be decoupled into separate times t allows us to tackle more assets at every given time with the quantum annealer. The sampled instantaneous portfolios are considered in order of decreasing Sharpe ratio. As customary in finance, a large Sharpe ratio is an indicator of a good-quality portfolio, since it implies that the return is large for the risk being assumed. We note, however, that other figures of merit could be used, such as the returns (the larger the better), or the volatility (the lower the better). In our case, though we choose to work with the Sharpe ratio because of the afore-mentioned reasons. We only retain candidates which fulfil the minimum holding constraint. www.nature.com/scientificreports/ By following iteratively, the algorithm implements a tree-like selection and efficiently rules out trajectories which do not meet the minimal 7-day holding period, as is illustrated in Fig. 2. Note that at any time t, there always exists at least one solution which meets the minimal 7-day holding period: the solution ω t = ω t−1 .
Alternatively, we could include this constraint by penalizing in Eq. (3) trajectories which do not meet the minimum holding period. The corresponding penalty term would then correlate different time steps between each other, such that the cost function would not have the separable form of Eq. (4). In such a case, the problem would become even more computationally intractable, and a global optimization strategy over all possible trading times would be necessary from the very beginning. It is therefore remarkable that our simple approach, based on quantum sampling and post-selection, produces such high-quality portfolios in an extremely efficient manner.  In Ref. 29 , authors described a method to cluster assets based on their time series' similarity. In the following, we apply this method, and discard all except for the best asset in each cluster (based on its historical Sharp ratio). The dimensional reduction strategy follows from applying first a Hodrick-Prescott smoothing, to extract data trends, and then computing the euclidean distance between different pairs of trends of assets. Thanks to this, we can perform a clustering of the assets into clusters of assets of similar behavior. This dimensional reduction strategy also allows us to significantly lower the portfolio's risk, by diversifying our investment among maximally uncorrelated assets. Moreover, our target investors typically select their investment risk category. To meet this risk requirement, we performed a pre-selection among available assets. We computed each assets' historical volatility, and discarded any asset which significantly exceeded the agreed risk level. This allows us to construct portfolios for investors with different risk profiles: high-risk, medium-risk, low-risk, and so forth.
In practice, we determined the optimal number of clusters N c by studying each clusters' variance. In the studied data, we found that the mean clusters' variance decreased exponentially with the number of clusters. For N c > 7 , the clusters' variance did not significantly decrease, indicating that 7 clusters capture almost all of the system's variations. We therefore set N c = 7 in our calculations.

Results
Following the protocol detailed in previous sections, based on looking at optimal Sharpe ratios of portfolios computed with the D-Wave quantum annealer, we found the optimal investment trajectory between May the 31st 2019 and May the 31st 2020 among a total of 50 international assets and indices. Daily transactions were allowed, but investments had to meet a minimum 7 day holding period. We allowed investments to be split into a maximum of K = 5 bundles, and set the bit depth to N q = 2.
We found the optimal portfolios by minimizing Eq. (3) on the D-Wave 2000Q processor. Computed investment trajectories for five risk packages (minimum risk, 5% , 10% , 15% , and 20% annual risk) are plotted in Fig. 3. Also plotted are 1000 random investment trajectories which meet the minimum 7 day holding period. These give us some intuition about the space of allowed portfolios. We see that allowed portfolios tend to exist between −20 and 20% annualized profit, and between 10 and 17% risk (as measured by the annualized volatility). It is clearly visible that the investment trajectories obtained using our method based on quantum optimization, dimensional reduction and post-selection, tend to gather much closer to the efficient frontier in portfolio space than those obtained by investing at random.
In Fig. 4, we show an example investment trajectory. This corresponds to the 15% risk package (purple dot in Fig. 3), which yields ≈ 60% return on investment. We can see from this figure that the N c = 7 most profitable assets which meet the 15% risk limit are Bankia Fonduxo Universal FI, Bankia Futuro Sostenible Universal FI, Bankia Indice Eurostoxx Universal FI, Amundi Fds Pioneer US Bd A2 EUR C, Fidelity Global Technology A-Acc-EUR, Nordea 1 − Chinese Equity Y EUR, and Franklin India A(acc)EUR. As can be read off from Fig. 4, the optimal investment on May the 31st 2019 meeting these constraints is: 1/5th in Fidelity Global Technology A-Acc-EUR, 1/5th in Amundi Fds Pioneer US Bd A2 EUR C, 1/5th in Franklin India A(acc)EUR, and 2/5th in Nordea 1 − Chinese Equity Y EUR. The everyday composition of the trading trajectory can be read form the chart, and it always fulfils by construction the minimum holding period.
Finally, we would like to discuss about the computational running time and quality of the computed portfolios, in comparison with other classical and quantum optimization strategies. The running time to compute the optimal portfolios in Fig. 3 using our strategy of combining D-Wave and post-selection was of just a few minutes,   29 . Compared to classical solvers, our method is way faster than brute-force classical search, which scales exponentially fast in the number of qubits and quickly hits astronomical numbers (the problem is clearly intractable in this sense for 50 assets). Compared to other solvers, D-Wave quantum annealing is known to be faster and allow larger portfolios than solvers such as Gekko, see Table IV in Ref. 29 , which equally applies in our case.

Conclusion
Here we have proposed a hybrid quantum-classical algorithm for dynamic portfolio optimization with minimal holding period. Our algorithm is based on quantum sampling of near-optimal portfolios at each trading step, combined with dimensional reduction and post-selection. We tested our idea using quantum optimization with the D-Wave 2000Q processor, and real data from 50 assets during 1 year. Our study shows that the method is remarkably efficient and produces in few minutes results close to the optimal efficient frontier in portfolio space, much better than typical random portfolios. We also showed how our approach can be used to produce trajectories corresponding to different risk profiles, in accordance with typical products offered by financial institutions to investors. Our results are a clear example of how the combination of quantum and classical techniques can provide novel and valuable methods to deal with real-life financial data in today's NISQ quantum processors. Moreover, we believe that the use of distributed quantum computer architectures, such as those discussed in Ref. 3 , could accelerate even more the implementation of this type of quantum computing algorithms to solve real-life financial mathematical problems.
Received: 19 December 2020; Accepted: 1 September 2021 Figure 4. The optimal investment trajectory between May the 31st 2019 and May the 31st 2020 among considered assets for an investor wanting to take ≤ 15% risk. As can be seen from Fig. 3, this portfolio provided ≈ 60% annualized return on investment (purple dot in the Fig. 3). Data for this calculation was obtained from daily prices of international assets and indices as explained in the text, which are public and can be obtained from e.g. Bloomberg, Yahoo Finance, and/or Morningstar databases.