Asymptotically Optimal Function Estimation by Minimum Complexity Criteria

Andrew R. Barron, Yuhong Yang, and Bin Yu

Abstract – The familiar minimum description length criterion for parametric model selection chooses a model that optimizes −log(likelihood) +(m/2) log n where m is the dimension of the model and n is the sample size. This criterion has been interpreted by Rissanen and others as an approximation to the length of a code for a sample, and this codelength is asymptotically optimal for samples drawn from densities in smooth finite-dimensional families. Sequences of finite-dimensional families may also be used to approximate smooth densities in nonparametric function classes. Here we show that for such smooth densities an optimal description-length criterion takes the form −log(likelihood) + const·m, for sufficiently large dimension families. The more parsimonious description of the data is achieved by restricting the candidate parameter values in each family to a grid for which the density satisfies smoothness constraints. Conditions are given such that both the statistical risk and the coding redundancy are optimal in an asymptotic minimax sense simultaneously for multiple function classes. Adaptive selection of the model size provides convergence at the optimal rate without prior knowledge of which function class contains the target function. The analysis provides a minimum description length reconciliation of the information criteria of Rissenen, Schwarz, and Akaike.

Index Terms – model selection, minimimum description-length criteria, information criteria, resolvability of function estimates, universal data compression, optimal convergence rates, adaptive inference.

I Introduction

We consider a class of model selection and distribution estimation procedures that have a minimum complexity or minimum description length (MDL) interpretation. Let {p_{m}(x_{1}, x_{2}, ..., x_{n}|θ) : θ ∈ Θ_{m} ⊂ R^{d_{m}}} be a sequence of parametric families of joint probability mass functions or joint probability density functions, where m is the model index taking values in a countable index set M and d_{m} is the dimension of parameter space of model m. These parametric families serve as models for the description of the behavior of a data sequence. In particular, each candidate probability distribution determines a uniquely decodable code. (For real-valued variables, the distributions determine codes for data represented to arbitrarily fine precision.) It is assumed that the data sequence behaves as random sample from an unknown probability distribution. The target function that we wish to estimate is the density function of this unknown distribution or other functions derived from the distribution for purposes of classification, prediction, or regression. Though analogs of what we present may hold for certain types of dependence models, it will be assumed throughout this paper that X_{1}, X_{2}, ..., X_{n} are independent random variables, that each model takes the product form p_{m}(x_{1}, x_{2}, ..., x_{n}|θ) = ∏_{i=1}^{n} p_{m}(x_{i}|θ), and that all the density functions are taken with respect to a given reference measure ν(dx) (e.g., counting measure or Lebesque measure) that is known to dominate the target distribution P (dx) = p(x)ν(dx).

As described in Rissanen [], there are several methods for coding data based on a parametric family, including two-stage codes and the mixture and predictive formulations of stochastic complexity. Here we focus on two-stage codes and the corresponding minimum description length criterion. The criterion takes a similar form to Akaike’s information criterion (AIC) [], and the Bayesian information criterion (BIC) of Schwarz [], but has a very different formulation. Here we present the description-length formulation in sufficient generality for our purposes.

Each parametric family may be defined for a continuum of parameter vectors θ, but to permit finite description we take Θ_{m} to be a countable subset of the parameter space in R^{d_{m}}. In some cases Θ_{m} = ∪_{ℓ}Θ_{m,ℓ} is a union of sets for various choices of bits of precision ℓ to which the coordinates of the parameter vector are

For submission to the IEEE Transactions of Information Theory, 1994

