Online Learning for Distribution-Free Prediction
We develop an online learning method for prediction, which is important in problems with large and/or streaming data sets. We formulate the learning approach using a covariance-fitting methodology, and show that the resulting predictor has desirable computational and distribution-free properties: It is implemented online with a runtime that scales linearly in the number of samples; has a constant memory requirement; avoids local minima problems; and prunes away redundant feature dimensions without relying on restrictive assumptions on the data distribution. In conjunction with the split conformal approach, it also produces distribution-free prediction confidence intervals in a computationally efficient manner. The method is demonstrated on both real and synthetic datasets.READ FULL TEXT VIEW PDF
We study three notions of uncertainty quantification—calibration,
We develop a method for assessing counterfactual predictions with multip...
Successful application of machine learning models to real-world predicti...
We introduce an efficient algorithmic framework for model selection in o...
Dyadic Data Prediction (DDP) is an important problem in many research ar...
An ever increasing volume of data is nowadays becoming available in a
Topology optimization by optimally distributing materials in a given dom...
Online Learning for Distribution-Free Prediction
Prediction is a classical problem in statistics, signal processing, system identification and machine learning[1, 2, 3, 4]. The prediction of stochastic processes was pioneered in the temporal and spatial domains, cf. [5, 6, 7, 8], but the fundamental ideas were generalized to arbitrary domains. In general, the problem can be formulated as predicting the output of a process, , for a given input test point after observing a dataset of input-output pairs
In this paper, we are interested in learning predictor functions in scenarios where is very large or is increasing. We consider two common classes of predictors:
The first class consists of predictors in linear regression (Lr) form. That is, given a regressor function , the predictor
is expressed as a linear combination of its elements. This class includes ridge regression,Lasso, and elastic net [9, 10, 11]. The regressor can be understood as a set of input features and a standard choice for it is simply .
The second class comprises predictors in linear combiner (Lc) form. That is, the predictor is expressed as a linear combination of all observed samples , using a model of the process. This class includes kernel smoothing, Gaussian process regression, and local polynomial regression [12, 13, 14, 15, 16]. A standard model choice is a squared-exponential kernel or covariance function for .
The weights in the Lr and Lc classes are denoted and , respectively. For a given prediction method in either class, the weights are functions of the data
. The set of possible weights is large and must be constrained in a balanced manner to avoid a high variance and bias of, which correspond to overfitting and underfitting, respectively.
In most prediction methods, the constraints on the set of weights are controlled using a hyperparameter, which we denoteand which is learned from data. There exist several methods for doing so, including cross-validation, maximum likelihood, and weighted least-squares [1, 17, 18]. These learning methods are, however, nonconvex and not readily scalable to large , which is the target case of this paper. Moreover, in the process of quantifying the prediction uncertainty, using e.g. the bootstrap  or conformal approaches [20, 21], the learning methods compound the complexity and render such uncertainty quantification intractable for large .
In this paper, we consider a class of predictors that can equivalently be expressed in either linear regression or combiner form. In this class, constrains the weights for each dimension of individually. Thus irrelevant features can be suppressed and overfitting mitigated [22, 23, 4]. To learn the hyperparameters , we employ a covariance-fitting methodology [18, 24, 25] by generalizing a fitting criterion used in [26, 27, 28]. We extend this learning approach to a predictive setting, which results in a predictor with the following attributes:
computable online in linear runtime,
implementable with constant memory,
does not suffer from local minima problems,
enables tractable distribution-free confidence intervals,
prunes away irrelevant feature dimensions.
These facts render the predictor particularly suitable for scenarios with large and/or growing number of data points. It can be viewed as an online, distribution-free alternative to the automatic relevance determination approach [22, 23, 4]. Our contributions include generalizing the covariance-fitting methodology to non-zero mean structures, providing connections between linear regression and combiner-type predictors, and an analysis of prediction performance when the distributional form of the data is unknown.
The remainder of the paper is organized as follows. In Section 2, we introduce the problem of learning hyperparameters. In Section 3, we highlight desirable constraints on the weights in an Lr setting, whereas Section 4 highlights these constraints in an Lc setting. The hyperparameters constrain the weights in different ways in each setting. In Section 5, the covariance-fitting based learning method is introduced and applied to an Lc predictor. The resulting computational and distribution-free properties are derived. Finally, in Section 6, the proposed online learning approach is compared with the offline cross-validation approach on a series of real and synthetic datasets.
In the interest of reproducible research, we have made the code for the proposed method available at https://github.com/dzachariah/online-learning.
Notation: is the elementwise Hadamard product. The operation stacks all
into a single column vector, whiledenotes the th column of . The sample mean is written as . The number of nonzero elements in a vector is denoted as . The Kronecker delta is denoted .
Abbreviations: Independent and identically distributed (i.i.d.).
Let denote a predictor with hyperparameters . For linear regression and combiner-type predictors, the choice of constrains the set of weights. For any input-output pair , the risk of the predictor is taken to be the mean-squared error 
The optimal choice of
is therefore the hyperparameter that minimizes the unknown risk. A common distribution-free learning approach is cross-validation, which estimates the risk for a fixed choice of. The risk estimate is formed by first dividing the training data into subsets and then predicting the output in one subset using data from the remaining subsets[1, ch. 7]:
where is the number of samples in subset and denotes the predictor using training data from all subsets except . In this approach, the hyperparameter is learned by finding that minimizes . For one-dimensional hyperparameters, a numerical search for the minimum of (2) is feasible when is small; otherwise it may well be impractical.
An alternative learning approach is to assume that the distribution of belongs to a family that is parameterized by
. Then the hyperparameter can be learned by fitting a covariance model of the data to the empirical moments of the training data, cf.. If additional assumptions are made about the distributional form, it is possible to formulate a complete probabilistic model of the data and learn using the asymptotically efficient maximum likelihood approach, cf. [30, 17].
The aforementioned statistical learning methods are, however, in general nonconvex and may therefore give rise to multiple minima problems [31, 32]. This becomes problematic when is multidimensional as the methods may require a careful choice of initialization and numerical search techniques [33, ch. 5]. In addition, the distributional assumptions employed in the maximum likelihood approach may lack robustness to model misspecifications .
More importantly for the data scenarios considered in this paper, these learning methods cannot readily be implemented online. That is, for a given dataset, the computational complexity of learning does not scale well with and the process must be repeated each time increases. Consequently, for each new , the predictor function must be computed afresh.
Our main goal is to develop a learning approach that is implementable online, obviates local minima problems, and has desirable distribution-free properties. Before we introduce the proposed learning approach, we introduce the hyperparameters in the context of Lr and Lc predictors, respectively.
Lr predictors are written in the following form
where is a -dimensional regression function and are weights. The empirical risk of a predictor can then be written as
When the empirical risk is minimized under a constraint on , irrelevant dimensions of can be suppressed which mitigates overfitting [1, 35]. For notational convenience we write the regressor matrix as
where denotes the th column.
Let denote corresponding optimal Lr predictor. This sparse predictor does not rely on any assumptions on the data, and is therefore robust to model misspecifications. We collect the prediction errors of the training data in the vector
where . In the subsequent analysis, we also make use of
Solving (5) is a computationally intractable problem, even for moderately sized . Thus will be taken as a reference against which any alternative can be compared by defining the prediction divergence :
where is a hyperparameter that shrinks all weights toward zero in a sparsifying fashion. This Lr predictor has desirable risk-minimizing properties and can be computed online in linear runtime for any fixed [40, 41]. The corresponding divergence can be bounded per the following result.
If the hyperparameter satisfies
then the divergence of the Lasso-based predictor from the optimal sparse predictor is bounded by
See Appendix A.1 ∎
If the hyperparameter in (8) satisfies (9), then (10) ensures that the maximum divergence of Lasso-based predictor from the optimal Lr predictor remains bounded. Then redundant dimensions of are pruned away, which is a desirable property. The hyperparameter in (8) is typically learned offline using (2), which is a limitation for the scenarios of large and/or increasing considered herein. Moreover, the individual weights in (3) are constrained uniformly in the Lasso approach. To provide a more flexible predictor that constrains the weights individually, we consider an Lc approach next.
For the special case in which the prediction errors of are i.i.d. zero mean Gaussian with variance , and the regressors are normalized such that , the inequality (9
) is satisfied with high probability whenfor some positive constant , cf. . However, this is an infeasible choice since is unknown and must be estimated.
Lc predictors are written in the following form
where contains the observed samples and are weights. Given any , we will find an unbiased predictor of that minimizes the conditional risk. That is, we seek the solution to
and denotes all inputs in . The weights are then constrained by assuming a model of given .
Here we specify a simple model that does not rely on the distributional form of but merely constrains its moments :
The conditional mean of is parameterized as
where is a given function and are unknown coefficents. In the simplest case we have an unknown constant mean, i.e., , by setting . Alternatively, one may consider an unknown affine function.
To enable an computationally efficient online implementation, we consider a model with a truncated sum of terms in (14). We write
and the hyperparameter as
for notational simplicity.
The result follows from an elementwise application of (14) to the training data:
The weights of the optimal Lc predictor in (11) are given by
and is a projector onto .
See Appendix B. ∎
Consequently, the hyperparameters in constrain the optimal weights via the model (13) and (14). If the model were correctly specified, we could in principle search for a risk-minimizing . This would, however, be intractable in an online setting.
See Appendix C. ∎
The Lc predictor has an Lr intepretation: when (14) is formed using the Fourier or Laplace operator basis, the corresponding regressor (18) will approximate any isotropic covariance function by varying . This includes covariance functions belonging to the important Matérn class, of which the popular squared-exponential function is a special case [30, ch. 6.5]. Cf. [42, 43, 45] for regressor functions with good approximation properties. Similarly, by partitioning an arbitrary regressor into the form (18) yields an Lc interpretation of an Lr predictor. The regressor matrix is then . For instance, a standard regressor is the affine function , which can be interpreted as a constant mean with a covariance function that is quadratic in .
The elements of constrain each of the weights (19) individually, unlike the uniform approach in (8). Thus the hyperparameter in the optimal Lc predictor determines the relevance of individual features in , similar to the automatic relevance determination framework in [23, 46, 47]. While this enables a more flexible predictor, the -regularization term in (19) does not yield the sparsity property of the Lasso-based predictor. More importantly, the learning approach (2) is intractable when is multidimensional. We show how to get around this issue in the next section.
We consider a covariance-fitting approach for learning in the flexible Lc predictor. We show that using the learned hyperparameter results in a predictor with desirable, computational and distribution-free, properties.
The normalized sample covariance matrix of the training data can be written as
using a weighted norm that penalizes correlated residuals, cf. [18, 24, 25]. To match the magnitude of the normalized sample covariance we subject the parameters in (20) to the normalization constraint
The learning problem defined by (20) is convex in and will therefore not suffer from local minima issues.
See Appendix D. ∎
The Lc predictor with a learned parameter has the following Lr form:
and the elements of are set to
See Appendix E. ∎
Thus the covariance-based learning approach endows the predictor with a sparsity property, via the weighted -regularization term in (23), that prunes away redundant feature dimensions. Similar to (10), we can bound the associated prediction divergence.
If all elements (24) are positive and satisfy
then the divergence of the Spice predictor from the optimal sparse predictor is bounded by
See Appendix A.2. ∎
The above results are valid even when the model (13) and (14) is misspecified. We have therefore shown a distribution-free property of the Spice predictor, which will prune away redundant feature dimensions automatically. It parallels (10) but is not dependent on any hyperparameters. Thus we can view the covariance-fitting approach as a vehicle for constructing a general Lr predictor with weights on the form (23).
In addition, the Spice predictor has computational properties that are appealing in cases with large and/or increasing , where an online implemention is desirable.
The Spice predictor (22) can be updated at each new data point in an online manner. The computations are based on
which have fixed dimensions and are readily updated sequentially. The pseudocode for the online learning algorithm is summarized in Algorithm 1. It is initialized by setting the weights to .
See Appendix F. ∎
The total runtime of the algorithm is , where is the number of cycles per sample, and its memory requirement is .
The complexity of the loop in Algorithm 1 is proportional to for each sample. Regarding the memory, the number of stored variables is constant and the largest variable is the matrix . ∎
In the given algorithm, the individual weights are updated in a fixed cyclic order. Other possible updating orders are described in a general context, in . As , however, the algorithm converges to the global minimizer (23) at each sample , irrespective of the order in which are updated [51, 52].
Let denote the convex cost function in (23) at sample . The termination point of the cyclic algorithm for provides the starting point for minimizing the subsequent cost function . For finite , the termination point will deviate from the global minimizer (23). In practice, however, we found that can be chosen as a small integer and that even yields good prediction results when . Furthermore, we observed that the results are robust with respect to the orders in which are updated.
In summary, the covariance-fitting methodology above has the following main attributes:
avoids local minima problems in the learning process,
results in a predictor that can be implemented online
maximum divergence from the optimal sparse predictor can be evaluated.
Its computational and distribution-free properties also make it possible to combine this approach with the split conformal method in , which provides computationally efficient uncertainty measures for a predictor under minimal assumptions.
Suppose the input-output data consist of i.i.d. realizations from an unknown distribution
For a generic point , we would like to construct a confidence interval for the predictor (22) with a targeted coverage. That is, find a finite interval
that covers the predicted output with a probability that reaches a prespecified level .
For simplicity, assume that is an even number and randomly split into two equally-sized datasets and . For a given targeted coverage level , the split conformal interval is constructed using the following three steps :
Train the Spice predictor using .
Predict the outputs in and compute the residuals
Sort the residuals and let denote the th smallest , where .
Setting as above in (28), yields an interval that covers the predicted output with a probability
Thus the targeted level can be ensured. In addition, when the residuals have a continuous distribution, the probability is also bounded from above by .
See [21, sec. 2.2]. ∎
, provide credibility intervals but, in contrast to the proposed approach, lack coverage guarantees. In fact, even when the said models are correctly specified, their uncertainty is systematically underestimated after learning the hyperparameters.
The Spice-based split conformal prediction interval in (28) provides a computationally efficient, distribution-free prediction and inference methodology.
In this section, we compare the online Spice approach developed above with the well-established offline -fold cross-validation approach for learning predictors that use a given regressor function . To balance the bias and variance of (2), as well as the associated computational burden, we follow the recommended choice of folds [1, ch. 7]. Finding the global minimizer of (2) in high dimensions is intractable, and in the following examples we restrict the discussion to the learning of predictors based on the scalable Ridge and Lasso regression methods, since these require only one hyperparameter.
Offline cross-validation was performed by evaluating (2) on a grid of 10 hyperparameter values and selecting the best value. Each evaluation requires retraining the predictor with a new hyperparameter value and the entire process is computationally intensive.
For the prediction problems below, we apply regression functions of the form (18):
using a constant mean . For we use either a linear function, , or the Laplace operator basis due to its attractive approximation properties . For the latter choice, suppose is -dimensional and belongs to . Then the elements of are defined by
where are the indices for dimension . We have that has dimension . The rectangular domain can easily be translated to any arbitrary point. When is large, we may apply the basis to each dimension separately. Then the resulting has dimension .
In this experiment we study the learning of predictors under two challenging conditions: heavy-tailed noise and colinear regressors. The input is of dimension . The dataset was generated using an conditional Student-t distribution with mean
and variance . The input were generated using an i.i.d. degenerate zero-mean Gaussian variable with covariance matrix , where the numerical rank of and the variances are normalized by setting .
For the predictors, we let be linear and thus . We ran Monte Carlo simulations to evaluate the performances of the predictors. In the first set of experiments we estimate the risk . To clarify the comparison between the predictors, the risk is normalized by the noise variance and presented in decibel scale (dB) in Table 1. As expected for this data generating process, the sparse predictors outperform Ridge. For Spice and Lasso, the difference is notable when but is less significant as more samples are obtained.
Next, we evaluate the inferential properties by repeating the above experiments with samples. The dataset is randomly partitioned into two sets and , each of size , to produce confidence intervals as in (28). We target the coverage level , and report the average confidence interval length as well as average coverage of the interval in Table 2. Note that the probability that is nearly exactly equal to the targeted level without relying on any distributional assumptions. Thus the reported confidence intervals are accurate. Furthermore, the average interval lengths are significantly smaller for the sparse predictors compared to Ridge. The reported intervals for Spice and Lasso are similar in length, with the former being slightly smaller. The interval lengths can be compared to the dynamic range of , which has a length of approximately 60.
|50||7.74 (0.90)||21.04 (0.90)||8.13 (0.90)|
|100||6.33 (0.90)||9.83 (0.90)||6.40 (0.90)|
|200||5.48 (0.90)||8.02 (0.90)||5.56 (0.90)|
The runtime for the offline learning approach using Ridge or Lasso, is which is similar to the runtime for the online Spice method , where in this example. The average runtimes are reported in Table 3. While all three methods scale linearly in , using a cross-validated Lasso predictor is slower in this implementation.
The ozone density determines the transmission of ultraviolet radiation through the atmosphere which has an important impact on biochemical processes and health. For this reason, measuring the total column ozone has been of interest to scientists for decades. In 1978, the Nimbus-7 polar orbiting satellite was launched, equipped with a total ozone mapping spectrometer. The satellite was sun synchronous, and due to the rotation of the Earth, its scan covered the entire globe in a 24 hour period with a certain spatial resolution . For illustrative purposes, we consider a set of spatial samples of ozone density from the satellite, measured in Dobson units (DU) and recorded on October 1st, 1988, cf. . The spatial coordinates were transformed from longitude and latitude using the area-preserving Mollweide map projection .
For , we use the Laplace operator basis with so that . The boundaries were set slightly larger than those given by the Mollweide projection: and , where is the radius of the Earth.
In the first experiment, the training is performed using samples and the ozone density is predicted on a fine spatial scale. Note that this dataset is nearly three orders of magnitude larger than that used in the previous example, which makes it too time consuming to implement the offline cross-validation method due to its computational requirement. Therefore we only evaluate the Spice predictor here. Fig. 1 illustrates both the training samples and the predicted ozone density
. It can be seen that the satellite data is not uniformly sampled and, moreover, it contains significant gaps. The predictions in these gapped areas appear to interpolate certain nontrivial patterns.
In the second experiment, we evaluate the inferential properties of the Spice predictor by learning from a small random subset of the data. We use samples, or approximately of the data. The resulting predictions exhibit discernible patterns as in Fig. 2.
The remaining samples are used for validating the predictor. Its risk is estimated by the out-of-sample prediction errors,
which we translate to Dobson units by taking the square-root. The result was a root-risk of DU. In addition, the confidence interval with had a length of DU and an empirical coverage of 0.90. Both performance metrics compare well with the dynamic range of the data, which spans [179.40, 542.00] or DU. Fig. 4 also shows that the empirical distribution of the prediction errors is symmetric. These results illustrate the ability of the proposed method to learn, predict and infer in real, large-scale datasets.
In this paper we considered the problem of online learning for prediction problems with large and/or streaming data sets. Starting from a flexible Lc predictor, we formulated a convex covariance-fitting methodology for learning its hyperparameters. The hyperparameters constrain individual weights of the predictor, similar to the automatic relevance determination framework.
It was shown that using the learned hyperparameters results in a predictor with desirable computational and distribution-free properties. We denote it as the Spice predictor. It was implemented online with a runtime that scales linearly in the number of samples, its memory requirement is constant, it avoids local minima issues, and prunes away redundant feature dimensions without relying on assumed properties of the data distributions. In conjunction with the split conformal approach, it also produces distribution-free prediction confidence intervals. Finally, the Spice predictor performance was demonstrated on both real and synthetic datasets, and compared with the offline cross-validation approach.
In future work, we will investigate input-dependent confidence intervals, using the locally-weighted split conformal approach.
We expand the empirical risk of by
Using Hölder’s inequality along with the triangle inequality yields
Thus, the prediction divergence is bounded by
For notational simplicity, we write
since . By inserting and into the cost function of (23), we have that
Multipling the inequality by and rearranging, yields
where the second inequality follows from using in (33).
Given (25), we have that
By re-arranging and dividing by , we obtain the following inequality
Therefore (35) can be bounded by
Note that in the special case when the prediction errors of are i.i.d. , and the regressors are normalized as , we have that
with probability greater than where is a positive constant, cf. [39, ch. 6.2]. In this setting is a consistent estimate of . Therefore (25) is satisfied with high probability in this case if the elements of are multiplied by a factor , so that
for some .
The conditional risk can be decomposed as
using (15). Here is a constant and the bias vanishes due to the unbiasedness constraint on . This constraint, in turn, can be expressed as
Note that a weak assumption in our case is , especially when and . Hence there exist that satisfy the equality above.
Thus the problem (12) can be reformulated using the convex Lagrangian
where is the vector of Lagrange multipliers.
A stationary point of the Lagrangian in (38) with respect to and satisfies
The first equality can be expressed as . We can insert the first equality into the second and solve for the Lagrange multipliers
Then we have
Noting that the orthogonal projector is given by concludes the proof.
Let be partitioned as in (18). Then
in combination with (16) allows us to identify the weights:
Next, define the problem
for which the minimizing , namely,
By expanding the criterion (20) we obtain
Next, we define an auxiliary variable that satisfies
The derivation follows in three steps. First we show that dropping the constraint (21) yields the simpler unconstrained problem in (20). Next, we prove that the fitted hyperparameters in the simplified problem yields the same predictor. Finally, we show how the corresponding Lr form arises as a consequence.
By dropping the constraint (21) we may consider the problem in the following form:
where is the unnormalized sample covariance matrix.