Minimum description length revisited

1. Introduction

The Minimum Description Length (MDL) Principle^1,2,3,4 is a theory of inductive inference that can be applied to general problems in statistics, machine learning and pattern recognition. Broadly speaking, it states that the best explanation for a given set of data is provided by the shortest description of that data. In 2007, one of us published the book The Minimum Description Length Principle (Ref. 4 referred to as G07 from now on), giving a detailed account of most works in the MDL area that had been done until then. During the last 10 years, several new practical MDL methods have been designed, and there have been exciting theoretical developments as well. It therefore seemed time to present an up-to-date combined introduction and review.

Why read this overview?

While the MDL idea has been shown to be very powerful in theory, and there have been a fair number of successful practical implementations, massive deployment has been hindered by two issues: first, in order to apply MDL, one needs to have basic knowledge of both statistics and information theory. To remedy this situation, here we present, for the first time, the MDL Principle without resorting to information theory: all the material can be understood without any knowledge of data compression, which should make it a much easier read for statisticians and machine learning researchers novel to MDL. A second issue is that many classical MDL procedures are either computationally highly intensive (for example, MDL variable selection as in Example 4 below) and hence less suited for our big data age, or they seem to require somewhat arbitrary restrictions of parameter spaces (e.g., NML with $v\equiv 1$ as in Sec. 2). Yet, over the last 10 years, there have been exciting developments — some of them very recent — which mostly resolve these issues. Incorporating these developments, MDL can be seen as a powerful extension of both penalized likelihood and Bayesian approaches, in which penalization functions and prior distributions are replaced by more general luckiness functions, average-case methodology is replaced by a more robust worst-case approach, and in which methods classically viewed as highly distinct, such as AIC versus BIC and cross-validation versus Bayes can, to some extent, be viewed from a unified perspective; as such, this paper should also be of interest to researchers working on the foundations of statistics and machine learning.

History of the field, recent advances and overview of this paper

MDL was introduced in 1978 by Jorma Rissanen in his paper Modeling by the Shortest Data Description. The paper coined the term MDL and introduced and analyzed the two-part code for parametric models. The two-part code is the simplest instance of a universal code or, equivalently, universal probability distribution, the cornerstone concept of MDL theory. MDL theory was greatly extended in the 1980s, when Rissanen published a sequence of ground-breaking papers at a remarkable pace, several of which introduced new types of universal distributions. It came to full blossom in the 1990s, with further major contributions from, primarily, Jorma Rissanen, Andrew Barron and Bin Yu, culminating in their overview paper³ and the collection⁵ with additional chapters by other essential contributors such as Kenji Yamanishi. The book G07 provides a more exhaustive treatment of this early work, including discussion of important precursors$/$alternatives to MDL such as MML,⁶ “ideal”, Kolmogorov complexity-based MDL⁷ and Solomonoff’s theory of induction.⁸ Universal distributions are still central to MDL. We introduce them in a concise yet self-contained way, including substantial underlying motivation, in Sec. 2, incorporating the extensions to and new insights into these basic building blocks that have been gathered over the last 10 years. These include more general formulations of arguably the most fundamental universal code, the Normalized Maximum Likelihood (NML) Distribution, including faster ways to calculate it as well. We devote a separate section to new universal codes, with quite pleasant properties for practical use, most notably the switch distribution (Sec. 3.1), which can be used for model selection combining almost the best of AIC and BIC; and the Reverse Information Projection (RIPr)-universal code (Sec. 3.3) specially geared to hypothesis testing with composite null hypotheses, leading to several advantages over classical Neyman–Pearson tests. In Sec. 4 we review recent developments on fast calculation of NML-type distributions for model selection for graphical models (Bayesian networks and the like), leading to methods which appear to be more robust in practice than the standard Bayesian ones. Recent extensions of MDL theory and practical implementations to latent variable and irregular models are treated in Sec. 5. Then, in Sec. 6 we review developments relating to consistency and convergence properties of MDL methods. First, while originally MDL estimation was formulated solely in terms of discretized estimators (reflecting the fact that coding always requires discretization), it has gradually become clear that a much larger class of estimators (including maximum likelihood for “simple” models, and, in some circumstances, the Lasso — see Example 4) can be viewed from an MDL perspective, and this becomes clearest if one investigates asymptotic convergence theorems relating to MDL. Second, it was found that MDL (and Bayes), without modification, can behave sub-optimally under misspecification, i.e., when all models under consideration are wrong, but some are useful — see Sec. 6.3. Third, very recently, it was shown how some of the surprising phenomena underlying the deep learning revolution in machine learning can be explained from an MDL-related perspective; we briefly review these developments in Sec. 6.4. Finally, we note that G07 presented many explicit open problems, most of which have been resolved — we mention throughout the text whenever a new development solved an old open problem, deferring some of the most technical issues to the Appendix.

Notational preliminaries

We shall mainly be concerned with statistical models (families of probability distributions) of the form ${ℳ}=\{p_{\theta }:\theta \in \mathrm{\Theta }\}$ parametrized by some set $\mathrm{\Theta }$ which is usually but not always a subset of Euclidean space; and families of models$\{{{ℳ}}_{\gamma }:\gamma \in \mathrm{\Gamma }\}$, where each ${{ℳ}}_{\gamma }=\{p_{\theta }:\theta \in {\mathrm{\Theta }}_{\gamma }\}$ is a statistical model, used to model the data $z^n:=(z_1,\dots ,z_n)$ with each $z_i\in {𝒵}$, for some outcome space ${𝒵}$. Each $p_{\theta }$ represents a probability density function (pdf) or probability mass function, defined on sequences of arbitrary length. With slight abuse of notation we also denote the corresponding probability distribution by $p_{\theta }$$($rather than the more common $P_{\theta })$. In the simple case that the data are i.i.d. according to each $p_{\theta }$ under consideration, we have $p_{\theta }(z^n)={\prod }_{i=1}^n{p}_{\theta }(z_i)$.

We denote the maximum likelihood (ML) estimator given the model ${ℳ}=\{p_{\theta }:\theta \in \mathrm{\Theta }\}$ by ${\widehat{\theta }}_{\mathrm{\text{ML}}}$, whenever it exists and is unique; the ML estimator relative to model ${{ℳ}}_{\gamma }$ is denoted by ${\widehat{\theta }}_{\mathrm{\text{ML}}|\gamma }$. We shall, purely for simplicity, generally assume its existence and uniqueness, although nearly all results can be generalized to the case where it does not. We use $\breve{\theta }$ to denote more general estimators, and ${\widehat{\theta }}_v$ to denote what we call the MDL estimator with luckiness function v, see (5).

2. The Fundamental Concept: Universal Modeling

MDL is best explained by starting with one of its prime applications, model comparison — we will generalize to prediction and estimation later, in Secs. 2.3 and 2.4. Assume then that we are given a finite or countably infinite collection of statistical models ${{ℳ}}_1,{{ℳ}}_2,\dots ,$ each consisting of a set of probability distributions. The fundamental idea of MDL is to associate each ${{ℳ}}_{\gamma }$ with a single distribution ${\bar{p}}_{\gamma }$, often called a universal distribution relative to ${{ℳ}}_{\gamma }$. We call the minus-log-likelihood $-log{\bar{p}}_{\gamma }(Z^n)$ the code length of data$Z^n$ under the universal code ${\bar{p}}_{\gamma }$. This terminology, and how MDL is related to coding (lossless compression of data), is briefly reviewed in Secs. 2.3 and 2.4; but a crucial observation at this point is that the main MDL ideas can be understood abstractly, without resorting to the code length interpretation. We also equip the model indices $\mathrm{\Gamma }:=\{1,2,\dots ,{\gamma }_{max}\}$$($where we allow $|\mathrm{\Gamma }|={\gamma }_{max}=\infty )$ with a distribution, say $\pi$; if the number of models to be compared is small (e.g., bounded independently of $n$ or at most a small polynomial in $n$), we can take $\pi$ to be uniform distribution — for large (exponential in $n$) and infinite $\mathrm{\Gamma }$, see Sec. 2.3 and Example 4. We then take, as our best explanation of the given data $z^n$, the model ${{ℳ}}_{\gamma }$ minimizing

or, equivalently, we maximize $\pi (\gamma ){\bar{p}}_{\gamma }(z^n)$; when $\pi$ is uniform this simply amounts to picking the $\gamma$ maximizing ${\bar{p}}_{\gamma }(z^n)$. (1) will later be generalized to $\pi$ that are not distributions but rather more general “luckiness functions” — see Sec. 2.3.

(1) The Bayesian universal distribution

The reader may recognize this as being formally equivalent to the standard Bayesian way of model selection, the Bayes factor method⁹ as long as the $\gamma$ are defined as Bayesian marginal distributions, i.e., for each $\gamma$, we set ${\bar{p}}_{\gamma }=p_{w_{\gamma }}^{\mathrm{\text{bayes}}}$, where

for some prior probability density $w_{\gamma }$ on the parameters in ${\mathrm{\Theta }}_{\gamma }$, which has to be supplied by the user. When $w_{\gamma }$ is clear from the context, we shall write ${\bar{p}}_{\gamma }^{{b}{a}{y}{e}{s}}$ rather than $p_{w_{\gamma }}^{{b}{a}{y}{e}{s}}$. Using Bayesian marginal distributions ${\bar{p}}^{{b}{a}{y}{e}{s}}$ is indeed one possible way to instantiate MDL model selection, but it is not the only way: MDL can also be based on other distributions such as ${\bar{p}}^{{N}{M}{L}}=p_v^{{N}{M}{L}}$ (depending on a function $v$), ${\bar{p}}^{{p}{r}{e}{q}}=p_{\breve{\theta }}^{{p}{r}{e}{q}}$ (depending on an estimator $\breve{\theta }$) and others; in general we add a bar to such distributions if the “parameter” $w,v$ or $\breve{\theta }$ is clear from the context. Before we continue with these other instantiations of ${\bar{p}}_{\gamma }$ we proceed with an example.

Example 1 (Bernoulli). Let ${ℳ}=\{p_{\theta }:\theta \in [0,1]\}$ represent the Bernoulli model, extended to $n$ outcomes by independence. We then have for each $z^n\in {\{0,1\}}^n$ that $p_{\theta }(z^n)={\theta }^{n_1}{(1-\theta )}^{n_0}$, where $n_1={\sum }_{i=1}^n{z}_i$ and $n_0=n-n_1$. Most standard prior distributions one encounters in the literature are beta priors, for which $w(\theta )\propto {\theta }^{\alpha }{(1-\theta )}^{\beta }$, so that $p_w^{{b}{a}{y}{e}{s}}(z^n)\propto \int {\theta }^{n_1+\alpha }{(1-\theta )}^{n_0+\beta }d\theta$. Note that $p_w^{{b}{a}{y}{e}{s}}$ is not itself an element of the Bernoulli model. One could use $p_w^{{b}{a}{y}{e}{s}}$ to compare the Bernoulli model, via (1), to, for example, a first-order Markov model, with Bayesian marginal likelihoods defined analogously. We shall say a lot more about the choice of prior below.

Example 2 (Gauss and general improper priors). A second example is the Gaussian location family ${{ℳ}}_{{g}{a}{u}{s}{s}}$ with fixed variance (say 1), in which ${𝒵}=ℝ$, and $p_{\theta }(z^n)\propto exp({\sum }_{i=1}^n{(z_i-\theta )}^2/2)$. A standard prior for such a model is the uniform prior, $w(\theta )=1$, which is improper (it does not integrate, hence does not define a probability distribution). Improper priors cannot be directly used in (2), and hence they cannot be directly used for model comparison as in (1) either. Still, we can use them in an indirect manner, as long as we are guaranteed that, for all ${{ℳ}}_{\gamma }$ under consideration, after some initial number of $m$ observations, the Bayesian posterior $w_{\gamma }(\theta |z^m)$ is proper. We can then replace $p_{w_{\gamma }}^{{b}{a}{y}{e}{s}}(z^n)$ in (2) by $p_{w_{\gamma }}^{{b}{a}{y}{e}{s}}(z_{m+1},\dots ,z_n|z^m):=\int p_{\theta }(z_{m+1},\dots ,z_n)w_{\gamma }(\theta |z^m)d\theta$. We extend all these conditional universal distributions to distributions on ${{𝒵}}^n$ by defining $p_{w_{\gamma }}^{{b}{a}{y}{e}{s}}(z_1,\dots ,z_n):=p_{w_{\gamma }}^{{b}{a}{y}{e}{s}}(z_{m+1},\dots ,z_n|z^m)p_0(z^m)$ for some distribution $p_0$ on ${{𝒵}}^m$ that is taken to be the same for all models ${{ℳ}}_{\gamma }$ under consideration. We can now use (1) again for model selection based on $p_{w_{\gamma }}^{{b}{a}{y}{e}{s}}(z_1,\dots ,z_n)$, where we note that the choice of $p_0$ plays no role in the minimization, which is equivalent to minimizing $-log\pi (\gamma )-log{p}_{w_{\gamma }}^{{b}{a}{y}{e}{s}}(z_{m+1},\dots ,z_n|z^m)$.

Now comes the crux of the story, which makes MDL, in the end, quite different from Bayes: defining the ${\bar{p}}_{\gamma }$ as in (2) is just one particular way to define an MDL universal distribution — but it is by no means the only one. There are several other ways, and some of them are sometimes preferable to the Bayesian choice. Here we list the most important ones.

(2) NML or Shtarkov¹⁰ distribution, and MDL estimators

This is perhaps the most fundamental universal distribution, leading also to the definition of an MDL estimator. In its general form, the NML distribution and “MDL estimators” depend on a function $v:\mathrm{\Theta }\to {ℝ}_0^{+}$. The definition is then given by

which is defined whenever the normalizing integral is finite. The logarithm of this integral is called the model complexity and is thus given by

$$\mathrm{\text{COMP}}({ℳ};v):=log\int \underset{\theta \in \mathrm{\Theta }}{max}(p_{\theta }(z^n)v(\theta ))d{z}^n\stackrel{(\mathrm{\text{if}}v\mathrm{\text{ constant}})}{=}log\int p_{{\widehat{\theta }}_{\mathrm{\text{ML}}}(z^n)}(z^n)d{z}^n.$$(4)

Here the integral is replaced by a sum for discrete data, and $max$ is replaced by $sup$ if necessary. This means that any function $v:\mathrm{\Theta }\to {ℝ}_0^{+}$ such that (4) is finite is allowed; we call any such $v$ a luckiness function, a terminology we explain later. Note that $v$ is not necessarily a probability density — it does not have to be integrable. For any luckiness function $v$, we define the MDL estimator based on v as