Andrew Barron and Yuhong Yang are with Yale University, Department of Statistics, P.O. Box 208290, New Haven, CT 06520. Bin Yu is with the University of California, Department of Statistics, Berkeley, CA represented. These subsets also may be divided according to the size and shape of suitable geometric sets in which appropriate parameter values may reside. The manner in which the parameter space is carved up will lead to natural choices for the description of parameter vectors θ. For each model m, we assign codelengths L(θ|m) for θ ∈ Θ_{m} satisfying the Kraft inequality ∑_{θinΘ_{m}} 2^{−L(θ|m)} ≤ 1. Likewise codelengths L(m) are assigned for m in M, satisfying ∑_{m} 2^{−L(m)} ≤ 1. The first stage of a code consists of the description of a model and parameters, using L(m)+L(θ|m) bits. The second stage is the Shannon code for the data given the parameterized model. For probability mass functions p_{m}(x^{n}|θ), the Shannon codelength is log 1/p_{m}(x^{n}|θ) rounded up to an integer. In general for probability densities we consider log 1/p_{m}(x^{n}|θ) as an idealized codelength. (This idealization is justified in part, when the data to be described is discretized into small cells of equal measure, by replacing the probability of a cell with the product of the density times the measure of the cell.) For each model m, the parameter estimate θ̂ is defined to achieve the minimum description length

