On the benefits of maximum likelihood estimation for Regression and Forecasting
Authors: Pranjal Awasthi
Abhimanyu Das
Rajat Sen
Ananda Theertha Suresh
Venue: ICLR 2022 Poster (2022)
Keywords: Forecasting
Time-Series
Regression
MLE

Abstract: We advocate for a practical Maximum Likelihood Estimation (MLE) approach towards designing loss functions for regression and forecasting, as an alternative to the typical approach of direct empirical risk minimization on a specific target metric. The MLE approach is better suited to capture inductive biases such as prior domain knowledge in datasets, and can output post-hoc estimators at inference time that can optimize different types of target metrics. We present theoretical results to demonstrate that our approach is competitive with any estimator for the target metric under some general conditions. In two example practical settings, Poisson and Pareto regression, we show that our competitive results can be used to prove that the MLE approach has better excess risk bounds than directly minimizing the target metric. We also demonstrate empirically that our method instantiated with a well-designed general purpose mixture likelihood family can obtain superior performance for a variety of tasks across time-series forecasting and regression datasets with different data distributions.

Published as a conference paper at ICLR 2022
ON THE BENEFITS OF MAXIMUM LIKELIHOOD ESTIMATION FOR REGRESSION AND FORECASTING
Pranjal Awasthi, Abhimanyu Das, Rajat Sen & Ananda Theertha Suresh Google Research
(pranjalawasthi, abhidas, senrajat, theertha) ) @goog com
ABSTRACT
We advocate for a practical Maximum Likelihood Estimation (MLE) approach towards designing loss functions for regression and forecasting, as an alternative to the typical approach of direct empirical risk minimization on a specific target metric. The MLE approach better suited to capture inductive biases such as prior domain knowledge in datasets, and can output post-hoc estimators at inference time that can optimize different types of target metrics. We present theoretical results to demonstrate that our approach is competitive with any estimator for the target metric under some general conditions. In two example practical settings, Poisson show that our competitive results can be used to prove
and Pareto regression, that the MLE approach has better excess risk bounds than directly minimizing the target metric. We also demonstrate empirically that our method instantiated with a well-designed general purpose mixture likelihood family can obtain superior performance for a variety of tasks across time-series forecasting and regression datasets with different data distributions.
INTRODUCTION
The task of fitting a regression model for a response variable y against a covariate vector x Rd is ubiquitous in supervised learning in both linear and non-linear settings (Lathuilière et al., 2020; Mohri et al., 2018) as well non-i.i.d settings like multi-varian forecasting (Salinas et al., 2020;
Wang et al., 2019). The end goal in regression and forecasting problems is often to use the resulting model to obtain good performance in terms of some target metric of interest on the population level (usually measured deviation are examples of common target metrics.
previously unseen test set). The mean-squared error the mean absolute
In this paper, our focus is on the choice of loss function that is used to train such models, which is an important question that is often overlooked, especially in the deep neural networks context where the emphasis is more the choice of network architecture (Lathuilière al., 2020).
Perhaps the most common method used by practitioners for choosing the loss function for a particular regression model is directly use the target metric of interest as the loss function for empirical risk minimization (ERM) over function class on the training set. We denote this approach for choosing a loss function as Target Metric Optimization (TMO). This especially more common with the advent of powerful general purpose function optimizers like deep networks and has also been rigorously analyzed for simpler function classes (Mohri et al., 2018).
Target Metric Optimization seems like a reasonable approach - the practitioner knows about the target metric of interest for prediction using the model, it seems intuitive to optimize for the same objective on the training data. Prior work (both theoretical and applied) has both advocated for and argued against TMO for regression problems. Many prominent works on regression (Goldberger et al., 1964; Lecue & Mendelson, 2016) use the TMO approach, though most of them assume that the data well behaved (e.g. sub-Gaussian noise). terms of applications, many recent works on time-series forecasting (Wu et al., 2020; Oreshkin et al., 2019; Sen et al., 2019) also use the TMO approach directly on the target metric. On the other hand, the robust regression literature has long advocated for not using the target metric directly for ERM in the case of contamination
heavy tailed response/covariate behaviour (Huber, 1992; Hsu & Sabato, 2016; Zhu & Zhou, 2021; Lugosi & Mendelson, 2019a; Audibert et al., 2011: Brownlees et al., 2015) account of its suboptimal

---

Published as a conference paper at ICLR 2022
High-probabiliin risk bounds. However, noted in (Prasad et al., 2020). many of these methods are either not practical (Lugosi & Mendelson, 2019a; Brownlees et al., 2015) or have sub-optimal empirical performance (Hsu & Sabato. 2016). Even more practical methods such as (Prasad et al., 2020) would lead to sufficiently more computational overhead over standard TMO.
Another well known approach for designing loss function is Maximum Likelihood Estimation (MLE). Here one assumes that the conditional distribution of given æ belongs to a family of distributions p[y|2; 2;; 0) parameterized by O (McCullagh & Nelder, 2019). Then one can choose the negative log likelihood as the loss function to optimize using the training set, to obtain an estimate mile. This approach is sometimes used in the forecasting literature (Salinas al., 2020; Davis & Wu, 2009) where the choice of a likelihood can encode prior knowledge about the data. For instance a negative binomial distribution can be used to model count data. During inference, given new instance x', one can output the statistic from p(y x mlee) that optimizes the target metric, as the prediction value (Gneiting, 201 1). MLE also seems like a reasonable approach for loss function design it is folklore that the MLE asymptotically optimal for parameter estimation, in terms of having the smallest asymptotic variance among all estimators (Heyde, 1978; Rao, 1963), when the likelihood is well-specified. However, much less known about finite-sample, fixed-dimension analysis of MLE, which is the typical regime of interest for the regression problems we consider in this paper. An important practical advantage for MLE is that model training is agnostic to the choice of the target metric- the same trained model can output estimators for different target metries inference time. Perhaps the biggest argument against the MIF approach is the requirement of knowing the likelihood distribution family. We address both these topics in Section 5.
Both TMO and MLE can be viewed as offering different approaches to selecting the loss function for a given regression model. In this paper, we argue that in several settings, both from practical and theoretical perspective, MLE might be better approach than TMO. This result might not be immediately obvious apriori while MLE does benefit from prior knowledge of the distribution, TMO also benefits from prior knowledge of the target metric at training time.
Our main contributions are as follows:
Competitiveness of MLE: In Section 3, we prove that under some general conditions on the family of distributions and a property of interest, the MLE approach is competitive with any estimator for the property. We show that this result can be applied to fixed design regression problems in order to prove that MLE can be competitive (up logarithmic terms) with any estimator in terms of excess square loss risk, under some assumptions.
Example Applications: In Section 1.1, we apply our general theorem to prove an excess square loss bound for the the MIE estimator for Poisson regression with the identity link (Nelder & Wedderburn, 1972; Lawless, 1987). We show that these bounds can be better than those of the TMO estimator, which in this case is least-square, regression. Then in Section 4.2, we show similar application in the context of Pareto regression ie y follows a Pareto distribution. We show that MLE can be competitive with robust estimators like the one in (Hsu & Sabato, 2016) and therefore can be better than TMO (least-squares).
We propose the of general purpose mixture likelihood family (see Section 5) that can capture a wide variety of prior knowledge across datasets, including zero-inflated or bursty data, count data, sub- Gaos-Gassian continuous data as well heavy tailed data, through different choices of (learnt) mixture weights. Then we empirically show that the MLE approach with this likelihood can outperform ERM for many different commonly used target metrics like WAPE, MAPE and RMSE (Hyndman & Koehler, 2006) for two popular forecasting and two regression datasets. Moreover the MLE approach is also shown have better probabilistic forecasts (measured by quantile losses (Wen et al., 2017)) than quantile regression (Koenker & Bassett Jr, 1978; Gasthaus et al., 2019; Wen et al., 2017) which is the TMO approach in this case.
Empirical Results:
PRIOR WORK ON MLE
Maximum likelihood estimators (MLE) have been studied extensively in statistics starting with the work of Wald (1949); Redner (1981), who showed the maximum likelihood estimates are asymptotically consistent for parametric families. Fahrmeir & Kaufmann (1985) showed asymptotic
2

---

