Information-Corrected Estimation: A Generalization Error Reducing Parameter Estimation Method

Abstract

Modern computational models in supervised machine learning are often highly parameterized universal approximators. As such, the value of the parameters is unimportant, and only the out of sample performance is considered. On the other hand much of the literature on model estimation assumes that the parameters themselves have intrinsic value, and thus is concerned with bias and variance of parameter estimates, which may not have any simple relationship to out of sample model performance. Therefore, within supervised machine learning, heavy use is made of ridge regression (i.e., L2 regularization), which requires the the estimation of hyperparameters and can be rendered ineffective by certain model parameterizations. We introduce an objective function which we refer to as Information-Corrected Estimation (ICE) that reduces KL divergence based generalization error for supervised machine learning. ICE attempts to directly maximize a corrected likelihood function as an estimator of the KL divergence. Such an approach is proven, theoretically, to be effective for a wide class of models, with only mild regularity restrictions. Under finite sample sizes, this corrected estimation procedure is shown experimentally to lead to significant reduction in generalization error compared to maximum likelihood estimation and L2 regularization.

1. Introduction

Kullback and Leibler [1] showed that minimizing a divergence ${\rho }_{KL}(f,g_{\mathit{\theta }})$ between the truth, f, and a parametric model density, $g_{\mathit{\theta }}$, is necessary and sufficient for making accurate predictions about data using the model defined by $\mathit{\theta }$. Recent work [2] on Berk–Nash equilibria has shown the central role that KL divergence plays in game theoretic choice models such as multi-armed bandits and stochastic multi-party games. KL divergence thus plays a leading role in machine learning and neuroscience, with several inferential approaches developed in the information theory literature. Such approaches for minimizing KL divergence employ a range of methods, including data partitioning, Bayesian indirect inference and M-estimation [3,4,5]. These approaches are quite distinct from the standard penalized loss minimization framework and, as such, are non-trivial to combine with supervised learning methods such as neural networks.

It is well known that maximum likelihood estimation (MLE) introduces an asymptotic bias in the KL divergence minimizer which is problematic for both model estimation and model selection. For many models, where the parameters $\mathit{\theta }$ are themselves important, this may be investigated as parameter bias and parameter variance. However, for models common in modern machine learning, the parameters themselves do not have any easily interpreted meaning. For these models, the parameters themselves are irrelevant and only the accuracy (in terms of KL divergence) of the model predictions matter. Within the information theory literature, this has often been referred to simply as bias (e.g., $b\left(G\right)$ from [6]). To distinguish it from parameter bias, one might refer to it as “prediction bias” or “generalization error”. Generalization error is the more common terminology (see, for example, Equation 1.1.6 [7]) and will be used here.

Before the widespread use of machine learning, most models had interpretable parameters, and thus there is a large literature focused on reducing parameter bias. For instance, the jackknife [8] (leave-one-out cross-validation) estimator is an early example. More relevant to this paper is the approach of Firth [9] and later Kosmidis [10,11]. More recently, Pagui, Salvan, and Sartori [12] proposed a parameter bias reducing estimation methodology. An extensive review of the literature around this point can be found in [13]. Unfortunately, these approaches do not consider the impact on KL divergence-based generalization error and thus are not applicable to the field of machine learning where the parameters themselves are devoid of meaning. Heskes [14] shows that classifiers do have a notion of bias-variance decomposition for generalization error, but it is not computable from parameter bias and parameter variance. Therefore, parameter bias reducing formulations are not useful within machine learning unless it can be shown that they also reduce generalization error.

In fact, to seat the approach taken in this paper to generalization error, we recall much earlier and seminal work at the intersection of statistics and information theory. Akaike [15], and later Takeuchi [16], proposed information criteria (AIC and TIC, respectively) for model selection designed explicitly to reduce generalization error. Konishi and Kitagawa [6] extended the approach of Takeuchi to cases where MLE was not used to fit the underlying model, but still restricted themselves to the question of model selection. Stone [17] proved that Akaike’s Information Criterion (AIC) is asymptotically equivalent to jackknifing when the estimator is finite. Takeuchi himself showed that TIC is an extension of AIC with fewer restrictions, and thus it too is equivalent to jackknifing whenever AIC would be valid.

For highly parameterized models, as are common in machine learning, model selection such as this is of limited utility. The parameter count may necessarily be very large, and thus none of the models fit using MLE may be acceptable. Then, merely choosing among them is unlikely to produce acceptable results. Within this field, typically $L_2$ or similar regularization is used to reduce generalization error. See Section 11.5.2 [18], for a typical example. For a more recent innovation, refer to [19]. Note that regularization schemes such as this often increase parameter bias while decreasing generalization error. Golub, Heath, and Wahba [20] showed that $L_2$ regularization is asymptotically equivalent to cross-validation for linear models, subject to certain assumptions. For nonlinear models, it has long been known that $L_2$ regularization is not always valid, and it is trivial to construct example models (See Section 4.1 for one such example) where this approach is always harmful in expectation.

Therefore, it is important to develop a method to reduce generalization error in model estimation analogous to the way that $L_2$ regularization would commonly be used for a highly parameterized model, but having applicability for a wider family of models, especially those for which $L_2$ regularization is not applicable. It is not the goal of this paper to perform a wide survey of generalization error reducing approaches, but we will rather propose an additional approach, investigate its properties, and show that it has superior performance when compared against $L_2$ regularization, which is currently the dominant generalization error reducing estimation procedure within the field of machine learning.

To this end, this paper introduces a generalization error reducing estimation approach referred to as Information Corrected Estimation (ICE). This estimator is proven to have a generalization error of only $O\left(n^{-\frac{3}{2}}\right)$ instead of $O\left(n^{-1}\right)$ as is the case for MLE, and is shown to be valid within a neighborhood around the MLE parameter estimate. Optimizing over this ICE objective function instead of the negative log likelihood thus produces parameters with superior out of sample performance.

Takeuchi’s TIC and Firth’s approach have never seen widespread use due to the computational and numerical issues that arise from the computation of this adjustment [21], and the ICE estimator in its raw form would have similar problems. Therefore, this paper also proposes an efficient approximation of this correction term, and shows through numerical experiments that the approximation is effective at improving model performance across a range of models.

2. Preliminaries

Let us assume that we have data ${\mathit{x}}_n:=\{x_1,\cdots ,x_n\}$ generated from an unknown joint density function $f\left(\mathit{x}\right)$ of $X_n:=\{X_1,\cdots ,X_n\}$. Where necessary, we define $Z_n$ to denote a second sample drawn from $f\left(\mathit{x}\right)$, independent of $X_n$, and ${\mathit{x}}_n^{\prime }$ is the observed realization of $Z_n$. We consider a model ${\mathcal{M}}_p$ given by a parametric family of densities ${\mathcal{M}}_p:=\left\{g\right(·\left|\mathit{\theta }\right)|\mathit{\theta }\in \mathrm{\Theta }\subseteq {\mathbb{R}}^p\}$, for some compact Euclidean parameter space $\mathrm{\Theta }$, which is misspecified and hence excludes the truth f. Henceforth, the distribution over x identified by $\mathit{\theta }$ may be referred to as $g_{\mathit{\theta }}\left(\mathit{x}\right):=g\left(\mathit{x}\right|\mathit{\theta })$ where it is notationally convenient to do so.

