Online Learning for Prediction via Covariance Fitting: Computation, Performance and Robustness
Authors: Muhammad Osama
Dave Zachariah
Peter Stoica
Thomas B. Schön
Venue: Accepted by TMLR ()

Abstract: We consider the problem of online prediction using linear smoothers that are functions of a nominal covariance model with unknown parameters. The model parameters are often learned using cross-validation or maximum-likelihood techniques. But when training data arrives in a streaming fashion, the implementation of such techniques can only be done in an approximate manner. Even if this limitation could be overcome, there appears to be no clear-cut results on the statistical properties of the resulting predictor. 

Here we consider a covariance-fitting method to learn the model parameters, which was initially developed for spectral estimation. We first show that the use of this approach results in a computationally efficient online learning method in which the resulting predictor can be updated sequentially. We then  prove that, with high probability, its out-of-sample error approaches the optimal level at a root-$n$ rate, where $n$ is the number of data samples. This is so even if the nominal covariance model is misspecified. Moreover, we show that the resulting predictor enjoys two robustness properties. First, it corresponds to a predictor that minimizes the out-of-sample error with respect to the least favourable distribution within a given Wasserstein distance from the empirical distribution. Second, it is robust against errors in the covariate training data. We illustrate the performance of the proposed method in a numerical experiment.

Published in Transactions on Machine Learning Research (01/2023)
Online Learning for Prediction via Covariance Fitting: Computation, Performance and Robustness
Muhammad Osama Deportment of Information Technology,
muhammad.osama su.
Uppsala University
Dave Zachariah
dave.zacharial waw.chharraaaa....oo wu.se
Department of Information Technology, Uppsala University
Petre Stoica
pt su.se
Department of Information Technology, Uppsala University
Thomas B. Schön
thomas schon Oit. wu.se
Department of Information. Technology, Uppsala University
Reviewed OpenReview: https: opte:ree//penreviee.u net/ forum? id=nAr9PhyEbQ
Abstract
We consider the problem online prediction using linear smoothers that are functions of a nominal covariance model with unknown parameters. The model parameters are often learned using cross-validation maximum-likelihood techniques. But when training data arrives in streaming fashion, the implementation of such techniques only be done in an approximate manner. Even if this linitation could be overcome, there appears to be no clear-cut results on the statistical properties of the resulting predictor.
Here consider covariance-fitiing method to learn the model parameters, which was initially developed for spectral estimation. We first show that the use of this approach results in a computationally efficient online learning method in which the resulting predictor can be updated sequentially. We then prove that, with high probability, error approaches the optimal level a root-n rate, where the number of data samples. This is even if the nominal covariance model misspecified. Moreover, show that the resulting predictor enjoys two robustness properties. First, it corresponds to predictor that minimizes the out-of-sample error with respect to the least favourable distribution within a given Wasserstein distance from the empirical distribution. Second, it is robust against errors in the covariate training data. We illustrate the performance of the proposed method in a numerical experiment.
out-of-sample
1 Introduction
We consider scenarios in which we observe stream of randomly distributed data
D, (x11,)1), ..., (Xn: Mn)}
1,2,3,

---

