Year: 2021

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.

heading: 1 INTRODUCTION
text: 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 as non-i.i.d settings like multi-variate 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 on a previously unseen test set). The mean-squared error or the mean absolute deviation are examples of common target metrics.
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 on the choice of network architecture (Lathuilière et al., 2020).
Perhaps the most common method used by practitioners for choosing the loss function for a particular regression model is to directly use the target metric of interest as the loss function for empirical risk minimization (ERM) over a function class on the training set. We denote this approach for choosing a loss function as Target Metric Optimization (TMO). This is 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 - if 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 is well behaved (e.g. sub-Gaussian noise). In 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 or heavy tailed response/covariate behaviour (Huber, 1992; Hsu & Sabato, 2016; Zhu & Zhou, 2021; Lugosi & Mendelson, 2019a; Audibert et al., 2011; Brownlees et al., 2015) on account of its suboptimal high-probability risk bounds. However, as 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 a loss function is Maximum Likelihood Estimation (MLE). Here one assumes that the conditional distribution of y given x belongs to a family of distributions p(y|x;θ) parameterized by θ ∈ Θ (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 θ̂mle. This approach is sometimes used in the forecasting literature (Salinas et 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 a new instance x′, one can output the statistic from p(y|x′; θ̂mle) that optimizes the target metric, as the prediction value (Gneiting, 2011). MLE also seems like a reasonable approach for loss function design - it is folklore that the MLE is 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 is 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 metrics at inference time. Perhaps the biggest argument against the MLE 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 a practical and theoretical perspective, MLE might be a 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 to logarithmic terms) with any estimator in terms of excess square loss risk, under some assumptions.
Example Applications: In Section 4.1, we apply our general theorem to prove an excess square loss bound for the the MLE 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-squares regression. Then in Section 4.2, we show a similar application in the context of Pareto regression i.e y|x 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).
Empirical Results: We propose the use of a 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-Gaussian continuous data as well as 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 to 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.

heading: 2 PRIOR WORK ON MLE
text: 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 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, we 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-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.e when the dimension grows with the number of samples. In 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.

heading: 3 COMPETITIVENESS OF MLE
text: 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 on the MLE in some fixed design regression settings that can be better than that of TMO. We will first setup some notation.
Notation: Given a positive semi-definite symmetric matrix M , ‖x‖M := xTMx is the matrix norm of the vector x. λ(M) denotes an eigenvalue of a symmetric square matrix M ; specifically λmax(M) and λmin(M) denote the maximum and minimum eigenvalues respectively. The letter f is used to denote general probability densities. We use p to denote the conditional probability density of the response given the covariate. ‖·‖1 will be overloaded to denote the `1 norm between two probability distributions i.e ‖p− p′‖1 := ∫ |p(z)− p′(z)|dz. DKL(p1; p2) will be used to denote the KL-divergence between the two distributions. If Z is a set equipped with a norm ‖·‖, then N ( ,Z) will denote an -net i.e any point z ∈ Z has a corresponding point z′ ∈ N ( ,Z) s.t ‖z − z′‖ ≤ . Bdr denotes the ball centered at the origin with radius r and Sd−1r denotes its surface. We define [n] := {1, 2, · · · , n}. | · | denotes the size of the enclosed set. General Competitiveness: We first consider a general family of distributions F over the space Z . For a sample z ∼ f (for z ∈ Z and f ∈ F), the MLE distribution is defined as fz = argmaxf∈F f(z). We are interested in estimating a property π : F → W of these distributions from 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.
Definition 1. The tuple (F , π), where F is a set of distributions and π : F → W a property of those distributions, is said to be (T, , δ1, δ2)-approximable, if there exists a set of distributions F̃ ⊆ F s.t |F̃ | ≤ T and for every f ∈ F , there exists a f̃ such that ‖f − f̃‖1 ≤ δ1 and Prz∼f
(∥∥∥π(fz)− π(f̃z)∥∥∥ ≥ ) ≤ δ2, where f̃z = argmaxf̃∈F̃ f̃(z) andW has a norm ‖·‖. The above definition states that the set of distributions F has a finite δ-cover, F̃ in terms of the `1 distance. Moreover the cover is such that solving MLE on 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 is 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 ∈ F and z ∼ f , Pr(‖π(f)− π̂(z)‖ ≥ ) ≤ δ. Let Ff be a subset of F that contains f such that with probability at least 1− δ, fz ∈ Ff and (Ff , π)
is (T, , δ1, δ2)-approximable. Then the MLE satisfies the following bound,
Pr(‖π(f)− π(fz)‖ ≥ 3 ) ≤ (T + 3)δ + δ1 + δ2.
We provide the proof of Theorem 1 in Appendix A. We also provide a simpler version of this result for finite distribution families as Theorem 3 in Appendix A 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 a fixed design matrix X ∈ Rn×d where n is the number of samples and d the feature dimension. We will work in a setting where n d. The target vector is a random vector given by yn ∈ Rn. Let yi be the i-th coordinate of yn and xi denote the i-th row of the design matrix. We assume that the target is generated from the conditional distribution given xi such that,
yi ∼ p(·|xi;θ∗), θ∗ ∈ Θ.
We are interested in optimizing a target metric `(·, ·) given an instance of the random vector yn. The final objective is to optimize,
min h∈H Eyi∼p(·|xi;θ∗)
[ 1
n n∑ i=1 `(yi, h(xi))
] ,
whereH is a class of functions. In this context, we are interested in comparing two approaches. TMO (see (Mohri et al., 2018)). This is standard empirical risk minimization on the target metric where given an instance of the random vector yn one outputs the estimator ĥ = minh∈H 1 n ∑n i=1 `(yi, h(xi)).
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.,
θ̂mle := argmin θ∈Θ L(yn; θ), where L(yn; θ) := n∑ i=1 − log p(yi|xi;θ)
Then during inference given a sample xi the predictor is defined as, h̃(xi) := argminŷ Ey∈p(·|xi;θ̂mle)[`(y, ŷ)] or in other words we output the statistic from the MLE distribution that optimizes the loss function of interest. For instance if ` is the square loss, then h̃(xi) will be the mean of the conditional distribution p(y|xi; θ̂mle). We will prove a general result using Theorem 1 when the target metric ` is the square loss andH is a linear function class. Moreover, the true distribution p(·|xi;θ∗) is such that E[yi] = 〈θ∗,xi〉 for all i ∈ [n] i.e we are in the linear realizable setting. In this case the quantity of interest is the excess square loss risk given by,
E(θ) := 1 n n∑ i=1 Eyn‖yn −Xθ‖22 − 1 n n∑ i=1 Eyn‖yn −Xθ∗‖22 = ‖θ − θ ∗‖2Σ, (1)
where Σ := ( ∑ i xix T i )/n is the normalized covariance matrix, and θ
∗ 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 θ ∈ Θ ⊆ Bdw and |N ( ,Θ ∩ Bdw)| ≤ |N ( ,Bdw)| for a small enough . Further the following conditions hold,
1. DKL(p(yi;θ); p(yi;θ′)) ≤ L‖θ − θ′‖2. 2. The negative log-likelihood L(yn;θ) as a function of θ is α-strongly convex and β-smooth, w.p.
at least 1− δ. Further suppose there exists an estimator θest such that E(θest) ≤ (c1 + c2 log(1/δ))η/n, where c1, c2 are problem dependent quantities and η > 0. Then the MLE estimator also satisfies,
E(θ̂mle) = O
(( c1 + c2d ( log n+ log(wLλmax(Σ)) + log(β/α) + log 1 δ ))η n ) w.p at least 1− δ.
We provide the proof in Appendix C. The proof involves proving the conditions in Theorem 1 and bounding the size of the cover T .
In order to better understand Theorem 2, let us consider a typical case where there exists a possibly complicated estimator such that E(θest) = O((d + log(1/δ))/n). In this case the above theorem implies that MLE will be competitive with this estimator up to a log n 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.

