Maximum Regularized Likelihood Estimators:
A General Prediction Theory and Applications

Maximum regularized likelihood estimators (MRLEs) are widely used in generalized linear regression, tensor response regression, and graphical modeling with high-dimensional data. It is thus of major interest to develop theory for this class of estimators.

Our specific goal is a general finite sample theory for prediction. Existing results are typically derived on a case-by-case basis. Moreover, many of these results also invoke restricted eigenvalues-type conditions (Bühlmann & van de Geer, 2011, Section 6). Such conditions are not only stringent and unverifiable in practice but also unsuitable for prediction. For example, restricted eigenvalue conditions in regression limit the correlations among the covariates. However, although correlations can affect the identifiability of the parameters, for prediction, even perfectly collinear covariates do not necessarily have a negative impact; in contrast, collinearity can even be beneficial (Hebiri & Lederer, 2013; Dalalyan et al., 2017). We are thus interested in a theory that does not involve additional assumptions and provides bounds for a general class of MRLEs. Besides its abstract value, such a general theory can also provide support for specific examples of MRLEs, such as the recently introduced approaches to tensor regression (Zhou et al., 2013; Li et al., 2013; Sun & Li, 2016), whose prediction properties have not been fully grasped.

In this paper, we establish a general oracle inequality in terms of the Kullback-Leibler divergence. Oracle inequalities are a standard way to formulate finite sample bounds in high-dimensional statistics. The Kullback-Leibler divergence is a standard way to quantify prediction accuracies; it applies to any model and yet specializes to well-established and interpretable notions of prediction performance. Our proofs invoke only the convexity of the parametrization and the definiteness and positive homogeneity of the regularizers. This makes the result applicable to a variety of parametric and non-parametric models and allows for a broad class of convex and non-convex regularizers.

The remainder of this paper is organized as follows. We introduce the framework and the general result in Section 2. We then provide examples in Section 3. We finally conclude with a brief discussion in Section 4. All proofs are deferred to the Appendices in the supplementary material: proofs for the main result to Appendix A, proofs for the examples to Appendix B, and proofs for the bounds of the empirical process terms in Appendix C. In addition, Appendix D contains notation and properties of tensors.

There are two types of oracle inequalities in the literature: so-called “fast rate bounds” and so-called “slow rate bounds.” “Fast rate bounds” are proportional to the square of the regularization parameter. Many representatives of this type of bounds are found in the literature, such as Bunea et al. (2007b); Raskutti et al. (2015) for regression, Ravikumar et al. (2011) for graphical models, and more generally, Bühlmann & van de Geer (2011); van de Geer (2016) and references therein. For example, the corresponding bounds for lasso prediction are of the form $s\log p/(w^{2}n),$ where $s$ is the number of non-zero elements in the true regression vector, $p$ is the number of parameters, $w$ is the restricted eigenvalue, and $n$ is the number of observations. These bounds are typically considered fast, because they can match minimax rates, see Verzelen (2012) and references therein. However, they rely on sparsity, and more importantly, they invoke restricted eigenvalue-type conditions or concern computationally challenging estimators instead (Bunea et al., 2007a; Dalalyan & Tsybakov, 2007; Rigollet & Tsybakov, 2011; Dalalyan & Tsybakov, 2012a; Dalalyan & Tsybakov, 2012b). Moreover, these eigenvalue-type assumptions are unverifiable and often unrealistic in practice, and even if they hold, the additional factors (such as $s$ and $1/w^{2}$ for lasso) can be large.