Published as a conference paper at ICLR 2022
normality for MLE for generalized linear models. It is also known that under some regularity assumptions, MLE achieves the Cramer- Rao lower bound asymptotically (Lehmann & Casella, 2006). However, note that none of these asymptotic results directly yield finite sample guarantees.
Finite sample guarantees have been shown for certain problem scenarios. Geer van de Geer (2000); Zhang (2006) provided uniform convergence bounds in Hellinger distance for maximum likelihood estimation. These ideas were recently used by Foster & Krishnamurthy (2021) to provide algorithms for contextual bandits. There are other works which study MLE for non-parametric distributions e.g., Dümbgen & Rufibach (2009) showed convergence rates for log- log-concave distributions. There has been some works (Sur & Candès, 2019; Bean et al., 2013; Donoho & Montanari, 2016; El Karoui, 2018) that show that MLE can be sub-optimal for high dimensional regression i. when the dimension grows with the number of samples. this work we focus on the setting where the dimension does not scale with the number of samples.
Our MLE results differ from the above work as we provide finite sample competitiveness guarantees. Instead of showing that the maximum likelihood estimator converges in some distance metric, we show that under some mild assumptions it can work as well as other estimators. Hence, our methods are orthogonal to known well established results in statistics.
Perhaps the closest to our work is the competitiveness result of Acharya et al. (2017), who showed that MLE is competitive when the size of the output alphabet is bounded and applied to profile maximum likelihood estimation. In contrast, our work applies for unbounded output alphabets and can provide stronger guarantees in many scenarios.
3 COMPETITIVANESS OF MLE
In this section, we will show that under some reasonable assumptions on the likelihood family, the MLE is competitive with any estimator in terms of estimating any property of a distribution from the family. We will then show that this result can be applied to derive bounds the MLE in some fixed design regression settings that can be better than that of TMO. We will first setup some notation.
Given positive semi-definite symmetric matrix M, ||~|| == xx. M is the matrix norm of the vector x. X(M) denotes an eigenvalue of a symmetric square matrix M; specifically Amax (M) and Amin (M) denote the maximum and minimum eigenvalues respectively. The letter is used to denote general probability densities. We use p to denote the conditional probability density of the response given the covariate. |1-11, will be overloaded to denote the E norm between two probability distributions i.e ||p - p'll1 : Ip(=) (p(=))p'(<>) (zz) (pI P2) will be used denote the Ki--divergence between the two distributions. If 2 is a set equipped with a norm ||-|], then N(e, 2) will denote an e-net i. any point 2 E ZZ has a corresponding point z' (e,, Z) s.t || - 2'1l e. BC denotes the ball centered at the origin with radius r and S denotes its surface. We define
Notation:
[n] := 21, ,n.. |- denotes the size of the enclosed set.
General Competitiveness: We first consider a general family of distributions F over the space 2. For a sample 2 ~ f (for z E Z and f E F). the MLE distribution is defined as fz argmax f(z). We are interested in estimating property an observed sample. The following definition will be required to impose some joint conditions on the distribution family and the property being estimated, that are needed for our result.
F - W of these distributions from
Definition 1. The tuple (F,T), where F is a set of distributions and x : FF of those distributions, is said be (T,e, 82) apprrxxxxiee,eee, if there exists a set of distribu-
- W a property
tions F F s.t |FI S T and for every f E F, there exists f such that |lf - -fll1 S 51 and Pra~/ (1 (112)
2 < 82, where fa argmax f(z) and W has a norm I1-I|-
The above definition states that the set of distributions FF has a finite 6-cover, F in terms of the E distance. Moreover the cover is such that solving MLE the cover and applying the property on the result of the MLE is not too far from π applied on the MLE over the whole set F. This property satisfied trivially if F is finite. We note that it is also satisfied by some commonly used set of distributions and corresponding properties. Now we state our main result.
Theorem 1. Let à be an estimator such that for any f E F and z ~ f. (f( - #(2)]l 2 c) S Let Fy be subset of F that contains I such that with probability at least 1 - 8, fz E F, and (F,
3

---

Published as a conference paper at ICLR 2022
is (T,e,61.02) dap--p)xxxiii Then the MLE satisfies the following bound,
Pr([]=(f) =(f=)l 2 3e) - (T - 3)6 81 +o 82.
We provide the proof of Theorem in Appendix A. We also provide a simpler version of this result for finite distribution families as Theorem 3 in Appendix for the benefit of the reader.
Competitiveness in Fixed Design Regression: Theorem 1 can be used to show that MLE is competitive with respect to any estimator for square loss minimization in fixed design regression. We will first formally introduce the setting. Consider fixed design matrix X R where is the number of samples and d the feature dimension. We will work in setting where n > d. The target vector is random vector given by y" E R Let yi be the i-th coordinate of y" and æ; denote the i-th row of the design matrix. We assume that the target is generated from the conditional distribution given æ; such that,
31 p(, 0"), 0" E O.
We are interested in optimizing target metric ((.,)) -) given an instance of the random vector y". The final objective is to optimize,
1
min yu ~p(-122;00)) NEH
Uiu,, (ys, h(x,))
where H is class of functions. In this context, we interested in comparing two approaches. TMO (see (Mohri et al., 2018)). This is standard empirical risk minimization on the tar- get metric where given an instance of the random vector y" one outputs the estimator h. = miNĄCH Σ; fyy, h(x))).
MLE and post-hoc inference (see (Gneiting, 2011)). In this method one first solves for the parameter in the distribution family that best explains the empirical data by MLE i.e.
i mle argmin 0), where (ly":0) Σ log dogp(yx<,;;0) oce
Then during inference given sample the predictor is defined as T(::)) argming E (epp))) 0)] or other words we output the statistic from the MLE distribu- tion that optimizes the loss function of interest. For instance if is the square loss, then h(x;) will be the mean of the conditional distribution p(y pr; m.le).
We will prove general result using Theorem 1 when the target metric is the square loss and H is linear function class. Moreover, the true distribution p(- |x,; 0*) is such that E'y]] (0" æx)) for all i [n] ie we are the linear realizable setting.
In this case the quantity of interest is the excess square loss risk given by,
E()) = Σ E Il1" X01
Σ Ey 118"
x0" 1 = 1 1-01-- 101
(1)
where ; = -E, xx /n is the normalized covariance matrix, and 0 is the population minimizer of the target metric over the linear function class. Now we are ready to state the main result.
Theorem 2. Consider a fixed design regression setting where the likelihood family is parameterized
by eeR O C Bdd and (N(c, e B < [c B&) )I for a small enough c. Further the following conditions hold,
DKL (p(y; ply:: 0') < L 0'112-
The negative log- likelihood C(y*: 0) as a function of 0 is a-strongly convex and B-smooth, w.p. 5. -
at least
Further suppose there exists an estimator 0 such that {(0ma) < (ci + c2 log(1/6))"//o where C C2 are problem dependent quantities and 1) O. Then the MLE estimator also satisfies,
4

---

Published as a conference paper at ICLR 2022
((0mle) 0 (
cod(oon n loo LAmax (E)) + log(3/a) - log ())
n
w.p at least
We provide the proof in Appendix C. The proof involves proving the conditions in Theorem 1 and bounding the size of the cover
In order to better understand Theorem 2, let us consider a typical case where there exists possibly complicated estimator such that e0.ss = O((d + (og(1/6))/n)) implies that MLE will be competitive with this estimator up to log factor. In many cases the MLE might be much simpler to implement than the original estimator but would essentially match the same error bound. We now provide specific examples in subsequent sections.
this case the above theorem
4 APPLICATIONS OF COMPETITIVENESS RESULT
In this section we will specialize to two examples, Poisson regression and Pareto regression, where show that MLE can be better than TMO through the of our competitive result in Theorem 2.
4.1 POISSON REGRESSION
We work in the fixed design setting in Section 3 and assume that the conditional distribution of y/æ is Poisson i.e,
ply: k:ao,; 0*) u
where, f4; = ,0", > O,
(2)
k!
for all i E [n]. Poisson regression is popular model for studying count data regression which naturally appears many applications like demand forecasting (Lawless, 1987). Note that here we study the version of Poisson regression with the identity link function (Nelder & Wedderburn, 1972), while another popular variant is the one with exponential link function (McCullagh & Nelder, 2019). We choose the identity link function for fair comparison of the two approaches it is realizable for both the approaches under the linear function class i.e the globally optimal estimator in terms of population can be obtained by both approaches. The exponential link function would make the problem non-realizable under a linear function class for the TMO approach.
We make the following natural assumptions. Let Σ (E x;,:T)/ //n be the design covariance matrix before and M -Ei= pu " where u. 22:/| |a I|2- Let x and be the condition numbers of the matrices M and E respectively.
Assumption 1. The parameter space & and the design matrix X satisfy the following,
(AI) The parameter space θ {0 _R ||0112 S min (0, x0,z,)) 2 y> 0, vi [n]]. (A2) The design matrix is such that Xmiu (E) > 0 and ||>al| S for all [n].
(A3) Let Amin (M) 2
x² .n (d (olloo22) log(1/8)) and (Mm)
log (1)))) <
VnAmmoo (M) (M6/16, for a small 8 E (0, 1)
The above assumptions are fairly mild. For instance Amin is Ô(1/d) for random Gaussian covariance matrices (Bai& & Yin, 2008). The other part merely requires that Amin (M) = (u(((V(dan dd/aax (ddmmxxm)/// (M)))).
We interested in comparing MLE with TMO for the square loss which is just the least-squares estimator i.e 0 == argmingge 1 Ily" x0l| Note that apriori unclear as to which approach would be better in terms of the target metric because on one hand the MLE method knows the distribution family but on the other hand TMO explicitly geared towards minimizing the square loss during training.
n
Least squares analysis typically provided for regression under Gaussaa noise. By adapting existing techniques (Hsu et al., 2012), we show the following guarantee for Poisson regression with least square loss. We provide proof in Appendix D for completeness.
Note that the constants be further tightened in our analysis.

---

Published as a conference paper at ICLR 2022
Lemma 1. Let Pmaz w.p. at least 6,
max; P4z. The least squares estimator 0, satisfies the following loss bounds
(0 Hmaz (log 2 + d)) 10 (log 3 + +))²
if Umarx 2 (log(110) dloo6)/2
E(Om)=
otherwise
Now present our main result in this section which uses the competitiveness bound in Theorem coupled with the existence of a superior estimator compared to TMO, to show that the MLE estimator have better bound than TMO.
In Theorem (in Appendix F), under some mild assumptions on the covariates, we construct an estimator Bett with the following bound for the Poisson regression setting,
+ log(3)
E S |00 Amax (E)
The construction of the estimator is median- of of-means tournament based along the lines of (Lugosi & Mendelson, 201 30) and therefore the estimator might not be practical. However, this immediately gives the following bound on the MLE as a corollary of Theorem 2.
Corollary 1. Under assumption and the conditions of Theorem the MLE estimator for the Poisson regression setting satisfies w.p. at least
((mmle o (100) Amax (2)) (2)) dllogn + log(w (E)) + log X log !)
The bound in Corollary can be better than the bound for 0 in Lemma In the sub- sub-Gaussian region, the bound in Lemma scales linearly with Mmaa which can be prohibitively large even when few covariates have large norms. The bound for the MLE estimator Corollary has no such dependence. Further, in the sub-Exponential region the bound in Lemma scales (dd /n) while the bound in Corollary has (d/n)) dependency, up to log-factore In Appendix G, we also show that when the covariates are one-dimensional,
even sharper analysis is possible, that shows that the MLE estimator is always better than least squares in terms of excess risk. Appendix 1.8, we perform simulated experiment that shows linear growth of the test enror w.r.t Amax (2)), further showing the efficacy of our theoretical bounds.
4.2 PARETO REGRESSION
Now will provide example of heavy tailed regression setting where it is well-known that TMO for the square loss does not perform well (Lugosi & Mendelson, 2019a). We will consider the Pareto regression setting given by,
bm
b b-1 xx)) for yı 2 m
p(y kx = b+T m.
(4)
V
provided (0", a0:) > 7 for all i [m). Thus y: is Pareto given 2;; and Elys [as] - p4 0 æz)). We
will assume that > 4 such that 4 + e-moment exists for e >> 0. As in the Poisson setting, we choose this parametrization for fair comparison between TMO and MLE i. in the limit of infinite samples fr lies in the linear solution space for both TMO (least squares) and MLE.
As before, apply Theorem we need an estimator with a good risk bound. We use the estimator in Theorem of (Hsu & Sabato, 2016), which in the fixed design pareto regression setting yields,
(1 O d dog E(
||0012 2
6(b-2
n
Note that the above estimator might not be easily implementable, however this yields the following corollary of Theorem 2, which is bound the performance of the MIE estimator.
Corollary 2. Under assumptions of our Pareto regression setting, the MLE estimator satisfies w.p.
al least - 5,
d (1ogn + log ( + log bm RAmax + log
|10
9 (mmle) ) = 1 1+0 O
b(b 2)

---

Published as a conference paper at ICLR 2022
The proof is provided in Appendix H. It involves verifying the two conditions in Theorem 2 in the Pareto regression setting.
The above MLE guarantee is expected to be much better than what can be achieved by TMO which squarss It is well established in the literature (Hsu & Sabato, 2016; Lugosi & Mendelson, 2019a) that ERM on square loss cannot achieve (1log(1/6) dependency in heavy tailed regime; instead it can achieve only a O( 1/8) rate.
CHOICE OF LIKELIHOOD AND INFERENCE METHODS
In this section we discuss some practical considerations for MLE, such as adapting to a target metric of interest inference time, and the choice of the likelihood family.
Inference for different target metrics: In most practical settings, the trained regression model is used inference time predict the response variable some test set to minimize some target metric. For the MLE based approach, once the distribution parameters are learnt, this involves using an appropriate statistic of the learnt distribution inference time (see Section 3). For mean squared error and mean absolute error, the estimator corresponds to the mean and median of the distribution, but for several other commonly used loss metries in the forecasting domain such as Mean Absolute Percentage Error (MAPE) and Relative Error (RE) (Gneiting, 201 1; Hyndman & Koehler, 2006), this estimator corresponds to a median of transformed distribution (Gneiting, 2011). Please see Appendix for more details. This ability to optimize the estimator at inference time for different target metrics using a single trained model is another advantage that MLE based approaches have over TMO models that are trained individually for specific target metrics.
Mixture Likelihood: An important practical question when performing MLE-based regression is to decide which distribution family to use for the response variable. The goal is to pick distribution family that can capture the inductive biases present in the data. It is well known that a misspecified distribution family for MLE might adversely affect generalization error of regression models (Heagerty & Kurland, 2001). At the same time, it is also desirable for the distribution family to be generic enough cater to diverse datasets with potentially different types of inductive biases, or even datasets for which no distributional assumptions can be made in advance.
simple approach that observe work particularly well in practice with regression models using deep neural networks is to assume the response variable comes from a mixture distribution, where each mixture component belongs along with the parameters of the distribution. More specifically, we consider mixture distribution of k components p(y|æ; OI...., O2. Wi......... wa) Σニ= WyPy (y/ æ; 02), where each P5 characterizes different distribution family, and the mixture weights wj and distribution parameters O, are learnt together. We would like to have mixture distribution that can handle different situations like sparse data, sub-Exponential and sub-Gau Gaussian tails, count data as well as heavy tailed data. Moreover it should be applicable to continuous valued datasets in general.
different distribution family and the mixture weights are learnt
We use three component mixture of (i) the constant 0 izero-inflation for dealing with bi-modal sparse data), (it) a continuous version of negative binomial where n and p are learnt and (idi) a Pareto distribution where the scale parameter learnt. We provide more details about the continuous version of negative binomial distribution in Appendix 1.2. Our experiments in Section 6 show that this mixture shows promising performance on variety of datasets.
This will increase the number of parameters and the resulting likelihood might require non-convex optimization. However, we empirically observed that with sufficiently par-parameterized networks and gradient-based optimizers, this is not a problem in practice (Fig. 4 shows a convergence curve).
6 EMPIRICAL RESULTS
We present empirical results time-series forecasting problems and two regression problems using neural networks. We will first describe our models and baselines. Our goal is to compare the MLE approach with the TMO approach for three target metrics per dataset.
7

---

Published as a conference paper at ICLR 2022
Model
Favorita
MS
MAPE
WAPE
RMSE
MAPE
WAPE
RMSE
TMO(MSE)
0.6121+000077
0.2891-+0.0023
175.3782-+0.223
0.50045+0.000
0.2839+000008
7.507++0.023
TXO(MAE)
0.3983:10.0012
0.2258+10.0006
161.4999:0.4748
0.4452:+0.0005
0.266:10.0001
7.0503:10.009
IMO(MAPE)
0.3199::0.0011
0.2528:10.000
0.3892:1..000
0.3143 0007
11.3799 1965
IMO(Huber)
0.432-00.0033
0.2366++:00018
164. 7.07-00777
0.4722-40.0007
0.269++0.0000 0002
7.093-00.113
MLE (2NBP)
0.3139::..001
0.2238:10.0009
16456521+::55885
0.3064++000001
0.2677++0.0002
7.2133++.00122
Table 1:
We provide the MAPE, WAPE and RMSE metries for all models the set of two time-series datasets, The confidence intervals provided are one standard etrer over 50 experiments, for cach cuity. TMC(<1oss>) refers to TMO using the <loss>. For the MLE row, we only train one model per dataset. The same mudel used oulpul a diflerent statistic each columm during inference. MAPE, we output the optimizer of MAPE given in Section I.5. For WAPE we output the median and for RMSE we output the mean.
Model
Bicycle Share
Turbine
MAPE
WAPE
RMSE
MAPE
WAPE
RMSE
IMO(MSE)
0.2503:10.0008
0.1421:0.0000
878.5815:t 1.3059
0.8884:10.01 18
0.3496::0.0041
IMO(MAE)
0.25941100001
0.1436:10.0003
1011157:114945
0.774=0.0054
0.008910000119
1.5789=000067
IMO(MAPE)
0.2382:+000012
0.1469::0.0012
899.9163+448811
0.8108++0.0099
0.8189=0.00 0.000
3.0573=0000019
(Huber))))
0.2536:10.0011
0.141411.0004
889.1173:11.9655
0,902=0.0128
0.3598
1.5992 0.0082
MLE (2NBP)
0.1969:10.0088
0.1235+00001
767.4368:+7.1274
0.9877440.0019
0.3379:000004
1.4547=0.0054
Table 2:
We provide the MAPE, WAFE and RMSE metrics for all the models on the test set of two regression datasets. The confidence intervals provided we only train one model per dataset. The same model is used to output a different statistic for each column during inference. For MAPE, we output the optimizer MAPE given in Section L.5, For WAPE we oulgul median and for RMSE we oulgut the mean.
standard
50 experiments, for each entry. (MC((lss))) refers TMO using the <losss. For the MLE row.
Common Experimental Protocol: Now we describe the common experimental protocol on all the datasets (we get into dataset related specifics and architectures subsequently). For a fair comparison the architecture is kept the same for TMO and MLE approaches. For each dataset, we tune the hyper-parameters for the TMO(SEE) objective. Then these hyper-parameters are held fixed for all models for that dataset i.e only the output layer and the loss function is modified. We provide all the details in Appendix I. Our code will be available here.
For the MLE approach, the output layer of the models map to the MLE parameters of the mixture distribution introduced in Section 5, through link functions. The MLE output has 6 parameters, three for mixture weights, two for negative binomial component and one for the scale parameter in Pareto. The choice of the link functions and more details are specified in Appendix I.3. The loss function used is the negative log-likeling ligllikel implemented in Tensorflow Probability (Dillon et al., 2017). Note that for the MLE approach only one model trained per dataset and during inference we output the statistic that optimizes the target metric in question. We refer to our MLE based models that employ the mixture likelihood from Section 5 as MLE (2NBP) ) loss, where ZNBP refers to the mixture components Zero, Negative- Negative-Binomial and Pareto.
For TMO, the output layer of the models map to 9 and we directly minimize the target metric in question. Note that this means we need to train separate model for every target metric. Thus we have one model each for target metries in ( ('MSE', 'MAE', 'MAPE'). Further we also train a model using the Huber loss In order to keep the number of parameters the same that of MLE, we add an additional 6 neurons to the TMO models.
6.1 EXPERIMENTS ON FORECASTING DATASETS
We perform our experiments on two well-known forecasting datasets used in Kaggle competitions.
1. The M5 dataset (M5, 2020) consists of time series data of product sales from 10 Walmart stores in three US states. The data consists of two different hierarchies: the product hierarchy and store location hierarchy. For simplicity, in our experiments we use only the product hierarchy consisting of 3K individual time-series and 1.8K time steps.
2. The Favorita dataset (Favorita, 2017) is a similar dataset, consisting of time series data from Corporación Favorita, South-American grocery store chain. As above, we use the product hierarchy, consisting of 4.5k individual time-series and 1.7k time steps.
The task is predict the values for the last 14 days all at once. The preceding 14 days are used for validation. We provide more details about the dataset generation for reproducibility in Appendix I.
TThe Huber loss is commonly used in robust regression (Huber, 1992; Lugasi & Mendelson, 2019a)
8

---

Published as a conference paper at ICLR 2022
Model
MAPE
WAPE
Model
PlOQL
p90QL
MLLE (NB)
0.3314-/0.0016
0.2021+00000
175.501+/1.1226
TNO (Quatille
0.1073+0.0002
0.062--0.0019
ALE (2NB)
0.3056-/0.0011
0.24534/ 0.002
170.0075//1 1.22
MLE (2NBF)
0.0788
0.0536_0.0007
MLR (ZNRP)
0.3139 0.22WW-..O D0000000000 16468214/-1.51 51885
Table 3:
The MLLE model predicts the empirical quantile of inter-
Table 4:
We perform ablation study on Favorita dataset, where progressively add the componenis of mixture distribu- tion. There are three MI E. models the progression: Negative Bino- mial (NB), Zero-Inflated Negative Binomial (ZNB) and finally ZNBP.
during inference. It is compared with Quantile regression (TMO hased). The results, averaged over 50 runs along with the correspond- ing confidence intervals are presented.
The base architecture for the baselines as well as our model is a seq-2-seq model (Sutskever et al., 2014). The encoders and decoders both are LSTM cells (Hochreiter & Schmidhuber, 1997). The architecture is illustrated in Figure and described in more detail in Appendix I.
We present experimental results in Table 1. On both the datasets the MLE model with the appropriate inference-time estimator for a metric is always better than TMO trained on the same target metric, except for WAPE in M5 where MLE's performance is only marginally worse. Note that the MLE model is always the best or second best performing model on all metrics, among all TMO models. For TMO the best performance is not always achieved for the same target metric. For instance, (MO(MAE) performs better than TMO(MSE) for the RMSE metric on the Favorita dataset. In Table we perform an ablation study the Favorita dataset, where we progressively add mixture components resulting in three MLE models: Negative Binomial. Iefoltnf Negative Binomial and finally ZNBP. This shows that each of the components add value in this dataset.
6.2 EXPERIMENTS ON REGRESSION DATASETS
We perform our experiments on two standard regression datasets,
1. The Bicyle Share dataset (Bicycle, 2017) has daily counts of the total number of rental bikes. The features include time features as well as weather conditions such as temperature, humidity, windspeed etc. A random 10% of the dataset is used test and the rest for training and validation. The dataset has a total of 730 samples.
2. The Gas Turbine dataset (Kaya et al., 2019) has 11 sensor measurements per example (hourly) from gas turbine in Turkey. We consider the level of NOx as our target variable and the rest as predictors. There are 36733 samples in total. We use the official train/test split. A randomly chosen 20% of the training set is used for validation. The response variable is continuous.
For all our models, the model architecture is fully connected DNN with one hidden layer that has 32 neurons. Note that for categorical variables, the input is first passed through an embedding layer (one embedding layer per feature), that is jointly trained. We provide further details like the shape of the embedding layers in Appendix I. The architecture is illustrated in Figure 2.
We present results in Table 2.
On the Bicycle Share dataset, the MLE (ZNBP) based model performs optimally in all metrics and often outperforms the TMO models by a large margin even though TMO is separate model per target metric. On the Gas Turbine dataset, the MIE based model is optimal for WAPE and RMSE, however it does not perform that well for the MAPE metric.
In Table 3. we compare the MLE based approach versus quantile regression (TMO based) on the Bicycle Share dataset, where the metric presented is the normalized quantile loss (Wang et al., 2019). We train the TMO model for the corresponding quantile loss directly and the predictions are evaluated on normalized quantile losses as shown in the table. The MLE based model is trained by minimizing the negative log- likelihood and then during inference we output the corresponding empirical quantile from the predicted distribution. MLE (ZNBP) ) outperforms TMO(Quantile) significantly.
Discussion: We compare the approaches of direct ERM on the target metric (TMO) and MLE followed by post-hoc inference time optimization for regression and forecasting problems. We prove a general competitiveness result for the MLE approach and also show theoretically that it can be better than TMO in the Poisson and Pareto regression settings, Our empirical results show that our proposed general purpose likelihood function employed in the MLE approach can uniformly perform well several tasks across four datasets. Even though this addresses some of the concerns about choosing the correct likelihood for a dataset. some limitations still remain for example concerns about the non-convexity of the likeliiil.o We provide a more in-depth discussion in Appendix J.
9

---

Published as a conference paper at ICLR 2022
REFERENCES
Martín Abadi, Paul Barham, Jianmin Chen, Zhifeng Chen, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Geoffrey Irving, Michael Isard, et Tensorflow: system for large-scale machine learning. In 12th USENIX symposium on operating systems design and implementation (OSDI 16), pp. 265-283, 2016.
Jayadev Acharya, Hirakendu Das, Alon Orlitsky, and Ananda Theertha Suresh. A unified maximum likelihood approach for estimating symmetric properties of discrete distributions. In International Conference on Machine Learning, pp. 11-21. PMLR, 2017.
Jean- Yves Audibert, Olivier Catoni, al. Robust linear least squares regression. The Annals of Statistics, 39(5):2766-2744 2011.
Zhi-Dong Bai and Yong- fong-Qua Yin. Limit of the smallest eigenvalue of a large dimensional sample covariance matrix. In Advances In Statisties, pp. 108-1221 World Scientific, 2008
Derek Bean, Peter Bickel, Noureddine El Karoui, and Bin Yu. Optimal m-estimation in high- dimensional regression. Proceedings of the National Academy of Sciences, 11063)):4563-14560 2013.
Bicycle. Bicycle share dataset.
https://www. kaggle wwt.coo//coo///o cottactpratct/d/////ooo
sike-share-daily-eee 2017
Christian Brownlees, Emilien Joly, and Gábor Lugosi. Empirical risk minimization for heavy-tailed
losses. The Annals of Statistics, 43(6), Dec 2015. ISSN 0090-5364. doi: 10.10
URL http:t/: / /x dow.dooooooooo org. ... 1214 1115AA111130
Yunshun Chen, Aaron TL Lun, and Gordon K Smyth. From reads to genes to pathways: differential expression analysis rna-seq experiments using rsubread and the edger quasi-likelihood pipeline. F1000Research, 5, 2016.
Richard A Davis and Rongning Wu. A negative binomial model for time series of counts. Biometrika, 96(3):33--74, 2009.
Joshua Dillon, Ian Langmore, Dustin Tran, Eugene Brevdo, Srinivas Vasudevan, Dave Moore, Brian Patton, Alex Alemi, Matt Hoffman, and Rif A Saurous. Tensorflow distributions. arXiv preprint arXiv: 2017.
David Donoho and Andrea Montanari. High dimensional robust m-estimation: Asymptotic variance via approximate message passing. Probability Theory and Related Fields, 166(3):935-995)) 2016.
Lutz Dümbgen and Kaspar Rufibach. Maximum likelihood estimation of log-concave density and its distribution function: Basic properties and uniform consistency. Bernoulli, 15(1):00-68, 2009.
Noureddine Karoui. On the impact of predictor geometry on the performance on high- dimensional ridge-regularized generalized robust regression estimators. Probability Theory and Related Fields, 170 2018
Ludwig Fahrmeir and Heinz Kaufmann. Consistency and asymptotic normality of the maximum likelihood estimator in generalized linear models. The Annals of Statistics. 13(1):342-368, 1985.
Favorita. Favorita forecasting dataset. Eavoritta-grocer-saleee for.eessss.rrcceceneo.e.aao 2017.
https: http://wwwwwwg.cccooo/ com/c/
Dylan J Foster and Akshay Krishnamurthy. Efficient first-order contextual bandits: Prediction, allocation, and triangular discrimination. arXiv preprint arXiv:2107.02237, 2021.
Jan Gasthaus, Konstantinos Benidis, Yuyang Wang, Syama Sundar Rangapuram, David Salinas, Valentin Flunkert, and Tim Jamuschowski. Probabilistic forecasting with spline quantile function In The 22nd international conference on artificial intelligence and statistics, pp. 1901-1910. PMLR, 2019.
Sara A Geer and Sara van de Geer. Empirical Processes in M- M-estimation, volume 6. Cambridge university press, 2000.
10