$${\widehat{\theta }}_v:=arg\underset{\theta \in \mathrm{\Theta }}{max}{p}_{\theta }(z^n)v(\theta )=arg\underset{\theta \in \mathrm{\Theta }}{min}\{-log{p}_{\theta }-[-logv(\theta )]\}.$$(5)

The v-MDL estimator is a penalized ML estimator, which coincides with the Bayes MAP estimator based on prior v whenever v is a probability density. Although this has only become clear gradually over the last 10 years, estimators of form (5) are the prime way of using MDL for estimation; there is, however, a second, “improper” way for estimating distributions within MDL though, see Sec.2.4. In practice, we will choose v that are sufficiently smooth so that, if the number of parameters is small relative to n, ${\widehat{\theta }}_v$ will usually be almost indistinguishable from the ML estimator ${\widehat{\theta }}_{\mathrm{\text{ML}}}$. COMP indeed measures something one could call a “complexity” — this is easiest to see if $v=1$, for then, if ${ℳ}$ contains just a single distribution, we must have $\mathrm{\text{COMP}}({ℳ},v)=0$, and the more distributions we add to ${ℳ}$, the larger $\mathrm{\text{COMP}}({ℳ},v)$ gets — this is explored further in Sec. 2.2.

Now suppose we have a collection of models ${{ℳ}}_{\gamma }$ indexed by finite $\mathrm{\Gamma }$ and we have specified luckiness functions $v_{\gamma }$ on ${\mathrm{\Theta }}_{\gamma }$ for each $\gamma \in \mathrm{\Gamma }$, and we pick a uniform distribution $\pi$ on $\mathrm{\Gamma }$. As can be seen from the above, if we base our model choice on NML, we pick the model minimizing

over $\gamma$, where $\mathrm{\text{COMP}}({{ℳ}}_{\gamma };v_{\gamma })$ is given by

$$\mathrm{\text{COMP}}({{ℳ}}_{\gamma };v_{\gamma })=log\int \underset{\theta \in {\mathrm{\Theta }}_{\gamma }}{max}(p_{\theta }(z^n)v_{\gamma }(\theta ))d{z}^n=log\int p_{{\widehat{\theta }}_{v_{\gamma }}(z^n)}(z^n)v_{\gamma }({\widehat{\theta }}_{v_{\gamma }}(z^n))d{z}^n.$$(7)

Thus, by (6), MDL incorporates a trade-off between goodness of fit and model complexity as measured by $\mathrm{\text{COMP}}$. Although the $n$-fold integral inside $\mathrm{\text{COMP}}$ looks daunting Suzuki and Yamanishi¹¹ show that in many cases (e.g., normal, Weibull–Laplace models) it can be evaluated explicitly with appropriate choice of $v$.

Originally, the NML distribution was defined by Shtarkov¹⁰ for the special case with $v\equiv 1$, leading to the rightmost definition in (3), and hence the term NML (in the modern version, perhaps “normalized penalized ML” would be more apt). This is also the version that Rissanen¹² advocated as embodying the purest form of the MDL Principle. However, the integral in (3) is ill-defined for just about every parametric model defined on unbounded outcome spaces $($such as $\mathbb{N},ℝ$ or ${ℝ}^{+})$, including the simple normal location family. Using nonuniform $v$ allows one to deal with such cases in a principled manner after all, see Sec. 2.5. For finite outcome spaces though, $v\equiv 1$ usually “works”, and (3) is well defined, as we illustrate for the Bernoulli model (see Sec. 4 for more examples).

Example 3 (Continuation of Example 1). For the Bernoulli model, ${\widehat{\theta }}_{\mathrm{\text{ml}}}(z^n)=n_1/n$ and $\mathrm{\text{COMP}}({ℳ},v)$ as in (7) with $v\equiv 1$ can be rewritten as $log{\sum }_{n_1=0}^n\left(\begin{array}{c}n\hfill \\ n_1\hfill \end{array}\right){(n_1/n)}^{n_1}{(n_0/n)}^{n_0}$, which, as we shall see in Sec. 2.2, is within a constant of $(1/2)logn$. As reviewed in that sub-section, the resulting $p_v^{{n}{m}{l}}$ is asymptotically (essentially) indistinguishable from $p_{w\text{J}}^{{b}{a}{y}{e}{s}}$ where the latter is equipped with Jeffreys’ prior, defined as $w_{{J}}(\theta )\propto \sqrt{|I(\theta )|}={\theta }^{-1/2}{(1-\theta )}^{-1/2}$, with $I(\theta )$ being the Fisher information at $\theta$.

(3) The two-part (sub-)distribution¹

Here one first discretizes $\mathrm{\Theta }$ to some countable subset $\ddot{\mathrm{\Theta }}$ which one equips with a probability mass function $w$; in contrast to the $v$ above, this function must sum to 1. One then considers

which is just a special case of (3). But since

$$\int \underset{\ddot{\theta }\in \ddot{\mathrm{\Theta }}}{max}{p}_{\ddot{\theta }}(z^n)w(\ddot{\theta })d{z}^n\le \int \sum _{\ddot{\theta }\in \ddot{\mathrm{\Theta }}}{p}_{\ddot{\theta }}(z^n)w(\ddot{\theta })d{z}^n=\sum _{\theta \in \ddot{\mathrm{\Theta }}}w(\theta )\left(\int p_{\theta }(z^n)d{z}^n\right)=1,$$(9)

we can approximate $p_w^{{n}{m}{l}}$ by the sub-distribution $p_w^{2-{p}}(z^n):={max}_{\ddot{\theta }\in \ddot{\mathrm{\Theta }}}{p}_{\ddot{\theta }}(z^n)w(\ddot{\theta })$. This “distribution” adds or integrates to something smaller than 1. This can be incorporated into the general story by imagining that $p_w^{2-{p}}$ puts its remaining mass on a special outcome, say “$◇$”, which in reality will never occur (while sub-distributions are thus “allowed”, measures that add up to something larger than 1 have no place in MDL). The two-part distribution $p_w^{2-{p}}$ is historically the oldest universal distribution. The fact that it can be considered a special case of NML has only become fully clear very recently¹³; in that same paper, an even more general formulation of (3) is given that has all Bayesian, two-part and NML distributions as special cases. Despite its age, the two-part code is still important in practice, as we explain in Sec. 2.3.

(4) The prequential plug-in distribution^14,15

Here, one first takes any reasonable estimator $\breve{\theta }$ for the given model ${ℳ}$. One then defines

where for i.i.d. models, the probability inside the product simplifies to $p_{\breve{\theta }(z^{i-1})}(z_i)$. For the normal location family, one could simply use the ML estimator: $\breve{\theta }(z^m):={\widehat{\theta }}_{\mathrm{\text{ML}}}(z^m)={\sum }_{j=1}^m{z}_j/m$. With discrete data though, the ML estimator should be avoided, since then one of the factors in (10) could easily become 0, making the product 0, so that the model for which $p_{\breve{\theta }}^{{p}{r}{e}{q}}$ is defined can never “win” the model selection contest even if most other factors in the product (10) are close to 1. Instead, one can use a slightly “smoothed” ML estimate (a natural choice for $\breve{\theta }$ is to take an MDL estimator for some $v$ as in (5), but this is not required). For example, in the Bernoulli model, one might take $\breve{\theta }(z^m)=(m_1+(1/2))/(m+1)$, where $m_1={\sum }_{i=1}^m{z}_i$. With this particular choice, $p_{\breve{\theta }}^{{p}{r}{e}{q}}$ turns out to coincide exactly with $p_{w_{\text{J}}}^{{b}{a}{y}{e}{s}}$ with Jeffreys’ prior $w_{{J}}$. Such a precise correspondence between ${\bar{p}}^{{p}{r}{e}{q}}$ and ${\bar{p}}^{{b}{a}{y}{e}{s}}$ is a special property of the Bernoulli and multinomial models though; with other models, the two distributions can usually be made to behave similarly, but not identically. The rationale for using ${\bar{p}}^{{p}{r}{e}{q}}$ is described in Sec. 2.4. In Sec. 3.2.1 we will say a bit more about hybrids between prequential plug-in and Bayes (the flattened leader distribution) and between prequential and NML (sequential NML).

Except for the just mentioned “hybrids”, these first four universal distributions were all brought into MDL theory by Rissanen; they are extensively treated by G07, in which one chapter is devoted to each, and to which we refer for details. The following two are much more recent.

(5) The switch distribution $\bar{p}$^switch(Ref. 16)

In a particular type of nested model selection, this universal distribution behaves arguably better than the other ones. It will be treated in detail in Sec. 3.1.

(6) Universal distributions $\bar{p}$^ripr based on the Reverse Information Projection

These universal distributions¹⁷ lead to improved error bounds and optional stopping behavior in hypothesis testing and allow one to forge a connection with group-invariant Bayes factor methods; see Sec. 3.3.

2.1. Motivation

We first give a very high-level motivation that avoids direct use of data compression arguments. For readers interested in data compression, Sec. 2.3 does make a high-level connection, but for more extensive material we refer to G07. We do, in Sec. 2.4, give a more detailed motivation in predictive terms, and, in Sec. 6, we shall review mathematical results indicating that MDL methods are typically consistent and enjoy fast rates of convergence, providing an additional motivation in itself.

Consider then models ${{ℳ}}_{\gamma }$, where for simplicity we assume discrete data, and let ${\widehat{\theta }}_{\mathrm{\text{ML}}|\gamma }$ be the maximum likelihood estimator within ${{ℳ}}_{\gamma }$. Define “the fit of the model to the data” in the standard way, as $F_{\gamma }(z^n):=p_{{\widehat{\theta }}_{\mathrm{\text{ML}}|\gamma }(z^n)}(z^n)$, the likelihood assigned to the data by the best-fitting distribution within the model. Now if we enlarge the model ${{ℳ}}_{\gamma }$, i.e., by adding several distributions to it, $F_{\gamma }(z^n)$ can only increase; and if we make ${{ℳ}}_{\gamma }$ big enough such that for each $z^n$, it contains a distribution $p$ with $p(z^n)=1$, we can even have $F_{\gamma }(z^n)=1$ on all data. If we simply picked the $\gamma$ maximizing $F_{\gamma }(z^n)$, we would be prone to severe overfitting. For example, if models are nested, then, except for very special data, we would automatically pick the largest one.

As we have seen, a central MDL idea is to instead associate each model ${{ℳ}}_{\gamma }$ with a single corresponding distribution ${\bar{p}}_{\gamma }$, i.e., we set $F_{\gamma }(z^n):={\bar{p}}_{\gamma }(z^n)$. Then the total probability mass on all potential outcomes $z^n$ cannot be larger than 1, which makes it impossible to assign overly high fit $F_{\gamma }(z^n)$ to overly many data sequences: no matter what distribution ${\bar{p}}_{\gamma }$ we chose, we must now have ${\sum }_{z^n}{F}_{\gamma }(z^n)=1$, so a good fit on some $z^n$ necessarily implies a worse fit on others, and we will not select a model simply because it accidentally contained some distribution that fitted our data very well — thus, measuring fit by a distribution ${\bar{p}}_{\gamma }$ instead of $F_{\gamma }$ inherently prevents overfitting. This argument to measure fit relative to a model with a single $\bar{p}$ is similar to Bayesian Occam’s Razor arguments¹⁸ used to motivate the Bayes factor; the crucial difference is that we do not restrict ourselves to ${\bar{p}}_{\gamma }$ of the form (2); inspecting the “Bayesian” Occam argument, there is, indeed, nothing in there which forces us to use distributions of Bayesian form.

The next step is thus to decide which $\bar{p}$ are best associated with a given ${ℳ}$. To this end, we define the fitness ratio for data $z^n$ as

where $v:\mathrm{\Theta }\to {ℝ}_0^{+}$ is a nonnegative function. To get a feeling for (11), it is best to first focus on the case with $v(\theta )\equiv 1$; it then reduces to We next postulate that a good choice for $\bar{p}$ relative to the given model is one in which $\mathrm{\text{FR}}(\bar{p},n)$tends to be as large as possible. The rationale is that, overfitting having already been taken care of by picking some $\bar{p}$ that is a probability measure (integrates to 1), it makes sense to take a $\bar{p}$ whose fit to data (as measured in terms of likelihood) is proportional to the fit to data of the best-fitting distribution in ${ℳ}$: whenever some distribution in the model ${ℳ}$ fits the data $z^n$ well, the likelihood $\bar{p}(z^n)$ should be high as well. One way to make “$\mathrm{\text{FR}}$ tends to be large” precise is by requiring it to be as large as possible in the worst-case, i.e., we want to pick the $\bar{p}$ achieving where the maximum is over all probability distributions over samples of length $n$. It turns out that this maximin problem has a solution if and only if the complexity (4) is finite; and if it is fine, the unique solution is given by setting $\bar{p}={\bar{p}}^{{N}{M}{L}}$, with ${\bar{p}}^{{N}{M}{L}}$ given by (3). The NML distribution thus has a special status as the most robust choice of universal $\bar{p}$ — even though $\bar{p}$ is itself a probability distribution, it meaningfully assesses fit in the worst-case over all possible distributions, and its interpretation does not require one to assume that the model ${ℳ}$ is “true” in any sense. The nicest sub-case is the one with $v(\theta )\equiv 1$, since then all distributions within the model ${ℳ}$ are treated on exactly the same footing; no data or distribution is intrinsically preferred over any other one. Unfortunately, for most popular models with infinite ${𝒵}$, when taking $v(\theta )\equiv 1$, (13) usually has no solution since the integral $\int p_{{\widehat{\theta }}_{\mathrm{\text{ML}}}(z^n)}(z^n)d{z}^n$ diverges for such models, making the complexity (4) infinite. For all sufficiently “regular” models (curved exponential families, see below), this problem can invariably be solved by restricting $\mathrm{\Theta }$ to a bounded subset of its own interior — one can show that the complexity (4) is finite with $v\equiv 1$, and thus (13) has a solution given by (3) if ${\widehat{\theta }}_{\mathrm{\text{ml}}}$ is restricted to a suitably bounded set. Yet, restricting $\mathrm{\Theta }$ to a bounded subset of itself is not satisfactory, since it is unclear where exactly to put the boundaries. It is more natural to introduce a nonuniform $v$, which can invariably be chosen so that the complexity (4) is finite and thus (13) has a solution — more on choosing $v$ at the end of Sec. 2.4.

A few remarks concerning this high-level motivation of MDL procedures are in order.