Suppose that ${\mathit{\theta }}_0$ is the quasi-true parameter of model $\mathcal{M}$, and $\widehat{\mathit{\theta }}\left(X_n\right)$ is the random variable representing the MLE of ${\mathit{\theta }}_0$ fit on a dataset, ${\mathit{x}}_n$. The negative log-likelihood of $X_n$ under the distribution $g_{\mathit{\theta }}$ is

$$-\ell (\mathit{\theta },X_n):=-\frac{1}{n}\sum _{i=1}^{n}log{g}_{\mathit{\theta }}\left({\mathit{x}}_i\right),$$

(1)

where $-\ell (\mathit{\theta },X_n)$ is written including a $\frac{1}{n}$ to make the expectation of this quantity $O\left(1\right)$ and asymptotically independent of n. Similarly, the minus sign is incorporated because $-\ell (\mathit{\theta },X_n)$ is a strictly non-negative quantity if $g_{\mathit{\theta }}\left({\mathit{x}}_i\right)$ is a probability. The MLE, $\widehat{\mathit{\theta }}\left(X_n\right)$, minimizes the negative log likelihood of the data set with respect to the model:

$$\widehat{\mathit{\theta }}\left({\mathit{x}}_n\right):=\underset{\mathit{\theta }}{argmin}[-\ell (\mathit{\theta },{\mathit{x}}_n)].$$

(2)

The expectation of $-\ell (\mathit{\theta },X_n)$ is the cross entropy between f and $g_{\mathit{\theta }}$:

$$-\mathcal{L}\left(\mathit{\theta }\right):={\mathbb{E}}_{X_n}[-\ell (\mathit{\theta },X_n)].$$

(3)

Here, the expectation is a function only of $\mathit{\theta }$ and of the distribution f that generated the data $X_n$. As a function of the distribution f, this value is $O\left(1\right)$, but could be large for poorly conditioned f. The quasi-true parameter ${\mathit{\theta }}_0$ is

$${\mathit{\theta }}_0:=\underset{\mathit{\theta }}{argmin}[-\mathcal{L}\left(\mathit{\theta }\right)].$$

(4)

Generalization Error in KL Divergence Based Loss Functions

Kullback and Leibler [1] viewed “information” as discriminating the sample data drawn from one distribution against another, and defined the KL-divergence ${\rho }_{KL}$ between distributions in terms of the ability to make predictions about one by knowing the other. Here,

$${\rho }_{KL}(f,g_{\mathit{\theta }})=\int log\left[\frac{f\left(x\right)}{g_{\mathit{\theta }}\left(x\right)}\right]f\left(x\right)dx.$$

(5)

This value is in general unknowable, but given a sample $X_n$ from f, $-\ell (\mathit{\theta },X_n)$ will converge asymptotically to ${\rho }_{KL}(f,g_{\mathit{\theta }})$ plus an additive constant that depends only on f. The convergence relies on White’s regularity conditions [22].

A well known result by Stone [17] shows that the MLE is a biased estimator of the minimum KL-divergence:

$${\mathbb{E}}_{X_n}[-\ell (\widehat{\mathit{\theta }}\left(X_n\right),X_n)]<{\mathbb{E}}_{X_n}[-\ell ({\mathit{\theta }}_0,X_n)],$$

(6)

because it is evaluated on the data $X_n$ which was used to fit $\widehat{\mathit{\theta }}$. Cross-validation was developed as a model selection technique to select a model from a group that actually minimizes ${\mathbb{E}}_{X_n}[{\rho }_{KL}(g_{{\mathit{\theta }}_0},g_{\widehat{\mathit{\theta }}(X_n)})]$ and not merely ${\mathbb{E}}_{X_n}[-\ell (\widehat{\mathit{\theta }}(X_n),X_n)]$ in the limit of large n. Takeuchi [16] and Akaike [15] explicitly modeled this bias (generalization error) of an estimation procedure $\mathit{\theta }(X_n)$ as

$$b:={\mathbb{E}}_{X_n}\left[\ell (\mathit{\theta }(X_n),X_n)-{\mathbb{E}}_{X_n^{\prime }}[\ell (\mathit{\theta }(X_n),X_n^{\prime })]\right].$$

(7)

Our goal is to obtain an estimate, $b^{\ast }$, of the generalization error b without using the MLE. We will then add this term to the objective function to develop the estimator ${\mathit{\theta }}^{\ast }(X_n)$ so as to cancel the lower order terms of the generalization error. This estimator will then minimize ${\mathbb{E}}_{X_n}[{\rho }_{KL}(g_{{\mathit{\theta }}_0},g_{{\mathit{\theta }}^{\ast }(X_n)})]$ more effectively than MLE, and potentially would in turn produce improved predictions from the model fitted over finite training sets.

Remark 1.

We note that under MLE, $b=O(\frac{1}{n})$ [16]. Equivalently, one could say that a particular realization of the generalization error $\ell (\mathit{\theta }(X_n),X_n)-{\mathbb{E}}_{X_n^{\prime }}[\ell (\mathit{\theta }(X_n),X_n^{\prime })]$ is itself $O_p(\frac{1}{n})$. Here, $O_p(\frac{1}{n})$ is used to indicate that the quantity is a random variable with finite variance, whose mean is $O(\frac{1}{n})$.

3. Information Corrected Estimation (ICE)

We propose the following penalized likelihood function:

Definition 1

(ICE Objective).

$$-{\ell }^{\ast }(\mathit{\theta })=-\ell (\mathit{\theta })+\frac{1}{n}tr(I_{\mathit{\theta }}{J}_{\mathit{\theta }}^{-1}),$$

(8)

where $J_{\mathit{\theta }}$ is the negative expected Hessian

$$J_{\mathit{\theta }}:=-{\mathbb{E}}_X[{\partial }_{\mathit{\theta }}^{2}logg(X|\mathit{\theta })]=-\int f(x){\partial }_{\mathit{\theta }}^{2}logg(x|\mathit{\theta })dx,$$

(9)

and $I_{\mathit{\theta }}$ is the Fisher Information matrix

