A Fast Normalized Maximum Likelihood Algorithm for Multinomial Data

1 Introduction

The Minimum Description Length (MDL) principle developed by Rissanen [Rissanen, 1978; 1987; 1996] offers a wellfounded theoretical formalization of statistical modeling. The main idea of this principle is to represent a set of models (model class) by a single model imitating the behaviour of any model in the class. Such representative models are called universal. The universal model itself does not have to belong to the model class as often is the case.

From a computer science viewpoint, the fundamental idea of the MDL principle is compression of data. That is, given some sample data, the task is to find a description or code of the data such that this description uses less symbols than it takes to describe the data literally. Intuitively speaking, this approach can in principle be argued to produce the best possible model of the problem domain, since in order to be able to produce the most efficient coding of data, one must capture all the regularities present in the domain.

The MDL principle has gone through several evolutionary steps during the last two decades. For example, the early realization of the MDL principle, the two-part code MDL [Rissanen, 1978], takes the same form as the Bayesian BIC criterion [Schwarz, 1978], which has led some people to incorrectly believe that MDL and BIC are equivalent. The latest instantiation of the MDL is not directly related to BIC, but to the formalization described in [Rissanen, 1996]. Unlike Bayesian and many other approaches, the modern MDL principle does not assume that the chosen model class is correct. It even says that there is no such thing as a true model or model class, as acknowledged by many practitioners. The model class is only used as a technical device for constructing an efficient code. For discussions on the theoretical motivations behind the modern definition of the MDL see, e.g., [Rissanen, 1996; Merhav and Feder, 1998; Barron et al., 1998; Gr ¨ unwald, 1998; Rissanen, 1999; Xie and Barron, 2000; Rissanen, 2001].

The most important notion of the MDL principle is the Stochastic Complexity (SC), which is defined as the shortest description length of a given data relative to a model class M. The modern definition of SC is based on the Normalized Maximum Likelihood (NML) code [Shtarkov, 1987]. Unfortunately, with multinomial data this code involves a sum over all the possible data matrices of certain length. Computing this sum, usually called the regret, is obviously exponential. Therefore, practical applications of the NML have been quite rare,

In our previous work [Kontkanen et al., 2003; 2005], we presented a polynomial time (quadratic) method to compute the regret. In this paper we improve our previous results and show how mathematical techniques such as discrete Fourier transform and convolution can be used in regret computation. The idea of applying these techniques for computing a single regret termwas first suggested in [Koivisto, 2004], but as discussed in [Kontkanen et al., 2005], in order to apply NML to practical tasks such as clustering, a whole table of regret terms is needed. We will present here an efficient algorithm for this specific task. For a more detailed discussion of this work, see [Kontkanen and Myllym¨ aki, 2005].

2 NML for Multinomial Data

The most important notion of the MDL is the Stochastic Complexity (SC). Intuitively, stochastic complexity is defined as the shortest description length of a given data relative to a model class. To formalize things, let us start with a definition of a model class. Consider a set Θ ∈ R d , where d is a positive integer. A class of parametric distributions indexed by the elements of Θ is called a model class. That is, a model class M is defined as

Consider now a discrete data set (or matrix) x N = (x1, . . . , xN) of N outcomes, where each outcome xj is an element of the set X consisting of all the vectors of the form (a1, . . . , am), where each variable (or attribute) ai takes on values v ∈ {1, . . . , ni}. Given a model class M, the Normalized Maximum Likelihood (NML) distribution [Shtarkov, 1987] is defined as

where ˆ θ(x N ) denotes the maximum likelihood estimate of data x N , and R N M is given by

and the sum goes over all the possible data matrices of size N. The termR N M is called the regret. The definition (2) is intuitively very appealing: every data matrix is modeled using its own maximum likelihood (i.e., best fit) model, and then a penalty for the complexity of the model class M is added to normalize the distribution.

The stochastic complexity of a data set x N with respect to a model class M can now be defined as the negative logarithm of (2), i.e.,

As in [Kontkanen et al., 2005], in the sequel we focus on a multi-dimensional model class suitable for cluster analysis. The selected model class has also been successfully applied to mixture modeling [Kontkanen et al., 1996], case-based reasoning [Kontkanen et al., 1998], Naive Bayes classification [Gr¨ unwald et al., 1998; Kontkanen et al., 2000b] and data visualization [Kontkanen et al., 2000a].

Let us assume that we have m variables, (a1, . . . , am), and we also assume the existence of a special variable c (which can be chosen to be one of the variables in our data or it can be latent). Furthermore, given the value of c, the variables (a1, . . . , am) are assumed to be independent. The resulting model class is denoted by MT . Suppose the special variable c has K values and each ai has ni values. The NML distribution for the model classMT is now

where hk is the number of times c has value k in x N , fikv is the number of times ai has value v when c = k, and R N MT ,K is the regret term. In [Kontkanen et al., 2005] it was proven that an efficient way to compute the regret term is via the following recursive formula:

where k1 + k2 = K.

As discussed in [Kontkanen et al., 2005], in order to apply NML to the clustering problem, we need to compute a whole table of regret terms. This table consists of the terms R n MT ,k for n = 0, . . . , N and k = 1, . . . , K, where K is the maximum number of clusters.

The procedure of computing the regret table starts by filling the first column, i.e., the case k = 1, which is trivial (see [Kontkanen et al., 2005]). To compute the column k, for k = 2, . . . ,K, the recursive formula (7) can be used by choosing k1 = k - 1, k2 = 1. The time complexity of filling the whole table is O ( K · N 2 ) . For more details, see [Kontkanen et al., 2005; Kontkanen and Myllym¨ aki, 2005].

In practice, the quadratic dependency on the size of data limits the applicability of NML to small or moderate size data sets. In the next section, we will present a novel, significantly more efficient method for computing the regret table.

3 The Fast NML Algorithm

In this section we will derive a very efficient algorithm for the regret table computation. The new method is based on the Fast Fourier Transform algorithm. As mentioned in the previous section, the calculation of the first column of the regret table is trivial. Therefore, we only need to consider the case of calculating the column k given the first k-1 columns. Let us define two sequences a and b by

for n = 0, . . . ,N. Evaluating the convolution of a and b gives

where the last equality follows from the recursion formula (7). This derivation shows that the column k can be computed by first evaluating the convolution (11), and then multiplying each term by n!/n n .

The standard convolution theorem states that convolutions can be evaluated via the (discrete) Fourier transform, which in turn can be computed efficiently with the Fast Fourier Transform algorithm (see [Kontkanen and Myllym¨ aki, 2005] for details). It follows that the time complexity of computing the whole regret table drops to O(N logN · K). This is a major improvement over O ( N 2 · K ) obtained by the recursion method of Section 2.

4 Conclusion And FutureWork

The main result of this paper was a derivation of a novel algorithm for the regret table computation. The theoretical time complexity of this algorithm allows practical applications of NML in domains with very large datasets. With the earlier quadratic-time algorithms, this was not possible.

In the future, we plan to conduct an extensive set of empirical tests to see how well the theoretical advantage of the new algorithm transfers to practice. On the theoretical side, our goal is to extend the regret table computation to more complex cases like general graphical models. We will also research supervised versions of the stochastic complexity, designed for supervised prediction tasks such as classification.

Acknowledgements

This work was supported in part by the Academy of Finland under the projects Minos and Civi and by the National Technology Agency under the PMMA project. In addition, this work was supported in part by the IST Programme of the European Community, under the PASCAL Network of Excellence, IST-2002-506778. This publication only reflects the authors' views.