(1)	It is clear that, by requiring $\mathrm{\text{FR}}$ to add to 1, we will be less prone to overfitting than by setting it simply to $p_{{\widehat{\theta }}_{\mathrm{\text{ML}}}(z^n)}(z^n)$; whether the requirement to add (at most) to 1, making $\mathrm{\text{FR}}$ essentially a probability density function, is a clever way to avoid overfitting (leading to good results in practice) is not clear yet. For this, we need additional arguments, which we very briefly review. First, the sum-to-1 requirement is the only choice for which the procedure can be interpreted as selecting the model which minimizes code length of the data (the original interpretation of MDL); second, it is the only choice which has a predictive interpretation, which we review in Sec. 2.4 below; third, it is the only choice under which time-tested Bayesian methods fit into the picture; and fourth, with this choice we get desirable frequentist statistical properties such as consistency and convergence rates, see Sec. 6.
(2)	The motivation above only applies to the NML universal distributions. How about the other five types? Originally, in the pure MDL approach mainly due to Rissanen, the NML was viewed as the optimal choice per se; other $\bar{p}$ should be used only for pragmatic reasons, such as them being easier to calculate. One would then design them so as to be as close as possible to the NML distributions in terms of the fitness ratio they achieve. In the following sub-section we show that all six of them satisfy the same MDL$/$BIC asymptotics, meaning that their fitness ratio is never smaller than a constant factor of the NML one, either again in the worst-case over all $z^n$ or in some weaker expectation sense. Thus, they are all “kind of ok” in a rather weak sense, and in practice one would simply revert to the one that is closest to NML and still usable in practice; with the Bayesian ${\bar{p}}^{{b}{a}{y}{e}{s}}$, as we shall see, one can even get arbitrarily close to NML as $n$ gets larger. This classical story notwithstanding, it has become more and more apparent that in practice one sometimes wants or needs properties of model selection methods that are not guaranteed by NML — such as near-optimal predictions of future data or strong frequentist Type-I error guarantees. This translates itself into universal codes ${\bar{p}}^{{s}{w}{i}{t}{c}{h}}$ and ${\bar{p}}^{{r}{i}{p}{r}}$ that, for some special sequences, achieve much higher fitness ratio than ${\bar{p}}^{{n}{m}{l}}$, while for all sequences having only very slightly smaller fitness ratio. This more recent and pragmatic way of MDL is briefly reviewed in Secs. 3.1 and 3.3. This raises the question how we should define a universal distribution: what choices for ${\bar{p}}_{\gamma }$ are still “universal” (and define an MDL method) and what choices are not? Informally, every distribution ${\bar{p}}_{\gamma }$ that for no $z^n\in {{𝒵}}^n$ has ${\bar{p}}_{\gamma }(z^n)\ll {\bar{p}}_{\gamma }^{{n}{m}{l}}(z^n)$ is “universal” relative to ${{ℳ}}_{\gamma }$. For parametric models such as exponential families, the “$\ll$” is partially formalized by requiring that at the very least, they should satisfy (14) below (G07 is much more precise on this).
(3)	Third, we have not yet said how one should choose the “luckiness function” $v$ — and one needs to make a choice to apply MDL in practice. The interpretation of $v$ is closely tied to the predictive interpretation of MDL, and hence we postpone this issue to the end of Sec. 2.4.
(4)	Fourth, the motivation so far is incomplete — we still need to explain why and how to incorporate the distribution $\pi$ on model index $\gamma$. This is done in Sec. 2.3 below.

2.2. Asymptotic expansions

Now let $\bar{p}$ be defined relative to a single parametric model ${ℳ}$. It turns out that all universal codes we mentioned have in common that, for “sufficiently regular” $k$-dimensional parametric models, the log-likelihood for given data $z^n$ satisfies the following celebrated asymptotics, often called the MDL or BIC expansion: for all “sufficiently regular” data sequences $z_1,z_2,\dots$, there exists a constant $C\in ℝ$ independent of $n$ such that for all $n$,

For ${\bar{p}}^{{n}{m}{l}}$ and ${\bar{p}}^{{b}{a}{y}{e}{s}}$, this holds for any choice of luckiness function $v$ and prior $w$ that is continuous and strictly positive on the parameter space $\mathrm{\Theta }$. For ${\bar{p}}_w^{2-{p}}$, this holds for clever choices of the discretization $\ddot{\mathrm{\Theta }}$ and the probability mass function $w$; for $p_{\breve{\theta }}^{{p}{r}{e}{q}}$, this holds in a weaker expectation sense (see Sec. 3.2.1), as long as $\breve{\theta }$ is a suitably smoothed version of the ML estimator. Essentially, “sufficiently regular” parametric models are all exponential families (such as Bernoulli, multinomial, normal, gamma, beta, Poisson, etc.) and curved exponential families; corresponding results also hold for regression with (generalized) linear models. “Sufficiently regular data” are all sequences for which there is an INECCSI subset ${\mathrm{\Theta }}_0$ of the parameter space $\mathrm{\Theta }$ such that, for all large $n$, the ML estimator of the sequence lies within ${\mathrm{\Theta }}_0$. Here INECCSI stands for a set whose Interior is NonEmpty and whose Closure is a Compact Subset of the Interior of $\mathrm{\Theta }$. Essentially, this is any bounded subset of the same dimensionality as $\mathrm{\Theta }$ that does not touch the boundaries of $\mathrm{\Theta }$ itself; in the Bernoulli example, it would be any set of the form $[ϵ,1-ϵ]$ for $ϵ>0$. For all universal distributions considered except ${\bar{p}}^{{p}{r}{e}{q}}$, as long as appropriate priors$/$estimators$/$luckiness functions are used, (14) will hold uniformly for all sequences in any INECCSI subset ${\mathrm{\Theta }}_0$, but the constant $C$ may grow larger if we replace ${\mathrm{\Theta }}_0$ by a strictly larger INECCSI subset ${\mathrm{\Theta }}_0^{\prime }$ with ${\mathrm{\Theta }}_0⊊︀{\mathrm{\Theta }}_0^{\prime }⊊︀\mathrm{\Theta }$. (for ${\bar{p}}^{{p}{r}{e}{q}}$ see Sec. 3.2.1). For the first four universal distributions, the inequality is actually equality up to a constant — (14) also holds with $\le$ replaced by $\ge$, for a different constant. For the switch distribution ${\bar{p}}^{{s}{w}{i}{t}{c}{h}}$, however, the left-hand side will be significantly smaller for a small but important subset of possible data sequences. Finally, since (14) thus also holds with $\bar{p}={\bar{p}}^{{N}{M}{L}}$ and $\le$ replaced by $\ge$, exponentiating (14), we see that, if one restricts the minimum in (13) to all such “sufficiently regular” $z^n$, $\mathrm{\text{FR}}({\bar{p}}^u(z^n))$ is guaranteed to be within a constant (independent of $n$) factor of the optimal $\mathrm{\text{FR}}({\bar{p}}^{{N}{M}{L}},z^n)$, for $u\in \{{b}{a}{y}{e}{s},2-{p},{p}{r}{e}{q},{s}{w}{i}{t}{c}{h},{r}{i}{p}{r}\}$.

The NML$/$COMP expansion and the Jeffreys (Fisher information) integral

For the case that the model ${ℳ}=\{p_{\theta }:\theta \in \mathrm{\Theta }\}$ is a $k$-dimensional exponential family and $\bar{p}$ is the NML or Bayes distribution, we can be significantly more precise and evaluate the constant $C$ in (14) up to $o(1)$: we get, under some weak additional regularity conditions on ${ℳ}$ and $v$,

where $k$ is the dimension of the model, $|I(\theta )|$ is the determinant of the $k\times k$ Fisher information matrix at parameter $\theta$, the integral is over the parameter space $\mathrm{\Theta }$ and the remainder term $o(1)$ vanishes as the sample size grows unbounded. This was first shown (essentially) by Rissanen,¹² for the case that $\mathrm{\Theta }$ is restricted to an INECCSI subset of the full parameter space (so that ${\bar{p}}^{{N}{M}{L}}$ with $v\equiv 1$ is defined), and $v\equiv 1$. For this uniform $v$ case, Myung et al.¹⁹ gave a differential geometric interpretation of the Fisher information term, relating it to an intrinsic “volume” of the parameter space. The general result for nonuniform $v$, and without INECCSI restrictions, was very recently shown in a breakthrough paper by Suzuki and Yamanishi,¹¹ solving Open Problem 6 from G07.

Analogously to ${\bar{p}}^{{n}{m}{l}}$ (in fact much easier mathematically), we can expand ${\bar{p}}^{{b}{a}{y}{e}{s}}$ using a classical Laplace approximation; under the same conditions as before, with now the additional restriction that there exists an arbitrary INECCSI subset ${\mathrm{\Theta }}_0$ of $\mathrm{\Theta }$ such that for all large $n$, the data have ML estimator within ${\mathrm{\Theta }}_0$, we find that

From (15) and (16) we see that, if the generalized Jeffreys integral$\int v(\theta )\cdot \sqrt{|I(\theta )|}d\theta$ is finite (see the Appendix), then there is a special choice of prior $w$, the generalized Jeffreys’ prior, with $w(\theta )=v(\theta )\sqrt{|I(\theta )|}/\int v(\theta )\sqrt{|I(\theta )|}d\theta$, under which $-log{\bar{p}}^{{N}{M}{L}}(z^n)$ and $-log{\bar{p}}^{{b}{a}{y}{e}{s}}(z^n)$ do not just coincide up to $O(1)$, but become asymptotically indistinguishable. If $v\equiv 1$, this $w$ coincides with the well-known Jeffreys’ prior $w_{{J}}$ popular in Bayesian inference; the special case of this prior for the Bernoulli model was encountered in Example 1. Thus, the Bayesian universal distribution with the (generalized) Jeffreys’ prior can be a very good alternative of $p_v^{{n}{m}{l}}$.

2.3. Unifying model selection and estimation

Suppose we are given a countable collection of models $\{{{ℳ}}_{\gamma }:\gamma \in \mathrm{\Gamma }\}$. Recall that the basic idea above was to associate each individual model ${{ℳ}}_{\gamma }$ with a single distribution ${\bar{p}}_{\gamma }$. It seems reasonable to do the same at the level of “meta-parameters” $\gamma$: we set $\overline{{ℳ}}:=\{{\bar{p}}_{\gamma }:\gamma \in \mathrm{\Gamma }\}$ and in complete analogy to (3), we define the meta-universal distribution

for some nonnegative weight function $\pi$ on $\mathrm{\Gamma }$. It then makes sense to postulate that the best sub-model ${{ℳ}}_{\gamma }$ for the given data $z^n$ is given by the $\gamma$ achieving the maximum in (17). Note that for determining this maximum, the denominator in (17) plays no role.

Let us assume that all ${\bar{p}}_{\gamma }$ have already been defined. Then we can use any $\pi$ such that the overarching $p_{\pi }^{{N}{M}{L}}$ in (17) exists. We can now formulate a general MDL Principle for model selection: we start with a (potentially huge) set of candidate distributions ${{ℳ}}_{{f}{u}{l}{l}}$. We next carve up ${{ℳ}}_{{f}{u}{l}{l}}$ into interesting sub-models ${{ℳ}}_{\gamma }$ with $\gamma \in \mathrm{\Gamma }$, so that ${\bigcup }_{\gamma \in \mathrm{\Gamma }}{{ℳ}}_{\gamma }={{ℳ}}_{{f}{u}{l}{l}}$. We then associate each ${{ℳ}}_{\gamma }$ with a universal distribution ${\bar{p}}_{\gamma }$, and we equip $\overline{{ℳ}}$ as defined above with luckiness function $\pi$ [note that ${{ℳ}}_{{f}{u}{l}{l}}$, a countable union of (usually) uncountable sets, consists of all distributions under consideration, while $\overline{{ℳ}}$ is a countable set]. We then base the selection of a sub-model ${{ℳ}}_{\gamma }$ on (17). What we earlier called the “general MDL Principle” underneath (1) was the special case in which $\sum \pi (\gamma )=1$, i.e., $\pi$ is a probability mass function. Via (9) we see that for any such probability mass function $\pi$, the denominator in (17) is well-defined, hence $\pi$ is a valid luckiness function.

Now consider the special case in which every ${\bar{p}}_{\gamma }$ is chosen to be an NML distribution $p_{v_{\gamma }}^{{n}{m}{l}}$ for some luckiness functions $v_{\gamma }$. We take some function ${\pi }^{\prime }:\mathrm{\Gamma }\to {ℝ}_0^{+}$ (which we will relate to the $\pi$ above later on) and we set, for $\theta \in {\mathrm{\Theta }}_{\gamma }$, $v_{{f}{u}{l}{l}}(\gamma ,\theta ):={\pi }^{\prime }(\gamma )v_{\gamma }(\theta )$. We can use $p_{v_{{f}{u}{l}{l}}}^{{N}{M}{L}}$ for parameter estimation on the joint parameters $(\gamma ,\theta )$ just as we did earlier for parametric models, by using the MDL estimator $(\gamma ,{\widehat{{\theta }_{\gamma })}}_{v_{{f}{u}{l}{l}}}$ picking the $(\gamma ,{\theta }_{\gamma })$ minimizing, over $\gamma \in \mathrm{\Gamma },{\theta }_{\gamma }\in {\mathrm{\Theta }}_{\gamma }$,

where again $\mathrm{\text{COMP}}({{ℳ}}_{{f}{u}{l}{l}},v_{{f}{u}{l}{l}})$ plays no role in the minimization. This MDL estimator really combines model selection (estimation of $\gamma$) and parametric estimation (estimation of ${\theta }_{\gamma }$). If we now define $\pi (\gamma ):={\pi }^{\prime }(\gamma )/exp(\mathrm{\text{COMP}}({{ℳ}}_{\gamma },v_{\gamma }))$, we find that $p_{\pi }^{{N}{M}{L}}$ defined relative to model $\overline{{ℳ}}$ as in (17) is equal to $p_{v_{{f}{u}{l}{l}}}^{{n}{m}{l}}$ defined relative to the full union of models ${{ℳ}}_{{f}{u}{l}{l}}$, and the $\gamma$ achieving the maximum in (17) coincides with the $\gamma$ minimizing (18). This indicates that model selection and estimation is really the same thing with MDL: if we are given a single parametric model ${{ℳ}}_{\gamma }$ with luckiness $v_{\gamma }$, we pick the $\theta$ minimizing the first two terms in (18) for fixed $\gamma$; if we are interested in both $\theta$ and $\gamma$, we minimize over all terms; and if we are only interested in $\gamma$, we pick the $\gamma$ achieving the maximum in (17), which, by construction, will give us the same $\gamma$ as the joint minimization over (18).

Two-part versus one-part codes: The role of data compression