On the other hand, oracle inequalities for prediction have be derived without sparsity or restricted eigenvalue conditions for lasso-type estimators (Greenshtein & Ritov, 2004; Rigollet & Tsybakov, 2011; Massart & Meynet, 2011; Koltchinskii et al., 2011; Huang & Zhang, 2012; Chatterjee, 2013; Bühlmann, 2013; Chatterjee, 2014; Lederer et al., 2016; Dalalyan et al., 2017). For example, the corresponding bounds for lasso prediction are of the form $\sqrt{\log p/n}\,\!|\!|{\beta}^{*}|\!|_{1}$, where ${\beta}^{*}$ is the true regression vector. Such bounds are typically referred to as “slow rate bounds,” because on a high level, the rates are $1/\sqrt{n}$ rather than $1/n.$ However, there are no questionable assumptions involved, and for regression, it has even been shown that $1/\sqrt{n}$ is the optimal rate in the absence of further assumptions (Foygel & Srebro, 2011; Zhang et al., 2017; Dalalyan et al., 2017). Overall, this means that “fast rate bounds” are not necessarily fast and “slow rate bounds” are not necessarily slow. To correct the misleading nomenclature, Lederer et al. (2016) suggested replacing the term “fast rate bound” with “sparsity bound’ and “slow rate bound” with “penalty bound.”

Although some examples of MRLEs have been equipped with assumptionless bounds, many other examples still lack such guarantees (or any prediction guarantees altogether). More generally, a broadly applicable prediction theory for MRLEs is still in need.

The contribution of this work is two-fold. First, Theorem 2.1 provides a general prediction guarantee for MRLEs in terms of the Kullback-Leibler loss. Besides being the first assumptionless bound in such a broad setting, the result specializes correctly to known results, such as for lasso, where the corresponding rates have been shown to be optimal up to log-factors. Second, we show that applications of the general theorem to specific examples lead to new guarantees in tensor response regression, generalized linear tensor regression, and graphical modeling. The theory thus also establishes new insights into individual cases of MRLEs.

Our first example is tensor response regression. In a standard notation (see Kolda (2006) or our Appendix D for details), tensor response regression is based on models of the form

where $Y^{i}\in\mathbb{R}^{b_{2}\times\dots\times b_{p}}$ is a $(p-1)$th order tensor response, ${{\mathrm{\Lambda}}^{*}}\in{\mathcal{L}}\subset\mathbb{R}^{b_{1}\times\dots\times b_{p}}$ is a $p$th order tensor coefficient, ${\bf z}^{i}\in\mathbb{R}^{1\times b_{1}}$ is a fixed or random row-vector of covariates, and $E^{i}\in\mathbb{R}^{b_{2}\times\dots\times b_{p}}$ is random $(p-1)$th order tensor noise. The operation ${{\mathrm{\Lambda}}^{*}}\times_{1}{\bf z}^{i}$ denotes the mode-$1$ product of ${{\mathrm{\Lambda}}^{*}}$ and ${\bf z}^{i}$.

Our goal is to estimate the predictive structure of the above model. Assuming that the noise tensors $E^{i}$ are mutually independent, the MRLEs in (1) are of the form

For computational ease, ${\mathcal{L}}$ is usually chosen as a set of low-rank tensors (Rabusseau & Kadri, 2016; Sun & Li, 2016). In any case, if the conditional density $f_{{\mathrm{\Lambda}}}(Y^{i}\mid{\bf z}^{i})$ is parametrized such that ${\mathrm{\Lambda}}\mapsto\log f_{{\mathrm{\Lambda}}}-\mathbb{E}_{{{\mathrm{\Lambda}}^{*}}}\log f_{{\mathrm{\Lambda}}}$ is convex, we can derive statistical guarantees in conditional Kullback-Leibler loss from Theorem 2.1. Importantly, we do not impose any additional restriction on the covariates or the noise; for example, we allow the covariates to be correlated with the noise. Typical regularizers for third order tensors, for example, include the sparsity-inducing regularizer at the entry level $u({\mathrm{\Lambda}}):=\sum_{i_{1}=1}^{b_{1}}\dots\sum_{i_{3}=1}^{b_{3}}|{\mathrm{\Lambda}}_{i_{1}i_{2}i_{3}}|$, at the fiber level $u({\mathrm{\Lambda}}):=\sum_{i_{2}=1}^{b_{2}}\sum_{i_{3}=1}^{b_{3}}\,\!|\!|{\mathrm{\Lambda}}_{\cdot i_{2}i_{3}}|\!|_{2}$, and at the slice level $u({\mathrm{\Lambda}}):=\sum_{i_{3}=1}^{b_{3}}\,\!|\!|{\mathrm{\Lambda}}_{\cdot\cdot i_{3}}|\!|_{F}:=\sum_{i_{3}=1}^{b_{3}}\sqrt{\sum_{i_{1}=1}^{b_{1}}\sum_{i_{2}=1}^{b_{2}}|{\mathrm{\Lambda}}_{i_{1}i_{2}i_{3}}|^{2}}$ and the low-rank inducing regularizer $u({\mathrm{\Lambda}}):=\,\!|\!|{\mathrm{\Lambda}}|\!|_{*}$ with $\,\!|\!|\cdot|\!|_{*}$ the tensor nuclear norm defined in Raskutti et al. (2015). Our framework covers all these examples.