$$I_{\mathit{\theta }}:={\mathbb{E}}_X[{\partial }_{\mathit{\theta }}logg(X|\mathit{\theta }){\partial }_{{\mathit{\theta }}^T}logg(X|\mathit{\theta }).$$

(10)

with ${\widehat{I}}_{\mathit{\theta }}$, ${\widehat{J}}_{\mathit{\theta }}$ being their estimates over the data.

Let ${\mathit{\theta }}^{\ast }$ denote the minimizer of (8).

The trace term in Equation (8) will be familiar from Takeuchi [16]. However, Takeuchi showed only that this was the leading order of the bias for the MLE estimate $\widehat{\mathit{\theta }}$, and therefore the proof found there is not sufficient to justify a new estimator that will itself be the target of optimization, and is required to be valid away from $\widehat{\mathit{\theta }}$. As in Takeuchi, because I and J are unknowable, we will substitute their approximations computed from the training data, $\widehat{I_{\mathit{\theta }}}$ and $\widehat{J_{\mathit{\theta }}}$ during the actual computation of this objective. The numerical impact of this approximation will be examined in Section 4.2.1.

Remark 2.

Though AIC was developed before TIC, it is easily reproduced as a special case of TIC. Subject to certain conditions (guaranteed by the requirements of [15]), at least in expectation, $I_{\widehat{\mathit{\theta }}}=J_{\widehat{\mathit{\theta }}}$. Thus, the quantity within the TIC trace term, $I_{\widehat{\mathit{\theta }}}{J}_{\widehat{\mathit{\theta }}}^{-1}$, is the identity matrix. Therefore, its trace is equal to p, the parameter count of the model, recovering AIC. TIC itself can be derived using a proof that is similar to, though somewhat simpler than, the one we include in (A2), of which Takeuchi’s proof is a special case that is valid only at the MLE estimate $\widehat{\mathit{\theta }}$.

We also define ${\widehat{J}}^{\ast }$ to be the negative hessian of $-{\ell }^{\ast }(\mathit{\theta })$ rather than $-\ell (\mathit{\theta })$, and similarly for ${\widehat{I}}^{\ast }$, with expectations written as $J^{\ast }$ and $I^{\ast }$. Analogously, $-{\mathcal{L}}^{\ast }(\mathit{\theta })$ is the expectation of $-{\ell }^{\ast }(\mathit{\theta })$ and ${\mathit{\theta }}^{\ast }$ is the minimizer of $-{\ell }^{\ast }(\mathit{\theta })$, while ${\mathit{\theta }}_0^{\ast }$ is the minimizer of $-{\mathcal{L}}^{\ast }(\mathit{\theta })$.

We refer to the estimation of ${\mathit{\theta }}^{\ast }$, by minimization of this corrected likelihood function as Information-Corrected Estimation (ICE). As the terminology suggests, we depart from the corrective approach used in Information Criterion, by directly minimizing the bias corrected likelihood function. Note that unlike $L_2$ regularization, the correction term is parameter-free and thus would not require cross validation to estimate a hyperparameter such as the $\lambda$ used by $L_2$.

General properties of this estimator are proved, and a set of regularity conditions are provided such that the estimator is asymptotically normal, and produces a bias that is $O_p(n^{-3/2})$ instead of the usual $O_p(n^{-1})$. Though this adds only a half-order to the bias correction, for most problems with reasonably large n, any increase in order is likely to greatly reduce bias. Experimental results demonstrate superior properties of ICE for linear models compared to MLE with and without $L_2$ regularization.

Remark 3.

For models satisfying White’s regularity conditions (See [22]), it is known that $J_{{\mathit{\theta }}_0}$ is positive definite (thus non-singular) and continuous, and also that $I_{{\mathit{\theta }}_0}$ is continuous with respect to θ. Therefore, $\frac{1}{n}tr(I_{{\mathit{\theta }}_0}{J}_{{\mathit{\theta }}_0}^{-1})$ would always be well defined in an open region around ${\mathit{\theta }}_0$. Similarly, the solution ${\mathit{\theta }}^{\ast }$ would be expected to have the same properties, and hence (for large enough n) the estimate $\frac{1}{n}tr({\widehat{I}}_{{\mathit{\theta }}^{\ast }}{\widehat{J}}_{{\mathit{\theta }}^{\ast }}^{-1})$ would be well defined when computed using the estimates ${\widehat{I}}_{{\mathit{\theta }}^{\ast }}$ and ${\widehat{J}}_{{\mathit{\theta }}^{\ast }}$.

Remark 4.

N.B: Though $-{\ell }^{\ast }(\mathit{\theta })$ is an estimator of $\mathcal{L}(\mathit{\theta })$ accurate to within $O(n^{-\frac{3}{2}})$, that does not mean that $\mathcal{L}({\mathit{\theta }}^{\ast })$ is reduced by any particular amount relative to $\mathcal{L}(\widehat{\mathit{\theta }})$. We expect that using this corrected objective will always (if it can be calculated accurately) generate some improvement by virtue of more accurately representing the true performance of the model out of sample, but there is no proof that this level of improvement has any particular form or asymptotic behavior.

Our approach preserves the linear complexity of training with respect to n. However, the computation of ${\widehat{J}}_{{\mathit{\theta }}^{\ast }}^{-1}$ at each iteration of the numerical solver requires the inversion of a symmetric positive definite matrix with a complexity of $O(p^3)$. Hence the approach is not suitable for high dimensional datasets without adjustment. See Section 5 for optimized approximations that are viable for larger parameter counts. Further exploration of large models based on this approach are beyond the scope of the present work.

Remark 5.

It is clear from inspection that if $-\ell (\mathit{\theta })$ is strictly convex, then so too is $-{\ell }^{\ast }(\mathit{\theta })$ for large enough n.

We first provide a proof of asymptotic convergence of ${\mathit{\theta }}^{\ast }$ under certain regularity conditions. With this convergence result in place, we then show that minimizing (8) leads to an $O(n^{-3/2})$ bias term, an improvement over the $O(\frac{1}{n})$ term produced by MLE.

Local Behavior of the ICE Objective

Suppose the following conditions hold:

• $\mathcal{M}$ satisfies White’s regularity conditions A1–A6 (see Appendix A.1 or [22]).
• ${\mathit{\theta }}_0$ is a global minimum of $-\mathcal{L}(\mathit{\theta })$ in the compact space $\mathrm{\Theta }$ defined in A2.
• There exists a $\epsilon >0$ such that $-\mathcal{L}({\mathit{\theta }}_0)<-\mathcal{L}({\mathit{\theta }}_1)-\epsilon$ for all other local minima ${\mathit{\theta }}_1$.
• For $k=0,1,2,3,4,5$ the derivative ${\partial }_{\mathit{\theta }}^k\mathcal{L}(\mathit{\theta })$ exists, is continuous, and bounded on an open set around ${\mathit{\theta }}_0$.
• For $k=0,1,2,3,4,5$, the variance $\mathbb{V}[{\partial }_{\mathit{\theta }}^k\ell (\mathit{\theta },X_n)]\to 0$ as $n\to \infty$ on an open set around ${\mathit{\theta }}_0$.

Then for sufficiently large n there exists a compact subset $U\subset \mathrm{\Theta }$ containing ${\mathit{\theta }}_0,\widehat{\mathit{\theta }},$ such that:

• For $k=0,1,2,3$ the derivative ${\partial }_{\mathit{\theta }}^k{\ell }^{\ast }(\mathit{\theta },{\mathit{x}}_n)$ exists, is continuous, and bounded on U, almost surely.
• For $k=0,1,2,3$, $\mathbb{V}[{\partial }_{\mathit{\theta }}^k{\ell }^{\ast }(\mathit{\theta },X_n)]\to 0$ as $n\to \infty$ on U, almost surely.
• ${\mathit{\theta }}^{\ast }\in U$ as $n\to \infty$ almost surely.
• $\sqrt{n}({\mathit{\theta }}^{\ast }-{\mathit{\theta }}_0^{\ast })\to N(0,{({\widehat{J}}_{{\mathit{\theta }}_0^{\ast }}^{\ast })}^{-1}{\widehat{I}}_{{\mathit{\theta }}_0^{\ast }}^{\ast }{\left({\widehat{J}}_{{\mathit{\theta }}_0^{\ast }}^{\ast }\right)}^{-1})$ almost surely.
• $-\mathcal{L}\left(\widehat{{\mathit{\theta }}^{\ast }}\right)=-{\ell }^{\ast }({\mathit{\theta }}^{\ast }(X_n),X_n)+O_p(n^{-3/2})$ almost surely.

Items (1–3) follow from Lemma A1 (see Appendix A.2). These are additional regularity conditions that are prerequisites for later theorems.

Item (4) follows from Theorem A1 in Appendix A.3. This states that the estimate ${\mathit{\theta }}^{\ast }$ is asymtotically normal in a way that is analogous to classical asymptotic normality results for MLE. It is only true almost surely because results (1–3) upon which it relies are only true almost surely.

Item (5) follows from Theorem A2 in Appendix A.4. This item establishes the superior accuracy of the ICE objective compared to the MLE objective function in predicting out of sample errors. Like item (4) this is only true almost surely because intermediate results on which it relies are only true almost surely.

The reduction in generalization error seen arises from the optimization over the superior ICE objective function, analogous to the way that $L_2$ regularization is used for this purpose.

Remark 6.

The regularity conditions described here are only slightly more strict than the conditions described by White [22]. In particular, models having three continuous derivatives as required by White, but not 5 as needed here are thought to be very rare. Requirement (2) is just the definition of ${\mathit{\theta }}_0$, which White labels differently, and requirement (3) excludes a pathological corner case, the further study of which is beyond the scope of this paper.

Remark 7.

Note that as $-\ell (\mathit{\theta },{\mathit{x}}_n)$ is convex in the neighborhood of ${\mathit{\theta }}_0$, so too is $-{\ell }^{\ast }(\mathit{\theta })$ for large enough n because $-{\ell }^{\ast }(\mathit{\theta })\to -\ell (\mathit{\theta })$. Thus it can be concluded that the local behavior of $-{\ell }^{\ast }$ in the neighborhood of ${\mathit{\theta }}_0$ is not appreciably worse than the behavior of $-\ell$ if the problem is not too ill conditioned.

4. Direct Computation Results

The following experiments have been designed to compare MLE, MLE with $L_2$ regularization, and ICE for regression. Each experiment involves simulation of training and test sets and is implemented in R. See the attached code to run each experiment.

Each of these experiments has been performed using the raw formula for ICE provided in Equation (8) with minimal adjustments. All gradients are computed using R’s default finite difference approach. This means that for a model with p parameters, the objective function is dominated by the inversion of J, which costs $O(p^3)$ time and $O(p^2)$ space. The use of finite difference gradients further increases the time complexity to $O(p^4)$, compounding the problem. This approach is therefore viable for small models with few parameters, but not realistic for larger models. Optimizations to overcome this limitation will be considered in upcoming Section 5. The use of finite difference derivatives was not found to produce appreciable numerical differences in the final output, so analytic derivatives were not used for this analysis.

The code and results for this section is provided in [23]. Throughout this section, the following estimators will be compared.

$$\begin{array}{ccc}\mathrm{MLE}& \widehat{\mathit{\theta }}(Z_n)& :={argmin}_{\mathit{\theta }}[-\ell (\mathit{\theta })]\\ L_2\mathrm{regularization}& {\mathit{\theta }}_{L_2}^{\ast }(Z_n)& :={argmin}_{\mathit{\theta },\lambda }{[-\ell (\mathit{\theta })+\lambda \parallel \mathit{\theta }\parallel }_2^2]\\ \mathrm{ICE}& {\mathit{\theta }}_{ICE}^{\ast }(Z_n)& :={argmin}_{\mathit{\theta }}[-\ell (\mathit{\theta })+\frac{1}{n}tr({\widehat{I}}_{\mathit{\theta }}{\widehat{J}}_{\mathit{\theta }}^{-1})]\end{array}$$

4.1. Gaussian Error Model

We begin by considering the simplest case of univariate linear regression with Gaussian residuals. The advantage of this simple model is that the exact form of the correction term can be derived analytically and aids therefore in building intuition on its behavior. For such a toy model, $y\sim N(\mu ,{\sigma }^2)$ and, for simplicity, the following example will consider $\mu$ to be a constant, but it is equally applicable if $\mu =\mu (\mathit{x})$. Consider the parameters of the model to therefore be $\mathit{\theta }:=(\mu ,\sigma )$ with their optimal values being ${\mathit{\theta }}_0:=({\mu }_0,{\sigma }_0)$. The the probability density function is

$$g(y,\mathit{\theta })=\frac{1}{\sqrt{2\pi \sigma }}{e}^{-\frac{{(y-\mu )}^2}{2{\sigma }^2}}.$$

(11)

It is known a priori that $L_2$ regularization cannot improve this model, as if ${\mu }_0\ne 0$, any decrease in the magnitude of $\mu$ is likely to be systematically harmful. Similarly, a decrease in $\sigma$ below $\widehat{\sigma }$ results in a decrease in model distribution entropy, and hence would be generally making overfitting worse, and would generate a correspondingly higher KL-divergence than the MLE estimate. Consequently, we would expect any $\lambda$ computed through cross-validation to be statistically indistinguishable from zero, and $L_2$ regularization to be generally harmful whenever $\lambda \ne 0$.

Generalization Error Analysis

The Gaussian model described was generated with ${\mu }_0=0.2$, ${\sigma }_0=0.2$, and $dy=0.001$. For each of $n\in \{16,32,64,128,256,512,1024\}$, 500 independent simulations of the data $y_1,\cdots ,y_n$ were performed, and then the parameters were fit from that data. In each simulation, $\mathit{\theta }$ was computed using MLE, MLE with $L_2$ regularization, and ICE. The $\lambda$ parameter for $L_2$ regularization was computed using 2-way cross-validation on the available data, and as expected, none of the computed values of $\lambda$ were statistically different from zero.

For each estimate of $\mathit{\theta }$, the KL-divergence ${\rho }_{KL}(f,g_{\mathit{\theta }})$ was computed (using the known value of ${\mathit{\theta }}_0$), and the results were compared. The ICE parameter estimation method showed statistically significant improvement over MLE at the 5-sigma level out to $n=64$, and was improved by just under 1-sigma at $n=1024$.

The KL-divergence results graphed against n on a log-log scale are shown in Figure 1. Every value of n is normalized by the average KL divergence of the MLE methodology to improve legibility. The $L_2$ series is statistically indifferent from the MLE series at 2 standard deviations beyond $n=32$, and the two are not materially different for any n. The ICE series is at least $4.5$ standard deviations below the MLE series until $n=1024$.

Remark 8.

In addition to the series shown in Figure 1, a series was computed using the true value of J, estimated from a much larger sample n = 1024 from the underlying distribution, and this series was indistinguishable from the series computed using $\widehat{J}$ for every n, thus it was not graphed. This validates Takeuchi’s approach of approximating J with $\widehat{J}$ in this instance.

As expected, the difference in $\mu$ between ICE and MLE is not statistically significant (at three standard deviations) for any n, but the ICE computed value of $\sigma$ (shown in Figure 2) is considerably larger than the MLE estimate, especially for small values of n. This explains the greatly reduced KL-divergence noted in Figure 1.

Note that the difference in estimated $\sigma$ is always statistically significant when compared to the MLE value. This is because both MLE and ICE are fit on the same data, so ICE would always have a larger $\sigma$ than MLE regardless of the actual data chosen from the distribution f. This is the cause of the large z-scores shown, always exceeding 200. We know from elementary statistics that correlation between the mean and std. deviation causes the MLE estimate of $\widehat{\sigma }$ to be systematically low by a factor of $\frac{n-1}{n}$. Indeed, the ICE estimate of ${\sigma }^{\ast }$ is closely tracking ${\sigma }_0$ whereas $\widehat{\sigma }$ is closely tracking $\frac{{\sigma }_0(n-1)}{n}$ as expected. This is one example where reducing generalization error also reduces parameter bias as a side effect.

4.2. Friedman’s Test Case

We now extend the example from Section 4.1 to the case where $\mu$ is no longer constants. For this example, we chose a standard regression test set, which is nonlinear in the features, based on Section 4.3 of [24]:

$$y_i={\mu }_{\mathit{\theta }}(x_i)+{\epsilon }_i,\epsilon \sim N(0,{\sigma }^2),$$

(12)

where the Friedman model is

$${\mu }_{\mathit{\theta }}(x_i)={\mathit{\theta }}_{0}sin(\pi x_{(i,0)}{x}_{(i,1)})+{\mathit{\theta }}_1{(x_{(i,2)}-{\mathit{\theta }}_2)}^2+{\mathit{\theta }}_3{x}_{(i,3)}+{\mathit{\theta }}_4{x}_{(i,4)}.$$

(13)

The random features, $X_j$, are i.i.d. uniform random and the parameter values are fixed. The true parameter set, ${\mathit{\theta }}_0=(10.0,20.0,0.5,10.0,5.0,1.0)$, reserves the last parameter ($1.0$) for the value of $\sigma$.

Note that here $\sigma$ must be treated as an unknown parameter. To do otherwise implies that the modeler knows the amount of noise expected in the data. In the case of a known noise term, overfitting is impossible since overfitting arises when a model reduces the projected noise below its actual value, which can never arise when the noise level is known.

The model probability density $g(x,y|\mathit{\theta })$ of y is given by

$$g(x_i,y_i|\mathit{\theta })=\frac{1}{\sqrt{2\pi {\sigma }^2}}{e}^{-\frac{{({\mu }_i-y_i)}^2}{2{\sigma }^2}}.$$

(14)

Recall that in Section 4.1, the value of $\mu$ was considered to be a constant. This example is a natural extension of Section 4.1, and was chosen due to the well-explored difficulty of Friedman’s problem.

We simulate 500 batches of equally sized training sets of length $n\in \{16,32,64,128,256,$$512,1024\}$. The test set is always of length 1024 to ensure accuracy for the smaller values of n. The starting point of the optimization is generated by adding a random perturbation, $\delta \mathit{\theta }\sim N(0,0.1)$, to each parameter. As before, the KL-divergence is computed between the distribution represented by the parameters and the true distribution, and these values are compared between estimation methods.

For each test sample, the KL divergence is computed using numerical integration with a $dy$ increment of $0.01$ over the interval containing $\mu \pm 10\sigma$ for both the true and model distributions. The computed probabilities are verified to numerically sum to unity within an error of $\pm {10}^{-3}$.

In each simulation, $\mathit{\theta }$ is computed using MLE, MLE with $L_2$ regularization, and ICE. The $\lambda$ parameter for $L_2$ regularization is computed using 4-way cross validation on each batch of the training data.

As shown in Figure 3 and

Table 1. Comparison of the average KL divergence across 500 replications for several model estimators given a fitting set size of n. For estimators other than $\widehat{\mathit{\theta }}$, the values in parentheses denotes the t-statistic of the difference between this estimator and $\widehat{\mathit{\theta }}$, with negative values indicating that the listed estimator has a lower KL divergence.

n	${\mathit{\rho }}_{\mathit{KL}}(\mathit{f},{\mathit{g}}_{\widehat{\mathit{\theta }}})$	${\mathit{\rho }}_{\mathit{KL}}(\mathit{f},{\mathit{g}}_{{\mathit{\theta }}_{\mathit{L}2}})$	${\mathit{\rho }}_{\mathit{KL}}(\mathit{f},{\mathit{g}}_{{\mathit{\theta }}^{\ast }})$
16	$6.19\times {10}^{-1}$	$8.442\times {10}^{-1}$ (8.40)	$3.02\times {10}^{-1}$ (−13.54)
32	$1.74\times {10}^{-1}$	$1.74\times {10}^{-1}$ (0.16)	$1.37\times {10}^{-1}$ (−11.11)
64	$6.85\times {10}^{-2}$	$6.85\times {10}^{-2}$ (0.0)	$5.81\times {10}^{-2}$ (−12.65)
128	$3.82\times {10}^{-2}$	$3.82\times {10}^{-2}$ (0.0)	$3.53\times {10}^{-2}$ (−11.13)
256	$2.19\times {10}^{-2}$	$2.19\times {10}^{-2}$ (0.0)	$2.12\times {10}^{-2}$ (−7.84)
512	$1.53\times {10}^{-2}$	$1.53\times {10}^{-2}$ (0.0)	$1.51\times {10}^{-2}$ (−5.66)
1024	$1.25\times {10}^{-2}$	$1.25\times {10}^{-2}$ (0.0)	$1.24\times {10}^{-2}$ (−5.03)

, $L_2$ is not effective for any value of n, and is is completely inactivated for $n>32$. Where regularization is used (i.e., $\lambda \ne 0$), it generally underperforms MLE. ICE is effective across the entire data range, outperforming MLE for every n, and always by a statistically significant margin of at least 5 sigma.

4.2.1. Impact of $\widehat{J}$ Approximation

It was noted previously that Takeuchi used $\widehat{J}$ (and likewise, $\widehat{I}$) in place of the true value of J, and we do so here as well. Though there is no realistic way to avoid this approximation in the real world, and the optimized approach discussed in Section 5 has an entirely different set of approximations, the impact of this approximation will be briefly characterized here.

In

Table 2. Comparison of the average KL divergence across 500 replications for ICE estimators with and without approximation of J given a fitting set size of n. For estimators other than $\widehat{\mathit{\theta }}$, the values in parentheses denotes the t-statistic of the difference between this estimator and $\widehat{\mathit{\theta }}$, with negative values indicating that the listed estimator has a lower KL divergence.

n	MLE	ICE ($\widehat{\mathit{J}}$ )	ICE (J )
16	$6.19\times {10}^{-1}$	$3.02\times {10}^{-1}$ (−13.54)	$6.21\times {10}^{-1}$ (0.05)
32	$1.74\times {10}^{-1}$	$1.37\times {10}^{-1}$ (−11.11)	$1.51\times {10}^{-1}$ (−3.74)
64	$6.85\times {10}^{-2}$	$5.81\times {10}^{-2}$ (−12.65)	$4.21\times {10}^{-2}$ (−17.29)
128	$3.82\times {10}^{-2}$	$3.53\times {10}^{-2}$ (−11.13)	$2.99\times {10}^{-2}$ (−11.82)
256	$2.19\times {10}^{-2}$	$2.12\times {10}^{-2}$ (−7.84)	$1.99\times {10}^{-2}$ (−4.19)
512	$1.53\times {10}^{-2}$	$1.51\times {10}^{-2}$ (−5.66)	$1.48\times {10}^{-2}$ (−1.62)
1024	$1.25\times {10}^{-2}$	$1.24\times {10}^{-2}$ (−5.03)	$1.24\times {10}^{-2}$ (−0.72)

, we revisit Table 1, but now drop the $L_2$ regualarization column, and add a new column where the ICE objective is allowed to use a much better approximated value of J, in this case approximated from 1024 independently drawn samples regardless of n.

As can be seen from Table 2, using the true value of J is at most marginally helpful. In fact, for most values of n it displays slightly better average results, but slightly higher std. deviation of those results, and thus reduced T-statistics. Thus, we conclude that the Takeuchi’s approximation, replacing J with $\widehat{J}$ is reasonable. The same conclusion was reached in Section 4.1, see the remark there. We note also that the ICE estimator using $\widehat{J}$ exhibits substantially better performance for very low sample sizes, but further investigation of this phenomenon is beyond the scope of the current paper.

In

Table 3. Mean matrix norms of J, its approximations, and differences from these approximations across 500 replications.

n	$\parallel \mathit{J}\parallel$	$\parallel \widehat{\mathit{J}}\parallel$	$\parallel \mathit{D}\parallel$	$\parallel \widehat{\mathit{J}}-\mathit{J}\parallel$	$\parallel \mathit{D}-\mathit{J}\parallel$	$\frac{1}{\mathit{n}}\mathit{tr}({\mathit{IJ}}^{-1})$	$\frac{1}{\mathit{n}}\mathit{tr}(\widehat{\mathit{I}}{\widehat{\mathit{J}}}^{-1})$	$\frac{1}{\mathit{n}}\mathit{tr}(\widehat{\mathit{I}}{\mathit{D}}^{-1})$
16	1423.29	925,489.26	1738.30	926,165.93	2905.64	0.1452	0.0027	6.9484
32	875.03	89,414.11	96.49	89,786.19	793.43	0.0521	0.0256	0.0291
64	214.05	4820.55	89.55	4646.51	125.02	0.0202	0.0160	0.0162
128	200.86	200.68	85.00	27.65	115.86	0.0097	0.0086	0.0086
256	194.36	191.81	82.70	19.12	111.67	0.0048	0.0045	0.0045
512	191.20	190.10	82.76	14.18	108.44	0.0023	0.0023	0.0023
1024	188.91	187.39	81.89	11.58	107.02	0.0012	0.0012	0.0012

, we show the average matrix norms of J, $\widehat{J}$, and also of the diagonal of $\widehat{J}$, referred to as the matrix D. The matrix D will be examined further in Section 5, and is included here for completeness. We also show the norms of several matrix differences.

We note that the ICE objective values themselves exhibit much lower variation than the matrix norms show in Table 3. In particular though the matrix D is not actually converging to J as n increases, we see from the correction term it generates that this difference does not appear to have a material impact for larger n. We thus conclude that the major eigenvectors of $(\widehat{J}-J)$ and $(D-J)$ are very nearly orthogonal to the gradient vectors used to construct $\widehat{I}$ for large n.

It is not clear from examining the trace terms in Table 3 that D is a worse approximation of J than $\widehat{J}$ is, even for small n where the impact of the ICE approach is most significant. A more complete investigation of the spectrum of these matrices is beyond the scope of the present work.

4.3. Multivariate Logistic Regression

The previous experiment is based on a well-known test case. In this second experiment, we assess the general performance of ICE under (i) varying dimensionality of the true data distribution, (ii) increasing misspecification, and (iii) increasing training set sizes. To achieve this goal, we generate a more exhaustive set of data from a more complex data generation process.

4.3.1. Data Generation Process

The synthetic data are designed to exhibit a number of characteristics needed to broadly evaluate the efficacy of ICE. First, the regressors should be sufficiently correlated so as to ensure that model selection is representative of typical datasets. However, we avoid multi-collinearity by ensuring the smallest eigenvalue is above a certain threshold. We additionally control the condition number of the covariance matrix $\mathrm{\Sigma }$ by randomly generating a symmetric positive definite covariance matrix $\mathrm{\Sigma }\in {\mathbb{R}}^p$ using the eigen-decomposition

$$\mathrm{\Sigma }=UD{U}^T,$$

(15)

where U is an orthogonal random matrix with elements $U_{ij}\sim N(0,1)$ and D is diagonal matrix of positive eigenvalues. The eigenvalues are uniformly distributed over the interval $[a,b]$ so that the condition number of $\mathrm{\Sigma }$ is $b/a$ and the eigenvalues are kept distinct. Here, a is chosen to be $1\times {10}^{-4}$ and b is chosen to be 0.1.

Using a Cholesky decomposition $\mathrm{\Sigma }=\mathrm{\Gamma }{\mathrm{\Gamma }}^T$ and the random mean vector $\mu \sim N(0,1)$, we generate correlated gaussian vectors of dimension p with the properties

$$X_i=\mu +{\mathrm{\Gamma }}_{ij}{Z}_j,Z_j\sim N(0,1),\forall j\in 1,\cdots ,p.$$

(16)

The data $({\mathit{x}}_n,y_n)$ are generated under a logistic regression

$$p(y=1|\mathit{x},{\mathit{\theta }}_0)=f(\mathit{x}|{\mathit{\theta }}_0)=\frac{1}{1+e^{-\mathit{x}{\mathit{\theta }}_0}}.$$

(17)

A key challenge in assessing the efficacy of bias reduction is to avoid generating excessively low entropy distributions. In such cases, bias reduction will have marginal effect as the parameters are all nearly zero. To avoid such scenarios, the intercept parameter of the true model is adjusted a-posterior until the following conditions are met:

• $c<{\mathbb{E}}_Z[p(Y=1|\mathit{x},{\mathit{\theta }}_0)]<d$
• $-\mathcal{L}({\mathit{\theta }}_0)>ϵ$

where $c=0.35$, $d=0.65$, and $ϵ=0.2$. If these conditions can not be met, then the replication is discarded.

4.3.2. Model Performance Comparison

As in prior sections, KL divergence is computed between the estimated model and the true model for each of the estimation methods. The T-statistics of the difference with the corresponding MLE KL divergence are computed, with negative T-statistics showing that an approach is performing better than the MLE approach. For $L_2$ regularization in this section, the value of $\lambda$ is computed via cross-validation, using two folds, on the provided fitting set.

Table 4. Comparison of the KL divergence for the different estimation approaches applied to mis-specified data. The values in parentheses denote the t-statistic relative to MLE. For $p=\{5,10,20\}$ there are $m=\{2,4,8\}$ non-explanatory variables added.

p	n	${\mathit{\rho }}_{\mathit{KL}}(\mathit{f},{\mathit{g}}_{\widehat{\mathit{\theta }}})$	${\mathit{\rho }}_{\mathit{KL}}(\mathit{f},{\mathit{g}}_{{\mathit{\theta }}_{{\mathit{L}}_2}})$	${\mathit{\rho }}_{\mathit{KL}}(\mathit{f},{\mathit{g}}_{{\mathit{\theta }}^{\ast }})$
5	500	$4.79\times {10}^{-3}$	$3.01\times {10}^{-3}$ (−13.43)	$4.56\times {10}^{-3}$ (−13.53)
5	1000	$2.64\times {10}^{-3}$	$1.76\times {10}^{-3}$ (−10.94)	$2.57\times {10}^{-3}$ (−12.80)
5	2000	$1.29\times {10}^{-3}$	$1.09\times {10}^{-3}$ (−6.15)	$1.27\times {10}^{-3}$ (−7.69)
5	5000	$5.09\times {10}^{-4}$	$4.60\times {10}^{-4}$ (−4.69)	$5.07\times {10}^{-4}$ (−6.19)
10	500	$9.79\times {10}^{-3}$	$9.85\times {10}^{-3}$ (0.16)	$9.18\times {10}^{-3}$ (−6.27)
10	1000	$5.05\times {10}^{-3}$	$5.13\times {10}^{-3}$ (0.51)	$4.83\times {10}^{-3}$ (−4.90)
10	2000	$2.50\times {10}^{-3}$	$3.05\times {10}^{-3}$ (5.99)	$2.56\times {10}^{-3}$ (1.70)
10	5000	$1.06\times {10}^{-3}$	$1.49\times {10}^{-3}$ (7.72)	$1.04\times {10}^{-3}$ (−0.86)
20	500	$2.18\times {10}^{-2}$	$2.16\times {10}^{-2}$(−0.29)	$1.95\times {10}^{-2}$ (−8.71)
20	1000	$1.13\times {10}^{-2}$	$1.24\times {10}^{-2}$ (3.79)	$1.10\times {10}^{-2}$ (−1.95)
20	2000	$6.86\times {10}^{-3}$	$7.52\times {10}^{-3}$ (4.47)	$6.72\times {10}^{-3}$ (−1.67)
20	5000	$3.57\times {10}^{-3}$	$4.24\times {10}^{-3}$ (6.56)	$3.59\times {10}^{-3}$ (0.45)

compares the KL divergences ${\rho }_{KL}$ from the true distribution to the model distributions produced using various estimation approaches applied to misspecified data. Here m denotes the number of regressors that are not predictive, i.e., ${\mathit{\theta }}_0$ contains m zeros. The experiment is replicated 300 times using the data generation process described above and the test set is fixed at 100,000 observations.

We observe that the t-statistic for ${\mathit{\theta }}^{\ast }$ is most significant for relatively small sample sizes, particularly $n=500$. For these small sizes, the improvement over MLE is greater, though noisier. There is uniform decay in improvement over $\widehat{\mathit{\theta }}$ as n grows, until for $p=10$ and $p=20$ the largest sizes are no longer statistically significant. This is expected, as both the MLE and ICE estimates are converging towards the true value of ${\mathit{\theta }}_0$, and for large enough sample sizes the ICE correction would be dominated by numerical error, particularly the ill conditioning of J.

The $L_2$ estimate improves for small values of p, but then becomes progressively worse for large values of p. We observe that for dimensionality above $p=5$, the $L_2$ regularization described here is no longer effective in reducing the KL-divergence. For low values of p the value of $\mathit{\theta }\mathit{x}$ has comparatively low variance, and thus the logistic function is reasonably locally approximated as linear. For higher p this approximation is less realistic and the performance of $L_2$ regularization degrades.

For the ICE estimates, larger values of p show fluctuations that are often not statistically significant. It is apparent that larger p is increasing the variance of the ICE divergences, probably due to numerical errors and ill conditioning. Larger values of n reduce the absolute size of the divergence improvement whereas larger values of p seem to increase it.

Note that though the t-statistics are degrading for large n, the absolute magnitude of the differences is asymptotically small. For these sizes, the results are insignificant, but more importantly, immaterial.

4.3.3. Convergence Analysis for Large n

For 10 randomly chosen example problems, under which the model coefficients are now fixed, the convergence behavior for large n, the training set size, is explored. Note that the test set remains fixed at 100,000 observations for each problem.

Table 5. Comparison of the KL divergence under the MLE $\widehat{\mathit{\theta }}$, $L_2$ regularization and ICE regularization ${\mathit{\theta }}_{ICE}^{\ast }$ against a large sample size for the case when $p=10$ and $m=4$.

n	$\mathcal{L}({\mathit{\theta }}_0)$	${\mathit{\rho }}_{\mathit{KL}}(\mathit{f},{\mathit{g}}_{\widehat{\mathit{\theta }}})$	${\mathit{\rho }}_{\mathit{KL}}(\mathit{f},{\mathit{g}}_{{\mathit{\theta }}_{{\mathit{L}}_2}})$	${\mathit{\rho }}_{\mathit{KL}}(\mathit{f},{\mathit{g}}_{{\mathit{\theta }}^{\ast }})$
500	$0.5439$	$9.28\times {10}^{-3}$	$7.92\times {10}^{-3}$	$7.74\times {10}^{-3}$
1000	$0.5439$	$5.50\times {10}^{-3}$	$5.86\times {10}^{-3}$	$4.81\times {10}^{-3}$
2000	$0.5439$	$2.65\times {10}^{-3}$	$3.67\times {10}^{-3}$	$2.65\times {10}^{-3}$
5000	$0.5439$	$1.85\times {10}^{-3}$	$2.72\times {10}^{-3}$	$1.35\times {10}^{-3}$
10,000	$0.5439$	$5.75\times {10}^{-4}$	$1.57\times {10}^{-3}$	$9.32\times {10}^{-4}$
20,000	$0.5439$	$5.84\times {10}^{-4}$	$8.10\times {10}^{-4}$	$6.11\times {10}^{-4}$
50,000	$0.5439$	$3.83\times {10}^{-4}$	$4.09\times {10}^{-4}$	$3.64\times {10}^{-4}$
100,000	$0.5439$	$1.67\times {10}^{-4}$	$1.15\times {10}^{-3}$	$1.86\times {10}^{-4}$

compares the KL divergence (averaged over all 10 problems) under MLE ($\widehat{\mathit{\theta }}$), $L_2$ regularization, and ICE for progressively larger sample sizes. The divergences ${\rho }_{KL}(f,g_{\widehat{\mathit{\theta }}})$ and ${\rho }_{KL}(f,g_{{\mathit{\theta }}_{ICE}^{\ast }})$ converge to zero as $n\to \infty$, as does ${\rho }_{KL}(f,g_{{\mathit{\theta }}_{L_2}^{\ast }})$.

Generally the ${\mathit{\theta }}_{ICE}^{\ast }$ estimates are seen to converge slightly faster than the $\widehat{\mathit{\theta }}$ estimates. The regularization in ${\mathit{\theta }}_{L_2}^{\ast }$ is observed to be beneficial for very small sample sizes, but then becomes marginally detrimental for large n.

5. Optimized Computation Results

For any model satisfying White’s Regularity Criteria, it is known that the matrix J is positive definite near the MLE optimum $\widehat{\mathit{\theta }}$. This implies that J is diagonally dominated, and indeed considering just its diagonal elements D, it is known that $tr(I{D}^{-1})>0$. Indeed $tr(I{D}^{-1})$ differs strongly from $tr(I{J}^{-1})$ most strongly for models with strong regressor interactions. Therefore, using finite difference gradients, consider the following approximations for the ICE objective function:

• ${\mathit{\theta }}^{\ast }$: J is computed directly, ICE is implemented as written.
• ${\mathit{\theta }}_2^{\ast }$: J is taken to be constant w.r.t. $\mathit{\theta }$: ($J_{\mathit{\theta }}=J_{\widehat{\mathit{\theta }}}$).
• ${\mathit{\theta }}_3^{\ast }$: J is taken to be diagonal: ($J=D$).
• ${\mathit{\theta }}_4^{\ast }$: J is taken to be the identity: ($J=\mathcal{I}$).

Clearly, we expect that ${\mathit{\theta }}_4^{\ast }$ above is the least accurate approximation, and items ${\mathit{\theta }}_2^{\ast }$ and ${\mathit{\theta }}_3^{\ast }$ have varying levels of accuracy depending on the problem at hand. The cost comparison of these approaches is shown in

Table 6. The asymptotic computational cost (per iteration) of various proposed approximations as a function of parameter count p. Cost is amortized when $J_{\mathit{\theta }}=J_{\widehat{\mathit{\theta }}}$ assuming that $n\approx p$. Note that a typical model will cost $O(p)$ in time and space for both the objective function and its gradients.

Approximation	Objective Cost (Space)	Objective Cost (Time)	Gradient Cost (Space)	Gradient Cost (Time)
Direct Computation	$O(p^2)$	$O(p^3)$	$O(p^2)$	$O(p^4)$
$J_{\mathit{\theta }}=J_{\widehat{\mathit{\theta }}}$	$O(p^2)$	$O(p^2)$	$O(p^2)$	$O(p^3)$
$J=D$	$O(p)$	$O(p)$	$O(p)$	$O(p^2)$
$J=\mathcal{I}$	$O(p)$	$O(p)$	$O(p)$	$O(p^2)$

.

Remark 9.

When computing gradients for use in a solver, often approximation error will have only a marginal impact on the final result, though it may increase the number of iterations needed for convergence. Broyden’s method [25] is a typical example of this approach in action. Efficient approximations of $[{\partial }_{\mathit{\theta }}\widehat{J}]$ might similarly have only a minor effect on accuracy and iteration count. The construction of approximate analytical derivatives is beyond the scope of the present work.

These approximations were computed and compared for the Friedman (see Section 4.2) problem, and the results are shown in Table 6 below.

From

Table 7. Comparison of the average KL divergence across 200 replications for MLE and several variants of ICE given a fitting set size of n. For estimators other than $\widehat{\mathit{\theta }}$, the values in parentheses denotes the t-statistic of the difference between this estimator and $\widehat{\mathit{\theta }}$, with negative values indicating that the listed estimator has a lower KL divergence.

n	${\mathit{\rho }}_{\mathit{KL}}(\mathit{f},{\mathit{g}}_{\widehat{\mathit{\theta }}})$	${\mathit{\rho }}_{\mathit{KL}}(\mathit{f},{\mathit{g}}_{{\mathit{\theta }}^{\ast }})$	${\mathit{\rho }}_{\mathit{KL}}(\mathit{f},{\mathit{g}}_{{\mathit{\theta }}_2^{\ast }})$	${\mathit{\rho }}_{\mathit{KL}}(\mathit{f},{\mathit{g}}_{{\mathit{\theta }}_3^{\ast }})$	${\mathit{\rho }}_{\mathit{KL}}(\mathit{f},{\mathit{g}}_{{\mathit{\theta }}_4^{\ast }})$
8	$1.22\times {10}^{+1}$	$4.55\times {10}^{+0} ($−$4.89)$	$5.70\times {10}^{+0} (-5.26)$	$3.83\times {10}^{+0} (-5.22)$	$1.28\times {10}^{+0} (-4.67)$
16	$6.68\times {10}^{-1}$	$2.99\times {10}^{-1} (-8.13)$	$3.53\times {10}^{-1} (-10.56)$	$3.36\times {10}^{-1} (-8.30)$	$5.47\times {10}^{-1} (-2.11)$
32	$1.45\times {10}^{-1}$	$1.14\times {10}^{-1} (-6.90)$	$1.04\times {10}^{-1} (-8.18)$	$1.08\times {10}^{-1} (-10.16)$	$3.63\times {10}^{-1} (19.60)$
64	$5.93\times {10}^{-2}$	$4.80\times {10}^{-2} (-10.42)$	$4.70\times {10}^{-2} (-6.95)$	$4.81\times {10}^{-2} (-9.81)$	$2.38\times {10}^{-1} (39.08)$
128	$2.48\times {10}^{-2}$	$2.24\times {10}^{-2} (-6.38)$	$2.33\times {10}^{-2} (-2.26)$	$2.26\times {10}^{-2} (-6.00)$	$1.62\times {10}^{-1} (61.85)$
256	$1.21\times {10}^{-2}$	$1.16\times {10}^{-2} (-4.28)$	$1.20\times {10}^{-2} (-0.68)$	$1.16\times {10}^{-2} (-4.11)$	$1.01\times {10}^{-1} (68.82)$
512	$6.26\times {10}^{-3}$	$6.10\times {10}^{-3} (-2.41)$	$6.16\times {10}^{-3} (-0.84)$	$6.10\times {10}^{-3} (-2.39)$	$5.42\times {10}^{-2} (63.84)$
1024	$3.05\times {10}^{-3}$	$3.00\times {10}^{-3} (-2.66)$	$3.04\times {10}^{-3} (-0.37)$	$2.99\times {10}^{-3} (-2.73)$	$2.61\times {10}^{-2} (59.78)$

, it is apparent that approach ${\mathit{\theta }}_4^{\ast }$, taking $J=\mathcal{I}$ is not effective. This is not surprising as the actual J matrix has dramatic differences in scale between regressors. Approximation ${\mathit{\theta }}_3^{\ast }$, taking $J=D$ is accurate enough that it cannot be statistically distinguished from the direct computation of ICE by the test above. Approximation ${\mathit{\theta }}_2^{\ast }$ tends to underperform approximation (3).

Therefore, we propose taking $J=D$ as a more numerically stable approximation of the ICE objective.

6. Conclusions

Takeuchi [16] is believed to be the first to have proposed using an objective function similar to ICE in order to reduce generalization error, though it was applied via model selection. Firth [9] introduced a similar term to reduce parameter bias in model fitting, as opposed to model selection, though he derived it only for exponential model families and did not consider its effect on generalization error. It is not known why this approach did not find widespread use, but one may infer that the $O(p^4)$ computational cost and instability was enough to keep it from wider adoption.

In this paper, we reintroduce the objective function of [16] and provide a more general proof of its widespread applicability. We then show that efficient implementations costing only $O(p)$ are possible. Under finite sample sizes, this bias correction term is shown experimentally in several models to lead to significant reduction in bias compared to maximum likelihood estimation with and without $L_2$ regularization. ICE offers many advantages over $L_2$ penalized maximum likelihood estimation: (i) it’s suitable for most nonlinear models, (ii) it’s provably asymptotically convergent; and (iii) does not rely on any parameters which would need to be provided by the operator or deduced through cross-validation.