In the oldest (1978) version of MDL, only two-part codes on countable sets were used: the minimum over $\theta \in {\mathrm{\Theta }}_{\gamma }$ was taken over a discretized grid ${\ddot{\mathrm{\Theta }}}_{\gamma }$ and $v_{\gamma }$ was a probability mass function over this grid; then for all $\gamma$, $\mathrm{\text{COMP}}({{ℳ}}_{\gamma },v_{\gamma })\le 0$ and $\mathrm{\text{COMP}}({ℳ},v)\le 0$ [see (9)] and they were both approximated by 0. From the Kraft inequality²⁰ we see where the name “two-part code” comes from: this inequality says that for every probability mass function $\pi$ on a countable set ${𝒜}$, there exists a lossless code such that for all $a\in {𝒜}$, the number of bits needed to encode $a$, is given by $-log\pi (a)$. Thus the resulting method can be interpreted as picking the $(\ddot{\theta },\ddot{\gamma })$ minimizing the two-stage code length of the data, where first the parameters $(\theta ,\gamma )$ are encoded using $-log\pi (\gamma )-log{v}_{\gamma }(\theta )$ bits, and then $z^n$ is encoded “with the help of $\gamma$”, using $-log{p}_{{\theta }_{\gamma }}(z^n)$ bits [in fact, the encoding of $(\gamma ,\theta )$ itself has two sub-stages here so we really have a two-part code where the first part itself has two parts as well].

The discretization involved in using a probability mass function$/$code for continuous-valued $\theta$ makes things (unnecessarily, as was gradually discovered over the last 30 years) very complicated in general. Also, if one combines the choice of $\theta$ with the choice of $\gamma$, the approximation of $\mathrm{\text{COMP}}({{ℳ}}_{\gamma },v_{\gamma })$ as 0 introduces some sub-optimalities. Thus, one would like to code the data in a way that avoids these two issues. It turns out that this can be achieved by replacing two-part by one-part codes for the data, namely, to use codes with length $-log{p}_v^{{N}{M}{L}}(z^n)$: assuming for simplicity that data are discrete, the same Kraft inequality implies that there must also be a code, directly defined on $z^n$, which achieves code length for each $z^n$ given by $-log{\bar{p}}^{{N}{M}{L}}(z^n)$. Thus, even though for general luckiness functions $v$ this code length cannot be decomposed into two sub-code lengths, it remains a valid code length and the name MDL for the resulting procedure remains, we feel, justified. In the past, it was sometimes thought by some MDL fans that two-part codes on countable sets would somehow lead to inherently better estimates ${\widehat{\theta }}_w$ than estimators ${\widehat{\theta }}_v$ for general luckiness functions as in (5). However, after 30 years it turned out there is nothing either conceptually or mathematically that indicates the need for two-part codes and countable sets: for any luckiness function $v$, the resulting procedure has a code length interpretation, and Grünwald and Mehta¹³ show that all consistency and convergence results that hold for two-part estimators also hold for general MDL estimators (Sec. 6.1) — thus invalidating the conjecture in Open Problem 13 of G07 that postulated a special status for luckiness functions $v$ that are probability mass functions on countable sets. For the luckiness function $\pi$ on the discrete structure $\mathrm{\Gamma }$, however, it is quite reasonable to choose a probability mass function: no probability mass is wasted [since the denominator in (17) plays no role in choosing $\gamma$], and designing $\pi$ by thinking about the code lengths $-log\pi (\gamma )$ comes very naturally, as the following example illustrates.

Example 4 (Variable selection: $L_1$-versus $L_0$-penalties). Suppose that each data point $Z_i=(X_i,Y_i)$ where $Y_i$ denotes the variable to be predicted and $X_i=(X_{i1},\dots ,X_{im})\in {ℝ}^m$ is a vector of covariates or “features” that may or may not help for predicting $Y_i$. We consider a linear model ${ℳ}=\{p_{\beta }:\beta \in {ℝ}^m\}$, decomposed into sub-models ${{ℳ}}_{\gamma }$ expressing,

where ${ϵ}_1,{ϵ}_2,\dots$ represent zero-mean i.i.d. $N(0,{\sigma }^2)$ denotes the normally distributed noise and $\gamma =({\gamma }_1,\dots ,{\gamma }_m)\in {\{0,1\}}^m$ is a binary vector indicating which variables are helpful for predicting the $Z_i$. Thus, if $\gamma$ has $k$ zero components, then ${{ℳ}}_{\gamma }$ is effectively an ($m-k$)-dimensional model. Our task is to learn, from the data, the vector ${\gamma }^{\ast }$ indicating which variables are truly relevant, and$/$or such that predictions of new $Y$ given new $X$ based on ${{ℳ}}_{{\gamma }^{\ast }}$ are as good as possible.

In light of the above, a straightforward way to use MDL here is to pick the $\gamma$ minimizing

where we refer to Fig. 14.2 in G07 for an explanation why we can condition on $x^n$ here. $v_{\gamma }$ in (19) is an appropriately chosen luckiness function and $\pi$ is really a probability mass function, such that $L_{\pi }(\gamma ):=-log\pi (\gamma )$ can be interpreted as the number of bits needed to encode $\gamma$ using a particular code. In terms of coding, a natural choice for such a code would be to first encode the number of nonzero components $k_{\gamma }$ in $\gamma$ using a uniform code (that assigns equal code length to all possibilities). Since $0\le k_{\gamma }\le m$, there are $m+1$ possibilities, so this takes ${log}_2(m+1)$ bits. In a second stage, one encodes the location of these components. There are $\left(\begin{array}{c}m\hfill \\ k_{\gamma }\hfill \end{array}\right)$ possibilities here, so this takes ${log}_2\left(\begin{array}{c}m\hfill \\ k_{\gamma }\hfill \end{array}\right)$ bits using a uniform code. All in all, one needs “nits” (bits re-expressed in terms of natural logarithm $log$) to encode $\gamma$. Then (20) can be written as $-log\pi (\gamma )$, with $\pi (\gamma )=1/((m+1)\cdot (\begin{array}{c}m\hfill \\ k_{\gamma }\hfill \end{array}\left)\right)$, which, as predicted by the Kraft inequality, sums to 1 over $\gamma \in {\{0,1\}}^m$.

As to the left part of the code length (19), if the variance ${\sigma }^2$ is known, a natural luckiness function $v_{\gamma }$ to use is a ($m-k_{\gamma }$)-dimensional Gaussian with mean 0 and variance ${\sigma }^2\mathrm{\Sigma }$ for some (usually diagonal) covariance matrix ${\mathrm{\Sigma }}_{\gamma }$. This gives (see Chap. 12 of G07)

where $X=(X_1,\dots ,X_n)$ and $|\cdot |$ stands for determinant. We thus end up finding the $\gamma$ minimizing the sum of (21) and (20). If, as here, the noise is normal with fixed variance, then for this choice of luckiness function, $p_{v_{\gamma }}^{{n}{m}{l}}(y^n|x^n)$ actually coincides with $p_{w_{\gamma }}^{{b}{a}{y}{e}{s}}(y^n|x^n)$ for a particular prior $w_{\gamma }$, thus one has a Bayesian interpretation as well (Bartlett et al.²¹ show that such a precise correspondence between NML and Bayes only holds in quite special cases, see the Appendix). If the variance ${\sigma }^2$ is unknown, one can treat it as a nuisance parameter and equip it with the improper Haar prior, leading to a modification of the formula above; see Example 6. Even if the noise is not known to be normally distributed, one can often still use the above method — pretending the noise to be normally distributed and accepting that one is misspecified — by varying the learning rate, as briefly explored in Sec. 6.3.

Note that the code$/$prior we used here induces sparsity: if there exists a $\gamma$ with mostly zero components that already fits the data quite well, we will tend to select it, since, for $k\ll n$, $log\left(\begin{array}{c}n\hfill \\ k\hfill \end{array}\right)$ increases approximately linearly in $k$. That does not mean that we necessarily believe that the “truth” is sparse — it just expresses that we hope that we can already make reasonably good predictions with a small number of features.

An alternative, and very popular, approach to this problem is the celebrated Lasso,²² in which we consider only the full model ${ℳ}={{ℳ}}_{(1,1,\dots ,1)}$ and we pick the $\beta \in {ℝ}^m$ minimizing

for some regularization parameter $\lambda >0$ [the factor ${\sigma }^2$ plays no role in the minimization; it is incorporated only to facilitate comparison with (21)]. It is known this will tend to select $\beta$ with many zero components, thus also inducing sparsity, and it can be implemented computationally much more efficiently than the two-part approach sketched above, effectively replacing $L_0$-penalties by $L_1$-penalties. With our “modern” view of MDL, this can be thought of as a form of MDL too, where we simply impose the luckiness function $v(\beta )=exp(-(\lambda /{\sigma }^2)\times {\sum }_{j=1}^m|{\beta }_j|)$ and use the estimator ${\widehat{\theta }}_v$ given by (5). The luckiness function $v$ depends on $\lambda$; the optimal choice of $\lambda$ is then once again related to the optimal learning rate; see Sec. 6.3. Finally, we note that there exists a third MDL approach one can use here: one starts out with an NML approach similar to (19) but then performs a continuous relaxation of the resulting optimization problem; the resulting “relaxed” NML criterion is then once again tractable and similar to an $L_1$-optimization problem such as (22); this approach has been described by Miyaguchi and Yamanishi²³ who extend the idea to group Lasso and other settings.

2.4. Log-loss prediction and universal distributions

Now consider the simple case again with a finite set of models $\{{{ℳ}}_{\gamma }:\gamma \in \mathrm{\Gamma }\}$ where $\mathrm{\Gamma }$ is small compared to $n$ and we use the uniform prior $\pi$, picking the $\gamma$ maximizing ${\bar{p}}_{\gamma }$. It was the fundamental insight of Rissanen¹⁴ and Dawid¹⁵ that such model choice by maximizing ${\bar{p}}_{\gamma }(z^n)$ for a single distribution ${\bar{p}}_{\gamma }$ can be motivated in a different way as well — in essence, it selects the model with the best predictive performance on unseen data. This approach shows that MDL is quite similar in spirit to cross-validation, the main difference with leave-one-out cross-validation being that the cross in cross-validation is replaced by a forward and that the loss function used to measure prediction error is restricted to be the logarithmic score, also commonly known as log-loss (which, however, is often used in cross-validation as well).

Formally, the log-loss of predicting a single outcome $z\in {𝒵}$ with a distribution $p$ is defined as $-logp(z)$: the larger the probability density, the smaller the loss. If one predicts a sequence of $n$ outcomes $z^n=(z_1,\dots ,z_n)$ with $n$ predictions $p_1,p_2,\dots ,p_n$, then the cumulative log-loss is defined as the sum of the individual losses: ${\sum }_{i=1}^n-log{p}_i(z_i)$.

Now, if we adopt a probabilistic world view and represent our beliefs about $z^n$ by a probability distribution $\bar{p}$, then the obvious way to make sequential predictions is to set $p_i:=\bar{p}(Z_i=\cdot |z^{i-1})$, so that $-log{p}_i(z_i)=-log\bar{p}(z_i|z^{i-1})$. For arbitrary probability distributions, we have, by the formula for conditional probability: for all $z^n\in {{𝒵}}^n$, $p(z^n)={\prod }_{i=1}^{n}p(z_i|z^{i-1})$. Taking logarithms gives

In other words, for every possible sequence, the cumulative log-loss obtained by sequentially predicting$z_i$ based on the previously observed data$z^{i-1}$ is equal to the minus-log-likelihood.

Conversely, if we are given an arbitrary sequential prediction strategy$\bar{s}$ which when input with a sequence $z^{i-1}$ of arbitrary length $i-1$ outputs a prediction for the next outcome $z_i$ in the form of a probability distribution ${\bar{s}}_{z^{i-1}}$ on ${𝒵}$, we can define$\bar{p}(z_i|z^{i-1}):={\bar{s}}_{z^{i-1}}(z_i)$ and then further define $\bar{p}(z^n):={\prod }_{i=1}^n\bar{p}(z_i|z^{i-1})$. A simple calculation shows that we must have ${\int }_{z^n\in {{𝒵}}^n}\bar{p}(z^n)d{z}^n=1$, so we have constructed a probability distribution $\bar{p}$ which once again satisfies (23). The fundamental insight here is that, when the log-loss is used, every probability distribution defines a sequential prediction strategy and — perhaps more surprisingly — vice versa, every sequential prediction strategy defines a probability distribution, such that on all sequences of outcomes, the minus-log-likelihood is equal to the cumulative loss.

Example 5. Consider again the Bernoulli model ${ℳ}=\{p_{\theta }:\theta \in [0,1]\}$. Each element $p_{\theta }\in {ℳ}$ defines a prediction strategy which, no matter what happened in the past, predicts that the probability that the next outcome $Z_i=1$ is equal to $\theta$. It incurs cumulative loss, on sequence $z^n$ with $n_1$ ones and $n_0=n-n_1$ zeros, given by $n_1(-log\theta )+n_0(-log(1-\theta ))$. The Bayesian universal distribution $p_{w_{\text{J}}}^{{b}{a}{y}{e}{s}}$ with Jeffreys’ prior that we considered in Example 1 satisfies, as was already mentioned, $p_{w_{\text{J}}}^{{b}{a}{y}{e}{s}}(Z_{m+1}=1|z^m)=(m_1+(1/2))/(m+1)$, so it “learns” the probability of 1 based on past data and does not treat the data as i.i.d. any more. The asymptotic expansion (16) then shows that its cumulative loss is of order

We may now ask: given a parametric model ${ℳ}$, what distribution (i.e., prediction strategy) in ${ℳ}$ leads to the best predictions of data $z_1,\dots ,z_n$? For simplicity we will assume that data are i.i.d. according to all $p_{\theta }\in {ℳ}$. We have to distinguish between the best sequential prediction strategy with hindsight and the best prediction strategy that can be formulated before actually seeing the data. The former is given by the $p_{\theta }\in {ℳ}$ achieving

i.e., the best predictions with hindsight are given by the ML distribution ${\widehat{\theta }}_{\mathrm{\text{ML}}}(z^n)$. However, ${\widehat{\theta }}_{\mathrm{\text{ML}}}(z^n)$ is only knowable after seeing all the data $z^n$, whereas in reality, we have, at each time $i$, to make a prediction $\bar{p}(Z_i|z^{i-1})$ relying only on the previously seen data $z^{i-1}$. We might thus aim for a prediction strategy (distribution) $\bar{p}$ which will tend to have a small regret (additional prediction error)

$${R}{E}{G}(\bar{p},z^n)=\sum _{i=1}^n-log\bar{p}(z_i|z^{i-1})-\left[\underset{\theta \in \mathrm{\Theta }}{min}\sum _{i=1}^n-log{p}_{\theta }(z_i|z^{i-1})\right]=-log\bar{p}(z^n)+log{p}_{{\widehat{\theta }}_{\mathrm{\text{ML}}}(z^n)}(z^n).$$(24)