For illustration, we consider tensor response regression with zero-mean array normal noise (Akdemir & Gupta, 2011; Hoff, 2011), the most widely-used representative of the above model class. The conditional Lebesgue density of $Y^{i}$ given by Hoff (2011) is

where $\,\!|\!|\cdot|\!|^{2}$ is the array norm, $b:=\prod_{j=2}^{p}{b_{j}}$, $\Sigma_{k}=A_{k}A_{k}^{\top}\in\mathbb{R}^{b_{k}\times b_{k}}$ with non-singular real matrix $A_{k}\in\mathbb{R}^{b_{k}\times b_{k}}$, $\Sigma^{-1/2}=\{A_{2}^{-1},\dots,A_{p}^{-1}\}$, and $\times$ denotes the tensor product. One can check readily that $\log f_{{\mathrm{\Lambda}}}(Y^{i}\mid{\bf z}^{i})-\mathbb{E}_{{\mathrm{\Lambda}}^{*}}\log f_{{\mathrm{\Lambda}}}(Y^{i}\mid{\bf z}^{i})$ is linear in ${\mathrm{\Lambda}}$. Hence our theory applies; in particular, Theorem 2.1 specializes to array normal models as follows.

For an instantiation of the bound, assume $\sum_{i=1}^{n}({\bf z}^{i})_{j}^{2}/n=1,$ $j\in\{1,\dots,b_{1}\},$ and $(\Sigma^{-1}_{k})_{i_{k}i_{k}}=h_{k}^{2},$ $h_{k}>0,$ $k\in\{2,\dots,p\},$ $i_{k}\in\{1,\dots,b_{k}\}.$ Consider the sparsity-inducing regularizer at the entry level $u({\mathrm{\Lambda}}):=\sum_{i_{1}=1}^{b_{1}}\dots\sum_{i_{p}=1}^{b_{p}}|{\mathrm{\Lambda}}_{{}_{i_{1},\dots,i_{p}}}|.$ Then, the tuning parameter $r$ can be calibrated such that with probability at least $1-2\exp(-t^{2}),$ it holds that

This concrete bound follows from Lemma C1. That lemma and its proof — as well as all other technical derivations — are deferred to the Appendix.

While the above results for tensor regression are novel, assumptionless bounds for simple regression with lasso-type estimators such as lasso (Tibshirani, 1996), group lasso (Yuan & Lin, 2006), sparse group lasso (Simon et al., 2013), and slope estimator (Bogdan et al., 2015) have been derived before. Simple linear regression is thus an ideal test case to confirm that our results specialize correctly. To this end, we first observe that for $p=2$ and $b_{2}=1$, tensor response regression with array normal noise reduces to ordinary linear regression of the form

where $y^{i}\in\mathbb{R}$ is a scalar response, ${\bf z}^{i}\in\mathbb{R}^{1\times b_{1}}$ is a row-vector of covariates, ${\beta}^{*}\in\mathbb{R}^{b_{1}}$ is the regression vector, and $\varepsilon^{i}\in\mathbb{R}$ is noise distributed as $\mathcal{N}(0,\sigma^{2})$. For $u({\beta}):=\,\!|\!|{\beta}|\!|_{1}:=\sum_{i=1}^{b_{1}}|{\beta}_{i}|$, the MRLE (1) becomes

which, setting $r=r^{\prime}/(2\sigma^{2})$, can be written in the standard lasso-form

If $r^{\prime}\geq 2\,\!|\!|\sum_{i=1}^{n}({\bf z}^{i})^{\top}\varepsilon^{i}|\!|_{\infty}$, Lemma 3.1 now implies the bound