log 1/p_{m}(x^{n}|θ̂) + L(θ̂|m) = min_{θ∈Θ_{m}}{log 1/p_{m}(x^{n}|θ) + L(θ|m)  (1.1)

The minimum description length criterion selects the model m̂ and estimates the density p̂ minimizing the total description length of the data (including terms for the description of the model, the description of the estimated parameter vector, and the description of the data given the model and estimated parameter). Thus

m̂ = argmin_{m}CRIT_{m} and p̂(x) = pm̂(x|θ̂)  (1.2)

where

CRIT_{m} = log 1/p_{m}(x^{n}|θ̂) + L(θ̂|m) + L(m).  (1.3)

The minimized description length is

C(x^{n}) = log 1/pm̂(x^{n}|θ̂) + L(θ̂|m̂) + L(m̂).  (1.4)

It is evident that q(x^{n}) = 2^{−C(x^{n})} defines a subprobability density function. In the discrete case this means that Kraft’s inequality ∑_{x^{n}} 2^{C(x^{n})} ≤ 1 is satisfied in accordance with the interpretation of C(x^{n}) as the length of a uniquely decodable code for x^{n}. The universality of the codelength C(x^{n}) is studied by examining the redundancy

red_{n}(C, p) = (1/n)E(C(X^{n}) − log1/p(X^{n})) = (1/n)Elogp(X^{n})/q(X^{n})  (1.5)

which has the interpretation of as the extra average codelength required to describe X^{n} without knowledge of p compared to the optimal codelength that would be achievable with knowledge of p, normalized by division by the sample size n. (Following [,] the quantity (1.5) retains a redundancy interpretation in the continuous sample space case: it is the least upper bound on the redundancy of codes based on the subprobability q, where the infimum is taken over arbitrary discretizations of the sample space.)

For various class of distributions p there are bounds on the best achievable redundancy uniformly over a class. A minimax setting for such analysis is developed in Davisson [] and Davisson and Garcia [] and optimal minimax redundancy rates are known in both parametric and nonparametric settings (see for instance C&B [], S&Y []). The minimax results are strengthened by Rissanen [,] who shows that for p in a smooth parametric family of dimension d one cannot have asymptotically smaller redundancy than the minimax rate (d/2)((log n)/n) except for a negligible set of parameter values. This minimax optimal parametric rate is achieved by the minimum description length criterion using parameters described to precision of order 1/\sqrt{n} using (1/2) log n bits per parameter. The corresponding MDL criterion assigns a parameter description penalty of order (d/2) log n. As we shall see other parameter precisions are needed to yield suitable estimates when p is not in any of the finite dimensional families. Upper bounds on the redundancy rate of two-stage codes are conveniently analyzed in terms of the index of resolvability as defined in general in [].

In accordance with [], one may also think of q(x^{n}) as a universal (sub)probability assignment, where the quality of this assignment is assessed by the relative entropy rate (1/n)Elogp(X^{n})/q(X^{n}). Perhaps more relevant to the question of the suitability of the estimator p̂(x) = pm̂(x|θ̂) is an examination of the distance between p̂ and p. In particular we consider the expected squared Hellinger distance Ed^{2}(p, p̂) where d^{2}(p, p̂) = ∫ ((p(x))^{1/2} − (p̂(x))^{1/2})^{2}ν(dx). This is the distance that for suitably well behaved densities is locally equivalent to the relative entropy D(p||p̂) = ∫ p(x) log p(x)/p̂(x)ν(dx) ([],[]).

As shown in [], both the redundancy red_{n}(p) (or equivalently the reative entropy rate (1/n)Elogp(X^{n})/q(X^{n})) and the expected squared Hellinger distance Ed^{2}(p, p̂) are upper bounded in terms of the index of resolvability. In the present setting the index of resolvability takes the form

R_{n}(p) = min_{m} min_{θ∈Θ_{m}}\{D(p||p_{m,θ}) + \frac{L(θ|m) + L(m)}{n}.  (1.6)

This index of resolvability captures the tradeoff between the accuracy of approximations to p and the complexity relative to sample size. The advantage for us is that for various function classes it is much easier to bound R_{n}(p) than it is to directly analyze the redundancy or the statistical risk. The redundancy is related to the resolvability by [,prop. ]: for all densities p and all n,

Red_{n}(p) ≤ R_{n}(p).  (1.7)

The statistical risk is related to the resolvability by Theorem 4 of []: there is a constant c such that for all p and all n

Ed^{2}(p, p̂) ≤ cR_{n}(p).  (1.8)

This conclusion holds under the condition that for some 0 < α < 1, ∑_{m} 2^{−αL(m)} ≤ 1 and ∑_{θ∈Θ_{m}} 2^{−αL(θ|m)} ≤

1. The constant c depends only on α. As discussed in [], if we are given codelengths L′(m) and L′(θ|m) that statisfy Kraft’s inequality then we can satisfy the more restrictive condition by dividing the codelengths by α < 1. This modification will increase the redundancy by no more than a factor of 1/α. In this paper we will only be interested in identifying the optimal rates up to an undetermined multiplicative constant. Consequently, the incorportation of the modification involving α has no adverse affect on the rates we obtain. For simplicity we will ignore this modification in what follows.

The aim of the present paper is to establish a single estimator that is simultaneously minimax rate optimal for each of the usual nonparametric function classes and for each of the parametric models. The idea is that the model selection criteria can be used to not only select an appropriate sequence of approximating families, but also to adapt the precision to which the parameters are represented, and to adapt geometric sets in which suitable parameter values occur. It is seen in certain cases that the best precision to which parameters are represented does not need to be of order 1/\sqrt{n}. Instead a precision of order (1/d_{m})^{γ} and parameters restricted to certain constraint sets leads to minimax optimal rates for densities in certain function classes, where 0 < γ < 1 and the shape and size of the constraint set depends on the class. Moreover, the description length of the parameters that satisfy these constraints will be of order d_{m} (times a constant) instead of (d_{m}/2) log n. However, to achieve the goal of simultaneous asymptotic minimax optimality we can not preselect the size of the precision nor can we preconstrain a norm on the parameters. (Indeed, the parameter estimates would be inconsistent in the parametric case if the precision did not become increasingly fine as a function of n). Instead, the criterion optimizes the location and precision of the parameters based on the data. If the data happen to behave as a sample from a density in one of the nonparametric classes (but not in any of the finite-dimensional families) the criterion may select larger size models and select a precision that leads to a parameter description length of order d_{m}. If the data happen to behave as a sample from one of the finite-dimensional models then the criterion may ultimately select the right fixed size model and it may select a precision of order 1/\sqrt{n} to provide parameter estimates that (if not fully efficient) are at least accurate at the right rate. In the latter case the penalty for the description of the estimated parameters must turn out to be of order (d_{m}/2) log n for most of the densities in a model (in accordance with the theory of Rissanen), but, to provide the flexibility to handle the nonparametric cases, it is not forced in advance to be of that fixed order.

We examine in this paper two standard types of nonparametric settings in which the phenomen occurs. The first concerns function classes of ellipsoidal type. Typical examples involve L^{2} norm constraints on derivatives of the target function. In these cases series expansion models are appropriate and suitable parameter vectors are restricted to certain ellipsoids. In the absence of knowledge of the order of derivatives and the magnitude of the L^{2} norm, an MDL criterion is used to select the radius and shape parameters of the ellipse as well as to optimize the precision. Built in to the criterion is the ability to retain an optimal rate of description in the case that the target function has a finite expansion.

The second type of nonparametric setting we examine involves function classes that demand more local control of the estimators. A typical example is the class of functions of bounded variation. In this case piecewise constant (histogram-type) estimators are appropriate provided we build in the flexibility to adjust the locations of the jumps. In the absence of knowledge of the total variation of the function, an MDL criterion is used that adapts a suitable function variation while estimating the location and heights of the jumps. This case is complicated by the need to have greater precision in the location of the jumps when the jump is large or when there are a large number of jumps in a small interval. Coding ideas are used to suitably resolve these issues. As in the previous case the convergence rate is minimax optimal and the penalty term for the description of the parameters is of the order d_{m}, i.e., there is no logarithmic factor.

An unknown probability distribution P governing a (finely discretized) random sample X_{1}, ..., X_{n} is estimated by minimizing the total length of a two-stage description: that is,

P̂_{n} = arg min_{Q∈Γ}{L(Q) + log 1/ ∏_{i=1}^{n} Q(X_{i})}  (1.1)

where Γ is a class of candidate distributions, L(Q) are codelengths for Q ∈ Γ satisfying the Kraft inequality ∑_{Q} 2^{−L(Q)} ≤ 1, and for each Q the likelihood term log 1/ ∏_{i=1}^{n} Q(X_{i}) is the length of the Shannon code for the sample (ignoring a necessary rounding to integer codelengths).

The index of resolvability [4] of the distribution P (relative to the class Γ and the codelengths L(Q)) defined by

R_{n}(P) = min_{Q∈Γ}\{\frac{L(Q)}{n} + D(P ||Q)\}

characterizes the best tradeoff between complexity and accuracy of approximation, where D(P ||Q) is the relative entropy (Kullback-Leibler information divergence). This index may be interpreted as the minimum of the redundancy of the two-stage codes. From Barron and Cover [4] it is known that the resolvability bounds the expected squared Hellinger distance between P and P̂_{n}. Here we will use the resolvability to verify that suitably chosen two-stage codelengths can achieve asymptotically optimal rates for data compression, density estimation, and regression, using sequences of parametric approximations to the unknown function.

We allow for parametric families of probability densities {p_{m}(x|θ) : θ ∈ Θ_{m} ⊂ R^{d_{m}}}, with index m in a countable index set M, where d_{m} is the dimension of the parameter space for model m. (In several of the cases we study we have one model of each dimension d_{m} = m, for m = 1, 2, ....) The description of a candidate density proceeds in stages. First we describe the model index m using L(m) bits. (e.g., if M = {1, 2, ...} we may set L(m) = log^{∗} m where log^{∗} is defined as in [] as a valid codelength for the integers). Then we describe a finite subset Θ_{m,k} indexed abstractly by k ∈ K where K is a countable set (possibly depending on m). This usually amounts to the description of a precision to which parameters are represented and the description of the size and shape of suitable geometric sets in which appropriate parameter values may reside. (For instance, in one of the important cases that follows, the coordinates of the index k = (ℓ, r, s) prescribes a precision, radius, and shape of a grid in an ellipse in R^{m}.) Typically Θ_{m,k} is dense in Θ_{m}. Finally we describe a candidate θ ∈ Θ_{m,k}. In general we allow arbitrary codelengths L(m), L(k|m), and L(θ|m, k) satisfying the Kraft inequalities: ∑_{m} 2^{−L(m)} ≤ 1, ∑_{k} 2^{−L(k|m)} ≤ 1, and ∑_{θ∈Θ_{m,k}} 2^{−L(θ|m,k)} ≤ 1. In particular, it is permitted to take this last description length to be equal to the log cardinality of Θ_{m,k}, that is, L(θ|m, k) = log CardΘ_{m,k}.

The minimum description-length principle leads to a selected model m̂, a selected set Θ_{m̂,k̂}, a parameter estimate θ̂ in Θ_{m̂,k̂}, and a density estimate p̂ = pm̂(·|θ̂) achieving

(m̂, k̂, θ̂) = arg min_{m,k,θ∈Θ_{m,k}} {log 1/p_{m}(X^{n}|θ) + L(θ|m, k) + L(k|m) + L(m)}  (1.3)

. Here p_{m}(X^{n}|θ) = ∏_{i=1}^{n} p_{m}(X_{i}|θ). As explained in [4] when a term is added to account for an arbitrarily high precision to which the data X_{i} are observed there is an interpretation of log 1/p_{m}(X^{n}|θ) as the length of a Shannon code for the data given the model and parameters. The other terms account for the bits of description of the model and parameters. If we let [X^{n}] = ([X_{1}], ..., [X_{n}]) denote the cells containing the observations in an arbitrary discretization of the sample space, then the minimum description-length principle yields a universal code of length

codelength[X^{n}] = min_{m,k,θ∈Θ_{m,k}} {log 1/P_{m}([X^{n}]|θ) + L(θ|m, k) + L(k|m) + L(m)} , (1.4)

where we are ignoring a necessary rounding to integer codelengths. The redundancy of the code when X^{n} is destributed according to P as the extra codelength beyond the entropy necessitated by lack of knowledge of P. It is defined by

redundancy_{n}(code, P) = (1/n)E{codelength[X^{n}] − log1/P [X^{n}]}.(1.5)

One of the results in [4] shows that the index of resolvability bounds the redundancy of the given code, that is, redundancy_{n}(code, P) ≤ R_{n}(P), uniformly over all possible discretizations of the sample space.

In cases popularized by Jorma Rissanen, Γ consists of the union of a sequence of finite dimensional smooth parametric families of densities p_{m,θ} with parameter vector θ restricted to a grid spaced at width proportional to 1/\sqrt{n}, and the dominate term in the codelength L(p_{m,θ}) is (m/2) log n where d_{m} is the dimension of the parameter vector for the mth family. These encodings achieve the minimum order redundancy for distributions in such families (Rissanen [6,7]).

On the other hand, when the sequence of families provides approximation to the density in the limit of large dimensions (as in the case of series expansion of the log-density [5]), the usual minimum description length criterion leads to a redundancy and a statistical risk of order ((log n)/n)^{2s/(2s+1)} compared to the optimal convergence rate of order (1/n)^{2s/(2s+1)}, as discussed in [4]. Here s is a number of bounded derivatives of the log-density.

It is shown that more parsimonious two-stage descriptions are possible by taking advantage of possible smoothness of the density function. The idea is to use for each family p_{m,θ} a union of constrained grids of candidate parameter values, one for each possible order of smoothness. This can lead to a considerably smaller penalty than the familiar (m/2) log n, while retaining an accurate description length interpretation. In particular, we consider cases involving orthogonal polynomials and trigonometric series. For each possible order of smoothness s and integer r > 0, consider the ellipse of parameter values {θ : ∑_{k=1}^{m} θ_{k}^{2}k^{2s} ≤ r^{2}}. If we take a grid spaced at width proportional to (1/m)^{s+1/2} in each coordinate, then the logarithm of the number of grid points in the ellipse is of order m (instead of (m/2) log n).

A suitable set of encodings and selection criterion is defined which includes data-based selection of the order of smoothness s, r as well as the size of the model m. The penalty term on the selection criterion becomes order m instead of (m/2) log n. When large size models are needed to get a good approximation, a better overall description length is obtained than with naive use of the minimum description length criterion.

Examination of the index of resolvability shows that the optimal convergence rate is now achieved for densities in the nonparametric smoothness classes. In particular, suppose the logarithm of the density function p, on a given interval, has s square-integrable derivatives with ∫ (D^{s}logp)^{2} ≤ r. This optimum rate of order (1/n)^{2s/(2s+1)} is achieved with the best sequence of model sizes m, where the two terms in the resolvability bound are of the form

O (\frac{m}{n}) + O (\frac{1}{m})^{2s} .

Here O(m/n) bounds the codelength divided by the sample size, for the best candidate densities in the family, and O(1/m)^{2s} is the order of the relative entropy distance between the density and its best approximation in the m’th family. The parameter values in an approximation achieving this rate are in an ellipse of the type considered above and the specified grid widths are of the right size to achieve the same order of best approximation while restricting the parameters to the grid. The desired consequence for the risk of the minimum complexity density estimator follows from the resolvability bound, since E(d^{2}(p, p̂_{n})) = O(R_{n}(p)) from [4], where d(p, p̂_{n}) is the Hellinger distance. We note that the optimal rate is obtained without prior knowledge of the best smoothness class s or the norm r. Using techniques in [3], the results are adapted to give similar conclusions for minimum complexity estimators in nonparametric regression settings. For orthogonal polynomials and trigonometric series the ellipsoids with ∑_{k} θ_{k}^{2}k^{2s} ≤ r^{2} correspond to L_{2} bounds on derivatives of the function. The theory extends to give optimal convergence rates for minimum complexity estimators of functions with orthogonal expansions satisfying arbitrary computable ellipsoidal constraints.

Similar results can be given for models with parameters satisfying L_{1} constraints as in the case of piecewise constant approximation and histogram density estimation in the case of a functions of bounded derivative. In this setting, two-stage codes are given that achieve asymptotic performance comparable to that achieved by the predictive stochastic complexity criterion in Rissanen, Speed and Yu [8] and in Speed and Yu [9]. Once again, the penalty term for the description of the histograms that satisfy suitable constraints is of order m rather than (m/2) log n. It is under investigation whether the optimal convergence rate can also be given for the class of densities with bounded variation using a two-stage code. Here it becomes critical to adapt the histogram bins. The minimum description length criterion is used to select the bin locations as well as the number of bins. The index of resolvability is used to bound the risk of the resulting estimator.

A code bit can be incorporated in the criterion to allow for data-based selection between the traditional MDL codelengths and the new codelengths. In this way one criterion is obtained that achieves the asymptotic optimality properties in both parametric and nonparametric cases. For smooth densities that are not in any of the given finite dimensional families, the criterion switches between prefering the (m/2) log n penalty and a penalty of order m, where the switch occurs at a model size depending on the data. In the parametric case, the selected model size stays bounded, and the traditional MDL penalty is prefered for all large n: indeed, the grid widths shrink as n grows to get a better likelihood value in the two-stage codelength.

The results have implications for the longstanding statistical debate between the relative merits of Akaike’s AIC criterion (with penalty of m) [1] and Schwartz’s BIC criterion (with penalty of (m/2) log n) [2]. It is revealed that both AIC and BIC can correspond to information-theoretic codelength in certain cases. The minimum description length principle leads to one criterion that retains the nonparametric advantages of AIC and the parametric advantages of BIC, depending on which leads to an overall shorter codelength for the data.

References

[1] H. Akaike (1973). “Information theory and an extension of the maximum likelihood principle.” In Proc. 2nd International Symposium on Information Theory P. N. Petrov and F. Csaki, Eds. Budapest: Akademia Kiado, pp. 267-281.

[2] G. Schwarz (1978). “Estimating the dimension of a model.” Annals of Statistics vol. 6, pp. 461–464.

[3] A. R. Barron (1991). “Complexity regularization with applications to artificial neural networks.” In Nonparametric Functional Estimation, G. Roussas, Ed. Boston, MA and Dordrecht, the Netherlands: Kluwer Academic Publishers. pp. 561–576.

[4] A. R. Barron and T. M. Cover (1991). “Minimum complexity density estimation.” IEEE Transactions on Information Theory, vol. 37, pp. 1034–1054.

[5] A. R. Barron and C.-H. Sheu (1991). “Approximation of density functions by sequences of exponential families.” Annals of Statistics, vol. 19, pp. 1347–1369.

[6] J. Rissanen (1983). “A universal prior for integers and estimation by minimum description length.” Annals of Statistics, vol. 11, pp. 416-431.

[7] J. Rissanen (1984). “Universal coding, information, prediction, and estimation.” IEEE Transactions on Information Theory, vol. 30, pp. 629-636.

[8] J. Rissanen, T. Speed, and B. Yu (1992). “Density estimation by stochastic complexity.” IEEE Transactions on Information Theory, vol. 38, pp. 315–323.

[9] T. Speed and B. Yu (1993). “Data compression and histograms.” Probability theory and related fields. To appear.