But what does “tend” mean here? One strict way to implement the idea is to require (24) to be small in the worst case — one looks for the distribution $\bar{p}$ achieving where the minimum is over all probability distributions on ${{𝒵}}^n$. But comparing (24) and (25) with (13), using that $-log$ is strictly decreasing, we see that the $\bar{p}$ achieving (25) is just the NML distribution with $v\equiv 1$, which was already our “favorite” distribution to use in MDL model comparison any way! And, just like before, if (25) has no solution, we may add a $-logv$ luckiness term to (24) so as to regularize the problem, and then the optimal prediction strategy will be given by $p_v^{{N}{M}{L}}$. We also see that with $v\equiv 1$, $\mathrm{\text{COMP}}({ℳ},v)$ is equal to the minimax regret (25); and with nonuniform $v$, $\mathrm{\text{COMP}}({ℳ})$ will become equal to the minimax luckiness regret, i.e., (24) with a $-logv$ term added.

We now also see where the idea to use the prequential plug-in distribution ${\bar{p}}^{{p}{r}{e}{q}}$ instead of ${\bar{p}}^{{n}{m}{l}}$ comes from: if calculating $p_v^{{n}{m}{l}}$ is too difficult, or if the horizon$n$ (which is needed to calculate $p_v^{{n}{m}{l}}$) is unknown, we might simply pick any estimator $\breve{\theta }$ which we think is “reasonable” and replace our prediction $p_v^{{n}{m}{l}}(Z_i|z^{i-1})$ by $p_{\breve{\theta }(z^{i-1})}(Z_i)$ — if the estimator was chosen cleverly, we can expect the resulting cumulative regret to be small. Reconsidering (14), we see that all the universal distributions, viewed as prediction strategies, with the right choice of luckiness functions, priors and$/$or estimates, can be made to achieve a logarithmic (in $n$) worst-case regret — since the cumulative log-loss achieved by the best predictor in hindsight usually grows linearly in $n$, a logarithmic regret is quite satisfactory. Returning to our Bernoulli example, we see that the cumulative log-loss obtained by ${\widehat{\theta }}_{\mathrm{\text{ML}}}$, the best with hindsight, is equal to $nH(n_1/n)=nH({\widehat{\theta }}_{\mathrm{\text{ML}}}(z^n))$, where $H(\theta )$ is the binary entropy, $H(\theta ):=-\theta log\theta -(1-\theta )log(1-\theta )$. Note that, in line with the above discussion, $nH({\widehat{\theta }}_{\mathrm{\text{ML}}})$ is linear in $n$ unless ${\widehat{\theta }}_{\mathrm{\text{ML}}}$ tends to 0 or 1, but the regret of ${\bar{p}}^{{b}{a}{y}{e}{s}}$ with Jeffreys’ prior is logarithmic in $n$.

We thus get a novel interpretation of MDL: it associates each model ${{ℳ}}_{\gamma }$ with a sequential prediction strategy ${\bar{p}}_{\gamma }$ that is designed to achieve small regret compared to the hindsight-optimal prediction strategy within ${{ℳ}}_{\gamma }$; it then picks the model for which the corresponding prediction strategy achieves the smallest cumulative loss on the data.

Related works

Dawid (see Ref. 15 and many subsequent works) suggests to use this prequential model choice also with respect to loss functions other than the logarithmic loss; minimax-optimal cumulative prediction strategies without making stochastic assumptions about the data, with log-loss but (mainly) with other loss functions, are one of the main topics in machine learning theory; see for example Ref. 24; but there they are generally not used for model comparison or selection.

Why the logarithmic score?

Why does it make sense to minimize cumulative log-loss? Outside of the MDL world, the log-loss is often used for two reasons: first, it is (essentially) the only local proper scoring rule.²⁵ Second, it has an interpretation in terms of money: for every sequential prediction strategy, there is a corresponding “sequential investment” strategy such that the smaller the cumulative log-loss, the larger the monetary gains made with this strategy (“Kelly Gambling”^17,20).

Within the MDL field however, the use of the log-loss comes from the Kraft inequality, which directly relates it to lossless data compression. As we already saw before Example 4, for any sequential prediction strategy, i.e., every distribution $p$ on sequences of length $n$, there is a lossless code $C$ such that, for all sequences of length $n$,

Conversely, for any code $C$, there is a corresponding distribution $p$ such that the above holds (see Chap. 3 of G07 for a very extensive explanation). Thus, the original MDL idea to “take the model that compresses the data most” is first made more formal by replacing it by “associate each model with a code that compresses well whenever some distribution in this model compresses well”, and this turns out to be equivalent to “associate each model ${{ℳ}}_{\gamma }$ with a distribution ${\bar{p}}_{\gamma }$ that assigns high likelihood whenever some distribution in the model assigns high likelihood”.

MDL prediction and “improper” estimation.

As is clear from the prequential interpretation of MDL given above, once a universal distribution $\bar{p}$ has been fixed, one can use it to predict $Z_i$ given $z^{i-1}$ by $\bar{p}(Z_i|z^{i-1})$. At least for i.i.d. data, we can estimate the underlying “true” distribution $p^{\ast }$ based on such predictions directly, by simply interpreting $\bar{p}(Z_i|z^{i-1})$ as an estimate of $p^{\ast }$! This is different from the previous form of MDL estimation described in Sec. 2.3, which was based on MDL (penalized ML) estimators ${\widehat{\theta }}_v$. Note that this standard MDL estimator is “in-model” or proper (to use machine learning terminology²⁶), whereas $\bar{p}(Z_i|z^{i-1})$ is out-model or improper: in general, there may be no $p\in {ℳ}$ such that $\bar{p}(\cdot |z^{i-1})=p$. For example, with Bayes universal distributions, $p_{w_{Z_i|z^{i-1}}}^{{b}{a}{y}{e}{s}}$ will be a mixture of distributions in ${ℳ}$ rather than a single element; see G07 for more discussion.

2.5. The luckiness function

The choice of a luckiness function is somewhat akin to the choice of a prior in Bayesian statistics, yet — as explained at length in Chap. 17 of G07 — there are very important differences, both technically (luckiness functions are not always integrable) and philosophically. Basically, a luckiness function just determines for what type of data one will be “lucky” [$v({\widehat{\theta }}_v(z^n))$ large] and get small cumulative regret based on small samples (and presumably, good model selection results as well), and for what data one will be less lucky and get good results only when the dataset grows much larger — $v$ may thus be chosen for purely pragmatic reasons. For example, as in Example 4 (see the italicized text there), if one assigns a large value $\pi (\gamma )$ to some model ${{ℳ}}_{\gamma }$ within a large collection of models $\{{{ℳ}}_{\gamma }:\gamma \in \mathrm{\Gamma }\}$ where $\gamma$ is sparse, one may do this because one hopes that this sub-model will already lead to reasonable predictions of future data, even though one feels that, at the same time, when more data becomes available, a model ${{ℳ}}_{{\gamma }^{\prime }}$ with a much larger number of nonzero parameters may at some point almost certainly become better (Example 4). Such an interpretation is not possible with a Bayesian prior $\pi$, where a large value of $\pi (\gamma )$ indicates a strong belief that ${{ℳ}}_{\gamma }$ is true, or at least, that predictions based on acting as if it will be true will be optimal — a Bayesian with high prior on $\gamma$ considers ${{ℳ}}_{\gamma }$ likely rather than just useful — a distinction worked out in detail by Grünwald.²⁷ Nevertheless, just like a Bayesian prior, the luckiness function has to be chosen by the user$/$statistician, and often contains a subjective element. Still, in contrast to Bayesian priors, since we invariably take a worst-case log-loss stance in MDL, there often is a uniquely preferable choice of luckiness function $v$ for parametric models ${ℳ}$. First, if $\mathrm{\text{COMP}}({ℳ},v)<\infty$ with uniform $v$ and no clear prior knowledge or preference is available, then uniform $v$ is usually preferable over other $v$, since it achieves the worst-case optimal prediction performance. Second, if $\mathrm{\text{COMP}}({{ℳ}}_{\gamma },v)=\infty$ with uniform $v$ for some $\gamma \in \mathrm{\Gamma }$ we can often still set the first few, say $m$, outcomes aside and pick a luckiness function $v_{\gamma }(\theta ):=p_{\theta }(z_1,\dots ,z_m)$ for $\theta \in {\mathrm{\Theta }}_{\gamma }$. The corresponding estimator (for fixed $\gamma$) ${\widehat{\theta }}_{v_{\gamma }}$ based on the data $z_{m+1},\dots ,z_n$ as given by (5) will then be equal to the ML estimator based on the full data, ${\widehat{\theta }}_{v_{\gamma }}(z_{m+1}^n)={\widehat{\theta }}_{\mathrm{\text{ML}}|\gamma }(z^n)$, and by choosing $m$ large enough, one can often get that $\mathrm{\text{COMP}}({{ℳ}}_{\gamma },v_{\gamma })$ is finite after all for all $\gamma \in \mathrm{\Gamma }$ and one may then compare models by picking the $\gamma$ maximizing $p_{v_{\gamma }}^{{n}{m}{l}}(z_{m+1},\dots ,z_n)$. Thus, the luckiness function is now determined by the first few “start-up” data points, and one uses NML based on the remaining data points with an estimator that coincides with the ML estimator based on the full data. G07 argues why this data-driven luckiness function is the best default choice available; note that it is analogous to our use of improper Bayes priors as described in Example 2.

3. Novel Universal Distributions

3.1. The switch distribution and the AIC–BIC dilemma

The AIC–BIC dilemma (see for example Ref. 28 and the many citations therein) is a classic conundrum in the area of statistical model selection: if one compares a finite number of models, the two standard benchmark methods, with (different) asymptotic justifications, are AIC and BIC. Suppose one first selects a model using AIC or BIC. One then predicts a future data point based on, e.g., maximum likelihood estimation, or by adopting the Bayes$/$MDL predictive distribution, within the chosen model. If one compares a finite number of models, then AIC tends to select the one which is optimal for prediction (compared to BIC, the predictions converge faster, by a factor of order $logn$, to the optimal ones). On the other hand, BIC is consistent: with probability 1, it selects the smallest model containing the true distribution, for all large $n$; the probability that AIC selects an overly large model does not go to 0 for large $n$. Both the predictive-optimality and the consistency property are desirable, but, like AIC and BIC, common methods all fail on one of the two. For example, MDL, with each of the four classical distributions, and Bayes factor model selection will behave like BIC for large $n$ and be consistent but prediction-sub-optimal; for any fixed $k$, leave-$k$-out and $k$-fold cross-validation will tend to behave like AIC and have the reverse behavior. Yang²⁸ shows that, in general, this dilemma cannot be solved: every consistent method has to be slightly prediction-sub-optimal in some situations; he also shows that prediction by model averaging cannot solve this dilemma either.

Nevertheless, as first shown by van Erven et al.¹⁶ (who thereby solved Open Problems 8 and 17 of G17), one can design universal distributions that “almost” get the best of both worlds: basing MDL model selection on them using (1) one gets a criterion which is strongly consistent while at the same time losing only an exceedingly small-order $loglogn$ factor in terms of prediction quality compared to the AIC-type methods. Although it can be applied to arbitrarily large model collections, the idea of this so-called switch distribution ${\bar{p}}^{{s}{w}{i}{t}{c}{h}}$ is best explained by considering the simplest case with just two nested models ${{ℳ}}_0\subset {{ℳ}}_1$: one starts with two standard universal distributions (say, Bayesian or luckiness-NML) ${\bar{p}}_0$ for ${{ℳ}}_0$ and ${\bar{p}}_1$ for ${{ℳ}}_1$. For every $i>0$, ${\bar{p}}_1$ defines a conditional distribution ${\bar{p}}_1(z_i,\dots ,z_n|z^{i-1})$. One now picks a “prior” distribution $\pi$ on the integers [typically one that decreases polynomially, e.g., $\pi (i)=1/(i(i+1))$], and one defines a new universal distribution for ${{ℳ}}_1$ by

This distribution is best understood from the prequential interpretation of MDL (Sec. 2.4). It will satisfy

$$\begin{array}{l}\sum _{i=1}^n-log{\bar{p}}^{{s}{w}{i}{t}{c}{h}}(z_i|z^{i-1})=-log{\bar{p}}^{{s}{w}{i}{t}{c}{h}}(z^n)\\ =-log\sum _{i=1}^n\pi (i){\bar{p}}_1(z_i,\dots ,z_n|z^{i-1}){\bar{p}}_0(z^{i-1})\\ \le \underset{i\in \{1,\dots ,n\}}{min}-log\pi (i){\bar{p}}_1(z_i,\dots ,z_n|z^{i-1})\cdot {\bar{p}}_0(z^{i-1})\\ \le \underset{i}{min}-log{\bar{p}}_1(z_i,\dots ,z_n|z^{i-1})\cdot {\bar{p}}_0(z^{i-1})\\ +2logn.\end{array}$$

In other words, the cumulative log-loss achieved by ${\bar{p}}^{{s}{w}{i}{t}{c}{h}}$ is “almost” (within an order $logn$ term) as small as that of the strategy that first predicts by ${\bar{p}}_0$ and then switches from ${\bar{p}}_0$ to ${\bar{p}}_1$ at the switching point $i$ that is optimal with hindsight. By clever choice of the prior $\pi$, one can get the extra term down to order $loglogn$. In cases where the data are actually sampled from a distribution in ${{ℳ}}_1$ that is “close” (defined suitably) to ${{ℳ}}_0$, the predictions based on ${\bar{p}}^{{s}{w}{i}{t}{c}{h}}$ will, with high probability, be substantially better than those based on ${\bar{p}}_1$ — a dramatic example (that makes very clear why this happens) is given in the first figure of Ref. 29. If the data come from a distribution that is “far” from ${{ℳ}}_0$, they will tend to be worse than those based on ${\bar{p}}_1$ by a negligible amount. Working out the math shows that associating ${{ℳ}}_1$ with ${\bar{p}}^{{s}{w}{i}{t}{c}{h}}$ and ${{ℳ}}_0$ with ${\bar{p}}_0$ indeed gives a strongly consistent model selection criterion that is almost [to within an $O(loglogn)$ factor] prediction-optimal, thus almost solving the AIC–BIC dilemma. van Erven et al.²⁹ describe in great detail why the standard NML or Bayesian universal model ${\bar{p}}_1$ does not lead to the optimal cumulative log-loss if the data come from a distribution close to, but not in, ${{ℳ}}_0$.

In case the number of models on the list is larger than two or even infinite, one has to associate each model with a separate switch distribution. The technique for doing so is described by van Erven et al.²⁹ who also give an efficient implementation and prove consistency and prediction-optimality of the switch distribution in a weak, cumulative sense for both finite and infinite numbers of models. van der Pas and Grünwald³⁰ mathematically show the “almost” prediction-optimality for a finite number of models.

3.2. Hybrids between NML Bayes and prequential plug-in