---

Published as a conference paper at ICLR 2022
Tilmann Gneiting. Making and evaluating point forecasts. Journal of the American Staristical Association, 106(494):746-762, 2011.
Arthur Stanley Goldberger et al. Econometric theory. Econometric theory., 1964.
lan Goodfellow, Yoshua Bengio, Aaron Courville, and Yoshua Bengio. Deep learning, volume 1. MIT Press, 2016.
Patrick J Heagerty and Brenda F Kurland. Misspecified maximum likelihood estimates and gener- alised linear mixed models. Biometrika, 88(4):973-985, 2001.
C.C. Heyde. On an optimal asymptotic property of the maximum likelihood estimator of parameter from stochastic process. Stochastic Processes and their Applications, 8(1): 1-9, 1978. ISSN 0304-4149, URL httt::///ww..siennddeeeeee. cow/science//////ccccccccmccmccc
Sepp Hochreiter and Jürgen Schmidhuber. Long short-term memory. Neural compulation, 9(8): 1735-1780, 1997.
Daniel Hsu and Sivan Sabato. Loss minimization and parameter estimation with heavy tails. Journal of Machine Learning Research, 17(18):1-40, 2016. URL http: http://mmrrrrrrr..o wtw...o///ppprr////////////// papprs/v17/ 14-273. 14-273.html.
Daniel Hsu, Sham M Kakade, and Tong Zhang. Random design analysis of ridge regression. In Conference on learning theory, pp. 9-1. JMLR Workshop and Conference Proceedings, 2012.
Peter J Huber. Robust estimation of a location parameter. In Breakthroughs in statistics, pp. 492-518. Springer, 1992.
Rob John Hyndman and Ann B Koehler. Another look at measures of forecast accuracy. International Journal of Forecasting, 22(4):679-688, 2006.
Heysem Kaya, PINAR TÜFEKCÍ, and Erdine Uzun. Predicting co and no x emissions from gas turbines: novel data and benchmark pems. Turkish Journal of Electrical Engineering & Computer Sciences, 2766):77334 2019.
Bernhard Klar. Bounds on tail probabilities of discrete distributions. Probability in the Engineering and Informational Sciences, 14(2):161-77,, 2000.
Roger Koenker and Gilbert Bassett Jr. Regression quantiles. Econometrica: journal of the Economet- ric Society, pp. 33-50, 1978.
Stéphane Lathuilière, Pablo Mesejo, Xavier Alameda-Pinede and Radu Horaud. A comprehensive analysis of deep regression. IEEE Transactions on Puttern Analysis and Machine Intelligence, 42 (9)20055-081, 2020. doi: 10. 109/TPAMM.2219.991022
Jerald Lawless. Negative binomial and mixed poisson regression. The Canadian Journal of Statistics/La Revue Canadienne de Statistique, pp. 209-225, 1987.
G. Lecue and S. Mendelson. Learning subgaussian classes: Upper and minimax bounds. In S. Boucheron and N. Vayatis, editors, Topics in Learning Theory. Societe Mathematique de France, 2016.
Erich Lehmann and George Casella. Theory of point estimation. Springer Science Business Media, 2006.
Gábor Lugosi and Shahar Mendelson. Mean estimation and regression under heavy-tailed distribu- tions: A survey. Foundations of Computational Mathematics, 19(5): 145-11990 2019a.
Gábor Lugosi and Shahar Mendelson. Sub-gaussian estimators of the mean of random vector. The annals of statistics, 47(2):783-794, 2019h.
M5.
M5
forecasting dataset.
com/
m5 m Sos-foorcaatingg-e intpreccaatin-g-ccccraae 2020.