This result equals the classical penalty bound for lasso prediction, see Hebiri & Lederer (2013, Equation 3) for example. It has been shown that these bounds are essentially optimal in the absence of other assumptions (Foygel & Srebro, 2011; Zhang et al., 2017; Dalalyan et al., 2017). Along the same lines, one can also show that our bounds specify correctly for the other lasso-type estimators mentioned above (Lederer et al., 2016, Section 3), and similarly, for trace regression (Koltchinskii et al., 2011, Theorem 1).

Our second example is generalized linear tensor regression. The corresponding models consist of two components (Zhou et al., 2013; Li et al., 2013): an exponential family distribution and a link function. The exponential family distribution reads

where $y^{i}\in\mathbb{R}$ is a scalar response, $\theta^{i}\in\mathbb{R}$ is the natural parameter, $\alpha>0$ is the overdispersion factor, and $b,c$ are known real-valued functions. The link function $g:\mathbb{R}\mapsto\mathbb{R},$ assumed strictly increasing, provides a linear connection between the mean functions $\mathbb{E}(y^{i}\mid\theta^{i})$ and tensor predictors ${\bf z}^{i}\in\mathbb{R}^{b_{1}\times\dots\times b_{p}}$ according to

where ${{\mathrm{\Lambda}}^{*}}\in{\mathcal{L}}\subset\mathbb{R}^{b_{1}\times\dots\times b_{p}}$ and $\langle\cdot,\,\cdot\rangle$ is the tensor inner product. One can check that $b^{\prime}(\theta^{i})=\mathbb{E}(y^{i}\mid\theta^{i})$. With canonical link $g:=(b^{\prime})^{-1}$, it holds that $\theta^{i}=\langle{{\mathrm{\Lambda}}^{*}},\,{\bf z}^{i}\rangle$ and the distribution of the response $y^{i}$ conditioned on ${\bf z}^{i}$ has density

Further, by introducing basis functions in the mean models, it is straightforward to extend this parametric setting to non-parametric frameworks. In sum, generalized linear tensor regression provides very flexible model classes for scalar responses.

We can now turn to the corresponding MRLEs. Given $n$ independent observations $(y^{i},{\bf z}^{i})$ and considering the canonical link, the MRLEs in (1) become

Similarly to tensor response regression, optimization over the full set $\mathbb{R}^{b_{1}\times\dots\times b_{p}}$ is computationally challenging due to high-dimensionality of the problem. Thus, $\mathcal{L}$ is typically chosen considerably smaller, with the belief that the true parameter has some additional structure. For example, a choice proposed in Zhou et al. (2013) is ${\mathcal{L}}:=\{{\mathrm{\Lambda}}\in\mathbb{R}^{b_{1}\times\dots\times b_{p}}\mid{\mathrm{\Lambda}}=\sum_{i=1}^{m}\beta_{1}^{(i)}\circ\dots\circ\beta_{p}^{(i)}\}$, where $m$ is a fixed integer, $\beta_{j}^{(i)}\in\mathbb{R}^{b_{j}}$, and $\circ$ denotes the outer product.

Let us now apply Theorem 2.1 to equip MRLEs in generalized linear tensor regression with theoretical guarantees. For this, note that the log-parametrization here is again linear, so that the main theorem indeed applies and yields the following results.

As a special case, Lemma 3.2 applies to ordinary logistic regression, where $p=1$, the canonical link is $g(x)=\log\big(x/(1-x)\big)$, and the MRLEs with a general regularizer are in the form of

Lemma 3.2 then implies the bound

for $r\geq\tilde{u}\big(\sum_{i=1}^{n}\big(y^{i}-e^{\langle{{\mathrm{\Lambda}}^{*}},\,{\bf z}^{i}\rangle}/\big(1+e^{\langle{{\mathrm{\Lambda}}^{*}},\,{\bf z}^{i}\rangle}\big)\big){\bf z}^{i}\big)$. For $\ell_{1}$-regularization, the tuning parameter $r$ can be calibrated such that with probability at least $1-2\exp(-t^{2}),$ it holds that