A Problem for NML: Unknown horizon Bayesian universal distributions with fixed priors have the property that the probability assigned to any initial sequence $z^{n^{\prime }}$, where $n^{\prime }<n$, is independent of the total length of the sequence. For other universal models, such as NML, this is not always the case. Take for example the Bernoulli model extended to sequences by independence: For sequence length $n=2$, the normalizing term in the NML equals $1+{(1/2)}^2+{(1/2)}^2+1=5/2$. For sequence length $n=3$, the normalizing term equals $2+6\times (1/3){(2/3)}^2=78/27$. For $n=2$, the NML probability of the sequence 00 is $1/(5/2)=0.4$. However, for the sequence length $n=3$, the probability of the initial sequence 00 is obtained as the sum of the probabilities of the sequences 000 and 001, which becomes $1/(78/27)+(4/27)/(78/27)\approx 0.397<0.4$. As shown by Ref. 21 (see also Ref. 31), there do exist cases in which NML is, like Bayes, horizon-independent, but these are very rare — see the Appendix.

The above discrepancy between the initial sequence probabilities for different sequence lengths may be a problem in situations where we need to obtain predictions without necessarily knowing the total length of the sequence, or the horizon. Another related issue is that even if the total sequence length was given, it can be computationally expensive to obtain marginal and conditional probabilities along the initial sequences. One possible solution would be to restrict to Bayesian universal distributions. However, while these solve the horizon issue, they are (a) still often computationally inefficient and (b) they lack NML-style worst-case regret interpretations. This has spurned research into universal codes that can be calculated without knowing the horizon in advance and that behave better as regards to (a) or (b), which we now review.

3.2.1. Prequential plug-in and the $(k/2)logn$ formula

The most straightforward way to deal with issue (a) is to use the prequential ${\bar{p}}^{{p}{r}{e}{q}}$ which, by construction, is horizon-independent. However, for the prequential ${\bar{p}}^{{p}{r}{e}{q}}$ (10) the BIC asymptotics (14) only hold in expectation if the data are sampled from one of the distributions in the model ${ℳ}$. This makes the result much weaker than for the other five universal distributions considered, for which the asymptotics hold for every individual sequence in some large set, i.e., without making any stochastic assumptions at all. One might thus wonder what happens for general data. Extending the earlier works by Takeuchi and Barron³² and Kotłowski et al.³³ shows that, if data are sampled from a distribution $p$, and $p_{\tilde{\theta }}$ is the distribution in ${ℳ}$ that is closest in Kullback–Leibler (KL) divergence to $p$, then (14) holds in expectation, with a correction term involving the variances of both distributions; for one-dimensional models, we get

a formula that can be extended to multidimensional models and individual sequence settings. Solving Open Problem 2 from G07, Grünwald and Kotłowski³⁴ show that, essentially, there exists no “in-model” estimator that can achieve the standard asymptotics in general; a correction such as (26) is always needed, whatever estimator one tries. Here an “in-model” estimator (or “proper” estimator, see end of Sec. 2.4) is an estimator that always outputs a distribution inside the model ${ℳ}$; the ML and Bayes MAP estimators are in-model, but the Bayes predictive distribution is not in-model, since it is a mixture distribution over all distributions in ${ℳ}$.

Solving Open Problem 3 from G07, Kotłowski et al.³³ also provide a new universal distribution, in which for any given estimator $\breve{\theta }$, $p_{\breve{\theta }}^{{p}{r}{e}{q}}(Z_{i+1}|z^i)=p_{\breve{\theta }(z^i)}(Z_{i+1})$ is turned into a slightly “flattened” version $p_{\breve{\theta }}^{{p}{r}{e}{q}\ast }(Z_{i+1}|z^i)$, which is not in ${ℳ}$ any more (it is not an in-model estimator), but it does achieve the standard $(k/2)logn$ asymptotics without correction. For example, in case of the normal location family with fixed variance ${\sigma }^2$, it coincides with a Bayesian predictive distribution based on a standard conjugate prior, which in this case is a normal with mean $\breve{\mu }$ (the Bayes MAP estimate) but variance ${\sigma }^2+O(1/n)$. More generally, $p_{\breve{\theta }}^{{p}{r}{e}{q}\ast }(Z_{i+1}|z^i)$ becomes a hybrid between the estimator $\breve{\theta }$ and a Bayes predictive distribution, but it has the advantage over the latter that it can be calculated without performing an integral over the parameter space. It thus provides an alternative to the NML distribution that is horizon-free and that is often faster to compute than ${\bar{p}}^{{b}{a}{y}{e}{s}}$.

Roos and Rissanen³⁵ and Rissanen et al.³⁶ developed other prequential, horizon-free universal codes that are non-Bayesian, yet remain more closely to NML in spirit than $p_{\breve{\theta }}^{{p}{r}{e}{q}\ast }$. They work out the details for discrete models including Bernoulli as well as linear regression models. For Bernoulli models, the resulting universal code coincides with the so-called one-step lookahead model proposed earlier by Takimoto and Warmuth.³⁷ For linear regression models the asymptotic consistency of the resulting model selection criterion was studied by Määttä et al.³⁸ and Määttä and Roos.³⁹ Relatedly, Watanabe and Roos⁴⁰ show that no horizon-independent strategy can be asymptotically minimax in the multinomial case and propose simple Bayesian universal models with a horizon-dependent Dirichlet prior that achieve asymptotic minimaxity and simplify earlier proposals. Among the proposed priors is $\mathrm{\text{Dir}}(\alpha ,\dots ,\alpha )$ with $\alpha =1/2-ln2/2lnn$ which converges to the Jeffreys’ prior $\mathrm{\text{Dir}}(1/2,\dots ,1/2)$ but has a mild dependency on the horizon $n$.

3.3. Hypothesis testing: Universal distributions based on the reverse information projection

Suppose we compare just two models, ${{ℳ}}_0$ and ${{ℳ}}_1$, as explanations for data $z^n$, a situation similar to classical null hypothesis testing, the standard method for evaluating new treatments in the medical sciences and scientific hypotheses in most applied sciences such as psychology, biology and the like: we can think of ${{ℳ}}_0$ and ${{ℳ}}_1$ as two hypotheses, where, usually, ${{ℳ}}_0$ represents the status quo (“treatment not effective”, “coin unbiased”). In case^a${{ℳ}}_0=\{P_0\}$ represents a simple (singleton) hypothesis, there is a strong additional motivation for using MDL as a hypothesis testing method, and in particular, for quantifying the evidence against ${{ℳ}}_0$ in terms of

the code length or cumulative-log-loss difference (see Sec. 2.4) between encoding (or sequentially predicting) the data with $p_0$ and with ${\bar{p}}_1$. This additional motivation is given by the no hyper-compression inequality (G07), a mathematical result stating that, no matter how ${\bar{p}}_1$ is defined, as long as it is a probability distribution, we have for all $K>0$, and $0\le \alpha \le 1$, This expresses, in terms of sequential log-loss prediction (compression), that, if $P_0$ is true, then the probability that one can predict data better, by $K$ or more loss units, by predictions based on ${\bar{p}}_1$ rather than $p_0$, is exponentially small in $K$ — and this holds independently of the sequence length $n$. In terms of more classical quantities, it states that, no matter how we chose ${\bar{p}}_1$, if $P_0$ holds true then the likelihood ratio is a $p$-value. In fact it is a conservative $p$-value, giving usually somewhat less evidence against ${{ℳ}}_0$ than a standard $p$-value, for which the rightmost inequality in (27) is an equality. The inequality (27) goes in the right direction to retain the cornerstone of classical Neyman–Pearson testing: if one sets significance level $\alpha$ before seeing the data and one chooses ${{ℳ}}_1$ whenever $D(Z^n)\le -log\alpha$, i.e., $p_0(Z^n)/{\bar{p}}_1(Z^n)\le \alpha$, then the probability, under the null $P_0$, of making a false decision is bounded by $\alpha$. But the fact that the rightmost inequality in (27) is usually strict has pleasant practical repercussions: as explored by Ref. 17, the Type-I error guarantees are retained under optional continuation. This is the (common) practice to decide, on the basis of an initial sample $z^n$, whether or not to gather new data and do a second test. One may for example decide to gather new data if the result based on $z^n$ was hopeful yet not conclusive. This is highly problematic for standard, strict $p$-value-based hypothesis testing, but with MDL testing with a simple ${{ℳ}}_0$, one can simply multiply the likelihood ratios of the two (or more) tests performed, or equivalently, add the code length differences for each test performed. The resulting code length difference$/$likelihood ratio will still lead to valid Type-I error bounds.¹⁷

But, all this holds only for simple ${{ℳ}}_0$. Yet the tests most used in practice, such as the $t$-test and contingency table tests, all involve composite ${{ℳ}}_0=\{P_{\theta }:\theta \in {\mathrm{\Theta }}_0\}$. For composite ${{ℳ}}_0$, the no-hypercompression inequality (27) usually only holds for some $P_0\in {{ℳ}}_0$, but for Type-I error guarantees and the like we would want to have it hold for all $P_{\theta }$ with $\theta \in {\mathrm{\Theta }}_0$. That is, we would like to employ universal distributions ${\bar{p}}_1$ and ${\bar{p}}_0$ such that we have

In general, this will not hold for standard choices (NML, Bayes, prequential plug-in, etc.) of ${\bar{p}}_1$ and ${\bar{p}}_0$. However, Grünwald et al.¹⁷ show that, for any given (arbitrary) ${\bar{p}}_1$, one can, under very mild conditions, construct a ${\bar{p}}_0$ such that (28) holds, thereby solving Open Problems 9 and 19 of G07. This ${\bar{p}}_0$ is the Reverse Information Projection^41,42 of ${\bar{p}}_1$ onto ${{𝒫}}_{{b}{a}{y}{e}{s}}({{ℳ}}_0)$, where ${{𝒫}}_{{b}{a}{y}{e}{s}}$ is the set of densities ${\bar{p}}_0$ for $z^n$ that can be written as Bayes marginal distributions ${\bar{p}}_0^{{b}{a}{y}{e}{s}}(z^n)=\int p_{\theta }(z^n)w_0(\theta )d\theta$ for some prior $w_0$ on ${\mathrm{\Theta }}_0$ — for every prior $w_0$, ${{𝒫}}_{{b}{a}{y}{e}{s}}({{ℳ}}_0)$ contains a separate distribution on $Z^n$. The RIPr is defined as the density achieving ${min}_{{\bar{p}}_0\in {{𝒫}}_{{b}{a}{y}{e}{s}}({{ℳ}}_0)}\times D({\bar{p}}_1\parallel {\bar{p}}_0)$, where $D(\cdot \parallel \cdot )$ is the Kullback–Leibler divergence. Thus, one constructs a ${\bar{p}}_0$ with the desired no-hypercompression property, and at the same time, it will minimize KL divergence to ${\bar{p}}_1$, which implies that if data were sampled from ${\bar{p}}_1$, it would yield optimal log-loss predictions. This, in turn, implies that the ${\bar{p}}_0$ constructed this way will satisfy the standard asymptotics (14) as long as the ${\bar{p}}_1$ on which it is based does. Based on the likelihood ratio between ${\bar{p}}_1$ and its RIPr ${\bar{p}}_0$, one is also allowed to do optional continuation while retaining Type-I Error guarantees. Thus, even if one is an adherent of classical, frequentist testing theory, there are strong reasons for MDL-style testing based on the RIPr universal distribution. Grünwald et al.¹⁷ further extend the reasoning to give guidelines on how ${\bar{p}}_1$ can be chosen to get further good frequentist properties.

Example 6 (Right-Haar priors and the Bayesian $t$-test). In a series of papers (highlights include Refs. 43 and 44), Berger and collaborators established Bayes factor testing methods for composite ${{ℳ}}_0=\{P_{\theta }:\theta \in {\mathrm{\Theta }}_0\}$ where the only free parameters in ${\mathrm{\Theta }}_0$ are “nuisance” parameters that are shared by ${\mathrm{\Theta }}_1$ and are governed by a group structure. A prime example is the unknown variance in the $t$-test. Berger uses a special type of improper prior, the so-called right-Haar prior, which can be defined for every such type of nuisance parameter. While Bayes factors usually do not combine well with improper priors, the Bayes factors for group invariance parameters equipped with the right-Haar prior behave remarkably well. Grünwald et al.¹⁷ show that, even though the right-Haar priors are usually improper, they can also be understood from a purely MDL perspective: if ${\bar{p}}_1^{{b}{a}{y}{e}{s}}$ and ${\bar{p}}_0^{{b}{a}{y}{e}{s}}$ are equipped with the right-Haar prior on the nuisance parameters, and the prior on the additional parameters in ${\bar{p}}_1^{{b}{a}{y}{e}{s}}$ satisfies some additional requirements, then both ${\bar{p}}_1^{{b}{a}{y}{e}{s}}$ and ${\bar{p}}_0^{{b}{a}{y}{e}{s}}$ can be interpreted as sequential prediction strategies, and the log of the Bayes factor can be interpreted as the code length$/$cumulative-log-loss difference. Moreover, ${\bar{p}}_0^{{b}{a}{y}{e}{s}}$ is (essentially) the RIPr for ${\bar{p}}_1^{{b}{a}{y}{e}{s}}$ and the no-hypercompression inequality (28) that is so desirable from a frequentist perspective holds uniformly for all ${\theta }_0\in {\mathrm{\Theta }}_0$.

Let us consider the one-sample Bayesian $t$-test as an example. Here ${{ℳ}}_0=\{p_{0,\sigma }:\sigma >0\}$ is the set of all normal distributions with mean 0; the variance ${\sigma }^2$ is a free parameter. ${{ℳ}}_1=\{p_{\mu ,\sigma }:\mu \in ℝ,\sigma >0\}$ is the set of all normal distributions with as free parameters $\mu$ and $\sigma$. The question of interest is to establish whether $\mu =0$ or not; $\sigma$ is an unknown “nuisance” parameter — it determines the scale of the data but is not itself of intrinsic interest. In the Bayesian$t$-test one equips both ${{ℳ}}_0$ and ${{ℳ}}_1$ with the improper right-Haar prior, $w(\sigma )=1/\sigma$. To complete the definition of ${\bar{p}}_1^{{b}{a}{y}{e}{s}}$, ${{ℳ}}_1$ is equipped with a conditional prior density (given $\sigma$) on the effect size $\delta :=\mu /\sigma$. This second density has to be symmetric around zero and proper (this is what we called the “additional requirement on the prior on ${\mathrm{\Theta }}_1$”, instantiated to the case where the nuisance parameter is a variance). One now proceeds by testing using the Bayes factor ${\bar{p}}_0^{{b}{a}{y}{e}{s}}/{\bar{p}}_1^{{b}{a}{y}{e}{s}}$. In this special case, the procedure was already suggested by Jeffreys,⁴⁵ and the right-Haar prior coincides with Jeffreys’ prior on the variance. Berger et al. extend the method to general group-invariant parameter vectors such as the joint mean and variance in the two-sample $t$-test, testing a Weibull against a log-normal and many other scenarios.