---

Published as a conference paper at ICLR 2022
Davis J McCarthy, Yunshun Chen, and Gordon K Smyth. Differential expression analysis of multifactor ma-seq experiments with respect to biological variation. Nucleic acids research, 40 (10):4288-429)) 2012.
Peter McCullagh and John A Nelder. Generalized linear models. Routledge, 2019.
Mehryar Mohri, Afshin Rostamizadeh, and Ameet Talwalkar. Foundations of machine learning. MIT press, 2018.
John Ashworth Nelder and Robert WM Wedderburn. Generalized linear models. Journal of the Royal Statistical Society: Series (General), 135(3)3370-384, 1972.
Ryan O'Donnell. Analysis of boolean functions. Cambridge University Press, 2014.
Boris N Oreshkin, Dmitri Carpov, Nicolas Chapados, and Yoshua Bengio. N-beats: Neural basis expansion analysis for interpretable time series forecasting. arXiv preprint arXi::1905.10437 2019.
Adarsh Prasad, Arun Sai Suggala, Sivaraman Balakrishnan, and Pradeep Ravikumar. Robust estima- tion via robust gradient estimation. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 82(3):601-622 2020.
C.R. Rao. Criteria of estimation in large samples. Sankhya, 25. Ser A, 1963.
Richard Redner. Note on the consistency of the maximum likelihood estimate for nonidentifiable distributions. The Annals of Statistics, pp. 225-228, 1981.
Philippe Rigollet. High-dimensional statistics.
http:://oww.mitt.ee///orrr Becture-notes/MMIBB. Lectueerrccee//M11111993311..2a. 2015.
Alesandro Rinaldo. Sub-exponential concentration.
http://www. wtwp:sstat..cmm cmu. edu/ pdÉ,
2019.
Mark D Robinson and Gordon K Smyth. Small-sample estimation of negative binomial dispersion, with applications to sage data. Biostatistics, 9(2):321-332, 2008
David Salinas, Valentin Flunkert, Jan Gasthaus, and Tim Januschowski. Deepar: Probabilistic forecasting with autoregressive recurrent networks. International Journal of Forecasting, 36(3): 1181-1191, 2020.
Rajat Sen, Hsiang-Fu Fu Yu, and Inderjit Dhillon. Think globally, act locally: A deep neural network approach high-dimensional time series forecasting. arXiv preprint vrXiv:1905.03800 2019.
Niranjan Srinivas, Andreas Krause, Sham M Kakade, and Matthias Seeger. Gaussian process opti- mization in the bandit setting: No regret and experimental design. arXiv preprint arXiv:0912.3995,
2009.
Pragya Sur and Emmanuel Candès. A modern logistic regression. Proceedings of the National Academy of Sciences, 116(2)): (66(2)):45166 2019.
theory for high-diimessioal
Ilya Sutskever, Oriol Vinyals, and Quoe Le. Sequence to sequence learning with neural networks. arXiv preprint arXiv: 1409.335, 2014.
Abraham Wald. Note on the consistency of the maximum likelihood estimate. The Annals of
Mathematical Statistics, 20(4):595-601, 1949.
Yuyang Wang. Alex Smola, Danielle Maddix, Jan Gasthaus, Dean Foster, and Tim Januschowski. Deep factors for forecasting. In International Conference Machine Learning, pp. 6607-6617. PMLR, 2019.
Ruofeng Wen Wen, Kari Torkkola, and Balakrishnan Narayanaswamy. multi-horizon quantile recurrent forecaster. In NIPS Time Series Workshop, 2017.
12

---

Published as a conference paper at ICLR 2022
Zonghan Wu, Shirui Pan, Guodong Long, Jing Jiang, Xiaojun Chang, and Chenggi Zhang. Connecting the dots: Multivariate time series forecasting with graph neural networks. In Proceedings of the 26th ACM SIGKDD International Conference Knowledge Discovery & Data Mining, pp. 753-763, 2020.
Tong Zhang. From epsilon-entropy to kl-entropy: Analysis of minimum information complexity density estimation. The Annals of Statistics, 34(5):2180-22) 2006.
Ziwei Zhu and Wenjing Zhou, Taming heavy-tailed features by shrinkage. International Conference on Artificial Intelligence and Statistics, pp. 3268-3276. PMLR, 2021.
13

---