where $p^{i}:=\mathbb{E}_{{\mathrm{\Lambda}}^{*}}(y^{i})$. We refer to Lemma C2 in Appendix C for details. The bound complements results for (weighted) $\ell_{1}$-regularized logistic regression that have been derived under additional assumptions, see van de Geer (2008) and van de Geer (2016, Chapter 12.4).

Our third example is graphical modeling. We consider the exponential trace framework (Lederer, 2016), which encompasses standard types of graphical models, such as Gaussian graphical models (Yuan & Lin, 2007; Friedman et al., 2008), non-paranormal graphical models (Gu et al., 2015), and Ising models (Lenz, 1920; Brush, 1967); we refer to Lederer (2016, Section 2) for details.

Exponential trace models are based on densities of the form

with respect to some $\sigma$-finite measure ${{\color[rgb]{0,0,0}\nu}}$ on $\mathbb{R}^{p}$. Here, the matrix-valued parameter ${\mathrm{\Lambda}}$ encodes the dependence structure of the random vector $X\in\mathcal{X}\subset\mathbb{R}^{p}$, the matrix-valued function $T(\cdot)$ on $\mathcal{X}$ determines how the data enters the model, ${\color[rgb]{0,0,0}a}({\mathrm{\Lambda}})$ is the normalization, and ${\mathcal{L}}$ is a convex set of matrices ${\mathrm{\Lambda}}$ with finite normalization. Given independent observations $X^{1},\dots,X^{n}$ of $X$, the MRLEs in (1) are of the form

Because the function $\log f_{{\mathrm{\Lambda}}}-\mathbb{E}_{{{\mathrm{\Lambda}}^{*}}}\log f_{{\mathrm{\Lambda}}}$ is linear in ${\mathrm{\Lambda}}$, we can then apply Theorem 2.1 to derive the following bound.

As a special case, Lemma 3.3 applies to the graphical lasso for multivariate Gaussian data. Recall that the graphical lasso (Yuan & Lin, 2007; Friedman et al., 2008) is formulated as

where $S_{++}^{p}$ is the set of positive definite $p\times p$ matrices and $X^{1},\dots,X^{n}$ are i.i.d. samples from a centered Gaussian distribution with unknown covariance matrix $({{\mathrm{\Lambda}}^{*}})^{-1}$. The Kullback-Leibler loss of ${{\color[rgb]{0,0,0}{\widehat{{\mathrm{\Lambda}}}}}}$ reads

which is equivalent (up to the factor 1/2) to Stein’s loss of the centered multivariate Gaussian distribution with covariance matrix $({{\color[rgb]{0,0,0}{\widehat{{\mathrm{\Lambda}}}}}})^{-1}$ (James & Stein, 1961). Thus, if $r^{\prime}\geq\,\!|\!|\mbox{vec}\big(({{\mathrm{\Lambda}}^{*}})^{-1}-\frac{1}{n}\sum_{i=1}^{n}X^{i}(X^{i})^{\top}\big)|\!|_{\infty}$, Lemma 3.3 yields

Following Lemma C3 in Appendix C, the tuning parameter $r^{\prime}$ can be calibrated such that with probability at least $1-4\exp(-t^{2}),$ it holds that

for all $t$ such that $0<t<\sqrt{n/4-\log\big(p(p-1)\big)}$. Our theory thus establishes the rate $\sqrt{\log p/n}$ for the graphical lasso in the Kullback-Leibler loss. This prediction rate complements the known estimation rates $\sqrt{\log p/n}$ in vectorized-matrix $\ell_{\infty}$-norm and $\sqrt{\min\{s+p,d^{2}\}\log p/n}$ in spectral norm (Ravikumar et al., 2011), where $s$ is the number of non-zero elements in ${{\mathrm{\Lambda}}^{*}}$ and $d$ is the maximum node degree. However, those results additionally require the mutual incoherence condition. Our prediction rate also complements the known estimation rate $\sqrt{(s+p)\log p/n}$ in Frobenius norm and spectral norm (Rothman et al., 2008), which requires only mild assumptions on the population covariance matrix.