4. Graphical Models

Graphical models are a framework for representing multivariate probabilistic models in a way that encompasses a wide range of well-known model families, such as Markov chains, Markov random fields and Bayesian networks; for a comprehensive overview, see Ref. 46. A key property of a graphical model is parsimony, which can mean, for instance, a low-order Markov chain or more generally a sparse dependency graph that encodes conditional independence assumptions. Choosing the right level of parsimony in graphical models is an ideal problem for MDL model selection.

In Bayesian network model selection the prevailing paradigm is, unsurprisingly, the Bayesian one. Especially the works of Geiger and Heckerman⁴⁷ and Heckerman et al.⁴⁸ have been extremely influential. The main workhorse of this approach is the so-called Bayesian Dirichlet (BD) family of scores which is applicable in the discrete case where the variables being modeled are categorical. Given a data sample, such scores assign a goodness value to each model structure. Exhaustive search for the highest scoring structure is possible when the problem instance (characterized by the number of random variables) is of limited size, but heuristic search techniques such as variants of local search or “hill-climbing” can be used for larger problem.

Different forms of the BD score imply different Dirichlet priors (different hyper-parameters) for the local multinomial distributions that comprise the joint distribution. For example, in the commonly used BDeu score, the priors are determined by a single hyper-parameter, $\alpha$. For a variable $X_i$ with $r$ distinct values and parents ${\mathrm{\text{Pa}}}_i$ that can take $q$ possible combinations of values (configurations), the BDeu prior is $\mathrm{\text{Dir}}(\alpha /rq,\dots ,\alpha /rq)$. One of the main motivations for adopting this prior is that it leads to likelihood equivalence, i.e., it assigns equal scores to all network structures that encode the same conditional independence assumptions. In light of the fact that Bayesian model selection embodies a particular form$/$variation of MDL, these methods fit, at least to some extent, in the MDL framework as well. However, there also exist more “pure”, non-Bayesian MDL methods for model selection in Bayesian networks; we mention Refs. 49 and 50 as early representative examples. These early methods are almost invariably based on the two-part coding framework. More recently, several studies have proposed new model selection criteria that exploit the NML distribution. One approach is a continuous relaxation of NML-type complexities proposed by Miyaguchi et al.⁵¹ in which the model selection problem takes on a tractable Lasso-type $L_1$-minimization form (see also Example 4). In other approaches, NML [or usually, approximations (but not relaxations) thereof] are used directly for encoding parts of the model; we now describe these latter approaches in a bit more detail.

4.1. Factorized NML and variants

Silander et al.⁵² propose the factorized NML (fNML) score for Bayesian network model selection which was designed to be decomposable meaning that it can be expressed as a sum that includes a term for each variable in the network. This property facilitates efficient search among the super-exponential number of possible model structures; see, e.g., Ref. 48. The fNML score factors the joint likelihood not only in terms of the variables but also in terms of distinct configurations of the parent configurations. Each factor in the product is given by a multinomial NML probability, for which a linear-time algorithm by Kontkanen and Myllymäki⁵³ can be used.

A similar idea where a Bayesian network model selection criterion is constructed by piecing together multiple NML models under the multinomial model was proposed recently by Silander et al.⁵⁴ In the proposed quotient NML (qNML) score, the local scores corresponding to each variable in the network are defined as log-quotients of the form

where ${\mathrm{\text{NML}}}_{{f}{u}{l}{l}}$ refers to an NML distribution defined by using a fully connected network to model the variable $X_i$ and its parents ${\mathrm{\text{Pa}}}_i$ in the numerator and the same thing for the parent set ${\mathrm{\text{Pa}}}_i$ in the denominator. Technically, this amounts to collapsing the configurations of the variables into distinct values of a single categorical variable. Even though the resulting categorical variable may have a huge number of possible values, the linear time algorithm⁵³ or efficient approximations (see the next sub-section) can be used to implement the computations. A notable property of the qNML score is that, unlike the fNML score, it is likelihood-equivalent (see above).

Eggeling et al.⁵⁵ apply similar ideas to a different model class, namely parsimonious Markov chains. There too, the likelihood is decomposed into factors depending on the configurations of other variables, and each part in the partitioning is modeled independently using the multinomial NML formula. The authors demonstrate that the fNML-style criterion they propose leads to parsimonious models with good predictive accuracy for a wide range of different scenarios, whereas the corresponding Bayesian scores are sensitive to the choice of the prior hyper-parameters, which is important in the application where parsimonious Markov chains are used to model DNA binding sites.⁵⁶

In all these papers, both simulated and real-world data experiments suggest that the MDL-based criteria are quite robust with respect to the parameters in the underlying data source. In particular, the commonly used Bayesian methods (such as the BDeu criterion) that are being used as benchmarks are much more sensitive and fail when the assumed prior is a poor match to the data-generating model, whereas the MDL methods are invariably very close to the Bayesian methods with the prior adapted to fit the data. This poses interesting questions concerning the proper choice of priors in the Bayesian paradigm.

In fact, the prevalence of the Bayesian paradigm and the commonly used BD scores is challenged by two recent observations: First, Silander et al.⁵² show that the Dirichlet prior with hyper-parameters $(1/2,\dots ,1/2)$, which is the invariant Jeffreys’ prior for the multinomial model, but not likelihood-equivalent when used in the BD score, is very close to the fNML model and consequently, enjoys better robustness properties than the BDeu score which is the likelihood-equivalent BD score variant. Second, Suzuki⁵⁷ shows that the BDeu criterion is irregular, i.e., prone to extreme overfitting behavior in situations where a deterministic relationship between one variable and a set of other variables holds in the data sample. The MDL scores discussed above are regular in this respect and their robustness properties seem to be better than those of the BD scores, see Ref. 54.

4.2. Asymptotic expansions for graphical models

Asymptotic results concerning MDL-based criteria in graphical models are interesting for several reasons. For one, they lead to efficient scores that can be evaluated for thousands of different model structures. Second, asymptotic expansions can lead to insights about the relative complexity of different model structures.

Various asymptotic forms exist for the point-wise and the expected regret depending on the model class in question. For convenience we repeat the classical expansion of the NML (as well as the Bayesian marginal likelihood with Jeffreys’ prior) regret$/$model complexity that applies for regular model classes ${ℳ}=\{p_{\theta }:\theta \in \mathrm{\Theta }\}$ for which $\mathrm{\text{COMP}}({ℳ},v)$ is finite with uniform $v$ (see Sec. 2.2 above) :

where $k$ is the dimension of the model, $|I(\theta )|$ is the determinant of the Fisher information matrix at parameter $\theta$, the integral is over the parameter space $\mathrm{\Theta }$ and the remainder term $o(1)$ vanishes as the sample size tends to infinity.

For discrete data scenarios, by far the most interesting case is the multinomial model (extension of the Bernoulli distribution to an i.i.d. sequence of $r$-valued categorical random variables) since it is a building block of a number of MDL-criteria such as fNML and qNML (see above). There are many asymptotic expansions for the NML regret under the multinomial model. Probably the most useful is the one proposed by Szpankowski and Weinberger⁵⁸ :

where $n$ is the sample size, $\alpha =\frac{r}{n}$ and $C_{\alpha }=\frac{1}{2}+\frac{1}{2}\sqrt{1+\frac{4}{\alpha }}$. This simple formula is remarkably accurate over a wide range of finite values of $n$ and $r$ (see Ref. 54). Note that the leading term is proportional to $n$ (rather than $logn$ as usual) because the formula is derived for the regime $r=\mathrm{\Theta }(n)$ where the alphabet size grows proportionally to the sample size. If $r$ grows slower than $n$ or not at all, the leading term tends to the classical form (29), where the leading term is $\frac{k}{2}logn$. In practice, the approximation (30) is applicable for a wide range of $r/n$ ratios.

Roos,⁵⁹ and Zou and Roos⁶⁰ studied the second term in the expansion (29), namely the Fisher information integral, under Markov chains and Bayesian networks using Monte Carlo sampling techniques. This approach reveals systematic differences between the complexities of models even if they have the same number of parameters.

5. Latent Variable and Irregular Models

Although thus far we have highlighted exponential family and regression applications, NML and other universal distributions can of course be used for model selection and estimation in complete generality — and many practical applications are in fact based on highly irregular models. Often, “classical” two-part distributions (based on discretized models) are used, since NML distributions often pose computational difficulties. However, Yamanishi and collaborators have managed to come up with tractable approximations of NML-type distributions for some of the most important irregular (i.e., nonexponential family) models such as hierarchical latent variable models,⁶¹ and the related Gaussian mixture models.^62,63 Suzuki et al.⁶⁴ provide an NML approach to nonnegative matrix factorization. Two-part codes (and corresponding MDL estimators) for mixture families that come close to achieving the minimax regret were considered very recently by Miyamoto et al.⁶⁵

When it comes to asymptotic approximations for code lengths$/$log-likelihoods based on NML and other universal distributions — all approximations so far (in Sec. 2.2) were derived essentially assuming that the model under consideration is an exponential family. Extensions to curved exponential families and generalized linear models are relatively straightforward (see G07 for details). For more irregular models, Watanabe has proposed the widely applicable information criterion (WAIC) and the widely applicable Bayesian information criterion (WBIC), see Refs. 66 and 67, where the latter can be viewed as an asymptotic expansion of the log-likelihood based on a Bayesian universal distribution. It coincides with BIC when applied to regular models but is applicable even for singular (irregular) models. The asymptotic form of WBIC is

where ${\theta }_0$ is the parameter value minimizing the Kullback–Leibler divergence from the model to the true underlying distribution, and $\lambda >0$ is a rational number called the real log-canonical threshold (see Ref. 67), which can be interpreted as the effective number of parameters (times two).

6. Frequentist Convergence of MDL and Its Implications

Rissanen first formulated the MDL Principle as — indeed — a Principle: one can simply start by assuming, as an axiom, that modeling by data compression (or, equivalently, sequential predictive log-loss minimization) is the right thing to do. One can also take a more conventional, frequentist approach, and check whether MDL procedures behave desirably under standard frequentist assumptions. We now review the results that show that, in general, they do — thus providing a frequentist justification of MDL ideas: with some interesting caveats, MDL model selection is typically consistent (the smallest model containing the true distribution is eventually chosen, with probability one) and MDL prediction and estimation achieves good rates of convergence (the Hellinger distance between the estimated and the true density goes to zero, with high probability, quite fast). In this section we review the most important convergence results. In particular, Sec. 6.1 shows that the link between data compression and consistent estimation is in fact very strong; and Sec. 6.4 shows that, by taking MDL as a principle, one can get useful intuitions about deep questions concerning deep learning; and the intuitions can then, as a second step, be once again validated by frequentist results.

Thus, let us assume, as is standard in frequentist statistics, that data are drawn from a distribution in one of the models under ${{ℳ}}_{\gamma }$ under consideration. We consider consistency and convergence properties of the main MDL procedures in their main applications: model selection, prediction and estimation.

Model selection

For model selection between a finite number of models, all universal codes mentioned here are consistent in wide generality; for example, this has been explicitly proven if the data are i.i.d. and all models on the list are exponential families, but results for more complex models with dependent data have also been known for a long time; see G07 for an overview of results. If the collection of models is countably infinite, then results based on associating each ${{ℳ}}_{\gamma }$ with ${\bar{p}}_{\gamma }^{{b}{a}{y}{e}{s}}$ have also been known for a long time; such results typically hold for “almost all” (suitable defined) distributions in all ${{ℳ}}_{\gamma }$; again, see G07 for a discussion of the (nontrivial) “almost all” requirement. These countable-$\mathrm{\Gamma }$ consistency results were extended to the switch distribution by van Erven et al.²⁹

Prediction and “improper” estimation

As to sequential prediction (Sec. 2.4), the rate of convergence results are very easy to show (see Chap. 15 of G07), but these typically only demonstrate that the cumulative-log-loss prediction error of sequentially predicting with a universal distribution $\bar{p}$ behaves well as $n$ increases. Thus, since the sum of prediction errors is small, say (for parametric models) of order $logn$, for most$t$ the individual prediction error at the $t$th sample point must be of order $1/t$, since ${\sum }_{t=1}^{n}1/t-logt=O(1)$. Still, it remains an open question how to prove for individual $t$ what exactly the expected prediction error is at that specific $n$. Since one can view each prediction as an “improper” estimate (end of Sec. 2.4), the convergence rates of the resulting estimators, which estimate the underlying distribution based on a sample of size $t$ as $\bar{p}(Z_{t+1}|z^t)$, usually also behave well in a cumulative sense, but again it is very hard to say anything about individual $t$. The asymptotic expansions (15) and (16) imply that, for fixed parametric models ${{ℳ}}_{\gamma }$, ${\bar{p}}_{\gamma }^{{b}{a}{y}{e}{s}}$ and ${\bar{p}}_{\gamma }^{{n}{m}{l}}$ achieve optimal cumulative prediction and estimation errors. If, however, they are defined relative to a full model class ${ℳ}={\bigcup }_{\gamma \in \mathrm{\Gamma }}{{ℳ}}_{\gamma }$ consisting of at least two nested models, then they may fail to achieve optimal rates by a $logn$ factor. van Erven et al.²⁹ show that sequential prediction$/$estimation based on the switch distribution achieves the minimax-optimal cumulative prediction$/$estimation error rates even in such cases. van der Pas and Grünwald³⁰ show that, if only two models are compared, then the optimal obtainable rate for individual $n$ for any consistent procedure is achieved as well.

6.1. Frequentist convergence of MDL estimation

Very strong results exist concerning the convergence of MDL estimation based on an MDL estimator ${\widehat{\theta }}_v$ as given by (5). A first, classical result was already stated by the ground-breaking work,⁶⁸ establishing that consistency and good convergence rates can be obtained for the special case of a two-part-code estimator ${\widehat{\theta }}_w$ based on a probability mass function $w$, as long as $w$ satisfies ${\sum }_{\theta \in \ddot{\mathrm{\Theta }}}w{(\theta )}^{\eta }<\infty$ for some $\eta <1$ and $w$ puts sufficient prior mass in a KL neighborhood of the true $\theta$. These results were greatly extended by Zhang^69,70 and further, very recently, by Grünwald and Mehta.¹³ The latter consider ${\widehat{\theta }}_v$ for general $v$. Let ${ℳ}=\{P_{\theta }:\theta \in \mathrm{\Theta }\}$ be a statistical model and suppose data are i.i.d. $\sim p$ with $p=p_{{\theta }^{\ast }}\in {ℳ}$. They find that a sufficient condition for consistency is that $v({\theta }^{\ast })<\infty$ and that for some $\eta <1$, the following generalized model complexity