Published as a conference paper at ICLR 2022
A PROOF FOR COMPETITIVENESS OF MLE
Theorem 3. Let ft be an estimator for property such that for any f E F and a ~ f. (r(
A<))ll 2 c) S Then the MLE based estimator satisfies the following bound,
Pr([I=(f) x (f=)ll 2 2x) - (1FI+ 1)ō.
Proof of Theorem 3. By triangle inequality and by the properties,
[T(()) - (S=)]l & [[² +())] + [1 - T
Hence,
1 + +(())>((+)122 I
We now take expectation of both LHS and RHS of the above equation with respect to the distribution f. For the LHS, observe that
D 701 Pr(II=(s)- - 2 2E).
For the first term in the RHS
E + - Pr([I=()- #2)] 2 e) S 6,
by the assumption in the theorem. Combining the above three equations yield that
Pr([I (f=))l 2 2e) = B[[ 5.
We now prove that E[1 ((f)) < FF
E -
(())
(a)
2 (A()){(R(E)
Σ (18))18(0)==(112>>
Σ J [(5)
Σ Pr (1 xf)ll 2
znf
Σ 6 1F
where (a) uses the fact that fz is the MLE estimate over set F and hence fz(2) 2 f().. (b) follows from Fubini's is theorem.
Now we prove an extension of our last theorem that can deal with infinite likelihood families.
Proof of Theorem 1. Let F be the cover of Fy that satisfies assumptions in Definition 1. By triangle inequality and by the properties of F. with probability least 82
[T((f) (ff2)l < (f) - (2)]l + [1* - (f.a)l + [((( T(fa)l < [ (ʃ) # (=)] + (=) (f=))l +
We have used the fact that fa E F, w.p at least 1 5. Hence with probability least 85 5,
|x())x( S ||<()>>()))2 1
14

---

Published as a conference paper at ICLR 2022
We now take expectation of both LHS and RHS of the above equation with respect to the distribution f. For the LHS, observe that
=|1 = (f [ff=) 3c).
For the first term in the RHS
- #()] 2 S
by the assumption in the theorem. Combining the above three equations yield that
(f( (f=)]l 2 3c) I E|1 [|*(2) ((f=))= 26 62.
We now prove that E (1())) - Tб + 51- Let f be the distribution in Fy that at most 61 away from f. Then,
[()) [(2) (1())))))))((0))10
(9)
f(+) (H6) + If - fll
/ (C)HMMM +
2
[ +
Σ f(:)) + 81 jeF
Σ (10)) +
je
Σ Pr (1[ - #(2) 2 E) +
jEF znf
6) + T8 51
jEF
where the (a) follows from the definition of & distance between distributions, (b) follows from the properties of the cover, (c) uses the fact that fz is the MLE estimate over set F and hence f=(2) 2 f(z). (d) follows from Fubini's 's Theorem.
B OTHER GENERAL RESULTS FOR THE MLE
In this section, prove that the log likelihood of the MLE distribution at the observed data point is close to the log-like likelihood the ground truth at the data point, upto an additive factor that depends on the Shtarkov sum of the family.
Lemma 2. Let z ~ f. The maximum likelihood estimator satisfies the following inequality with probability least 8,
1 8' log f()) log fa(z) log f(z) log E1((s) log
where ExEZ f2(2) is the Shtarkov sum of the family F. Furthermore, if F is finite, then
logf(2) f() < log fz(2) < log f() + log F + log 1 8
Further is the distribution family is finite, then
log f2) < log fz) s log fz) + log + log 8
15

---

Published as a conference paper at ICLR 2022
Proof. Let K = EEE2 fz(2)
fa(=)
=
21
Kf((2)
(1)
fe(z)
v IE
K Kf(=)
E 2EZ
(=))=
K
where (a) follows from Markov's inequality. The first part of the lemma follows by taking the logarithm and substituting the value of K. For the second part, observe that is finite then,
E1((z) ΣES03 f(z) = EEs() = E1=15. ZEZ
JEZ
JEF
FEF
We now prove result of Shtarkov's sum, which will be useful in other scenarios. Lemma 3. Let w (=) for some property x.
[v((w) E1(e)) Ed(l.)
Proof.
Et(()) Σ (=(=)(=)= Σ Σ fu(2)= = fu(w).
2:0(=)=w
C GENERAL FIXED DESIGN RESULT
In this section, we will prove Theorem 2 which is an application of Theorem to the fixed design regression setting with the square loss the target metrie. The theorem holds under some reasonable assumptions. We will fist prove some intermediate lemmas.
Lemma 4. Let (()0) = II p(- [z;; 0) in the fixed design setting. If 1 0'112 S then,
(If(:))) - -):::)) S V2nL5.
under the assumptions of Theorem 2.
Proof. By Pinskers' inequality we have the following chain,
[[ 0) f(:00)) S (2
2
(KL (1; p): 0)))
/ 2nL8
Lemma' 5. Assume the conditions of Theorem 2 holds. Let @mle == argmax Cly", 0), in the Poisson regression setting. Then with probability at least we have,
4w (E)e (E)
E((m.e) E(0mle)) s
n
16

---

Published as a conference paper at ICLR 2022
Proof. We assume that the event in (2) of Theorem 2 holds.
3
where 0 is the closest point in (V(, () to mle- Therefore, by definition
L(y"; mlee) z L(y"; + (B/2)²².
By virtue of the strong strong-convente convexity have,
2 3
mmle 0 mll
-2
2 a
Now by triangle inequality we have,
IV mlomt E((mll))
-110mm -00
0 0"
Amax (2)
0 mle
0 mlo
s
Now the fact Ja bl < 22mx
vb)l vbl and that O C Bd conclude,
(8(mmle) Elmteo) < 4w (Σ)ē
Proof of Theorem 2. Consider the net (ó, 0). I = 8² /(2nL). then by Lemma the net forms ô cover on F. Note that under the conditions of Theorem (W 0) S (3w/6)d.
Next note that in the application of Theorem 1, we should set
(c) + log(1 1/))""
Thus we set é in place of ē in Lemma we would need the following to apply Theorem 1,
Aw.Amax 12 B
(cr + 02 log(1/8)))
n
Thus, & can be such that,
log < max [log(2L),, 0.5 log a
3 + looo Amaxx (E)) + log(n))
From above we can apply Theorem with the above e, 8, dz and T (3w/0)² Therefore, we get
IP (0mm 2 (cr + log(l1////// log(1/5)))) < (T++5)
(5)
n
We can set 8 - 6/(T + 5) to then conclude that w.p at least 1 6,
E(Ômle) = =0( (cr c2d(log(w) + log(1 11 /d))))
Substituting the bound on log 1 yields the result.
17

---

Published as a conference paper at ICLR 2022
D LEAST SQUARES FOR POISSON
We begin by restating the following general statement about OLS.
Lemma 6 (see Theorem 2.2 in (Rigollet, 2015)). Suppose rank (X X) and & E R be an orthonormal basis for the column-span of X. Then we have the following result,
4
(() S sup
(u E)
(6)
ueSi
where o"c.. e. Here, E = yf" X0 for the given response sample y".
Now are position to prove Lemma 1.
Proof of Lemma J. We are interested in tail bounds on the RHS of (6). We will first analyze tail bounds for fixed S We have the following chain,
E es H Elexp e := Elexp(s((, e)}l
E| [lexp = B Blexp(se;(5) 14))
==1
where =0"" az)) Now bounding each term separately we have,
E(exp(su:() - p4))] exp( exp(-su:p;) exp(x; - 1))
have,
II exp(ps(exp(su:) 1)) exp
Σ (=(exp(sus) -1)
M4 (su;))
Σ
n
Σ Pmna j! (s0s))
1
exp
exp
1488 +
j!
exp
Σ M4SU
Mmnar (cxp(s)
1)
where we have used the fact that ||l||p S ||ll/2 for all p 2 2. Thus have,
(u H < exp(umaz(exp(s1 - 1)).
This means that the RV (u, é) is sub-exponential with parameter D² = 2/ and 0.56. Thus t < 4fmaa then we have that for a fixed u.
12
P(uu, E) 2 t) S exp
(8)
4f6mmx
Thus by union bound over the c-net with 2, have that wp 1 5,
sup ē) < 88pmaz (log(1/0)))))) log 6).
u
Thus we get the risk bound wp
32/maa (log(1/8) + log 6).
E E(( <
(9)
For the sub-exponential region when t 2 Mmaz we get the bound,
E(θ 16 max (1og(116) log 6)² p
We get the bound by setting d.
18

---

Published as a conference paper at ICLR 2022
E COMPETITIVENESS FOR POISSON REGRESSION
Lemma 7. Let f(;0) = II p(- , 0) where p is defined in Eq. (2). If 0 S à then,
IIf(=: -(::001l & R 2nw 8.
Proof. By Pinskers' inequality we have the following chain,
II[;: 0) f(::0))1 -
2 (KL ()
2 (0, æ;) (0,)) (og((0',x,))) ((0,z;) (0', x),))
Notice that the absolute value of the derivative of a log x wrt æ, is upper bounded by a/r 1 if x 7 >> and also 2 7. Therefore, following from above we have,
0 x2
11) (0, - (0',
0)((:)) - ((:0))11 S 2 2
1=1
0'xx;
2
11
110 0'1121
7
1=1
2n
110
10-0-
-
7
thus concluding the proof.
Now we prove high probability bounds on the strong convexity and smoothness of the likelihood in the context of Poisson regression.
Lemma 8. Let be the condition number of the matrix M, where u, xx/||<0,122 Then with probability at least 26, we have
Amin (V)C(y": 2 n
21101 Amin (M),
(10)
provided ndmmin (M) 2 x7 (d + log(1/0)) and Amax (M)(dlog(24)) + log((18)) <
/ndmiul (M)/16.
Proof. We start with the expression of the Hessian of C(y";
yi
yuuu
L()).
0 xx.)²
|1011
where u, = =
For fixed unit vector a E sd-1 let us define,
uou u u.u
L(u, 0)
0
2
youk
where 1; (ui, 0/
19

---

Published as a conference paper at ICLR 2022
Following the definition of sub-exponential RV in (Rinaldo, 2019),
Σ yar E SE V² Σ vt. a 0.56 max 2/
(11)
where v, 2/4uz. This implies that wp atleast
- -28/C we have,
C Σ y:us Σ Maul 2 Σ mrof log
8
if.
C Σ mse (mmax)))))) (
Σ fiio log
S 2(
Note that the above condition is mild and is satisfied when Amin (Ee panne z log(C/0)/47² as min, E, Heu; 2 Amin (E, u.u,u)) Here, C is a problem dependent constant that will be chosen later.
Now consider an c-net in & norm over the surface unit sphere denoted by Me,, d). well known that [W(e,, d)] S (3/e)*. Now any S sq-1 can be written as 2 x + a where æ E (N(e, d) and B Therefore, have that
L 0) L(x, Amax (L(0))e
A similar argument above gives us,
1 max L(x 0) + L(u, 0) xEV(e,d)
Amax (L(0)) max
L(u, 0) S
uEB
2 uuE
Amax ( (L()))
S2
max (EN(e.d)
L(x,0)
Thus by an union bound over the net we have wp 1 - 25 and other conditions,
min L(z 0) 2
1
Amia
Σ uuu.
2e
Amax
Σ Emaanit
$ESt
|10112
10112
Me u,) log
C
2(1 + 2c), WEB
max
Σ
101
if the conditions above hold with C = (3/e)d. We can now set 1/(8+ * x) where is the condition
number of the matrix E; uu uT. Now by virtue of that fact that / Amaxx (E u, log (og((/o))) <
Amin (E, (uurr ) /16 have the result.
Now we will prove result on the smoothness of the negative log likelihood for the Poisson regression setting.
Lemma 9. With probability at least - 28, we have
2nk2 Amax (M),
Amax (V) <
(12)
72
if Amiin (M) 2 (d 6 + log(1/8))/4n73
Proof. We start by writing out the Hessian again but this time upper bounding it in semi-definite ordering.
VOC
ya T
yi 72
y/
R uuu
: L(
O
x
72
1=1
1=1
20

---

Published as a conference paper at ICLR 2022
For fixed unit vector s let define,
yr R uu u.a f
L(u, 0) =
Σ yili
22
where U, R2 (u; /22. Proceeding as in Lemma 8 conclude that,
Σ
Euo yar SE 1² Σ v 0.56 max v/
(13)
where vs 22,,z This implies that wp atleast
28/C we have,
Σ you. - Σ Me0il s 2/ Σ Miot log 8
if.
Enoo Hiou log
(C)
2( Σ u (max)
Note that the above condition is mild and is satisfied when Amin (E, m,u R log(C//))442 min, M, 2 Amin (E: p.u,u?) Here, C' is problem dependent constant that will be chosen later.
Now consider an 1/2-net in & norm the surface unit sphere denoted by (0.55, d). It well known that [W(e, d)] S (3/e)d. Now any z E sQ-1 can be written x x+ where æ e (0.5,, and u E 1B.55 Therefore, we have that,
Amax (L(0)) = max Liu, 0) S (EN(eed)
max L(x, + max L(u, 0)
uEB x-1
uE3!
Amax (L()) 2 max L(x, 0)
# =eN(c.d)
Thus set C = 6d and obtain wp least 1 -20
Amacx (L(B)) S n R2
Amax (M),
72
provided 4v/ 4 Amaax (M) log(C/8) S Amaxx (M)/VT.
Proof of Corollary I. We show that the conditions of Theorem 2 hold.
Creating Nets: First we need to create an e-net over the parameter space θ We start by creating an 7/2R. Then remove all centers O. if Ei s.t (0c, 2D,) < 7/2. This is a valid c-net over θ as all net partitions that are removed do not have any points lying in O In subsequent section, we will always follow this strategy create e-nets over subsets of e.
6-net the sphere with radius Now suppose,
Strong Convexity and Smoothness: From Lemma 8 and 9 have that,
9 w
R2
72 X-
with probability 1 - O(8).
KL Divergence: The L in Theorem 2 is bounded by 2wR2 /y according Lemma
Combining the above into Theorem 2 and using the estimator in Section F we get our result.
21

---

Published as a conference paper at ICLR 2022
F MEDIAN OF MEANS ESTIMATOR FOR POISSON
In this section we will design a median of means estimator for the Poisson regression model based on the estimator proposed in the work of Lugosi & Mendelson (2019b). Recall that have fixed design matrix X E I xxd with rows 21,... Zn distribution with mean μ4 = - at. We will further assume that the design matrix is chosen from an (L, 4) hyper-contractive distribution. Mathematically this implies that for any unit vector w
and for each i E [n], ys is drawn from a Poisson
E((- x u)* S L- L-E((E=ax .u) 2)?
For simplicity we will assume that L = 0(1). In the above definition note that E is the empirical expectation over teh fixed design. This is benign assumption and for instance would be satisfied if the design matrix is drawn from sub-Gaussian distribution. In the general case, the obtained bounds will scale with L. We have the following guarantee associated with Algorithm
Algorithm 1: Median of Mcans Estimator
Input: Samples S {x1 8/ ....... (Zns yn)) confidence parameter 8.
Step 1: Compute Σ = E [x Form S' = {2ं = y/2 Σ ーー ii........ z'x = Σ x
Step 2: Randomly partition S' into blocks of size n n/k each where k = 20 [log( 11
Step 3: Feed in the k blocks to the median-of-nn means estimator of Lugosi & Mendelson (2019b)
get V.
Step Return 0 -= シー Av.
Theorem 4. There an absolute constant c such that with probability at least 6, Algorithm outputs such that
d + ( log(?))
0 S e. 110"112 Ammar Am)re(E)
(14)
n
Proof. The proof is exactly along the lines of the proof of Theorem I in the work of (Lugosi & Mendelson, 2019b) that we repeat here for the sake of completeness since the original theorem is not explicitly stated for a fixed design setting. To begin with notice that
E ElyE श
E x)Σ =
(15)
E .x)x)) x
(16)
(17)
Hence if & is the output of the median of the means estimator in Step 3 of Algorithm 1, then
the least squares error of 0 is exactly | 112 For convenience define şı' = Σ and E' E[='>
0 Σt. Exactly in Lugosi & Mendelson (2019b) our goal is to show that μí' beats any other vector v in the median of means tournament is far away from μ¹" To quantify this define
Tr(E')
Amax (E))
max (400
AV100
For fixed vector of length r, and block By, 1 beats if
2k ( H²) - v+r>0.
n
VE
Let us denote by oiy E 10. random variable representing whether data point i is in block j or not. By Chebychev's inequality we get that with probability least 9/10,
k
k
( p²)
V1O
Elož,((a' w²) v)2]
(18)
71 B,
n
k
10,
np
E[((x', w²) - u)
(19)
22

---

Published as a conference paper at ICLR 2022
Here is the probability of point belonging to block k. Noting that np ((n/k) get that with probability at least 9/10,
k
kAmax (E))
(ac p²)) 2 S v10r
(20)
iEB
Applying binomial tail estimates we get that with probability at least 1 - e-k1180 /e' beats v at least 8/ 10 of the blocks. By applying the covering argument verbatim as in the proof of Theorem in Lugosi & Mendelson (2019b) we get that with probability least 1 8, ú will satisfy
| = S Amax (E') - log(?)
Finally, it remains to bound the spectrum of E'. We have
'' Ely² xxx E -500000 Σt
(21)
z [[(( x=)² + 1.0x))) x 2 nxx -ー
(22)
Ti + T2:
(23)
where
T E[(0" x)) 2 - xx TE -
(24)
TB = E x) Σ x T
-ー $.
(25)
To bound Ti, T2 we note that for any function m(x) we have
( -ー [=mm²x))]T.
(26)
Using the above inequality and the fact that the design matrix (O(1), 4) hyper-contractive we get that
V E 2) 0(1 |1 Amax (Σ)
(27)
v F x) 0(110 >mmxx (E)).
(28)
Combining the above we get that
E' 1 Ti + T2 3 Mmxx (2))).
(29)
(30)
G 1-D POISSON REGRESSION
When the covariates are one dimensional, sharper analysis can actually be performed to show that the MLE dominates TMO in all regimes in the Poisson regression setting considered above.
Lemma 10. There exists an absolute constant probability least 8 that,
2 2 such that for any 8 E [ 1), it holds with
E(G..
1 >01 A
4 4.10
log
1
(31)
11
2i=1
The bound above is also light i.e with constant probability it holds that
8(0m)) = n "(2= (yi - 0,m)² 35=55 52( =2(10)
10 Σ
(32)
m
=
B(
We provide the proofs in later in the section. Having established the bound for least squares estimator, next prove the following upper bound on the mean squared emror achieved by the ML.E.
23

---

Published as a conference paper at ICLR 2022
Theorem 5. There exists absolute constant N 2 such that for any å E [ 1), it holds with probability at least 8 that,
E(0mle) 1 Σ (y) mmle)² 1=1
4-10"1
2 log
(33)
7t
5
[(0mme)
It is easy see that the covariate dependent term in the bound on the mean squared error achieved by mle (defined in Eq. 32) is always better the corresponding term in the bound achieved by 0i, (defined in Eq. 33). To see this notice that,
B((k) (E== 1 (z:)) B(0mle)
1.
(from Cauchy-Schhie Schwarz inequality.)
(E x
Furthermore, in many cases the bound achieved by dmte can be significantly better than the one achieved by l/.. As an example consider skewed data distribution where vn of the 2,'s 's take large value of /n, while the remaining data points take value of nt, where is small constant. In this case we have that,
B(((s) B(0mle)
= (E 2 = 2((n).
EE= x7)²
Proof of Lemma 10. Notice that
((0a) (0)s )²)² xf) n
(34)
Hence, it is enough to bound the parameter distance, i.e., (6,, 0²) In order to do that we first notice that you; a sub-exponential random variable with parameters v², a where where v² = 2/1,,7 [4 0 zi, and = 0.56x; (Rinaldo, 2019). In other words,
y: E SE(,,,
(35)
Thus from the bound on a sum of independent sub-exponential variables have that,
Σ yi E SE(v² = Σ
v² a 0.56max...)
(36)
Thus we have that,
12 exp
if L otherwise
P Σ yiuz Empll 2t 5
20²
2 eepp
+
This means, that w.p at least - 28,
1 7
Σ yoat -
Σ Uiatl 2
Σ m log
1
2 > !oyl
Σ ( log
log
88)
0.56 max; 34
The condition above is satisfied under our assumptions x; and 8, thereby leading to the bound
(6 0*) < log(1/)
Σ x
The bound the MSE claimed in the lemma then follows.
24

---

Published as a conference paper at ICLR 2022
Now prove the lower bound. Again it is enough to show a lower bound on |0,s 1 Notice that
1 (y: 14)0;
0 1
(37)
Ei=1 x2
Define the random variable 1 M.),y We will show anti-concentration for Z by
computing the cyper-contracttvityy of the random variable. Recall that random varibale Z is r--HC (hyper-contractive) if [E| 1E 22 Next we have
B|Z²]
y p:)ayl
p.at.
(38)
-
(y
14(1 + 3/t + Σ {
==1
igj
A.
(39)
Ei paa
Hence 7 is 1--HC with 1/2
From anti-concentration of hyper-contractive random variables (O' (O'Donnell, 2014) we have
P(|Z) 2 2 VE|z²1) 2 2(7*))
(40)
Hence we get that
P
p4
2 2
Σ paii
> 0(1)²)
Σ maat
(41)
A
(42)
Next notice that since 14; 2 y for all (Assumption 1), we have that 1 + 3p: S (3 + t)p, This
implies that A s (3 )) Eh=1 maa?) Hence if y is constant then with probability at least
3+T [1) we have
(E E((m) = (0(
(43)
Σ x7
Proof of Theurem 5. Recall that Onto is defined as
m
O argmin dee
Σ fx; - y, log().;).
(44)
Setting the gradient of the objective to zero, we get the following closed form expression for 0 m
mmle
y1
(45)
1= 22
25

---

Published as a conference paper at ICLR 2022
Next, note that Z Σ Y: is poission random random variable with parameter μ
Eisr 0g From tail bounds for Poisson random variables (Klar, 2000) we have that for any > 0,
Z d
2e
(46)
Taking c/H)log()) ) get that with probability at least - 6,
Σ v -h ul (2.2) log
(47)
provided that log()) < (that holds for our choice 6, once is large enough). Hence we conclude that with probability least
6, the mean squared error of mle bounded by
Σ 12 1= 0"
1mle 0
Ei=1
Iz
22
1 1og(3))
0
11 i= Tt.
10 1
=0
(48)
EE= Z4
The bound on the mean squared error follows from the above.
H COMPETITIVENESS FOR PARETO REGRESSION
We verify the conditions of Theorem for the Pareto regression setting. Lemma 11. Let (); I[ p(110;;; 0) where p is defined in Eq. (4). If 1 0'112 S 8 then,
26n5R
II f(:001l S $V
Proof. By Pinskers' inequality we have the following chain,
[|/(;; 0) f(: 0'1l S (2Dxx(((s0);(((000))
2
0): 0)))
i=1
2b log m: - log mm|
1=1
2b (0, ax))
(0' xxi)
?
i=1
26n 0
0' 11,2R
The second inequality follows from the fact that log x is Lipschitz with parameter i > 1/L. The last inequality follows from Cauchy-Schwan and the norm bound on x; 's.
Now we prove the rest of the conditions.
Lemma 12. have the following smoothness bound,
< (b 1)
Xmax X
(V2L)
Amax Σ).
26

---

Published as a conference paper at ICLR 2022
Proof. We start by writing out the Hessian,
V2L b 16-11000 a 20,000)0001
= L(0)
(0,
Using the fact that (0, x.) 2 7 gives the result.
Lemma 13. We have the following strong convexity bound,
Amin V2L)
Amia (E).
T R
Proof. It follows from the expression of the Hessian in Lemma 12 and using the fact (0, xx)) S |0| |
Proof of Corollary 2. We show that the conditions of Theorem 2 hold.
Creating Nets: First we need to create an e-net over the parameter space e We start by creating an e-net the sphere with radius w. Now suppose,
7/2R. Then remove all centers 0g if 3i s.t (Oc, æ.)) < 7/2. This is a valid e-net over θ as all net partitions that are removed do not have any points lying in O In subsequent section, will always follow this strategy create e-nets subsets of e.
Strong Convexity and Smoothness: From Lemma 13 and 12 we have that,
8
w R2
a 72 5.
KL Divergence: The L in Theorem 2 is bounded by 26R/7 according to Lemma 11.
Combining the above into Theorem 2 and using the estimator in (Hsu & Sabato, 2016) we get our result.
1 MORE ON EXPERIMENTS
We provide more experimental details in this section.
I.1 METRICS
The metrics and loss functions used are as follows:
MSE: The metric is
RMSE is just the square-root of this metric.
MAE: The metric
n
yil
71
WAPE: The metric
Σ
1
0:
yi
FL
1/%|
-=
MAPE: The metric is
1
Σ
72
Di
27

---

Published as a conference paper at ICLR 2022
Quantile Loss: The reported metrics in Table 3 are the normalized quantile losses defined as,
2p(us Dr 2. + 2p(0: Y4 Up)..
(49)
-i=1 1
During training the unnormalized version is used for quantile regression.
Huber Loss: The loss is given by.
( 0)3 8ly 5 2 otherwise.
if ly 5] S 8
Lo(y, 0)
I.2
MORE DETAILS ABOUT THE MIXTURE DISTRIBUTION
We use mixture distribution between zero, continuous extension of negative binomial (NB) distribution and Pareto. The continuous extension of negative binomial is such that the p.d.f at is proportional to,
r'n + k) T(K + 1)[(n)
p(v)
(1-")"
given parameters n and That is, the p.m.f of a regular discrete NB distribution is written in terms of Gamma functions and then we extend that to non-integral points up to proportionality, such that the measure sums to 1. This definition of NB is the standard implementation in Tensorflow (Abadi et al., 2016; Dillon et al., 2017), in order to support both discrete and continuous data has been used in many regression datasets before, especially in the field of genomics where the collected data a discrete continuous mixture (Robinson Smyth, 2008; Chen et al., 2016; McCarthy et al., 2012).
Why not discrete-contimann mixture of zero, discrete NB and Pareto?
mixture these three distributions a discrete continuous mixture, whose CDF is well-defined, is possible to define likelihood in standard manner such that it has component proportional to the pdf of Pareto at all points and to this we add dirac delta functions proportional to magnitude of the negative-binomial binomial pmf at non non-negative integral points. We have extra mass at for the zero component. The sum the dirac masses along with the integral of the Pareto density sums one. However, such likelihood is not very useful in practice. For instance, if the data is mostly continuous but not heavy tailed, the log-likelihood would mostly have the Pareto component (because non- non-integral points have no contribution from the other components) which is not desirable outcome, as NB distribution can better model sub-Gaus Ga-Gaussiin and sub-Exponential data in terms of moments.
account
1.3 MAPPING OF OUTPUTS FOR ZNBP
For the ZNBP model we require an output dimension size of 6. The first three dimensions are mapped through softmax layer (Goodfellow et al., 2016) to mixture weights. The fourth dimension is mapped the 'n' in begatiae-binomial likelihood through the link function,
[x+1, 1/(1 - x)
((x)
otherwise
The fifth dimension is mapped 'p' of the negative- binomial dimension is mapped to the scale parameter for the Pareto component using the ((x) link function above.
negative-binomial though the sigmoid function. The last
1.4 HARDWARE
We use the Tesla V100 architecture GPU for experiments. We use Intel Xeon Silver 4216 16-Core, 32- Thread, 2.1 GHz (3.2 GHz Turbo) CPU and our post-hoc inference for MLE is parallelized over all the cores.
28

---

Published as a conference paper at ICLR 2022
1.5 MORE DETAILS ABOUT INFERENCE FOR MLE
We follow the approach of monte-carlo sampling. For each inference sample æ', we generate 10k samples from the learnt distribution p( [x:; 0 mle)). Then we compute the correct statistics. For MSE, RMSE the statistic is just the mean and for WAPE, MAE is the median. For a quantile, is the corresponding quantile from the empirical distribution. For any loss of the form,
y
fly. ()
the optimal statistic is the median from the distribution proportional to y ( (Gneiting, 2011), Note that the MAPE falls under the above with -1 and relative error corresponds to 1. The statistic can be computed by importance weighing the empirical samples.
1.6 MODELS, HYPERPARAMETERS AND TUNING
MLE
TMO S+L
paramb.y
LSTe
LOTH
96-h
W4-bel
30-1
24-h
R6-6+1
St-1
Encoder
Decoder
Figure I: Time-Serics Seq-2-Seq models. The MLE config is shown on the left and the TMO config is shown at the right. The main difference is the output dimension and the loss function. In order to keep the number of parameters the same, in the TMO model we add an extra layer of size 6 with Relu activation (shown as Adj. (adjustment))
For the time-series datasets, the model is a seq-2-seq model with one hidden layer LSTMs for both the encoder and the decoder. The hidden layer size is h 256. The output layer of the LSTM is connected to the final output layer through one hidden layer also having f neurons. Note that f was tuned the set [8, 16. 32, 64, 128, 256] and for all datasets and models 256 was chosen. In order to be keep the number of parameters exactly the same, in the TMO models we add an extra layer with ReLU with neurons before the output.
We tuned the learning rate for Adam optimizer in the range [le-5, le- [e-1] in log-scale. The batch-size was also tuned in [64, 128. 256, 512] and the Huber-d in [2 for in range(-8, 8)]. The learning rate was eventually chosen as 2.77e-3 for both datasets, as it was close to the optimal values selected for all baseline models. The batch-size was chosen be 512 and the Huber-6 was and 64 for M5 and Favorita respectively.
For the regression datasets the model a DNN with one hidden layer of size 32. For the Bicycle dataset the categorical features had there own embedding layer. The features [season, month, weekday, weathersit] had embedding layer sizes [2,4,4.2].
We tuned the learning rate for Adam optimizer in the range [le-5, le- 1] in log-scale. The batch-size was also tuned in [64, 128, 256, 512] and the Huber-d in [2' for in range(-8, 8)]. The learning rate was eventually chosen 3e-3 for the gas turbine dataset based on the best perforamce of the baseline models and 1.98e-3 for the gas turbine dataset. The batch-size was chosen to be 512 and the Huber-6 was 128 and 32 for Bicyle Share and Gas turbine respectively.
For the ZNBP model we also tune the a parameter in the pareto component between [3, 4, 5). The value of selected for all datasets, except for gas turbine where used a
29

---

Published as a conference paper at ICLR 2022
MLE
TMO
Link Function
Linear
Relu
Linear layer
Relu
Linear Layer
X
Figure 2: Fully connected network for regression models. The MLE config is shown the left and the TMO config shown at the right. The main difference output dimension and the loss function.
Datasct
M5
1879
256
1653
Bicycle
584
32
Gis Turbine 22039
32
Table 5: Note that here d refers to the dimension of the last layer of the architecture used in the respective datasets.
We used batched version of GP-UCB (Srinivas al., 2009) tune the hyper-parameters, We used Tensorflow (Abadi et al., 2016) to train our models. In Table 5, we provide the number of samples and the dimension of the last layer in each of our datasets.
1.7 DATASETS
For the Favorita and M5 dataset we used the product hierarchy over the item-level time series, Along with the item-level (leaf) time-series, we also add all the higher-level (parent) time-series from the product hierarchy (i.e. we add family and class level time-series for Favorita, and department and category level time series for M5). The time-series for parent time-series is obtained as the mean of the time-series of its children. This is closer to real forecasting setting in practice where one is interested in all levels of the hicrarchy. The metrics reported are over all the time-series (both parents and leaves) treated equally. The history length for predictions is set to 28.
For the M5 dataset the validation scores are computed using the predictions from time steps 1886 to 1899, and test scores steps 1900 to 1913. For the Favorita dataset the validation scores computed using the predictions from time steps 1660 to 1673, and test scores on steps 1674 to 1687.
The train test splits are mentioned in the main paper. For the Gas turbine dataset use the official train test split. For the Bicycle share data there is no official split, but we use a randomly chosen fixed 10% the test set for all experiments.
30

---

Published as a conference paper at ICLR 2022
1.8 ADDITIONAL EXPERIMENTS AND FIGURES
In order to show the dependency of Amax (E) in the bound in Corollary we perform a simulated experiment. We generate dataset with n 5000 and d 10 such that each coordinate of is distributed i.i.d from a uniform distribution between [0, UJ- In out experiment, we vary in {1,2. ... 9). Then y genearted from Poisson distribution with rate (0* .x) for a fixed 0", This varies Amax (2:) which is equal to U²/12 We train and validate on 2500 samples with early stopping and plot the squared loss achieved the test by the mle based estimator in Figure 3 versus Amax (E). We can clearly see a linear relationship which further validates our theoretical results.
In Figure 4, plot the average training loss function of training iterations for the MLE (ZN3P) model. We can see that the loss converges to minima.
Test Square Loss
60
40
20
4
A max (X)
Figure 3: Test squared error versus Amax (E). can clearly see linear relationship. Each point in the plot is averaged over 10 runs and we plot the standard error bars.
Mean Training Loss
6
5
4
2
0
200
400
Iterations 100)
Figure 4: We plot the average training loss for the MLE (ZNBP) model as function of training iterations. We observe that the training curve converges.
J EXTENDED DISCUSSION AND LIMITATIONS
We advocate for MLE with suitably chosen likelihood family followed by post-hoc inference tuned to the target metric interest, in favor of ERM directly with the target metric. On the theory side, prove general result that shows competitiveness of MLE with any estimator for target metric under fairly general conditions. Application of the bound in the case MSE for Poisson regression and Pareto regression is shown. We believe that general result is of independent interest and can be used as a tool to prove competitiveness of MLE for a wide variety of problems. Such applications can be an interesting direction of future work.
On the empirical side we show that a well designed mixture likelihood like the one from Section can adapt quite well to different datasets the mixture weights trained. As have mentioned before,
31

---

Published as a conference paper at ICLR 2022
the MLE oog-likelhhod loss in such cases can be non-convex which might lead to some limitations in terms of optimization. However, we observed that this is usually not a problem in practice and the solutions that can be reached by mini- SGD can be quite good in terms of performance.
In conclusion we would recommend the following protocol for a practitioner based on our theoretical and empirical observations:
the overall problem is convex for TMO but introducing a MLE loss makes the problem non- convex, then the gains from the MLE approach might be neutralized by the added hardness of the non-convexity introduced. An example of such a situation is TMO for minimizing square loss on linear function class, which is just least-squares linear regression, but introduction of mixture likelihood like the one in Section 4 makes the problem non-convex. In this case it might be better to stick with TMO at least proceed with caution with the MLE approach. Note that if the chosen MLE retains the convexity of the problem, for example Poisson MLE in Section 1, then we would still recommend going with the MLE approach.
However, in many practical scenarios when training using a deep network, the TMO approach is non-convex to begin with, even when the target metric itself something simple and convex like the square loss. In such case we would recommend the MLE approach with a likelihood class that can capture inductive biases about the dataset. This is because both TMO and MLE are non-convex and better to capitalize the potential gains from the MLE approach.
Finally, note that the user can always choose between TMO and even between different likelihood classes through cross-validation in a practical setting. If the practitioner would like to forgo the decision making in choosing the likelihood class, we recommend using a versatile likelihood like the mixture likelihood in Sections 5.
We do not anticipate this work to have any negative social or ethical impact.
32


Reviews:
[{"review_id": "mWSbEDbEHf", "reviewer": "Reviewer_Fcze", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "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": "", "main_review": "**Strengths:**\n\n(a) The paper tackles an interesting and challenging topic to compare two well-known inferential methods for regression, namely, the maximum likelihood estimation (MLE) and the estimation based on loss functions (TMO).\n\n(b) A new theorem, Theorem 2, is given to compare the performance of the MLE with that of an estimator in terms of a quantity which measures the correctness of estimators of regression coefficients.\n\n(c) Applying Theorem 2, detailed discussion is given to compare MLE and TMO for Poisson regression and Pareto regression.\n\n**Weaknesses:**\n\n(d) I wonder the mathematical correctness of the discussion in Section 5 related to the three component mixture distribution which consists of the point distribution at 0, the negative binomial distribution and Pareto distribution. The point distribution and negative binomial distribution are discrete distributions, while Pareto distribution is a continuous distribution. Therefore the likelihood function of the mixture of these distributions should not be evaluated based on the function $\\sum_{j=1}^3 w_j p_j(y| \\boldsymbol{x} ; \\boldsymbol{\\theta}_j)$, where $p_1$ and $p_2$ are the probability mass functions of the point distribution and the negative binomial distribution, respectively, and $p_3$ is the probability density function of the Pareto distribution. I am worried from the discussion in Section 5 that this approach is used for the maximum likelihood estimation of the mixture distribution. There should be more explanation about the maximum likelihood estimation of this mixture model.\n\n(e) If my worry in Comment (e) is correct, part of the results of the experiments in Section 6 regarding the three component mixture should be reconsidered.\n\n(f) In the paper, MLE and TMO are compared through the quantity given in Equation (1). However I am not sure how meaningful this comparison will be because MLE and TMO are generally the minimizers of loss functions which are not based on the square loss risk as in (1). Since the two estimators are obtained through different loss functions, I think it is difficult to compare the goodness of the estimators through the single quantity such as the quantity (1). Or is there any special meaning to use the quantity (1) to evaluate the MLE and TMO?\n\n**Other Comments:**\n\n(g) I am confused about the MLE for regression models used in the paper. In the ordinary MLE for regression, the parameters of a probability distribution and regression coefficients are estimated simultaneously via the maximization of the likelihood function. However it seems from the paragraph titled \"MLE and post-hoc inference\" on page 4 that the parameters of the distribution $\\boldsymbol{\\theta}$ and the link function $\\tilde{h}$ are estimated in different steps. Is this approach used for Theorem 2 and other results of the paper? This should be explained in the paper because the approach discussed there is different from the common approach for MLE.\n\n(h) p.2, Section 1, Competitiveness of MLE, l.2 up: MLE can competitive  ===> MLE can be competitive\n\n(i) p.3, Section 3, Notation, ll.1-2 up: the sphere centered at ... ===> the ball centered at ... (I think this expression is less confusing.)\n\n(j) p.3, Definition 1, l.3: The definition of $| \\cdot |$ in $|\\tilde{\\cal F}| \\leq T$ should be given.\n\n(k) p.3, Theorem 1, l.2 up: Is \"the MLE based estimator\" exactly the same as the MLE or something different?\n\n(l) p.6, Equation (4):  As a comment related to Comment (g), how is the parameter $b$ in the Pareto regression model (4) estimated?\n\n(m) p.8, Section 6, Common Experimental Protocol, 2nd paragraph, l.2 up: that employs the mixture ===> that employ the mixture\n\n(n) p.28, l.7 up: be keep the number ===> keep the number", "summary_of_the_paper": "This paper compares two inferential methods for regression models, the maximum likelihood estimation and the estimation based on loss functions. For that purpose, a quantity is proposed to measure the correctness of an estimator of a regression model. Then it is shown that, under certain conditions on the quantity of an estimator, the proposed quantity of the maximum likelihood estimator can be evaluated. This result is applied to two regression models based on the Poisson distribution and Pareto distribution. Choices of the probability distributions and target metrics are discussed for the maximum likelihood estimation. Experiments are given to compare the performance of the maximum likelihood estimator and the estimators based on some loss functions.", "summary_of_the_review": "The paper tackles a challenging topic to compare two well-known inferential methods for regression. Some new theoretical results are presented on this topic. However I have concerns about the mathematical correctness of the paper (see Comments (d) and (e)) and the motivation for the study (see Comment (f)).", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works."}, "extra_text": {}}, {"review_id": "aCkK65KHdrt", "reviewer": "Reviewer_1pkV", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "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": "", "main_review": "Strengths:\n\n- It is well known that the MLE is efficient under very mild assumptions, so although less intuitive compared with the TMO, it is important to promote its usage. One disadvantage of MLE is the likelihood is different to determine and it is often misspecified in practice. The paper discusses this point and recommend using a mixture of distributions. \n- The finite sample error bound is obtained in a more general context than Acharya et al. (2017).\n- Specific examples (Poisson regression and Pareto regression) are provided to illustrate the theoretical results.\n\nWeaknesses:\n\n- The notations and presentation are difficult to follow. Some key notations do not seem to be properly defined. Here are some examples. What is the meaning of $\\ell_{1}$ norm for probability distributions? $f$ stands for the probability density function or cumulative distribution function or the underlying probability measure? What is the meaning of $|\\tilde{F}|$? From Definition 1, it seems to be the size of the set, but this does not seem to be the case when I check the proof. Due to these issues, I was not able to fully understand and appropriately evaluate the significance of the theoretical results, and I have reflect this in my confidence level.\n- The assumptions in Theorem 1 seems to already require that the MLE has a tight bound. I am not clear how to interpret the main result.\n- I am not sure if it is true that the proof for discrete distributions can be easily extended to continuous distributions. This is usually not the case in existing MLE theory partly due to fact that the support for continuous distributions are not countable.\n- Why both $\\ell_{1}$ norm and KL-divergence are required.\n- For the two distributions in section 4, it is well known that the MLE is more efficient than the least squares. \n- Minor issues:\n  - \"w.p\" should be \"w.p.\" in multiple places\n  - \"s.t\" should be \"s.t.\"\n  - The result in (Bai & Yin, 2008) for random covariance matrices is for Gaussian distributed covariates only. \n  - It might be better to give the definition of MAPE, WAPE and RMSE in the main paper than in the appendix so that readers can follow the main paper without going back and forth. \n", "summary_of_the_paper": "This paper promotes the maximum likelihood estimation over target metric optimization. The main theoretical results is a finite sample error bound. Specific examples and numerical results are provide to further illustrate the investigation.", "summary_of_the_review": "This paper in about an interesting and important idea, but the notations and presentation need some improvement. \n", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "SZL-7j70hED", "reviewer": "Reviewer_BdKr", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "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": "", "main_review": "The theoretical novelty of the paper lies in the fact that it is the first work that determines the performance of a MLE in terms of the performance of existing estimators using finite samples, as opposed to asymptomatic convergence to the true value. Subsequently, it derives tighter bounds, than existing results, for mle poisson and pareto regression. \n\n(1) The results seem to hold in low-dimensional cases where (n -- samples> d- covariates). The practical utility of the paper would be stronger if it could be applied in the high-dimensional regime ( for example in corollary 2, the d^2 factor seems to be warning in this scenario).\n\n(2) Similarly, it would be interesting if the authors could compare also with ridge/lasso estimators in section (6).\n\n(3) The bound of corollary (2) is the tightest that can be achieved (one cannot find a different \\theta_est satisfying the conditions of theorem (2) that could yield better excess square loss risk for the MLE?) \n\n(4) Regarding the bound of lemma 1, and given that I am not an expert in the area, I cannot grasp its tightness (how large is \\lambda_max (\\Sigma) ) in practical settings. I think it would help if the authors could show experimentally ( for experiments that do not make use of the mixture of distributions suggested- but only the poisson and/or pareto distribution separately ) and demonstrate empirically the tightness of the bounds.\n\n(5) Moreover, I think it would be useful if the authors could summarize in a Table n and d for the datasets considered.\n\n(6) Given that the resulting loss  that stems from the proposed mixture is con-convex,  it would be useful if the authors could demonstrate convergence empirically by providing some learning curves in the appendix.\n\n(7) Finally, it is not clear to me why theorem 2 can be applied on the mixture (even when it holds separately for each component).", "summary_of_the_paper": "The paper provides theoretical results that favor MLE estimators, in terms of the excess square loss risk,  compared to empirical risk estimators under mild assumption. In particular, the paper devises an estimator for Poisson regression and employs an existing estimator for heavy-tail Pareto regression, and subsequently derives upper bounds on the excess square loss risk that are tighter compared to the least squares estimator. Finally, a mixture of distributions (with each component coming from a different distribution class) is proposed. The MLE for the proposed mixture outperforms empirical risk estimators under different objective losses. ", "summary_of_the_review": "The paper makes theoretical contributions on the performance of mle compared to empirical risk losses. \nHowever, the applicability of the derived results in real-world regimes could be improved. ", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper has been independently assessed by three expert reviewers. The results place it at the borderline of acceptance decision: while one of the reviewers gave it a straight accept evaluation, two others assessed it as marginally rejectable, even after discussion with the authors. All of the reviewers agreed that the theoretical results provided should help promote the use of MLE estimators over perhaps more prevalently used in current practice TMO, and that is the main contribution of this work. The reviewers were concerned with the clarity of the presentation and with a confusing notation used. Some of these issues have been addressed in the authors' responses. All things considered, I conclude that this work can be of some interest to the ICLR audience, and as such it can be assessed as marginally acceptable for this conference: "accept if needed". I will recommend it as such for consideration by the Senior Area Chair and the Program Committee.
Decision: Accept (Poster)