heading: 4 APPLICATIONS OF COMPETITIVENESS RESULT
text: In this section we will specialize to two examples, Poisson regression and Pareto regression, where we show that MLE can be better than TMO through the use of our competitive result in Theorem 2.

heading: 4.1 POISSON REGRESSION
text: We work in the fixed design setting in Section 3 and assume that the conditional distribution of y|x is Poisson i.e,
p(yi = k|xi;θ∗) = µki e −µi
k! where, µi = 〈θ∗,xi〉 > 0, (2)
for all i ∈ [n]. Poisson regression is a popular model for studying count data regression which naturally appears in 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 a fair comparison of the two approaches as 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 Σ = ( ∑n i=1 xix T i )/n be the design covariance
matrix as before andM = ( ∑n i=1 µiuiu T i )/n, where ui = xi/‖xi‖2. Let χ and ζ be the condition numbers of the matricesM and Σ respectively. Assumption 1. The parameter space Θ and the design matrixX satisfy the following,
• (A1) The parameter space Θ = {θ ∈ Rd : ‖θ‖2 ≤ w,min(‖θ‖ 2 2, 〈θ,xi〉) ≥ γ > 0, ∀i ∈ [n]}. • (A2) The design matrix is such that λmin(Σ) > 0 and ‖xi‖2 ≤ R for all i ∈ [n]. • (A3) Let λmin(M) ≥ R 2 4nγ2 (d log(24χ) + log(1/δ)) and √ λmax(M)(d log(24χ) + log(1/δ)) ≤√
nλmin(M)/16, for a small δ ∈ (0, 1) 1.
The above assumptions are fairly mild. For instance λmin is Ω̃(1/d) for random Gaussian covariance matrices (Bai & Yin, 2008). The other part merely requires that λmin(M) = Ω̃( √ dλmax(M)/n).
We are interested in comparing MLE with TMO for the square loss which is just the least-squares estimator i.e θ̂ ls := argminθ∈Θ 1 n‖y
n −Xθ‖22. Note that it is 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 is explicitly geared towards minimizing the square loss during training.
Least squares analysis is typically provided for regression under sub-Gaussian noise. By adapting existing techniques (Hsu et al., 2012), we show the following guarantee for Poisson regression with least square loss. We provide a proof in Appendix D for completeness.
1Note that the constants can be further tightened in our analysis.
Lemma 1. Let µmax = maxi µi. The least squares estimator θ̂ ls satisfies the following loss bounds w.p. at least 1− δ,
E(θ̂ ls) =
{ O ( µmax n ( log 1δ + d )) if µmax ≥ (log(1/δ) + d log 6)/2
O (
1 n ( log 1δ + d )2) otherwise
Now we present our main result in this section which uses the competitiveness bound in Theorem 2 coupled with the existence of a superior estimator compared to TMO, to show that the MLE estimator can have a better bound than TMO.
In Theorem 4 (in Appendix F), under some mild assumptions on the covariates, we construct an estimator θest with the following bound for the Poisson regression setting,
E(θest) ≤ c · ‖θ∗‖2 · λmax(Σ) ( d+ log( 1δ )
n
) . (3)
The construction of the estimator is median-of-means tournament based along the lines of (Lugosi & Mendelson, 2019a) 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 1 and the conditions of Theorem 4, the MLE estimator for the Poisson regression setting satisfies w.p. at least 1− δ,
E(θ̂mle) = O ( ‖θ∗‖2 · λmax(Σ) d(log n+ log(wRλmax(Σ)) + logχ+ log 1 δ )
n
) .
The bound in Corollary 1 can be better than the bound for θ̂ ls in Lemma 1. In the sub-Gaussian region, the bound in Lemma 1 scales linearly with µmax which can be prohibitively large even when a few covariates have large norms. The bound for the MLE estimator in Corollary 1 has no such dependence. Further, in the sub-Exponential region the bound in Lemma 1 scales as Õ(d2/n) while the bound in Corollary 1 has a Õ(d/n) dependency, up to log-factors. In Appendix G, we also show that when the covariates are one-dimensional, an even sharper analysis is possible, that shows that the MLE estimator is always better than least squares in terms of excess risk. In Appendix I.8, we perform a simulated experiment that shows a linear growth of the test error w.r.t λmax(Σ), further showing the efficacy of our theoretical bounds.

heading: 4.2 PARETO REGRESSION
text: Now we will provide an example of a 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,
p(yi|xi) = bmbi yb+1i , mi = b− 1 b 〈θ∗,xi〉 for yi ≥ mi (4)
provided 〈θ∗,xi〉 > γ for all i ∈ [n]. Thus yi is Pareto given xi and E[yi|xi] = µi := 〈θ∗,xi〉. We will assume that b > 4 such that 4 + -moment exists for > 0. As in the Poisson setting, we choose this parametrization for a fair comparison between TMO and MLE i.e in the limit of infinite samples θ∗ lies in the linear solution space for both TMO (least squares) and MLE.
As before, to apply Theorem 2 we need an estimator with a good risk bound. We use the estimator in Theorem 4 of (Hsu & Sabato, 2016), which in the fixed design pareto regression setting yields,
E(θest) = ( 1 +O ( d log 1δ n )) ‖θ∗‖2Σ b(b− 2) .
Note that the above estimator might not be easily implementable, however this yields the following corollary of Theorem 2, which is a bound on the performance of the MLE estimator. Corollary 2. Under assumptions of our Pareto regression setting, the MLE estimator satisfies w.p. at least 1− δ,
E(θ̂mle) = 1 +O
d2 ( log n+ log ζ + log bwRλmax(Σ)γ + log 1 δ ) n  ‖θ∗‖2Σ b(b− 2) .
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 is least-squares. It is well established in the literature (Hsu & Sabato, 2016; Lugosi & Mendelson, 2019a) that ERM on square loss cannot achieve a O(log(1/δ)) dependency in a heavy tailed regime; instead it can achieve only a O( √ 1/δ) rate.

heading: 5 CHOICE OF LIKELIHOOD AND INFERENCE METHODS
text: In this section we discuss some practical considerations for MLE, such as adapting to a target metric of interest at inference time, and the choice of the likelihood family.
Inference for different target metrics: In most practical settings, the trained regression model is used at inference time to predict the response variable on 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 at 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 metrics in the forecasting domain such as Mean Absolute Percentage Error (MAPE) and Relative Error (RE) (Gneiting, 2011; Hyndman & Koehler, 2006), this estimator corresponds to a median of a transformed distribution (Gneiting, 2011). Please see Appendix I 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 a 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 to cater to diverse datasets with potentially different types of inductive biases, or even datasets for which no distributional assumptions can be made in advance.
A simple approach that we observe to 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 to a different distribution family and the mixture weights are learnt along with the parameters of the distribution. More specifically, we consider a mixture distribution of k components p(y|x;θ1, . . . ,θk, w1, . . . , wk) = ∑k j=1 wjpj(y|x;θj), where each pj characterizes a different distribution family, and the mixture weights wj and distribution parameters θj are learnt together. We would like to have a mixture distribution that can handle different situations like sparse data, sub-Exponential and sub-Gaussian tails, count data as well as heavy tailed data. Moreover it should be applicable to continuous valued datasets in general.
We use a three component mixture of (i) the constant 0 (zero-inflation for dealing with bi-modal sparse data), (ii) a continuous version of negative binomial where n and p are learnt and (iii) a Pareto distribution where the scale parameter is learnt. We provide more details about the continuous version of negative binomial distribution in Appendix I.2. Our experiments in Section 6 show that this mixture shows promising performance on a 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 over-parameterized networks and gradient-based optimizers, this is not a problem in practice (Fig. 4 shows a convergence curve).

heading: 6 EMPIRICAL RESULTS
text: We present empirical results on two 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.
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(MSE) 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.
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-likelihood implemented in Tensorflow Probability (Dillon et al., 2017). Note that for the MLE approach only one model is 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(ZNBP)loss, where ZNBP refers to the mixture components Zero, Negative-Binomial and Pareto.
For TMO, the output layer of the models map to ŷ and we directly minimize the target metric in question. Note that this means we need to train a separate model for every target metric. Thus we have one model each for target metrics in {’MSE’, ’MAE’, ’MAPE’}. Further we also train a model using the Huber loss 2. In order to keep the number of parameters the same as that of MLE, we add an additional 6 neurons to the TMO models.

heading: 6.1 EXPERIMENTS ON FORECASTING DATASETS
text: 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, a 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 to 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.
2The Huber loss is commonly used in robust regression (Huber, 1992; Lugosi & Mendelson, 2019a)
Model p10QL p90QL
TMO (Quantile) 0.0973±0.0002 0.0628±0.0019 MLE(ZNBP) 0.0788±0.0008 0.0536±0.0007
Table 3: The MLE model predicts the empirical quantile of interest during inference. It is compared with Quantile regression (TMO based). The results, averaged over 50 runs along with the corresponding confidence intervals are presented.
Model MAPE WAPE RMSE
MLE(NB) 0.3314+/-0.0016 0.2521+/-0.002 175.501+/-1.1928
MLE(ZNB) 0.3186+/-0.0011 0.2453+/-0.002 170.0075+/-1.282
MLE(ZNBP) 0.3139±0.0011 0.2238±0.0009 164.6521+/-1.5185
Table 4: We perform an ablation study on the Favorita dataset, where we progressively add the components of our mixture distribution. There are three MLE models in the progression: Negative Binomial (NB), Zero-Inflated Negative Binomial (ZNB) and finally ZNBP.
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 1 and described in more detail in Appendix I.
We present our 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, TMO(MAE) performs better than TMO(MSE) for the RMSE metric on the Favorita dataset. In Table 4 we perform an ablation study on the Favorita dataset, where we progressively add mixture components resulting in three MLE models: Negative Binomial, Zero-Inflated Negative Binomial and finally ZNBP. This shows that each of the components add value in this dataset.

heading: 6.2 EXPERIMENTS ON REGRESSION DATASETS
text: 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 as 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 a 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 a 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 our 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 a separate model per target metric. On the Gas Turbine dataset, the MLE 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 on 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 log-likelihood. We provide a more in-depth discussion in Appendix J.

heading: A PROOF FOR COMPETITIVENESS OF MLE
text: Theorem 3. Let π̂ be an estimator for a property π, such that for any f ∈ F and z ∼ f , Pr(‖π(f)− π̂(z)‖ ≥ ) ≤ δ. Then the MLE based estimator satisfies the following bound,
Pr(‖π(f)− π(fz)‖ ≥ 2 ) ≤ (|F|+ 1)δ.
Proof of Theorem 3. By triangle inequality and by the properties,
‖π(f)− π(fz)‖ ≤ ‖π(f)− π̂(z)‖+ ‖π̂(z)− π(fz)‖.
Hence, 1‖π(f)−π(fz)‖≥2 ≤ 1‖π(f)−π̂(z)‖≥ + 1‖π̂(z)−π(fz)‖≥ .
We now take expectation of both LHS and RHS of the above equation with respect to the distribution f . For the LHS, observe that
E[1‖π(f)−π(fz)‖≥2 ] = Pr(‖π(f)− π(fz)‖ ≥ 2 ).
For the first term in the RHS
E[1‖π(f)−π̂(z)‖≥ ] ≤ Pr(‖π(f)− π̂(z)‖ ≥ ) ≤ δ,
by the assumption in the theorem. Combining the above three equations yield that
Pr(‖π(f)− π(fz)‖ ≥ 2 ) ≤ E[1‖π̂(z)−π(fz)‖≥ ] + δ.
We now prove that E[1‖π̂(z)−π(f̃z)‖≥ ] ≤ |F|δ. E[1‖π̂(z)−π(fz)‖≥ ] = ∫ z f(z)1‖π̂(z)−π(fz)‖≥
(a) ≤ ∫ z fz(z)1‖π̂(z)−π(fz)‖≥
≤ ∫ z ∑ f f(z)1‖π̂(z)−π(f)‖≥
(b) = ∑ f ∫ z f(z)1‖π̂(z)−π(f)‖≥
= ∑ f Pr z∼f (‖π̂(z)− π(f)‖ ≥ )
≤ ∑ f δ = |F|δ,
where (a) uses the fact that fz is the MLE estimate over set F and hence fz(z) ≥ f(z). (b) follows from Fubini’s 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 Ff that satisfies assumptions in Definition 1. By triangle inequality and by the properties of F̃ , with probability at least 1− δ2 − δ,
‖π(f)− π(fz)‖ ≤ ‖π(f)− π̂(z)‖+ ‖π̂(z)− π(f̃z)‖+ ‖π(f̃z)− π(fz)‖ ≤ ‖π(f)− π̂(z)‖+ ‖π̂(z)− π(f̃z)‖+
We have used the fact that fz ∈ Ff w.p at least 1− δ. Hence with probability at least 1− δ2 − δ,
1‖π(f)−π(fz)‖≥3 ≤ 1‖π(f)−π̂(z)‖≥ + 1‖π̂(z)−π(f̃z)‖≥ .
We now take expectation of both LHS and RHS of the above equation with respect to the distribution f . For the LHS, observe that
E[1‖π(f)−π(fz)‖≥3 ] = Pr(‖π(f)− π(fz)‖ ≥ 3 ). For the first term in the RHS
E[1‖π(f)−π̂(z)‖≥ ] ≤ Pr(‖π(f)− π̂(z)‖ ≥ ) ≤ δ, by the assumption in the theorem. Combining the above three equations yield that
Pr(‖π(f)− π(fz)‖ ≥ 3 ) ≤ E[1‖π̂(z)−π(f̃z)‖≥ ] + 2δ + δ2.
We now prove that E[1‖π̂(z)−π(f̃z)‖≥ ] ≤ Tδ + δ1. Let f̃ be the distribution in F̃f that is at most δ1 away from f . Then,
E[1‖π̂(z)−π(f̃z)‖≥ ] = ∫ z f(z)1‖π̂(z)−π(f̃z)‖≥
(a) ≤ ∫ z f̃(z)1‖π̂(z)−π(f̃z)‖≥ + ‖f − f̃‖1
(b) ≤ ∫ z f̃(z)1‖π̂(z)−π(f̃z)‖≥ + δ1
(c) ≤ ∫ z f̃z(z)1‖π̂(z)−π(f̃z)‖≥ + δ1
≤ ∫ z ∑ f̃∈F̃ f̃(z)1‖π̂(z)−π(f̃z)‖≥ + δ1
(d) = ∑ f̃∈F̃ ∫ z f̃(z)1‖π̂(z)−π(f̃)‖≥ + δ1
= ∑ f̃∈F̃ Pr z∼f̃ (‖π(f̃)− π̂(z)‖ ≥ ) + δ1
≤ ( ∑ f̃∈F̃ δ) + δ1 = Tδ + δ1,
where the (a) follows from the definition of `1 distance between distributions, (b) follows from the properties of the cover, (c) uses the fact that f̃z is the MLE estimate over set F̃ and hence f̃z(z) ≥ f̃(z). (d) follows from Fubini’s Theorem.

heading: B OTHER GENERAL RESULTS FOR THE MLE
text: In this section, we prove that the log likelihood of the MLE distribution at the observed data point is is close to the log-likelihood of 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 at least 1− δ,
log f(z) ≤ log fz(z) ≤ log f(z) + log (∑ z∈Z fz(z) ) + log 1 δ ,
where ∑ z∈Z fz(z) is the Shtarkov sum of the family F . Furthermore, if F is finite, then
log f(z) ≤ log fz(z) ≤ log f(z) + log |F|+ log 1
δ ,
Further is the distribution family is finite, then
log f(z) ≤ log fz(z) ≤ log f(z) + log|F|+ log 1
δ .
Proof. Let K = 1δ ∑ z∈Z fz(z).
Pr(f(z) ≤ fz(z)/K) = Pr ( fz(z) Kf(z) ≥ 1 )
(a) ≤ E [ fz(z)
Kf(z) ] = 1
K ∑ z∈Z fz(z) = δ,
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 if F is finite then,∑
z∈Z fz(z) ≤ ∑ z∈Z ∑ f∈F f(z) = ∑ f∈F ∑ z∈Z f(z) = ∑ f∈F 1 = |F|.
We now prove a result of Shtarkov’s sum, which will be useful in other scenarios.
Lemma 3. Let w = π(z) for some property π.∑ w fw(w) ≤ ∑ z fz(z).
Proof. ∑ z fz(z) ≥ ∑ z pg(z)(z) = ∑ w ∑ z:g(z)=w fw(z) = ∑ w fw(w).

heading: C GENERAL FIXED DESIGN RESULT
text: In this section, we will prove Theorem 2 which is an application of Theorem 1 to the fixed design regression setting with the square loss as the target metric. The theorem holds under some reasonable assumptions. We will fist prove some intermediate lemmas.
Lemma 4. Let f(·;θ) = ∏n i=1 p(·|xi;θ) in the fixed design setting. If ‖θ − θ′‖2 ≤ δ then,
‖f(·;θ)− f(·;θ′)‖1 ≤ √ 2nLδ.
under the assumptions of Theorem 2.
Proof. By Pinskers’ inequality we have the following chain, ‖f(·;θ)− f(·;θ′)‖1 ≤ √ 2DKL(f(·;θ); f(·;θ′))
≤ √√√√2 n∑ i=1 DKL(p(·|xi;θ); p(·|xi;θ′)) ≤ √ 2nLδ
Lemma 5. Assume the conditions of Theorem 2 holds. Let θ̃mle := argmaxθ∈N (̃,Θ) L(yn,θ), in the Poisson regression setting. Then with probability at least 1− δ we have,
|E(θ̂mle)− E(θ̃mle)| ≤ 4wλmax(Σ)̃
n
√ β
α .
Proof. We assume that the event in (2) of Theorem 2 holds.
|L(yn; θ̂mle)− L(yn;θc)| ≤ β
2 ̃2,
where θc is the closest point in N (̃,Θ) to θ̂mle. Therefore, by definition
L(yn; θ̃mle) ≤ L(yn; θ̂mle) + (β/2)̃2.
By virtue of the strong-convexity we have,
∥∥∥θ̂mle − θ̃mle∥∥∥2 2 ≤ β α ̃2.
Now by triangle inequality we have, | √ E(θ̂mle)− √ E(θ̃mle)| = | ∥∥∥θ̂mle − θ∗∥∥∥ Σ − ∥∥∥θ̃mle − θ∗∥∥∥ Σ |
≤ ∥∥∥θ̂mle − θ̃mle∥∥∥ Σ ≤ √ λmax(Σ)β α ̃
Now we can use the fact |a− b| ≤ 2 max( √ a, √ b)| √ a− √ b| and that Θ ⊆ Bdw to conclude,
|E(θ̂mle)− E(θ̃mle)| ≤ 4wλmax(Σ)̃ √ β
α .
Proof of Theorem 2. Consider the net N (φ, θ). If φ = δ2/(2nL), then by Lemma 4 the net forms a δ cover on F . Note that under the conditions of Theorem 2 |N (φ, θ)| ≤ (3w/φ)d. Next note that in the application of Theorem 1, we should set
= (c1 + c2 log(1/δ))
η
n .
Thus if we set φ in place of ̃ in Lemma 5 we would need the following to apply Theorem 1,
4wλmax(Σ)φ
√ β
α =
(c1 + c2 log(1/δ)) η
n .
Thus, φ can be such that,
log 1
φ ≤ max
( log(2nL), 0.5 log β
α + log(4wλmax(Σ)) + log(n)
) .
From above we can apply Theorem 1 with the above , δ1 = δ2 = δ and T = (3w/φ)d. Therefore, we get
P ( E(θ̂mle) ≥ 3 (c1 + c2 log(1/δ)) η
n
) ≤ (T + 5)δ. (5)
We can set δ = δ/(T + 5) to then conclude that w.p at least 1− δ, E(θ̂mle) = O ( (c1 + c2d(log(w) + log(1/φ))) η
n
) .
Substituting the bound on log 1φ yields the result.

heading: D LEAST SQUARES FOR POISSON
text: We begin by restating the following general statement about OLS. Lemma 6 (see Theorem 2.2 in (Rigollet, 2015)). Suppose r = rank(XTX) and φ ∈ Rn×r be an orthonormal basis for the column-span of X. Then we have the following result,
E(θ̂ ls) ≤ 4
n sup
u∈Sr−11
(uT ̃)2, (6)
where ̃ = φT . Here, = yn −Xθ∗ for the given response sample yn. Now are at a position to prove Lemma 1.
Proof of Lemma 1. We are interested in tail bounds on the RHS of (6). We will first analyze tail bounds for a fixed u ∈ Sr−11 . We have the following chain,
E[exp{s〈u, ̃〉}] = E[exp{s〈φu, 〉}] := E[exp{s〈v, 〉}]
= n∏ i=1 E[exp{svi i}] = n∏ i=1 E[exp(svi(yi − µi))]
where µi = 〈θ∗,xi〉. Now bounding each term separately we have, E[exp(svi(yi − µi))] = exp(−sviµi) exp(µi(exp(svi)− 1))
We have, n∏ i=1 exp(µi(exp(svi)− 1)) = exp ( n∑ i=1 µi(exp(svi)− 1) )
= exp  n∑ i=1  ∞∑ j=1 µi j! (svi) j  ≤ exp  n∑ i=1 µisvi + ∞∑ j=2 µmax 1 j! (svi) j  ≤ exp
(∑ i µisvi + µmax(exp(s)− s− 1) ) , (7)
where we have used the fact that ||v||p ≤ ||v||2 = 1 for all p ≥ 2. Thus we have,
E[exp{s〈u, ̃〉}] ≤ exp(µmax(exp(s)− s− 1)).
This means that the RV 〈u, ̃〉 is sub-exponential with parameter ν2 = 2µmax and α = 0.56. Thus if t ≤ 4µmax then we have that for a fixed u,
P(〈u, ̃〉} ≥ t) ≤ exp ( − t 2
4µmax
) . (8)
Thus by a union bound over the -net with = 2, we have that wp 1− δ, sup u 〈u, ̃〉 ≤ √ 8µmax(log(1/δ) + r log 6).
Thus we get the risk bound wp 1− δ,
E(θ̂ ls) ≤ 32µmax
n (log(1/δ) + r log 6). (9)
For the sub-exponential region when t ≥ µmax we get the bound,
E(θ̂ ls) ≤ 16
n max
( (log(1/δ) + r log 6)
2 , µ2max
) .
We get the bound by setting r = d.

heading: E COMPETITIVENESS FOR POISSON REGRESSION
text: Lemma 7. Let f(·;θ) = ∏n i=1 p(·|xi;θ) where p is defined in Eq. (2). If ‖θ − θ′‖2 ≤ δ then,
‖f(·;θ)− f(·;θ′)‖1 ≤ R √ 2nw
γ δ.
Proof. By Pinskers’ inequality we have the following chain, ‖f(·;θ)− f(·;θ′)‖1 ≤ √ 2DKL(f(·;θ); f(·;θ′))
≤ √√√√2 n∑ i=1 DKL(p(·|xi;θ); p(·|xi;θ′))
≤ √√√√2 n∑ i=1 (〈θ,xi〉(log(〈θ,xi〉)− log(〈θ′,xi〉))− (〈θ,xi〉 − 〈θ′,xi〉))
Notice that the absolute value of the derivative of a log x− x wrt x, is upper bounded by a/γ − 1 if x ≥ γ > 0 and also a ≥ γ. Therefore, following from above we have,
‖f(·;θ)− f(·;θ′)‖1 ≤ √√√√2 n∑ i=1 ( θTxi γ − 1 ) |〈θ,xi〉 − 〈θ′,xi〉|
≤ √√√√2 n∑ i=1 ( θTxi γ − 1 ) ‖θ − θ′‖2R
≤
√ 2nwR2
γ ‖θ − θ′‖2,
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 ui = xi/‖xi‖2. Then with probability at least 1− 2δ, we have
λmin ( ∇2θL(yn;θ) ) ≥ n
2‖θ‖2 λmin(M), (10)
provided nλmin(M) ≥ 14γ2 (d log(24χ) + log(1/δ)) and √ λmax(M)(d log(24χ) + log(1/δ)) ≤√
nλmin(M)/16.
Proof. We start with the expression of the Hessian of L(yn;θ),
∇2θL = n∑ i=1 yixix T i 〈θ,xi〉2 < n∑ i=1 yiuiu T i ‖θ‖2 := L(θ),
where ui = xi/‖xi‖2.
For a fixed unit vector u ∈ Sd−11 let us define,
L(u,θ) := n∑ i=1 yiu Tuiu T i u ‖θ‖2 := n∑ i=1 yivi,
where vi = 〈ui,u〉2/‖θ‖2.
Following the definition of sub-exponential RV in (Rinaldo, 2019),∑ i yivi ∈ SE ( ν2 = ∑ i ν2i , α = 0.56 max i vi ) , (11)
where νi = 2µiv2i . This implies that wp atleast 1− 2δ/C we have,
| ∑ i yivi − ∑ i µivi| ≤ 2
√∑ µiv2i log ( C
δ ) if, √∑
µiv2i log
( C
δ
) ≤ 2( ∑ µiv 2 i )/(max
i vi).
Note that the above condition is mild and is satisfied when λmin (∑ i µiuiu T i ) ≥ log(C/δ)/4γ2, as
minu ∑ i µivi ≥ λmin (∑ i µiuiu T i ) . Here, C is a problem dependent constant that will be chosen later.
Now consider an -net in `2 norm over the surface unit sphere denoted by N ( , d). It is well known that |N ( , d)| ≤ (3/ )d. Now any z ∈ Sd−11 can be written as z = x+ u where x ∈ N ( , d) and u ∈ Bd−1 . Therefore, we have that
L(z,θ) ≥ L(x,θ)− λmax(L(θ))
A similar argument as above gives us,
λmax(L(θ)) = max u∈Bd−11 L(u,θ) ≤ max x∈N ( ,d)
L(x,θ) + 1
2 max u∈Bd−11 L(u,θ)
=⇒ λmax(L(θ)) ≤ 2 max x∈N ( ,d) L(x,θ)
Thus by an union bound over the net we have wp 1− 2δ and other conditions,
min z∈Sd−11 L(z,θ) ≥ 1 ‖θ‖22 λmin (∑ i µiuiu T i ) − 2 ‖θ‖22 λmax (∑ i µiuiu T i )
− max u∈Sd−11 2(1 + 2 ) √√√√∑ i µi〈u,ui〉2 ‖θ‖4 log ( C δ ) .
if the conditions above hold with C = (3/ )d. We can now set = 1/(8∗χ) where χ is the condition number of the matrix ∑ i µiuiu T i . Now by virtue of that fact that √ λmax( ∑ µiuiuTi ) log(C/δ) ≤
λmin (∑ i µiuiu T i ) /16 we have the result.
Now we will prove a result on the smoothness of the negative log likelihood for the Poisson regression setting. Lemma 9. With probability at least 1− 2δ, we have
λmax ( ∇2θL(yn;θ) ) ≤ 2nR 2
γ2 λmax(M), (12)
if λmin(M) ≥ R2(d log 6 + log(1/δ))/4nγ2.
Proof. We start by writing out the Hessian again but this time upper bounding it in semi-definite ordering,
∇2θL = n∑ i=1 yixix T i 〈θ,xi〉2 4 n∑ i=1 yixix T i γ2 4 n∑ i=1 yiR 2uiu T i γ2 := L(θ)
For a fixed unit vector u ∈ Sd−11 let us define,
L(u,θ) := n∑ i=1 yiR 2uTuiu T i u γ2 := n∑ i=1 yivi,
where vi = R2〈ui,u〉2/γ2. Proceeding as in Lemma 8 we conclude that,∑ i yivi ∈ SE ( ν2 = ∑ i ν2i , α = 0.56 max i vi ) , (13)
where νi = 2µiv2i . This implies that wp atleast 1− 2δ/C we have,
| ∑ i yivi − ∑ i µivi| ≤ 2
√∑ µiv2i log ( C
δ ) if, √∑
µiv2i log
( C
δ
) ≤ 2( ∑ µiv 2 i )/(max
i vi).
Note that the above condition is mild and is satisfied when λmin (∑ i µiuiu T i ) ≥ R2 log(C/δ)/4γ2,
as minu ∑ i µivi ≥ λmin (∑ i µiuiu T i ) . Here, C is a problem dependent constant that will be chosen later.
Now consider an 1/2-net in `2 norm over the surface unit sphere denoted by N (0.5, d). It is well known that |N ( , d)| ≤ (3/ )d. Now any z ∈ Sd−11 can be written as z = x+uwherex ∈ N (0.5, d) and u ∈ Bd−10.5 . Therefore, we have that,
λmax(L(θ)) = max u∈Bd−11 L(u,θ) ≤ max x∈N ( ,d)
L(x,θ) + 1
2 max u∈Bd−11 L(u,θ)
=⇒ λmax(L(θ)) ≤ 2 max x∈N ( ,d) L(x,θ)
Thus we set C = 6d and obtain wp at least 1− 2δ,
λmax(L(θ)) ≤ 2 nR2
γ2 λmax(M),
provided 4 √ λmax(M) log(C/δ) ≤ λmax(M)/ √ n.
Proof of Corollary 1. We show that the conditions of Theorem 2 hold.
Creating Nets: First we need to create an -net over the parameter space Θ. We start by creating an -net over the sphere with radius w. Now suppose, < γ/2R. Then we remove all centers θc if ∃i s.t 〈θc,xi〉 < γ/2. This is a valid -net over Θ as all net partitions that are removed do not have any points lying in Θ. In subsequent section, we will always follow this strategy to create -nets over subsets of Θ.
Strong Convexity and Smoothness: From Lemma 8 and 9 we have that,
β α = w2R2 γ2 χ.
with probability 1−O(δ). KL Divergence: The L in Theorem 2 is bounded by 2wR2/γ according to Lemma 7.
Combining the above into Theorem 2 and using the estimator in Section F we get our result.

heading: F MEDIAN OF MEANS ESTIMATOR FOR POISSON
text: 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 we have a fixed design matrix X ∈ Rn×d with rows x1, . . . ,xn, and for each i ∈ [n], yi is drawn from a Poisson distribution with mean µi = θ∗ · xi. 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 u
E[ ( Σ− 1 2x · u )4 ] ≤ L · E[ ( Σ− 1 2x · u )2 ]2.
For simplicity we will assume that L = O(1). In the above definition note that E is the empirical expectation over teh fixed design. This is a benign assumption and for instance would be satisfied if the design matrix is drawn from a sub-Gaussian distribution. In the general case, the obtained bounds will scale with L. We have the following guarantee associated with Algorithm 1.
Algorithm 1: Median of Means Estimator Input: Samples S = {(x1, y1), . . . , (xn, yn)}, confidence parameter δ. Step 1: Compute Σ = E[xx>]. Form S′ = {x′1 = y1Σ− 1 2x1, . . . ,x ′ n = ynΣ
− 12xn}. Step 2: Randomly partition S′ into k blocks of size n/k each where k = 20dlog( 1δ )e. Step 3: Feed in the k blocks to the median-of-means estimator of Lugosi & Mendelson (2019b) to get v̂. Step 4: Return θ̂ = Σ− 12 v̂.
Theorem 4. There is an absolute constant c > 0 such that with probability at least 1−δ, Algorithm 1 outputs θ̂ such that
‖θ̂ − θ∗‖2Σ ≤ c · ‖θ∗‖2 · λmax(Σ) (d+ log( 1δ )
n
) . (14)
Proof. The proof is exactly along the lines of the proof of Theorem 1 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[yΣ− 1 2x] = E[(θ∗ · x)Σ− 12x)] (15)
= E[(θ∗ · x)Σ− 12x] (16)
= Σ 1 2θ∗. (17)
Hence if v̂ is the output of the median of the means estimator in Step 3 of Algorithm 1, then the least squares error of θ̂ is exactly ‖v̂ − Σ 12θ∗‖2. For convenience define µ′ = Σ 12θ∗ and Σ′ = E[x′x′>]− Σ 12θ∗θ∗>Σ 12 . Exactly as in Lugosi & Mendelson (2019b) our goal is to show that µ′ beats any other vector v in the median of means tournament if v is far away from µ′. To quantify this define
r = max ( 400
√ Tr(Σ′)
n , 4 √ 10
√ λmax(Σ′)
n
) .
For a fixed vector v of length r, and block Bj , µ′ beats v if
−2k n ∑ i∈Bj (x′i − µ′) · v + r > 0.
Let us denote by σi,j ∈ {0, 1} a random variable representing whether data point i is in block j or not. By Chebychev’s inequality we get that with probability at least 9/10,∣∣k
n ∑ i∈Bj (x′i − µ′) · v ∣∣ ≤ k n √ 10 √∑ i E[σ2i,j((x′i − µ′) · v)2] (18)
= k
n
√ 10 √ np 1
n ∑ i E[((x′i − µ′) · v)2]. (19)
Here p is the probability of a point belonging to block k. Noting that np = Θ(n/k) we get that with probability at least 9/10, ∣∣k
n ∑ i∈Bj (x′i − µ′) · v ∣∣ ≤ √10r√kλmax(Σ′) n . (20)
Applying binomial tail estimates we get that with probability at least 1− e−k/180, µ′ beats v on at least 8/10 of the blocks. By applying the covering argument verbatim as in the proof of Theorem 1 in Lugosi & Mendelson (2019b) we get that with probability at least 1− δ, v̂ will satisfy
‖v̂ − Σ 12θ∗‖2 ≤ c · λmax(Σ′) (d+ log( 1δ )
n
) .
Finally, it remains to bound the spectrum of Σ′. We have
Σ′ = E[y2Σ− 1 2xx>Σ− 1 2 ]− Σ 12θ∗θ∗>Σ 12 (21)
E[((θ · x)2 + θ · x)Σ− 12xx>Σ− 12 ] (22) T1 + T2, (23)
where
T1 = E[ ( θ∗ · x )2 Σ− 1 2xx>Σ− 1 2 ], (24)
T2 = E[ ( θ∗ · x ) Σ− 1 2xx>Σ− 1 2 ]. (25)
To bound T1, T2 we note that for any function m(x) we have
E[m(x)Σ− 1 2xx>Σ− 1 2 ] √ E[m2(x)]I. (26)
Using the above inequality and the fact that the design matrix is (O(1), 4) hyper-contractive we get that √
E[(θ∗ · x)2] = O(‖θ∗‖ √ λmax(Σ)) (27)√
E[(θ∗ · x)4] = O(‖θ∗‖2λmax(Σ)). (28) Combining the above we get that
Σ′ T1 + T2 (29) O(λmax(Σ))I. (30)
G 1-D POISSON REGRESSION
When the covariates are one dimensional, a 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 c ≥ 2 such that for any δ ∈ [ 1nc , 1), it holds with probability at least 1− δ that,
E(θ̂ls) = 1
n n∑ i=1 (yi − θ̂ls)2 ≤ 4 · |θ∗| n · ∑n i=1 |xi|3∑n i=1 x 2 i log ( 1 δ ) . (31)
The bound above is also tight i.e with constant probability it holds that
E(θ̂ls) = 1
n n∑ i=1 (yi − θ̂ls)2 = Ω ( |θ∗| n · ∑n i=1 |xi|3∑n i=1 x 2 i︸ ︷︷ ︸
:=B(θ̂ls)
) . (32)
We provide the proofs in later in the section. Having established the bound for least squares estimator, we next prove the following upper bound on the mean squared error achieved by the MLE.
Theorem 5. There exists an absolute constant c ≥ 2 such that for any δ ∈ [ 1nc , 1), it holds with probability at least 1− δ that,
E(θ̂mle) = 1
n n∑ i=1 (yi − θ̂mle)2 ≤ 4 · |θ∗| n · ( ∑n i=1 x 2 i∑n
i=1 |xi|︸ ︷︷ ︸ :=B(θ̂mle)
) log ( 2
δ
) . (33)
It is easy to 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 θ̂ls (defined in Eq. 33). To see this notice that,
B(θ̂ls) B(θ̂mle)
= ( ∑n i=1 |xi|3)( ∑n i=1 |xi|)
( ∑n i=1 x 2 i ) 2 ≥ 1. (from Cauchy-Schwarz inequality.)
Furthermore, in many cases the bound achieved by θ̂mle can be significantly better than the one achieved by θ̂ls. As an example consider a skewed data distribution where √ n of the xi’s take a large
value of √ n, while the remaining data points take a value of n , where is a small constant. In this case we have that,
B(θ̂ls) B(θ̂mle)
= ( ∑n i=1 x 3 i )( ∑n i=1 xi)
( ∑n i=1 x 2 i ) 2 = Ω(n ).
Proof of Lemma 10. Notice that
E(θ̂ls) = 1
n (θ̂ls − θ∗)2( n∑ i=1 x2i ). (34)
Hence, it is enough to bound the parameter distance, i.e., (θ̂ls − θ∗)2. In order to do that we first notice that yixi is a sub-exponential random variable with parameters ν2i , α where where ν 2 i = 2µix 2 i , µi = θ ∗xi, and αi = 0.56xi (Rinaldo, 2019). In other words,
yixi ∈ SE(ν2i , α), (35)
Thus from the bound on a sum of independent sub-exponential variables we have that,∑ i yixi ∈ SE(ν2 = ∑ i ν2i , α = 0.56 max i xi). (36)
Thus we have that,
P ( | ∑ i yixi − ∑ i µixi| ≥ t ) ≤ { 2 exp ( − t 2 2ν2 ) if t ≤ ν 2 α 2 exp ( − t2α ) otherwise
This means, that w.p at least 1− 2δ,
| ∑ i yixi − ∑ i µixi| ≤ 2 √ ( ∑ i µix2i ) log ( 1 δ ) if √
( ∑ i µix2i ) log ( 1 δ ) ≤ 2 ∑ µix 2 i 0.56 maxi xi
The condition above is satisfied under our assumptions on xi and δ, thereby leading to the bound
(θ̂ls − θ∗) ≤ 2 √ w log(1/δ) √∑ x3i∑
i x 2 i
.
The bound on the MSE claimed in the lemma then follows.
Now we prove the lower bound. Again it is enough to show a lower bound on |θ̂ls − θ∗|. Notice that θ̂ls − θ∗ = ∑n i=1(yi − µi)xi∑n
i=1 x 2 i
. (37)
Define the random variable Z = ∑n i=1(yi − µi)xi. We will show anti-concentration for Z by computing the hyper-contractivity of the random variable. Recall that a random varibale Z is η-HC (hyper-contractive) if E[Z4] ≤ η4E[Z2]2. Next we have
E[Z2] = E[ n∑ i=1 (yi − µi)xi]2 = n∑ i=1 µix 2 i . (38)
E[Z4] = E[ n∑ i=1 (yi − µi)xi]2
= n∑ i=1 µi(1 + 3µi)x 4 i + 2 ∑ i 6=j µiµjx 2 ix 2 j
:= ∆. (39)
Hence Z is η-HC with η4 =
(∑n i=1 µix 2 i )2 ∆ .
From anti-concentration of hyper-contractive random variables (O’Donnell, 2014) we have
P(|Z| ≥ 1 2
√ E[Z2]) ≥ Ω(η4). (40)
Hence we get that
P | n∑ i=1 (yi − µi)xi| ≥ 1 2 √√√√ n∑ i=1 µix2i  ≥ Ω(η4) ≥ (∑n i=1 µix 2 i )2 ∆
(41)
(42)
Next notice that since µi ≥ γ for all i (Assumption 1), we have that 1 + 3µi ≤ (3 + 1γ )µi. This implies that ∆ ≤ (3 + 1γ ) (∑n i=1 µix 2 i )2 . Hence if γ is a constant then with probability at least
1 3+ 1γ = Ω(1) we have
E(θ̂ls) = Ω ( 1 n · ( ∑n i=1 x 3 i )∑n
i=1 x 2 i
) . (43)
Proof of Theorem 5. Recall that θ̂mle is defined as
θ̂mle = argmin θ∈Θ n∑ i=1 θxi − yi log(θxi). (44)
Setting the gradient of the objective to zero, we get the following closed form expression for θ̂mle.
θ̂mle = ∑n i=1 yi∑n i=1 xi . (45)
Next, we note that Z = ∑n i=1 yi is a poission random random variable with parameter µ =∑n
i=1 θ ∗xi. From tail bounds for Poisson random variables (Klar, 2000) we have that for any > 0,
P[|Z − µ| > ] ≤ 2e− 2 µ+ . (46)
Taking = c √ µ log( 1δ ), we get that with probability at least 1− δ,
| n∑ i=1 yi − µ| ∈ (1 2 , 2 )√ µ log( 1 δ ), (47)
provided that log( 1δ ) < µ (that holds for our choice δ, once n is large enough). Hence we conclude that with probability at least 1− δ, the mean squared error of θ̂mle is bounded by
|θ̂mle − θ∗| = ∣∣∣∑ni=1 θ∗xi∑n i=1 xi − ∑n i=1 yi∑n i=1 xi ∣∣∣ = O (√µ log( 1δ )∑n i=1 xi )
= O (√|θ∗| log( 1δ )√∑n
i=1 xi
) . (48)
The bound on the mean squared error follows from the above.

heading: H COMPETITIVENESS FOR PARETO REGRESSION
text: We verify the conditions of Theorem 2 for the Pareto regression setting. Lemma 11. Let f(·;θ) = ∏n i=1 p(·|xi;θ) where p is defined in Eq. (4). If ‖θ − θ′‖2 ≤ δ then,
‖f(·;θ)− f(·;θ′)‖1 ≤
√ 2bnδR
γ
Proof. By Pinskers’ inequality we have the following chain, ‖f(·;θ)− f(·;θ′)‖1 ≤ √ 2DKL(f(·;θ); f(·;θ′))
= √√√√2 n∑ i=1 DKL(p(·|xi;θ); p(·|xi;θ′))
≤ √√√√ n∑ i=1 2b| logmi − logm′i|
≤ √√√√ n∑ i=1 2b|〈θ,xi〉 − 〈θ′,xi〉| γ
≤
√ 2bn‖θ − θ′‖2R
γ .
The second inequality follows from the fact that log x is Lipschitz with parameter L if x > 1/L. The last inequality follows from Cauchy-Schwarz and the norm bound on xi’s.
Now we prove the rest of the conditions. Lemma 12. We have the following smoothness bound,
λmax(∇2θL) ≤ (b− 1) γ2 λmax(Σ).
Proof. We start by writing out the Hessian,
∇2θL = n∑ i=1 (b− 1) 〈θ,xi〉2 xix T i := L(θ)
Using the fact that 〈θ,xi〉 ≥ γ gives us the result.
Lemma 13. We have the following strong convexity bound,
λmin(∇2θL) ≥ (b− 1) w2R2 λmin(Σ).
Proof. It follows from the expression of the Hessian in Lemma 12 and using the fact 〈θ,xi〉 ≤ ‖θ‖2‖xi‖2
Proof of Corollary 2. We show that the conditions of Theorem 2 hold.
Creating Nets: First we need to create an -net over the parameter space Θ. We start by creating an -net over the sphere with radius w. Now suppose, < γ/2R. Then we remove all centers θc if ∃i s.t 〈θc,xi〉 < γ/2. This is a valid -net over Θ as all net partitions that are removed do not have any points lying in Θ. In subsequent section, we will always follow this strategy to create -nets over subsets of Θ.
Strong Convexity and Smoothness: From Lemma 13 and 12 we have that,
β α = w2R2 γ2 ζ.
KL Divergence: The L in Theorem 2 is bounded by 2bR/γ according to Lemma 11.
Combining the above into Theorem 2 and using the estimator in (Hsu & Sabato, 2016) we get our result.

heading: I MORE ON EXPERIMENTS
text: We provide more experimental details in this section.
I.1 METRICS
The metrics and loss functions used are as follows:
MSE: The metric is 1
n n∑ i=1 (ŷi − yi)2.
RMSE is just the square-root of this metric.
MAE: The metric is 1
n n∑ i=1 |ŷi − yi|.
WAPE: The metric is ∑n i=1|ŷi − yi|∑n i=1 |yi| .
MAPE: The metric is 1
n ∑ i:yi 6=0 ∣∣∣∣1− ŷiyi ∣∣∣∣.
Quantile Loss: The reported metrics in Table 3 are the normalized quantile losses defined as, n∑ i=1 2ρ(yi − ŷi)Iyi≥ŷi + 2ρ(ŷi − yi)Iyi<ŷi∑n i=1 |yi| . (49)
During training the unnormalized version is used for quantile regression.
Huber Loss: The loss is given by,
Lδ(y, ŷ) =
{ 1 2 (y − ŷ)
2, if |y − ŷ| ≤ δ δ|y − ŷ| − δ 2
2 , otherwise.
I.2 MORE DETAILS ABOUT THE MIXTURE DISTRIBUTION
We use a mixture distribution between zero, a continuous extension of negative binomial (NB) distribution and Pareto. The continuous extension of negative binomial is such that the p.d.f at y is proportional to,
p(y) ∝ Γ(n+ k) Γ(k + 1)Γ(n) (1− p)npk
given parameters n and p. 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 . It has been used in many regression datasets before, especially in the field of genomics where the collected data is a discrete continuous mixture (Robinson & Smyth, 2008; Chen et al., 2016; McCarthy et al., 2012).
Why not use a discrete-continuous mixture of zero, discrete NB and Pareto?
A mixture of these three distributions is a discrete continuous mixture, whose CDF is well-defined. It is possible to define a likelihood in a standard manner such that it has a 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 pmf at non-negative integral points. We have an extra mass at zero to account for the zero component. The sum of the dirac masses along with the integral of the Pareto density sums to one. However, such a 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-integral points have no contribution from the other components) which is not a desirable outcome, as a NB distribution can better model sub-Gaussian and sub-Exponential data in terms of moments.
I.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 a softmax layer (Goodfellow et al., 2016) to mixture weights. The fourth dimension is mapped to the ’n’ in negative-binomial likelihood through the link function,
φ(x) = { x+ 1, if x > 0 1/(1− x) otherwise .
The fifth dimension is mapped to ’p’ of the negative-binomial though the sigmoid function. The last dimension is mapped to the scale parameter for the Pareto component using the φ(x) link function above.
I.4 HARDWARE
We use the Tesla V100 architecture GPU for our 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.
I.5 MORE DETAILS ABOUT INFERENCE FOR MLE
We follow the approach of monte-carlo sampling. For each inference sample x′, we generate 10k samples from the learnt distribution p(·|x′; θ̂mle). Then we compute the correct statistics. For MSE, RMSE the statistic is just the mean and for WAPE, MAE it is the median. For a quantile, it is the corresponding quantile from the empirical distribution. For any loss of the form,
`(y, ŷ) = ∣∣∣∣1− (yŷ )β∣∣∣∣
the optimal statistic is the median from the distribution proportional to yβp(y|x′; θ̂mle) (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.
I.6 MODELS, HYPERPARAMETERS AND TUNING
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 h neurons. Note that h was tuned in 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 6 neurons before the output.
We tuned the learning rate for Adam optimizer in the range [1e-5, 1e-1] in log-scale. The batch-size was also tuned in [64, 128, 256, 512] and the Huber-δ in [2i for i 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 to be 512 and the Huber-δ was 32 and 64 for M5 and Favorita respectively.
For the regression datasets the model is 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 [1e-5, 1e-1] in log-scale. The batch-size was also tuned in [64, 128, 256, 512] and the Huber-δ in [2i for i in range(-8, 8)]. The learning rate was eventually chosen as 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-δ was 128 and 32 for Bicyle Share and Gas turbine respectively.
For the ZNBP model we also tune the α parameter in the pareto component between [3, 4, 5]. The value of 3 was selected for all datasets, except for gas turbine where we used α = 5.
MLE TMO
We used a batched version of GP-UCB (Srinivas et al., 2009) to 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.
I.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 a parent time-series is obtained as the mean of the time-series of its children. This is closer to a real forecasting setting in practice where one is interested in all levels of the hierarchy. The metrics reported are over all the time-series (both parents and leaves) treated equally. The history length for our 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 on steps 1900 to 1913. For the Favorita dataset the validation scores are computed using the predictions from time steps 1660 to 1673, and test scores on steps 1674 to 1687.
The train test splits are as mentioned in the main paper. For the Gas turbine dataset we use the official train test split. For the Bicycle share data there is no official split, but we use a randomly chosen fixed 10% as the test set for all our experiments.
I.8 ADDITIONAL EXPERIMENTS AND FIGURES
In order to show the dependency of λmax(Σ) in the bound in Corollary 1 we perform a simulated experiment. We generate a dataset with n = 5000 and d = 10 such that each coordinate of x is distributed i.i.d from a uniform distribution between [0, U ]. In out experiment, we vary U in {1, 2, · · · , 9}. Then y is genearted from a Poisson distribution with rate 〈θ∗,x〉 for a fixed θ∗. This varies λmax(Σ) which is equal to U2/12. We train and validate on 2500 samples with early stopping and plot the squared loss achieved on the test by the θ̂mle based estimator in Figure 3 versus λmax(Σ). We can clearly see a linear relationship which further validates our theoretical results.
In Figure 4, we plot the average training loss as a function of training iterations for the MLE(ZNBP)model. We can see that the loss converges to a minima.

heading: J EXTENDED DISCUSSION AND LIMITATIONS
text: We advocate for MLE with a suitably chosen likelihood family followed by post-hoc inference tuned to the target metric of interest, in favor of ERM directly with the target metric. On the theory side, we prove a general result that shows competitiveness of MLE with any estimator for a target metric under fairly general conditions. Application of the bound in the case of MSE for Poisson regression and Pareto regression is shown. We believe that our 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 5 can adapt quite well to different datasets as the mixture weights are trained. As we have mentioned before,
the MLE log-likelihood 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-batch 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:
If the overall problem is convex for TMO but introducing a MLE loss makes the problem nonconvex, 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 a linear function class, which is just least-squares linear regression, but introduction of a mixture likelihood like the one in Section 4 makes the problem non-convex. In this case it might be better to stick with TMO or 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 4.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 is something simple and convex like the square loss. In such a 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 it is better to capitalize on 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.

Review summaries:
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.

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. 

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.