Published in Transactions on Machine Learning Research (01/2023)
Given covariate Xx+1 in space X, goal is to predict the outcome Mc++ in R. A large class of predictors (also known as linear smoothers) can be described as weighted combination of observed outcomes:
xx; Σ w/ X)ui.
(1)
where x denotes any test point and the weights [w; 2)) are parametorizod by A. The scnsitivity of such a predictor function noise in the training data is often characterized by how close the in-sample prediction 7(x; A] is to Yi and quantified by the sum of in-sample weights,
0 < dfn A
w. (x;;
(2)
also known as the 'effective' degrecs of freedom (Ruppert et al. 2003 Wasserman 2006 Hastic et al., 2009 These degrees of freedom are oflen luned lo avoid overfitting to the irreducible noise in the training data with the aim of achieving good out-of-sample performance. This includes learning the parameters À from Dn via distribution-free cross- validation distribution- maximum likelihood methods, which however can typically be implemented only approximately in the online scenario.
In this paper, we consider an alternative method using covariance- criterion first proposed in the context of spectral cstimation (toocaa et al. 2010a b). We show that this method
enables sequential computation of predictor with learned parameters,
approaches optimal out-of-sample performance at root-n rate,
enjoys two types of robustness properties.
For illustration the online learning method, we include a numerical experiment.
Notation: |Z]| w ttr(Z WZ] is weighted Frobenious norm of matrix Z using a positive definite weight matrix W. The element-wie Hadamard product between z and z' is denoted z @ z'.
2 Problem formulation
The linear smoother predictor 1 can be written compactly as
0xx = w A)y, where w (x; A) = [w(x]: A),... (x,;; A)] and [Bi...., "bt......... I.
(3)
We will investigate a class of model-based weights for which (3) can be computed sequentially from the stream D... Specifically, suppose is modeled as zero-mean stochastic process with a nominal covariance function parameterized as
Covlu. xx Ago(x, (x,xxx'' x') + Σ Ancia (x)xx)))
(4)
where x and x' are two arbitrary covariates, 5(x, x') the Kronecker delta function and (dx(x)) are real- valued features x. Here À the unknown set d nonnegative model covariance parameters. When belongs a vector space, then periodic Fourier-type features provide convenient choice due their
excellent covariance approximating properties (Ruppert et al. 2003 Rahimi & Recht 2007 Stein 2012
Hensman et al. 2017 Solin & Särkkä 2020). Using the notation above, set of optimal weights can be
written as
w(x; A) = C *
(5)
Bishop 2006
Rasmussen & Williams
2006
Stein
2012
where
T (x)
Φ T (xi)
Φ(x)=
(x))
TT((x.)) (xn)
2

---

Published in Transactions on Machine Learning Research (01/2023)
contain the features and
Cx ΦΛΦ T + A0 > 0
(6)
is a nominal covariance matrix where = diag(Ai, ... .4)) The predictor function above includes a varicty of penalized regression methods (see the references cited above). The degroes of freedom of (3) are controlled by À and we have that.
0< dA) = - c = -r
+ + S min(n, d)
Stoica & Stanasila,
1982 Rupport et al.
2003
Let An denote the parameters that are filled to the data stream Dn in some way. The concern of this paper is two-fold: (i) sequential computation of (x; An) and (ii) performance guarantees that hold even when the model class is not well-specified. That is, the nominal covariance model (4) using {Φ)(x)) may not match the unknown covariance structure of the data. Nevertheless, we aim to derive meaningful guarantees that remain valid also in this misspecified case.
Cross-validation or maximum likelihood methods are two popular approaches to fitting A. However, these methods have challenges the online learning setting. First, computing X, for the covariance model above is non-convex problem that be riddled with multiple local minima. Second, for each additional training data point (X++1, Ma+1:Bu++); the parameters Xx+1 will have to be refitted and therefore [(x; X,,) recomputed from scratch approximated (see below). Third, to the best of our knowledge, finile out-of-sample prediction performance guarantees for (x;; Xn) neither exist for eross-validation nor maximum likelihood An-
We will consider an alternative covariance-based fitting criterion for X, used in another context, viz. spectral estimation (Stoica et al. 2010aa b). The main contribution is to derive performance guarantees of a sequen- tially computable predictor, thereby establishing predictive properties of the covariance-iittnng methodology. Similar results do not appear to be available in the previous literature.
3 Other approaches to online model-based prediction
We identify two main lines of work that enable online fast implementations of model-based linear smoothers, which include Gaussian process regression and kriging methods: The first line considers spectral model approximations (e.g., (Rahimi & Recht 2007 Hensman et al. 2017; Solin & Särkkä 2020 ) which is covered by the class of covariance models in (4). These mothods also enable efficient online computation of ((x; A), but for a fixed set of model parameters A. The second line considers sparse variational approx- imations (c.g., (Titsias] 2009; Bui et al., 2017 Stanton et al.] 2021)). These methods can rocompute the predictor 7xx; A) efficiently when new data arrives, but again for fixed A.
Fitting An and computing Tx: An) in sequential manner requires recomputing past covariance quantities and is computationally prohibitive. For the covariance model above, it means recomputing 5 and the inverse of 6 at each new sample. Work on this problem appears to be scarce. [Stanton et al. (2021 consider seeking maximum likelihood estimates An by a type of gradient-based search that projects past quantities onto a diver-dimeesional space. However, neither a convergence analysis towards the sought maximum likelihood estimate nor any resulting predictive properties are provided. By contrast, Cai & Yuan (2012)
provide an asymplotic performance guarantee of linear smoother with learned regularization parameter, but their approach is restricted to the offline setting.
Our focus below is on deriving finite-sample performance guarantees for a sequentially computable predictor.
4 Learning via covariance fitting
The parameters À can be learned by fitting the model covariance matrix CA in (5) to the empirical covariance matrix Ty Specifically, we will use the following fitting criterion, known as SPICE (Stoica et al. 2010a b).
An arg min I|yy Cxlle
(8)
X
3

---

Published in Transactions on Machine Learning Research (01/2023)
where A, is function of the datastream Dr. Since this criterion is convex in X, we can be sure that global minimizer can be determined.
4.1 Sequential computation
The covariance-fitting criterion in can be connected an alternative convex problem and this connection
be utilized to perform parameter estimation online, shown in (Zachariah & Stoica 2015). We first leverage this connection lo compule T(x; Xn) in sequential manner, and subsequently we derive a series of predictive performance guarantees of the proposed method.
Theorem
The predictor function T(x; 1 can be updated from ÿ(x: 2n) in a constant runtime O (d²).
The total memory requirement of the method is also on the order of (dP).
Proof 1 We first note that the predictor has equivalent form
x x)
A
(9)
0(2)
For later use, we also note that is invoriant to a uniform rescaling of X, since
Ō()) Ō Va 0
(10)
The covariance-fiting criterion in 8 can be expanded into
Iyy -
+K,
(11)
where constant. Thus An is also minimizer of
V(A) = y 1y
Next, we follow Zachariah & Stoica 2015 rion,
Appendia A) and consider the following augmented conver crite-
1 ly 2012 0 A - 100+r(CC)) (r(Cx),
V' A)
(12)
Ao
and show that ils minimizzes produce the predictor ((x;An) Φ (x)(((An) The first argument of (12) is minémized by
θ = - - "
A& Cx
where the second equality follows from using the matria inversion lemma. Inserting the minimizer it back into 12), have concentrated cost function:
Let us consider the minimizing A. By rescaling the parameters by ||||->> >0, 0, we have that
a v'(((o)), a2) = =a (y
=V(A)
where a = || > 0. It follows that ad,, is minimizer of 12
since
((((a))), 1 v(2m) - V(2) = v'(a())))

---

ad++1 In Published in Transactions on Machine Learning Research (01/2023)
This scaled minimizer can be related to the predictor with fitted parameters: (xx; x (x))(2n)) Φ (x))(a)))) using 100.
We now change the order of the minimization of (' X) to arrive at an alternative way of computation. The mandmizing parameters are
y +0 112
k 0
(0)
00 n
(13)
k l....,d d
where Wh = / Ei (xx)). Inserling
13
inlo
120
yields the following equivalent conver cost function
1 1b 0
arg min
lly 40
1 +
(14)
0
where vb [4/1 Val". Let On denote the minimizer of [14) then Tf(x; X,,) & (x)0n.
Eq. [4) is weighted square-root LASSO problem
Belloni et al. 2011 that can be solved in runtime on
the order of (dd) using variables
b" Σ Φ
(x:))s
d²
Vr
(15)
of fixed dimension that are updated recursively. Thus the memory requirement is dominated by the storing of the d x d-matriz A".
The pseudocode for the method is provided the appendix.
4.2 Out-of-sample performance
We now turn to evaluating the out-of-sample performance of the predictor learned from the data stream D,,. Specifically, we consider the mean-squared error
MSE [W (xx Am.))
(16)
for the subsequent sample (Xn+1, 11+1) ) in the stream. The expectation is conditional on Dn, thus the MSE will depend on the particular realization of the stream since the learned predictor is a function of all past samples.
To provide performance reference, note that all predictors of the form with 5 belong to the following class of predictor functions
d
F _2{ f(x) = Σ Φz (x) 0 E Rd
(17)
k=1
i.e. a linear combination of all features {xx(x)} in the nominal model [4] We can now benchmark (x: Xn) against the minimal achievable error among all predictors in F, even when the nominal covariance model is misspecified. Specifically, we provide the following finite out-of-sample performance guarantee for the learned predictor.
Theorem Assume the outcome y and features Φ6(x) bounded, and that the features such that (x,)) > 0. If the data pairs (xi, v.) in the stream are drawn ii.dd then the out- -aa-ole error of ((:: An) is given by
1 n E (Unt1 - (n))) < min [1 {xx+)))) + K
+ bn
(18)
yEF
with probability of at least e, where constant and by bounded as O(n 3344.. That is, with high probability, the out-of-sample error approaches the minimum achievable error root-n rate. Note that the number of features d increases only the second term at logarithmic rate.
5

---

Published in Transactions on Machine Learning Research (01/2023)
Proof 2 For notational simplicity, let (x.y) denote the random (n + 1)th sample. Let ÿ be any predictor in F and eapress ils oul-of-sumple mean-squree error os:
(())
((x) y
R (v- (oV-5T(x0)) (x)0)² E
0 E0,
(19)
1
where E E[zz') J. Similarly, the in-sample error can be expressed as R,, (T) 0 Σ9, where Σ n (=
Zn The gap between in- and out-of of sample errors can be bounded as:
( R(O)
0 Σ Σ)θ
d+1 d+1d+1
0.110 Eig - Egy
==
(20)
- (11011 + 1)2 max ESg Eigl
Nerl, we boand (see also Greenshtein & Rálov
2004 ) Since y and ((() are bounded random variables,
we have that [21zz] S for some and using Hoeffding's inequality
Egg 2 of S 2 2exp Pr P1
no²
(21)
2B2
Combining this result with the union bound over all (d 1)²
variables in T, we have that
1102
Pr (d+ 1) 2exp
(22)
2B²
Consequently, we can replace ō by
2
2(d + 1)2
B In
(23)
E
in 20 that
RR (7) R(O) S (1]01| + 1)20.
holds for any predictor in F with a probability of at least E. Thus if 0 bounded, ||0]] S P for some P, then
R(T) -P 5 (.) < R + (P+1))
(24)
holds with a probabilily of at least. 1 e.
Let us now study two specific predictor functions in F: An optimal predictor that minimizes the out-of- sample error (x) = o (x)0", where 0* = [l Elo(x)y] a bounded vector. The learned prealctoor y² (x: An) = o (x),,, where Qn a minimizer of 14 This vector is also bounded because the minómizer of (I4) coincides with that of
arg min Ily 20012 1
for some value of S < optimal and learned predictors, denoted by V" and T for brevity.
< 30. Thus both 100 11 and ||0 are bounded by some P and 24 applies the
Since On, minimizes the crilerion in il follows that
/ (T") + 1/2 21 0nlhi
V (y())*) 1/2 0" 0'
Vn

---

Published in Transactions on Machine Learning Research (01/2023)
After rearranging, we have
v (Trt) (00)--RRnnn (nn()*) S Sn (1 00 0 1 e 0,1l1)
In 1/2 @ 11 S 1/,
(25)
where B = l/blloc Multiplying both sides of the equality by the positive quantity (v (Q*) + +V (v))))), we have
(y( R, Rn S (v Ra - 3/2
(26)
S (2 Rn (yy) + 3P )n 1/2 3P,
where the second inequalily follows from using 25
Finally, by definition R(y") S R(T") and we have that
R(T") S (n,(T")) + (P
s
(2 Rn(y") + 1/3 )n 11
S + +( (2 R(y*) + (P+ 1)90 +n 1122 )n 3/2
(27)
= -Rly') + + + / In 2(d+1} + +
with probobility of at least E, where 24 was used in the first and third inequality and the second inequality.
26 was used in
4.3 Distributional robustness
In the previous section, we showed that the out-of-sample error fxx Xx) approaches the minimum achiev- able MSE at root-n rate. We will now that this predictor also provides robustness against distributional uncertainty for finite
The feature vector Φ(x) X + Rd for any predictor in FF leads to distribution of the random variables (Φ. y) which we denote p( y). The out-of-sample MSE can then be written
E (un+1 ((x++1)) Ep [0 o 0111
(28)
Using n i.dd. samples (o., we can define a predictor in F that minimizes the MSE under the least favourable distribution among all plausible distributions that are consistent with the data, Such predictor is called 'distributionally robust', see, e.g.. Duchi Namkoong (2018). To formalize set of plausible distributions, first define the empirical distribution
1 Σ
5(Φ 6i y 9:)
(29)
Then we consider set of distributions
{p: D(p,, p) cn),
(30)
where D(p,, p) is some divergence measure. A distributionally robust predictor minimizes the MSE under the least-lavourable distribution in the set 330,, viz.
max p: D(P....))..
Ep [1o+1 T 0)
(31)
Several different divergence measures D(p p) have been considered in the literature, including Kullback- Leibler divergence, chi-square divergence, and so on. One popular divergence measure is the Wasserstein distance Blanchet et al. 2019), which is defined as
D(p,, p) = inf Ez [c(b, y, $. &.y'),,
(32)

---

Published in Transactions on Machine Learning Research (01/2023)
where c, D'.8/) ) is nonnegative cost function and π is joint distribution over (Ø. y, ,',3') whose marginals equal ( and p(o',y') respectively. Thus D(pn. p) can be interpreted as measuring the expected cost of moving probability mass from one distribution to the other.
Theorem Suppose the observed features are standardized so that Ei=1 oi(),,
Then ((:; Xn)
corresponds to a predictor that minimizes the oul-of-sample MSE under the least favourable distribution in the set {p: Pp S n ?), defined by Wosserstein distance
32 with cost function
|1 I'I
y リ=У.
c( v. 6',4')
(33)
00
otherwise.
Thus the predictor is robust against distributional uncertainties in the features Φ, which may be high- dimensional. Note that the size of the distribution set shrinks with n.
Proof 3 When the features are standardized, then vb 1 and 14 becomes
arg min
n lly
vn 1
(34)
Using Theorem 1 in Blanchet et al. 2019 gence bound En n-
34
corresponds predictor that minimizes
37 with diver-
4.4 In-sample robustness
When learning ū(x; X,) it is possible that the observed covariates themselves are subject to errors so that the dataset is:
D. = {(x1, (81,,)),.... (Kn; Mn)}
Then the true feature vector Φ; = ((x;) can be viewed as a perturbed version of the observed vector Φ. =Xx),, where the perturbation 6; Φ; di is unknown. This problem (aka. errors-in-variables) leads to yet another interpretation of the predictor (xx;
Theorem Consider the bounded set of possible in-sample perturbations:
Sn = { [ 8n Ep.. [ <nn nnnn- E....
=k=
Then [(x; An) corresponds to predictor that minémizes the in na-sae root-MsE under the least-favourde perturbations in Sn:
V VB.. [(y- 6))
max ES,
(35)
where ÿ O
Proof 4 The problem in
35 can be written as:
8
max ESS
ly - (Ø + A)0l22 where A
(36)
Vn
6
Let [A]k denote the 1th column of the matrix A.
We can then upper bound the error as
d
d
max (5 Vn
y
$0
A 1
< max Es
- lly 3
[|
1 11
k=1
k=1
1 Ily
Φ0112
+ max
|| || 10,,
(37)
(8:)ES
vn
Ily
+ n
Ep..
22 02
nn
8

---

Published in Transactions on Machine Learning Research (01/2023)
where the bound is attainable when
40
A}s
V R.
|ly - hy-00/1
But the bound is of the same form as the cost function in (4). Thus solving problem 14 implies the mini- mization of (3).. Theorem 1 in Xu et al. (2009) established connection between a square-root marimization problem in the form of 36) and the penalized form on the right-hand side of (37).
5 Numerical Experiment
In the previous sections we have showed several computational and thoorctical propertics of the predictor function xx; )₂) which we shall call the SPICR-predictor. In this section presenl. a numerical experiment for sake of illustration.
5.1 Setup
We observe stream of samples generated by the following (unknown) process
x Uniform( 10]²),
(38)
y/x ~ GP(O, + o oxx()'))),
where
1(x,'') 0² (1+ v3
exp (-v3 |x
v3
x'
x
with noise variance a and scale / 7. other words, x is a two-dimensional covariate drawn from a uniform distribution and y is drawn from a Gaussian process (GP) with zero mean and Matérn covariance function. A realization of the process above and training data points are shown in Figures Ia and le
We consider class F with d = 100 periodic feaatures {x)(x)) using the Laplacian eigenfunction basis (oolin & Särkkä 2020). Note that this corresponds to a misspecified covariance model 4 We are interested in the online learning of predictors in F, and use the least-squares (Ls) and ridge regression methods as baseline references. Both methods can be implemented in an online fashion, but the latter requires fixing a regularization parameter. Here we simply set this parameter to 0.1 based on visual inspection.
5.2 Out-of-sampil performance
For illustration, consider the predictions produced by the Ls, ridge regression and SPICE methods, see Figure [] As expected, the Ls provides poor results at these sample sizes. Ridge regression with a fixed regularization parameter and SPICE with adaptively learned parameters appear to perform similarly here. To evaluate their out-ol-sample errors, we compare the MSE against that of the oracle predictor based on the (oracle) Gaussian process repression (GPR) predictor in (38 Table 1 that the out-of-sample error of SPICE is lower than that of Ls and ridge regression, and that the chosen class F is capable of predicting the GP in (38] well.
Following the discussion effective degrees of freedom dfa in Ruppert et al. 2003 we also provide a comparison between Ls, SPICE and the oracle GPR predictors in Figure 2 While Ls attains the maximum dfn at n 100, SPICE moderates growth rate a data-adaptive and online manner. The degrees of freedom of the oracle predictor increases gracefully and remains below its maximum value, even when n increases beyond d.
5.3 Run-time
We report the runtimes of RIDGE (since Ls is virtually identical) and SPICE, and consider well-specified GPR predictor with covariance parameters learned using the maximum likelihood method reference. As

---

Published in Transactions on Machine Learning Research (01/2023)
10000
35
9000
15
(a)
(b)
(c)
(d)
10
10
10
1
1.5
(e)
(I)
(g)
(h)
Figure 1: Contour plots. First column shows a realization of y in (38 along with sampling patterns {x; ) E X' for n 100 and n 250 (top and bottom rows, respectively). Second, third and fourth columns show the contour plots the SPICE-predicte ridge regression and the Ls-predictor, respectively. All three predictors belong to F.
MSE/ MSE
Ls
RIDGE
SPICE
50
4.38 x 10
1.71
1.11
100
21.12
1.47
1.09
250
1.47
1.19
1.06
500
1.11
1.06
1.02
Table 1: Mcan-square error (MSE) for Ls and SPICE methods, normalized by MSE* of an oracle predictor which is given the unknown covariance function in (35). For given set of training data Da, we compute the averaged squared error over 250 test points. The mean of this error the MSE and was approximated using 100 different realizations of Dn.
120
of: degrees of freedom
LS
100
GP
80
SPICE
60
40
20
0
50
100
150
200
of data points
Figure 2: Plot of degrees of freedom dfn, against number of data points for Ls, SPICE and oracle GP predictors.
can be seen in Table 2 the computational complexity GPR considcrably higher than that of the online alternatives. RIDGE has a lower run-time than SPICE but both perform as O(n). A visualization of the trends for 500 is given in Figure 3
10

---

Published in Transactions on Machine Learning Research (01/2023)
Run times (ms)
71
GPR
RIDGE
SPICE
50
59
0.15
6.0
100
113.8
0.16
11.6
250
231.3
0.20
29.2
500
1507.2
0.22
58.6
Table 2: Run times for oracle GPR, RIDGE and SPICE. realization.
Note that the run times are computed for single
10
Scace
Buu tinc (6s)
10
100
10
10"
R: 806
data points
Figure 3: Plot of run time of single run for learning parameters and prediction against. versus of data points for RIDGE, SPICE and oracle GP predictors.
5.4 Robustness
In Theorems and established two different types of robustness results with respect the features ((x). We therefore study the predictive performance when the training and test distributions over (Φ) y) diverge.
Specifically, the covariates in the training data are drawn X|x N(x, ofI), where follows (38 This
results in a distribution shift over the features (aka. irros-in-variables) Figure ! when data drawn from (38). We see that the out-of-sample error for RIDGE and SPICE increases consistently with Og, while Ls produces very poor and inconsistent results. SPICE is notably more robust against this dist ribut ional shift and these results corroborated the derived robustness properties.
evaluates the test MSE
6 Conclusion
We considered the problem of learning model-based linear smoothers online. If the model parameters were fixed, the resulting predictor which includes Ganssian process regression and kriging methods - can readily be computed sequentially. Since the model parameters unknown, however, they must be learned from data, using, c.g., maximum likelihood or vals-valliotion methods. But implementing them when the data. arrives stream requires either recomputing the predictor which is computationally prohibitive or resorting to approximations. In either case, these approaches do not offer clear-cut results the statistical properties of the resulting predictor.
We applied a covariance-fitting method to learn the model parameters, which was initially developed for spectral estimation. We first used its computational properties to show that the resulting predictor can computed soquentially. We then derived finite out-of-sample performance guarantecs of the resulting predic- tor and showed that its error approaches the minimun achievable level at root-n rate. Finally, we established connections to the distributional robustness literature by showing that the predictor is robust against distri-
11

---

Published in Transactions on Machine Learning Research (01/2023)
3.5
Spice
3
1.5
0.1
0.2
0.3
0.4
0.5
0.5
0.7
0.8
0.9
Figure 4: Out-of-sample error when the distribution of training data diverges from that of the test data. MSE (normalized by minimal MSE*) versus perturbation noise level dr on training data (n = 100). Evaluation based on 500 Monte Carlo runs.
butional uncertainties and errors the covariate training data. The performance, computational complexity and robustness of the proposed method were illustrated in a numerical experiment.
Acknowledgments We thank the anonymous reviewers from TMLR for their constructive feedback. This research was partially supported by the Swedish Research Council under contracts 2018-05040 and 2021- 05022, Wallenberg AI, Autonomous Systems and Software Program (WASP) funded by Knut and Alice Wallenberg Foundation, and by Kjell och Märta Beijer Foundation.
References
Alexandre Belloni, Victor Chernozhukov, and Lie Wang. Square-root lasso: pivotal recovery of sparse signals via conic programming. Biometrika, 98(4):791-806, 2011.
Christopher Bishop. Pattern recognition and machine learning. Springer, 2006.
Jose Blanchet, Yang Kang, and Karthyek Murthy, Robust wasserstein profile inference and applications to machine learning. Journal of Applied Probability, 56(3):830-857, 2019.
Thang D Bui, Cuong Nguyen, and Richard E Turner. Streaming sparse gaussian process approximations. Advances in Newral Information Processing Systems, 30, 2017.
T Tony Cai and Ming Yuan. Minimax and adaptive prediction for functional linear regression. Journal of the American Statistical Association, 101(119):1001-111,, 2012.
John Duchi and Hongseok Namkoong. Learning models with uniform performance via distributionally robust optimization. arXiv preprint 0rx1m:1800.0888 110.087000000 2018.
Eitan Greenshtein and Ya' Acov Ritov. Persistence in high-dimensional linear predictor selection and the virtue of overparametrization. Bernoulli, 10(6):971-988, 2004.
Trevor Hastie, Robert Tibshirani, and Jerome H Friedman. The elements of statistical learning: data mining, inference, and prediction, volume 2. Springer, 2009.
James Hensman, Nicolas Durrande, and Arno Solin. Variational fourier features for gaussian processes. Mach. Learn. Res., 18(151-1), 2017.
Ali Rahimi and Benjamin Recht. Random features for large-scale kernel machines. In NIPS, volume 3, pp.
5. Citeseer, 2007.
12

---

and Published in Transactions on Machine Learning Research (01/2023)
Carl Edward Rasmussen and Chris Williams. Gaussian processes for machine learning, volume 1. MIT press Cambridge, 2006.
David Ruppert, Matt P Wand, and Raymond Carroll. Semiporametric regression. Number 12. Cambridge university press, 2003.
Arno Solin and Simo Särkkä. Hilbert space methods for reduced-rank gaussian process regression. Statistics and Computing, 30(2):419 446, 2020.
Samuel Stanton, Wesley Maddox, lan Delbridge, and Andrew Gordon Wilson, Kernel interpolation for scalable online gaussian processes. In International Conference on Artificial Intellignnce and Stalistics, pp. 3133 3141. PMLR, 2021.
Michael L Stein. Interpolation of spotial data: some theory for kriging. Springer Science & Business Media, 2012.
Petre Stoica and Octavian Stanasila. Some spectral properties of the matrix b/a. Bul. Inst. Politehnic Bucuresti, ser. Electro., 44(3):3-8, 1982.
Petre Stoica, Prabhu Babu, and Jian Li. New method of sparse parameter estimation separable models and its use for spectral analysis of irregularly sampled data. IEEE Transactions on Signal Processing, 59 (1):35-47, 47, 2010a.
Petre Stoica, Prabhu Babu, and Jian Li. Spice: A sparse covariance-based estimation method for array processing, IEEE Transactions on Signal Processing, 59(2):629 638, 2010b.
Michalis Titsias. Variational learning of inducing variables in sparse gaussian processes. In Artificial intel- ligence and statistics, pp. 567 574. PMLR, 2009.
Larry Wasserman. All of nonparametric statistics. Springer Science Business Media, 2006.
Huan Xu, Constantine Caramanis, and Shic Mannor. Robust regression and lasso. In Advances in Neural Information Processing Systems, pp. 1801 1808, 2009.
Dave Zachariah and Petre Stoica. Online hyperparameter-fice sparse estimation mothod. IEEE Transactions on Signal Processing, 63113)33348 3359, 2015.
A Appendix: Sequential computation
Here we provide a pseudocode for computing [xx; An) = (x))0, sequentially. This is accomplished via a cyclic minimization of the convex problem [4
First, define the kth column as C₁ [D) and Ya = -yE Ejjk CzOj- Then the cost. function in 14 can be equivalently expressed
((+)) (15x c20x11311/2 +1 + Ck.
(39)
and minimized cyclically, one coordinate h = 1, ..., al time. It was shown in Zachariah & Stoica 2015) that the minimizer of (39h is given by
8kTk if Vn-Tr Tym > VarBa 7
0=
(40)
0,
else
where
04 = 1
=== -=>|c
7 = [c
(41)
= sign(c,
7kk
Cr - 2
(42)
B
Br
13

---

Published in Transactions on Machine Learning Research (01/2023)
The key observation here is that the variables 41
all expressible using the recursively computable quantities in 115 (along with prior parameter iterate 0' which is initialized at 0).
To arrive at this conclusion, define the following global variables
2 = - O0' = | g Φ T
(43)
which can be expressed using 15 of the global variables:
Then we have that Yx + ca l. so that we can express
41 in terms
ag + 0 + 2g 3x Akeo
19 + %b 86 sign + 4kx0)
(44)
Finally we can express the cyclic minimization approach in recursive form as outlined in Algorithm 1 initializing 0' = 0 at = 1. Additional computational considerations are developed in Sec. 111 of (Zachariah
& Stoica 2015 Link to implementation provided here: bttps://github com/Muhammad Osama
A Algorithm SPicE-predictor
1: Input: (Xa, Vn) and x (and initial 0*)
Update recursive variables A", b" and e" in 15
3: Set global variables go = ct + An0' 20' b and g = b" A"0'
repeat
- d
Compute 0g, 3r. Th and sy
7:
Compute 0g in (40
Update global variables So == 90 + Ahe (0/ 0x)2 + 2(0% 0k)g and g - + x 0%)
8:
Update 8% := 0g
10: until number of iterations equals T.
11:
Output: xxx Φ (x)0
14


Reviews:
[{"review_id": "GVsK9ecGc0", "reviewer": "Reviewer_tGLV", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\n1. The research problem is interesting.\n2. Theoretical analysis, if technically correct, would be very good contributions.\n\nWeaknesses:\n\n1. The proposed iteration procedures appear to be a straightforward application of weighted square-root LASSO. This raises concerns about the novelty of the paper.\n2. There are undefined notations. We find assumptions introduced in the proofs without being claimed in the statement of the theorem.\n3. Some technical details are missing, for example, there is no description of the main proposed iteration algorithm.\n4. Some writing problems as specified below.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper studies the regression problem where the output values of the sampling model have conditional covariance, conditioned on input vector. The covariance also depends on a training parameter (which confuses me a bit). The hypothesis space is designed to be all the weighted sums of training labels, with weights being functions of covariates. A convex optimization-based online algorithm is proposed and analyzed. The theoretical analysis provides computational complexity bound, out-of-sample error bound, and robustness analysis.", "claims_and_evidence_yn": "", "claims_explanation": "", "audience_yn": "", "audience_explanation": "", "requested_changes": "1. I feel confused about the claim (in the abstract) that MLE cannot be made online. I think that is easy. For example, just a google search gives me the following papers:\n\n* Surace, S.C. and Pfister, J.P., 2018. Online maximum-likelihood estimation of the parameters of partially observed diffusion processes. IEEE transactions on automatic control, 64(7), pp.2814-2829.\n\n* Rohde, D. and Cappé, O., 2011, June. Online maximum-likelihood estimation for latent factor models. In 2011 IEEE Statistical Signal Processing Workshop (SSP) (pp. 565-568). IEEE.\n\n* Malik, S. and Enzner, G., 2010, March. Online maximum-likelihood learning of time-varying dynamical models in block-frequency-domain. In 2010 IEEE International Conference on Acoustics, Speech and Signal Processing (pp. 3822-3825). IEEE.\n\nI do not mean to ask the authors to cite them, but just that the above claim is confusing.\n\n2. I am not certain about the claim at the end of page 2 that \"using data-dependent parameters λ in (4) does not readily provide any out-of-sample prediction performance guarantees.\" In particular, we see some works in the literature, and such data-dependent approaches are referred to as \"adaptive algorithms\". For example,\n\n* Caponnetto, Andrea; Yao, Yuan. Cross-validation-based adaptation for regularization operators in learning theory. Anal. Appl. (Singap.) 8 (2010), no. 2, 161--183.\n\n* Cai, T. Tony; Yuan, Ming. Minimax and adaptive prediction for functional linear regression. J. Amer. Statist. Assoc. 107 (2012), no. 499, 1201--1216.\n\n3. The norm notation in (8) is not defined.\n\n4. The method of computing the features in (3) is not specified. For this reason, some part of the proof of Theorem 1 needs to be re-assessed after the revision of the paper.\n\n5. I double the sentence \"It can also be shown that the minimizer of the concentrated criterion V(θ(λ),λ) equals cλ◦, that is a rescaling of the covariance-based parameters (8)\" between Equations (10) and (11). In particular, I think such equivalence may fail when the norm of the vector y does not equal 1. I suggest that the authors provide more details on this point.\n\n6. In (12), the factor 1/sqrt(n) of the second term on the right-hand side should be sqrt(n).\n\n7. Please cite the specific results in Zachariah & Stoica (2015) for the claimed computational complexity.\n\n8. Link \"toappear\" in page 3 is broken.\n\n9. I feel confused about (3), which claims that the data is not iid. Is there any target function that the training process is planned to recover? Is the vector lambda unknown? Is there any mechanism (for example, procedure, or probability distribution) on data generation? Is there any assumption on the input space script-X? For example, is script-X a metric space? Is it a topological space? Or is it just a compact subset of some Euclidean space? Is (2) an assumption, or the consequence of some other assumptions?\n\n10. Please provide details on the iteration procedure.\n\n11. Right after (13), we see the description of \"an unknown distribution p(x, y)\". Is this distribution also used to generate the training data Dn?\n\n12. Right after (15), the authors mentioned the probability \"at least 1 − ε\". What is the associated probability distribution? Note that we already see means in (15), so both sides of the inequality should already be constants.\n\n13. Right after (17), a new assumption that \"φ(x)\" is bounded, is introduced. Please include this assumption in the theorem.\n\n14. For the paragraph before (22): how is the \"optimal predictor y*\"  defined? Please provide more details to explain the claim \"Since z and Σ are bounded, the parameters are also bounded such that k...... ≤ P for some P .\" Note that the boundedness of a minimization problem coefficient does not guarantee the boundedness of its solution. For example, min(1/x + ax) over (0,inf) leads to unbounded solution x*=1/sqrt(a) when the constant \"a\" takes values in a bounded set (0,1).", "broader_impact_concerns": "This paper studies a machine learning algorithm. I have not any ethical concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "BAUJ8sCAax", "reviewer": "Reviewer_cTRW", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "\nStrengths: The theoretical results and proofs are presented in a clear way, and the setting considered in the paper is general enough to be interested in the audience of TMLR.\n\nWeaknesses and Requested changes:\n\n1. This paper lacks a more detailed discussion on related work. Such an online prediction problem should have been studied thoroughly in recent years, but most references in the current paper are before 2010.\n2. Lack of experiment justification: This paper advocates the proposed algorithm from three aspects, computational efficiency, out-of-sample generalization, and robustness. However, the experimental section only compares the out-of-sample error with simple baseline methods, including LS and ridge regression. More empirical results should be provided to demonstrate the proposed method's computational efficiency (run time comparison) and robustness (corresponding to the two notions considered in the paper).\n\nMinor comments:\n\n1. Typo 3rd line after equation (2): “This is includes”\n2. Last line on page 3, “is provided at TOAPPEAR”\n3. The paper can be more self-contained by adding remarks to discuss the results used from Zachariah & Stoica (2015), Blanchet et al. (2019), and Xu et al. (2009).", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "Summary of contributions： This paper proposes a covariance-fitting method to solve an online learning problem together with three types of theoretical guarantees, including computational complexity, excess risk bound, and robustness. Specifically, the proposed online learning method can be updated sequentially, which is more computationally efficient than existing techniques. Moreover, its out-of-sample error approaches the minimum achievable level at a root-n rate, where n is the number of data samples. Finally, it minimizes the out-of-sample error under the least favorable distribution within a given Wasserstein ball from the empirical distribution. It is also robust against errors in the covariate training data.\n\n", "claims_and_evidence_yn": "", "claims_explanation": "", "audience_yn": "", "audience_explanation": "", "requested_changes": "See above", "broader_impact_concerns": "There is no such concern for this theoretical paper.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "iMZEuXvV6o", "reviewer": "Reviewer_ZrgT", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strength\n\n- Perturbation-robustness guarantees in both test and training times are studied.\n\n- Theoretical findings are demonstrated via numerical experiments.\n\nWeakness\n\n- The main reason that the update can be performed sequentially is because the sufficient statistics needed for solving the least square problem can be expressed as summations over samples. Given this, I do not think that it is surprising to expect the predictor can be updated in an online manner. \n\n- Generalization bounds given in the paper are standard applications of basic probabilistic inequalities and has nothing to do with online aspect of the problem. \n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The submission considers a linear smoother predictor for streaming data. Specifically, data $(x_i, y_i)$ comes in an online fashion, and the learner wants to learn/update the prediction of an input $\\hat{y}(x)$ as observing a new sample without learning from the scratch (i.e., iterating again over all previous samples). The proposed estimator fits the model covariance matrix $C_\\lambda$ using the criterion known as SPICE and plug it to construct the predictor. Generalization and perturbation-robustness guarantees are given for the overall prediction performance. \n\nAs far as I can tell, not much thing is new in this paper: the algorithm is already some well-known one in literature, and the analysis is seemingly an straight-forward application of standard techniques. \n", "claims_and_evidence_yn": "", "claims_explanation": "", "audience_yn": "", "audience_explanation": "", "requested_changes": "The paper could have been more interesting if you also study non-i.i.d. streaming setup. For example, some data could have been corrupted (not just perturbed in a small amount), or underlying distribution can slowly change, etc. In its current setup, I do not see any particular challenge by considering an online setup. ", "broader_impact_concerns": "No specific broader impact concerns. ", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: While one of the reviewer complains about the lack of novelty in the paper, the other two found the derived results are worth publication. I follow their suggestions here. 
While the contributions are indeed limited, the paper is a nice extension and continuation of previous works by some of the authors that may be of interest to the machine learning community.
Decision: Accept (poster)