is bounded and $o(n)$, i.e., it grows slower than linear (the slower it grows, the faster the MDL estimator converges to the true distribution in Hellinger distance). This condition strictly and significantly weakens the Barron–Cover requirement. The result holds without any further conditions; for example, ${ℳ}$ may be a countable union of parametric models or even a huge nonparametric model. Note that ${\mathrm{\text{COMP}}}_1({ℳ},v)$ is the model complexity that we have encountered before in (7). Ironically, for any $\eta <1$, slow growth $\left(o\right(n\left)\right)$ of ${\mathrm{\text{COMP}}}_{\eta }({ℳ})$ is sufficient for consistency of ${\widehat{\theta }}_v$, but for $\eta =1$, which would be more fully in line with the MDL ideas, it is not.

6.2. From MDL to Lasso

As illustrated in Example 4, when used for high-dimensional variable selection, the original MDL approach would be to use a mixed two-part$/$one-part code as in (1) with a $-log\pi (\gamma )$ term to account for the model index $\gamma \in \mathrm{\Gamma }$. In such settings, there may well be $p>n$ variables of interest, each of which may or may not be included in the model, so that the minimization over $\mathrm{\Gamma }$ requires trying out $2^p\gg 2^n$ choices — which is practically infeasible. For this reason, in practice people have strongly preferred the Lasso and related methods based on $L_1$-penalties, which take linear rather than exponential time in $p$ (note that the classic MDL essentially penalizes by an $L_0$-penalty). However, Barron and Luo,⁷¹ and Barron et al.⁷² showed that, under some conditions on the true distribution (such as Gaussian noise), the Lasso method can be re-interpreted in terms of code length minimization after all; see also Ref. 73, and, for further extensions, Refs. 74 and 75. For a different approach to unify model selection with very high-dimensional models with the luckiness NML, see Miyaguchi and Yamanishi.²³

Although some of the details may differ, it seems that most of these works are subsumed by the aforementioned result of Grünwald and Mehta¹³ who show that general penalized estimators can be re-interpreted as minimizing a one-part code length as long as ${\mathrm{\text{COMP}}}_1({ℳ},v)$ is bounded, and can be proven consistent under the (still quite weak) condition that ${\mathrm{\text{COMP}}}_{\eta }$ as in (32) is bounded for some $\eta <1$. Thus, the connection between MDL and general (including Lasso and other $L_1$-penalties, but also with entirely different penalties) penalization methods is substantially stronger than it seemed before these developments took place.

Supervised machine learning

Importantly, all the works mentioned here except Ref. 13 cannot show convergence under misspecification — for example, when applied to the Lasso, they would require an assumption of normal noise (corresponding to the squared error used in the Lasso fit, which is equivalent to the log-loss under a normal distribution for the noise). In practice though, the Lasso (with the squared error) is often used in cases in which one cannot assume normally distributed errors. Reference 13 contains results that can still be used in such cases [although the formula for ${\mathrm{\text{COMP}}}_{\eta }({ℳ},v)$ changes], based on ideas which we sketch in the following sub-section.

More generally, one of the major areas within machine learning is supervised learning in which one assumes that data $(X_1,Y_1),(X_2,Y_2),\dots$ are i.i.d. $\sim P_0$, with $X_i\in {𝒳}$ and $Y_i\in {𝒴}$, and one aims to use the data to learn a predictor function $f:{𝒳}\to {{𝒴}}^{\prime }$ that has small expected loss or risk, defined as $E_{(X,Y)\sim P}[\ell (f(X),Y)]$, where $\ell :{𝒴}\times {{𝒴}}^{\prime }\to ℝ$ is some loss function of interest and $f$ is a member of some “predictor model” ${ℱ}$. For example, the statistical notion of “regression with random design” corresponds, in machine learning, to a supervised learning problem with ${𝒴}={{𝒴}}^{\prime }=ℝ$ and $\ell (y^{\prime },y)={(y^{\prime }-y)}^2$. Early MDL convergence results do not cover this “supervised” situation: they are not equipped to handle either random design or loss functions beyond the log-loss. Some of the more recent works mentioned above are able to handle random design but not general loss functions (for example, for Lasso-type applications they require the noise to be normally distributed). Reference 13 seems to be the first that can fully handle supervised learning scenarios: the convergence results can be used with random design, and they can also be used with large classes of loss functions including squared error (without normality assumption) and zero$/$one-loss. This is achieved by associating predictors $f$ with densities $p_f(x,y)\propto exp(-\ell (f(x),y))$, so that the log-loss relative to density $p_f$ on data $(x,y)$ becomes linearly related to the loss of $f$ on $(x,y)$; the analysis then proceeds via analyzing convergence of MDL for the densities $\{p_f:f\in {ℱ}\}$ as a misspecified probability model.

6.3. Misspecification

As beautifully explained by Rissanen,⁷⁶ one of the main original motivations for MDL-type methods is that they have a clear interpretation independent of whether any of the models under consideration is “true” in the sense that it generates the data: one chooses a model minimizing a code length, i.e., a prediction error on unseen data, which is meaningful and presumably might give something useful irrespective of whether the model is true (Rissanen even argues that the whole notion of a “true model” is misguided). This model-free paradigm also leads one to define the NML distribution as minimizing prediction error in a stringent worst-case-over-all data sense [Eq. (13)] rather than a stochastic sense. Nevertheless, it is of interest to see what happens if one samples data from a distribution for which all models under consideration are wrong, but some are quite useful in the sense that they lead to pretty good predictions. Doing this leads to rather unpleasant surprises: as first noted by Grünwald and Langford,⁷⁷ MDL (and Bayesian inference) can become inconsistent: one can give examples of $\{{{ℳ}}_{\gamma }:\gamma \in \mathrm{\Gamma }\}$ with countably infinite $\mathrm{\Gamma }$ and a “true” data generating distribution $P_0$ such that, when data are sampled i.i.d. from $P_0$, MDL will tend to select a sub-optimal model for all large $n$ — while all sub-models ${{ℳ}}_{\gamma }$ are wrong, one of them, ${{ℳ}}_{\tilde{\gamma }}$ is optimal in several intuitive respects (closest in KL divergence to $P_0$, leading to best predictions under a number of loss functions), yet it will not be selected for large $n$. While the models considered by Grünwald and Langford⁷⁷ were quite artificial, Grünwald and van Ommen⁷⁸ showed that the same can happen in a more natural linear regression setting; moreover, they also showed that even if $\mathrm{\Gamma }$ is finite, although then eventually MDL will select the best sub-model, for even relatively large $n$ it may select arbitrarily bad sub-models. De Heide⁷⁹ shows that the problem also occurs with MDL and Bayesian regression with some real-world datasets.

It turns out that the root of the problem is related to the no-hypercompression property (27). If the collection of models ${ℳ}={\bigcup }_{\gamma \in \mathrm{\Gamma }}{{ℳ}}_{\gamma }$ contains the density $p_0$ of the “true” distribution $P_0$, then any distribution $p\in {\bigcup }_{\gamma \in \mathrm{\Gamma }}{{ℳ}}_{\gamma }$ will satisfy no-hypercompression relative to the true $p_0$ :

This property underlies the proof of all MDL consistency and rate-of-convergence results, such as those by Barron and Cover,⁶⁸ Zhang,⁶⁹ and Grünwald and Mehta.¹³ However, if the model class ${ℳ}$ does not contain the true $p_0$, then, in order to prove consistency, one needs (33) to hold with the $P_0$ outside the brackets unchanged, but the $p_0$ inside the brackets replaced by $\tilde{p}$, the distribution$/$density in ${ℳ}$ that is closest to $P_0$ in KL divergence (why it should be KL is explained at length by Grünwald and van Ommen⁷⁸). Unfortunately though, (33) does not necessarily hold with $p_0$ replaced by $\tilde{p}$. If it does not, MDL (and Bayesian methods, whose consistency relies on similar properties) may become inconsistent. Grünwald and van Ommen⁷⁸, based on earlier ideas in Refs. 80 and 81, propose a solution that works for Bayesian universal distributions: it replaces the likelihoods $p_{\theta }(z^n)$ for every $p=p_{\theta }$ with $p\in {ℳ}$ by the generalized likelihood $p_{\theta }^{\eta }(z^n)$ for some $\eta >0$; usually $\eta <1$ — this $\eta$ has the same mathematical function as the $\eta$ appearing in (32). It turns out that with such a modification, if $\eta$ is chosen small enough, a version of the no-hypercompression inequality (33) holds after all. References 78 and 81 also provide a method for learning $\eta$ from the data, the “Safe Bayesian” algorithm (note that $\eta$ cannot be learned from the data by standard MDL or Bayesian methods). The recent work of Grünwald and Mehta¹³ suggests that the modification of likelihoods by exponentiating with $\eta$ should work for general MDL methods as well.

6.4. PAC-MDL bounds and deep learning

One of the great mysteries of modern deep learning methods in machine learning is the following⁸²: deep learning is based on neural network models which can have many millions of parameters. Although typically run on very large training samples $z^n$, $n$ is usually still so small that the data can be fit perfectly, with zero error on the training set. Still, the trained models often perform very well on future test sets of data. How is this possible? At first sight this contradicts the tenet, shared by MDL and just about any other method of statistics, that good generalization requires the models to be “small” or “simple” [small $\mathrm{\text{COMP}}({ℳ})$ in MDL analyses, small VC dimension or small entropy numbers in statistical learning analyses] relative to the sample size. One of several explanations (which presumably all form a piece of the puzzle) is that the local minimum of the error function found by the training method is often very broad — if one moves around in parameter space near the minimum, the fit hardly changes. Hochreiter and Schmidhuber⁸³ already observed that describing weights in sharp minima requires high precision in order to not incur nontrivial excess error on the data, whereas flat minima can be described with substantially lower precision, thus forging a connection to the MDL idea; in fact related ideas already appear in Ref. 84. In these papers, the MDL Principle is used in a manner that is less direct than what was done thus far in this paper: we (and, usually, Barron and Rissanen) directly hunt for the shortest description of the data. In contrast, the aforementioned authors simply note that, no matter how a vector of parameters for a model was obtained, if, with the obtained vector of parameters, the data can be compressed substantially, for example by coding first the parameters and then the data with the help of the parameters, then, if we believe the MDL Principle, with these parameters the model (network) should generalize well to future data. In modern practice, neural networks are often trained with stochastic gradient descent (SGD), and it has been empirically found that networks that generalize well do tend to have parameters lying in very flat minima.

While this use of the MDL Principle seems less precise than what we reviewed earlier in this paper, it can once again be given a frequentist justification, and this justification is mathematically precise after all: the so-called PAC-Bayesian generalization bounds⁸⁵ show that the generalization performance of any classifier can be directly linked to a quantity that gets smaller as soon as one needs (a) less bits to describe the parameter and as soon as one needs (b) less bits to describe the data given the parameters; both the results and their proofs are very similar to the MDL convergence results by Barron and Cover,⁶⁸ Zhang,^69,70 and Grünwald and Mehta.¹³ Although in general, the formulation is not as straightforward as a simple sum of the two description lengths (a) and (b), the connections between both the two-part code length and the Bayesian code length are quite strong, as was already noticed by Blum and Langford.⁸⁶ In particular, for discrete $\mathrm{\Theta }$, such PAC-Bayes bounds contain a term $-log\pi (\theta )$ which can be interpreted as the number of bits needed to encode $\theta$ using the codes based on some distribution $\pi$; for general, uncountable $\mathrm{\Theta }$, this term gets replaced by a KL divergence term that can still be related to a code length via a so-called “bits back argument” pioneered by Hinton and van Camp.⁸⁴ Dziugaite and Roy,⁸⁷ and Zhou et al.,⁸² inspired by earlier work by Langford and Caruana,⁸⁸ indeed show that, for some real-world datasets, one can predict nontrivial generalization using deep neural nets by looking at the number of bits needed to describe the parameters and applying PAC-Bayesian bounds.

7. Concluding Remarks

We have given a self-contained introduction to MDL, incorporating and highlighting recent developments. Of necessity, we had to make a choice as to what to cover in detail, and there are many things we omitted. We would like to end with briefly mentioning three additional developments. First, there has always been the question about how MDL relates to other complexity notions such as those considered in the statistical learning theory literature²⁶: Vapnik–Chervonkis dimension, entropy numbers, Rademacher complexity and so on. A major step towards understanding the relation was made by Grünwald and Mehta¹³ who show that for probability models with members of the form $p_{\theta }(z)\propto exp(-\eta {\mathrm{\text{loss}}}_{\theta }(z))$, where ${l}{o}{s}{s}$ is an arbitrary bounded loss function, the NML complexity can be precisely bounded in terms of the Rademacher complexity defined relative to ${l}{o}{s}{s}$. Second, we should note that Rissanen’s own views and research agenda have steered in a direction somewhat different from the developments we describe: Rissanen⁸⁹ published Information and Complexity in Statistical Modeling, which proposes foundations of statistics in which no underlying “true model” is ever assumed to exist. As Rissanen writes, “even such a well-meaning statement as “all models are wrong, but some are useful”, is meaningless unless some model is ‘true’.” Rissanen expands MDL and NML ideas in the direction of the Kolmogorov structure function, taking the idea of distinguishable distributions underlying Ref. 19 as the fundamental; while presumably compatible with the developments we describe here, the emphasis of this work is quite different.

We end with a word about applications: since 2007, numerous applications of MDL and MDL-like techniques have been described in the literature; as discussed in Sec. 6.2, highly popular methods such as Lasso and Bayes factor methods can often be seen as “MDL-like”. Even as to specific “pure” MDL applications (such as based on NML and two-part codes), the number and scope of applications are simply too large to give a succinct representative overview. However, there is one particular area which we would like to mention specifically, since that area had hardly seen any MDL applications before 2007 whereas nowadays such applications are flourishing: this is the field of data mining. Some representative publications are Refs. 90–92. Most of this work centers on the use of two-part codes, but sometimes NML and other sophisticated universal distributions$/$codes are used as well.⁹³

Notes

^a In this sub-section we view, for notational convenience, the elements of ${{ℳ}}_j$ as probability distributions $P_{\theta }$ with densities or mass functions $p_{\theta }$.

Acknowledgments

We would like to thank Kenji Yamanishi, Matthijs van Leeuwen and Bin Yu for providing references to new MDL work, Jun-Ichi Takeuchi for prompting us to write this paper and Jorma Rissanen and Andrew Barron for many fruitful conversations over the last 20 years. This work is part of the research programme Safe Bayesian Inference with Project No. 617.001.651, which is financed by the Dutch Research Council (NWO).