On the complexity of logistic regression models

TL;DR

This paper explores the complexity of logistic regression models, revealing how input distribution and correlations influence model complexity, and proposes a new model selection criterion based on input entropy.

Contribution

It introduces a novel complexity measure for logistic models considering input distribution and correlations, and develops an entropy-based model selection method.

Findings

Input correlations reduce model complexity.

Entropy-based criterion outperforms standard methods in various conditions.

Complexity contribution from input alphabet size is minimal.

Abstract

We investigate the complexity of logistic regression models which is defined by counting the number of indistinguishable distributions that the model can represent (Balasubramanian, 1997). We find that the complexity of logistic models with binary inputs does not only depend on the number of parameters but also on the distribution of inputs in a non-trivial way which standard treatments of complexity do not address. In particular, we observe that correlations among inputs induce effective dependencies among parameters thus constraining the model and, consequently, reducing its complexity. We derive simple relations for the upper and lower bounds of the complexity. Furthermore, we show analytically that, defining the model parameters on a finite support rather than the entire axis, decreases the complexity in a manner that critically depends on the size of the domain. Based on our findings, we propose a novel model selection criterion which takes into account the entropy of the input distribution. We test our proposal on the problem of selecting the input variables of a logistic regression model in a Bayesian Model Selection framework. In our numerical tests, we find that, while the reconstruction errors of standard model selection approaches (AIC, BIC, $\ell_1$ regularization) strongly depend on the sparsity of the ground truth, the reconstruction error of our method is always close to the minimum in all conditions of sparsity, data size and strength of input correlations. Finally, we observe that, when considering categorical instead of binary inputs, in a simple and mathematically tractable case, the contribution of the alphabet size to the complexity is very small compared to that of parameter space dimension. We further explore the issue by analysing the dataset of the "13 keys to the White House" which is a method for forecasting the outcomes of US presidential elections.

Equations

$$p(y|\mathbf{x},\mathcal{M})=\frac{e^{y(b+\mathbf{w}\cdot \mathbf{x})}}{2\cosh (b+\mathbf{w}\cdot \mathbf{x})},$$

$$p(\hat{y}|\hat{\mathbf{x}},\mathcal{M})=\int d\mathbf{\theta }{e}^{T\ell (\mathbf{\theta })}p(\mathbf{\theta })$$

$$\ell (\mathbf{\theta })=\mathbf{\theta }\cdot \overline{y\mathbf{x}}-\overline{\log (2\cosh (\mathbf{\theta }\cdot \mathbf{x}))},$$

$$J(\mathbf{\theta })=\frac{\sqrt{\det F(\mathbf{\theta })}}{\int d\mathbf{\theta }\sqrt{\det F(\mathbf{\theta })}}$$

$$\log p(\hat{y}|\hat{\mathbf{x}},\mathcal{M})=T\ell ({\mathbf{\theta }}^{\ast })-\frac{n}{2}\log \frac{T}{2\pi }-\log \int d\mathbf{\theta }\sqrt{\det F(\mathbf{\theta })}+O(1/T)$$

$$C=\log \int d\mathbf{\theta }\sqrt{\det F(\mathbf{\theta })}$$

$$F_{i,j}(\mathbf{\theta })=\sum _{\mu }\nu ({\mathbf{x}}^{\mu }){\cosh }^{-2}(\mathbf{\theta }\cdot {\mathbf{x}}^{\mu })x_i^{\mu }{x}_j^{\mu }.$$

$$e^C=\int d\mathbf{\theta }\sqrt{\sum _{\mathbf{\alpha }}{\det }^2(X^{\mathbf{\alpha }})\prod _{i=1}^n\nu ({\mathbf{x}}^{{\alpha }_i}){\cosh }^{-2}(\mathbf{\theta }\cdot {\mathbf{x}}^{{\alpha }_i})}.$$

$$e^C\le {\pi }^n\sum _{\mathbf{\alpha }\in \Omega }\prod _{i=1}^n\sqrt{\nu ({\mathbf{x}}^{{\alpha }_i})},$$

$$e^C\sim \frac{2}{\sqrt{\pi n}}{\left(\frac{8\pi {\Theta }^2{e}^2}{n}\right)}^{n/4}.$$

$$\log r=-\frac{n}{2}-\frac{n}{2}\log \left(\frac{T{H}_n}{n{H}_N}\right)+\log n,$$

$$p(y|\mathbf{x})=\frac{e^{y\theta ({\sum }_{i=1}^n{x}_i-b)}}{2\cosh (\theta ({\sum }_{i=1}^n{x}_i-b))}$$

$$\ell (\theta )=\theta \overline{y(X_n-b)}-\overline{\log (2\cosh (\theta (X_n-b)))},$$

$$e^C\le \pi \sum _{X_n}\sqrt{\nu (X_n)}.$$

$$F_{i,j}(\mathbf{\theta })=\sum _{\mathbf{x}}\nu (\mathbf{x}){\cosh }^{-2}(\mathbf{\theta }\cdot \mathbf{x})x_i{x}_j.$$

$$\begin{array}[]{lll}F_{i,j}(\bm{\theta})&=&\displaystyle\sum_{\bm{x}_{/k}}\nu^{+}(\bm{x}_{/k})\cosh^{-2}(\theta_{k}+\bm{\theta}_{/k}\cdot\bm{x}_{/k})x_{i}x_{j}\,+\\ &+&\displaystyle\sum_{\bm{x}_{/k}}\nu^{-}(\bm{x}_{/k})\cosh^{-2}(\theta_{k}-\bm{\theta}_{/k}\cdot\bm{x}_{/k})x_{i}x_{j},\end{array}$$

$$\begin{array}[]{lll}F_{i,j}(\bm{\theta})&=&\displaystyle\sum_{\bm{x}_{/k}}\nu^{+}(\bm{x}_{/k})\cosh^{-2}(\theta_{k}+\bm{\theta}_{/k}\cdot\bm{x}_{/k})x_{i}x_{j}\,+\\ &+&\displaystyle\sum_{-\bm{x}_{/k}}\nu^{-}(-\bm{x}_{/k})\cosh^{-2}(\theta_{k}+\bm{\theta}_{/k}\cdot\bm{x}_{/k})x_{i}x_{j},\end{array}$$

$$F_{i,j}(\mathbf{\theta })=\sum _{{\mathbf{x}}_{/k}}\left({\nu }^{+}({\mathbf{x}}_{/k})+{\nu }^{-}(-{\mathbf{x}}_{/k})\right){\cosh }^{-2}({\theta }_k+{\mathbf{\theta }}_{/k}\cdot {\mathbf{x}}_{/k})x_i{x}_j.$$

$$\nu (x_k=1;{\mathbf{x}}_{/k})+\nu (x_k=-1;-{\mathbf{x}}_{/k})=\nu ({\mathbf{x}}_{/k})$$

$$F_{i,j}(\mathbf{\theta })=\sum _{{\mathbf{x}}_{/k}}\nu ({\mathbf{x}}_{/k}){\cosh }^{-2}({\theta }_k+{\mathbf{\theta }}_{/k}\cdot {\mathbf{x}}_{/k})x_i{x}_j.$$

$$\begin{array}[]{lll}F_{i,j}&=&2^{-n}\sum_{\hat{\bm{x}}}\left[\cosh^{-2}(\theta_{i}+\theta_{j}+\hat{\bm{\theta}}\cdot\bm{\hat{x}})-\cosh^{-2}(\theta_{i}-\theta_{j}+\hat{\bm{\theta}}\cdot\bm{\hat{x}})+\right.\\ &-&\left.\cosh^{-2}(-\theta_{i}+\theta_{j}+\hat{\bm{\theta}}\cdot\bm{\hat{x}})+\cosh^{-2}(-\theta_{i}-\theta_{j}+\hat{\bm{\theta}}\cdot\bm{\hat{x}})\right],\\ &=&2^{-n}\sum_{\hat{\bm{x}}}\left[\cosh^{-2}(\theta_{i}+\theta_{j}+\hat{\bm{\theta}}\cdot\bm{\hat{x}})-\cosh^{-2}(\theta_{i}-\theta_{j}+\hat{\bm{\theta}}\cdot\bm{\hat{x}})+\right.\\ &-&\left.\cosh^{-2}(\theta_{i}-\theta_{j}-\hat{\bm{\theta}}\cdot\bm{\hat{x}})+\cosh^{-2}(\theta_{i}+\theta_{j}-\hat{\bm{\theta}}\cdot\bm{\hat{x}})\right]\\ &=&2^{-n+1}\sum_{\hat{\bm{x}}}\left[\cosh^{-2}(\theta_{i}+\theta_{j}+\hat{\bm{\theta}}\cdot\bm{\hat{x}})-\cosh^{-2}(\theta_{i}-\theta_{j}+\hat{\bm{\theta}}\cdot\bm{\hat{x}})\right],\end{array}$$

$$\begin{array}[]{lll}F_{i,j}(\bm{\theta})&=&\frac{1}{2}\int_{-\infty}^{\infty}dz\,p(z)\left[\cosh^{-2}(\theta_{i}+\theta_{j}+z)-\cosh^{-2}(\theta_{i}-\theta_{j}+z)\right]=\\ &=&\frac{1}{2}\int_{-\infty}^{\infty}dz\,p(z-\theta_{i})\left[\cosh^{-2}(z+\theta_{j})-\cosh^{-2}(z-\theta_{j})\right]\end{array}$$

$$\lambda (r)={\int }_{-\infty }^{\infty }dz\frac{e^{-z^2/2r^2}}{\sqrt{2\pi r^2}}{\cosh }^{-2}(z).$$

$$e^C=\frac{2{\pi }^{n/2}}{\Gamma (n/2)}{\int }_0^{\sqrt{n}\Theta }dr{r}^{n-1}{\lambda }^{n/2}(r),$$

$$e^C\sim \frac{2}{\sqrt{\pi n}}{\left(\frac{8\pi {\Theta }^2{e}^2}{n}\right)}^{n/4}.$$

Full text

On the complexity of logistic regression models

**Nicola Bulso1, Matteo Marsili2, Yasser Roudi1

1** The Kavli Institute for Systems Neuroscience and Centre for Neural Computation, NTNU, Trondheim, Norway

2 The Abdus Salam International Centre for Theoretical Physics (ICTP), Trieste, Italy

Keywords: model complexity, logistic regression, Bayesian Model Selection, Minimum Description Length, Jeffreys prior

Abstract

We investigate the complexity of logistic regression models which is defined by counting the number of indistinguishable distributions that the model can represent (Balasubramanian,, 1997). We find that the complexity of logistic models with binary inputs does not only depend on the number of parameters but also on the distribution of inputs in a non-trivial way which standard treatments of complexity do not address. In particular, we observe that correlations among inputs induce effective dependencies among parameters thus constraining the model and, consequently, reducing its complexity. We derive simple relations for the upper and lower bounds of the complexity. Furthermore, we show analytically that, defining the model parameters on a finite support rather than the entire axis, decreases the complexity in a manner that critically depends on the size of the domain. Based on our findings, we propose a novel model selection criterion which takes into account the entropy of the input distribution. We test our proposal on the problem of selecting the input variables of a logistic regression model in a Bayesian Model Selection framework. In our numerical tests, we find that, while the reconstruction errors of standard model selection approaches (AIC, BIC, $\ell_{1}$ regularization) strongly depend on the sparsity of the ground truth, the reconstruction error of our method is always close to the minimum in all conditions of sparsity, data size and strength of input correlations. Finally, we observe that, when considering categorical instead of binary inputs, in a simple and mathematically tractable case, the contribution of the alphabet size to the complexity is very small compared to that of parameter space dimension. We further explore the issue by analysing the dataset of the “13 keys to the White House” which is a method for forecasting the outcomes of US presidential elections.

1 Introduction

In many studies and in different fields of research, a recurring task is to find the set of explanatory variables, among many candidates, which accurately predict the probability of occurrence of an event. This is the case of identifying risk factors for developing a particular disease, finding the set of neurons in a population which are driving the activity of another neuron or discovering biomarkers from genomic data, to cite a few examples (Truett et al.,, 1967; Saeys et al.,, 2007; Hertz et al.,, 2013). When the event we want to predict is binary, such as when estimating the probability of developing a certain disease, a widely employed statistical method is Logistic Regression (Cox,, 1958).

Logistic regression model is also interesting because it is the building block of more sophisticated architectures. For instance, under the so called pseudo-likelihood approximation, the difficult task of learning an Ising Model reduces to that of learning many logistic regression models (Besag,, 1972; Ravikumar et al.,, 2010; Aurell and Ekeberg,, 2012). Another example is that of the Kinetic Ising model (Marre et al.,, 2009; Roudi and Hertz,, 2011; Battistin et al.,, 2015). Both models are popular methods for analysing neural assemblies in computational neuroscience (Schneidman et al.,, 2006; Roudi et al.,, 2015).

Selecting a set of inputs or predictors111In the statistics and machine learning literature, the term “input” is often used interchangeably with “predictor”, “feature”, “covariate” or more classically “independent variable” (Friedman et al.,, 2001). Through the paper we will mainly use the terms “input” or “predictor”. among some candidates corresponds to comparing different models for the output variable. Thus we search for models with a good trade-off between goodness of fit and model complexity: very simple models might not be able to capture the relevant inputs which are driving the variable we are interested in, whereas complex models are more prone to overfitting (Friedman et al.,, 2001). Thus assessing the complexity of a model is an interesting question for model selection (see e.g. Bulso et al.,, 2016).

Classical approaches to Model Selection include Bayesian Information Criteria (BIC) (Schwarz,, 1978) and Akaike Information Criteria (AIC) (Akaike,, 1974). These criteria are derived from different perspectives as approximations of more general quantities in the limit of large samples. However, when the sample is not large enough, these methods might require careful interpretation and application. In this limit, prior information becomes important in the inference process in Bayesian Model Selection. Yet, such information might not be available a priori. In these cases it might be advantageous to use an uninformative prior, such as Jeffreys prior (Jeffreys,, 1946).

As shown in Balasubramanian, (1997), the choice of Jeffreys prior in Bayesian Model Selection corresponds to measuring model complexity from a geometric perspective, namely by counting the number of indistinguishable distributions that the model can represent. In this perspective, complex models are models which are able to describe a wide range of probability distributions. In Myung et al., (2000), the authors further show how this approach for characterizing model complexity relates to the Minimum Description Length (MDL) principle (Rissanen,, 1987, 1996).

Jeffreys prior is a commonly used prior in Bayesian analysis. It exhibits many nice features among which being uninformative and parametrization invariant (George E. P. Box,, 1973) (yet see also Kass and Wasserman,, 1996). For Generalized Linear Models, its theoretical properties has been investigated in Ibrahim and Laud, (1991). In particular, for binomial data it has been shown that the penalization related to the model complexity under Jeffreys prior depends not only on the number of explanatory variables, which determines the dimensionality of the model, but also on their distribution (Chen et al.,, 2008). Interestingly, Jeffreys prior arises in the large sample size limit also from other perspectives in recently proposed approaches to Model Selection (LaMont and Wiggins,, 2017; Mattingly et al.,, 2018). In particular, Mattingly et al., (2018) show, among other things, that an optimal prior obtained by maximizing the mutual information between the parameters and their expected data approaches the Jeffreys prior in the limit of abundant data.

In this paper, we address the question of how complex a given logistic model is. We systematically characterize how the complexity depends on the input distribution and prove how the information on model complexity can be used to devise a successful model selection criterion which improves over other well known criteria (BIC, AIC, $\ell_{1}$ regularization) whose performances instead strongly depend on the sparsity of the data generating model.

The paper is organised as follows. In section 2, we introduce the notation and define the measure of complexity that we use in the paper. Then we proceed to study, analytically and numerically, the complexity of logistic regression models with binary variables by varying the dimensionality of the model and the correlations among inputs. We also investigate the special case of models whose parameters are defined on a finite support (regularized models). In section 3, we exploit our results on the complexity to devise a novel model selection criterion and test it on simulated data. Finally in section 4, we investigate the special problem of degenerate models, namely models where parameters might be constrained to being equal. This is an interesting case in model selection since standard recipes based on counting free parameters are not able to distinguish between, for instance, a logistic model with only one parameter multiplying one input and a model with the same parameter multiplying many different inputs. Furthermore, we will see how this issue relates to the problem of evaluating the complexity of categorical variables. We report analytical and numerical results on the complexity along with an application concerning U.S. presidential elections forecasts.

2 The complexity of logistic regression models with binary inputs

We consider a binary variable $y$ distributed according to a logistic model

[TABLE]

with parameters $\bm{\theta}=(b,\bm{w})$ and $N$ binary inputs $\bm{x}=(x_{1},x_{2},...,x_{N})$. The model is equivalent to a single-layer neural network with $N$ binary input nodes, $\bm{x}$, connected to an output node, $y$, through the weights $\bm{w}$, and bias $b$. Now suppose that we are given $t=1,..,T$ measurements of the inputs, $\hat{\bm{x}}=\{\bm{x}^{(t)}\}_{t=1}^{T}$, and the output, $\hat{y}=\{y^{(t)}\}_{t=1}^{T}$, and that the output data has been generated according to the conditional probability distribution in equation 1 for a given choice of the values of parameters (which can also be zero). We are interested in inferring, from the data, the set of non-zero parameters (non-zero components of the vector $\bm{\theta}$). This is a problem of input selection where, given $N$ candidate inputs, we ask which of them are relevant for the variable we want to explain, i.e. which of them is expected to have any influence (non-zero connection) on the output. We proceed following a Bayesian model selection approach. In this framework each possible combination of non-zero parameters represents a model ${\cal M}$ (a set of candidate inputs) and, given $N$ inputs, there are $2^{N+1}$ possible models 222From now on with“models” we mean non degenerate models, i.e. models whose non-zero parameters are independent of each other, and come back to discuss degenerate models only in the last section.. The posterior probability of a model ${\cal M}$ given the data $p({\cal M}|\hat{y},\hat{\bm{x}})$, can be expressed by Bayes rule as $p({\cal M}|\hat{y},\hat{\bm{x}})\propto p(\hat{y}|\hat{\bm{x}},{\cal M})p({\cal M})$, where $p({\cal M})$ is the prior over models and

[TABLE]

is the probability of measuring $\hat{y}$ given $\hat{\bm{x}}$ and model ${\cal M}$. In the last equation, $p(\bm{\theta})$ is the prior over parameters and $\ell(\bm{\theta})$ is the normalised log-likelihood function

[TABLE]

where $x_{0}^{(t)}=1$ and $\theta_{0}=b$. The overbars in equation 3 stand for averages over all observations, namely $\overline{f(X)}=1/T\sum_{t}f(X^{(t)})$ for any function $f(X)$. Thus, the set of non-zero parameters that has generated the output $\hat{y}$ from the input $\hat{\bm{x}}$ can be inferred by searching for the model that maximizes $p({\cal M}|\hat{y},\hat{\bm{x}})$ among all candidates. When all models are a priori equally likely, this procedure is equivalent of finding the model that maximizes $p(\hat{y}|\hat{\bm{x}},{\cal M})$, which is given by equation 2.

The choice of the prior over parameters represents an important issue when the sample is not very large. In this paper we use Jeffreys prior. Our choice is motivated by the geometrical interpretation of model complexity that originates from it and by its relation with MDL (Balasubramanian,, 1997; Myung et al.,, 2000). Jeffreys prior is defined as follows

[TABLE]

where $F(\bm{\theta})$ is the Fisher Information whose elements are $F_{i,j}(\bm{\theta})=-\text{E}[\partial_{\theta_{i},\theta_{j}}^{2}\ell(\bm{\theta})]$, where $\text{E}[\cdot]$ is the expected value with respect to the model distribution $p(y|\bm{x},\bm{\theta},{\cal{M}})$. Since the elements of the Hessian matrix $H_{i,j}(\bm{\theta})=-\partial_{\theta_{i},\theta_{j}}^{2}\ell(\bm{\theta})$, calculated from equation 3, do not depend on the random variable $y$, it follows that the Fisher Information matrix is simply equal to the Hessian.

With this choice for the prior, the probability $p(\hat{y}|\hat{\bm{x}},{\cal M})$ in equation 2 can be approximated by the following expression (see Balasubramanian,, 1997, for instance):

[TABLE]

where $n\leq N$ is the number of non-zero parameters (non-zero components of the parameter vector $\bm{\theta}$) employed by model $\mathcal{M}$ and $\bm{\theta}^{\ast}$ is the maximum likelihood estimate.

In MDL theory the first three terms in equation 5, all together, represent the stochastic complexity of the string $\hat{y}$ given the model $\mathcal{M}$ and the inputs $\hat{x}$, which is the length of the shortest code that can be obtained by encoding the data with the parametric family ${\cal M}$ (Rissanen,, 1987). The penalty terms in MDL (the second and third term in equation 5), also known as the geometric complexity (Myung et al.,, 2000), are independent of the specific instance of the output data $\hat{y}$ and rather related to intrinsic and geometric properties of the parametric family. Yet it is worth noting that since the Fisher Information in equation 5 depends on the input data $\hat{\bm{x}}$, the complexity also depends on the input distribution, as we will discuss thoroughly in the paper. Finally, the second term in equation 5 corresponds to the well known Bayesian Information Criteria (Schwarz,, 1978).

In the following sections we will investigate the typical properties of the geometric complexity for logistic regression models with binary variables (equation 1). In particular, we will focus on the penalty term

[TABLE]

and we will study its dependence on the model size and input distribution. For logistic regression models, as defined by equation 1, the elements of the Fisher Information matrix can be written explicitly in the form

[TABLE]

where $\nu(\bm{x}^{\mu})$ is the frequency of observing the configuration $\bm{x}^{\mu}$ in the data, with $\mu=1,...,2^{n}$. Beyond the interesting theoretical aspects, the complexity represents a higher order term in the expansion of the posterior 5 and it is therefore relevant for model selection, particularly when the sample size is not very large. From now on we will refer to $C$ in equation 6 as the complexity. Notice that we will also use the term complexity to refer to $e^{C}$ and the difference will be clear from the context.

2.1 A worked example: n = 2

We start by considering simple cases which are analytically solvable. The simplest case is when only one parameter is different from zero. When the only non-zero parameter is the bias, i.e $\theta=b$, then the complexity is simply $e^{C}=\int d\theta\cosh^{-1}(\theta)=\pi$. The same result is retrieved also when the only non-zero parameter is a weight, i.e. $\theta=w$ and the associated input is observed in the configuration $x$ with frequency $\nu(x)$: $e^{C}=\int d\theta\sqrt{F(\theta)}=\int d\theta\sqrt{\sum_{x}\nu(x)\cosh^{-2}(\theta x)}=\int d\theta\cosh^{-1}(\theta)=\pi$.

We consider now a model in which only two weights are different from zero, i.e. $\bm{\theta}=(w_{1},w_{2})$. The Fisher information is now a two by two matrix whose determinant can be easily calculated: $\det F(\bm{\theta})=4\nu(1-\nu)\cosh^{-2}(\theta_{1}+\theta_{2})\cosh^{-2}(\theta_{1}-\theta_{2})$, where $\nu$ is the frequency that the two inputs, $\bm{x}=(x_{1},x_{2})$, have been observed with the same sign. The complexity is therefore given by the integral $\sqrt{4\nu(1-\nu)}\int d\bm{\theta}\cosh^{-1}(\theta_{1}+\theta_{2})\cosh^{-1}(\theta_{1}-\theta_{2})$ which can be solved by the change of variable, $\phi_{1}=\theta_{1}+\theta_{2}$ and $\phi_{2}=\theta_{1}-\theta_{2}$, with the result $e^{C}=\sqrt{\nu(1-\nu)}\pi^{2}$.

As it is clear from the above formula, the complexity varies in the range $0\leq e^{C}\leq\pi^{2}/2$ according to the degree of correlation between the two inputs. The complexity reaches its maximum, i.e. $\pi^{2}/2$, when the two inputs are observed with the same sign half of the times, $\nu=1/2$. It is worth noting that the maximum value for the complexity of the logistic regression model with $n=2$ is smaller than the complexity of a model for two independent outputs, each of which being explained by a logistic regression model with one parameter, which is $\pi^{2}$. In both cases the models employ two parameters. However, in the latter case, the parameters are completely independent of each other and the volume spanned in the space of distributions will be the maximum achievable given the absence of constraints. Following the same argument, for a generic number of parameters $n$, the complexity of a logistic regression model should always be smaller than $\pi^{n}$.

In case of strongly correlated inputs, namely when $\nu$ is close to one or zero, the complexity will reach its minimum. In particular, for a sample of size $T$, the smallest non-zero value for $\nu$ is $\nu=1/T$ and the largest value different from one is $\nu=1-1/T$. In the latter cases, the complexity attains its lower (non trivial) bound which is equal to $e^{C}=\pi^{2}\sqrt{(T-1)/T^{2}}$. Intuitively, correlations among inputs induce effective dependences among parameters, thus constraining the model into a reduced volume in the space of distributions. In other words, the stronger the correlations, the higher the constraints and, consequently, the lower the complexity. This represents a hard case for model selection. In fact, in the case of two strongly correlated inputs, the three models with at least one non-zero weight parameter multiplying the inputs will have approximatively the same maximum likelihood values. It follows that the selected model will strongly depend on the assumptions pertaining the model selection recipe employed.

In the limiting case in which $\nu=0$ or $\nu=1$, we get $e^{C}=0$ and we say that the model is redundant: the model is constrained to a lower dimensional manifold defined by the deterministic dependence between the parameters induced by equal ($x_{1}=x_{2}$) or specular inputs ($x_{1}=-x_{2}$). In this case, the $n$ by $T$ dimensional matrix of input data, $x_{i}^{(t)}$, also called the design matrix, is not full rank. In all cases which are not redundant, i.e. $0<\nu<1$, the matrix is full rank, the complexity is finite and the BIC term is always larger than $C$ in absolute value, consistent with the expansion of equation 5 and the geometrical interpretation of the penalization factors as ratio between volumes (Balasubramanian,, 1997).

In the following we derive and discuss some results on the complexity of logistic regression models as a function of the dimension $n$ of the model and of the shape of the input distribution extending the observations drawn from these simple cases.

2.2 Bounds on geometric complexity

In this section, we will derive and discuss analytical results for the upper and lower bound of the complexity. Additionally, we will support and further investigate the results with numerical simulations. Without loss of generality, we consider models with only weight parameters, namely $\bm{\theta}=\bm{w}$ as defined at the beginning of section 2. In fact, all results derived in this section easily extend also to models with a bias parameter, i.e. $\bm{\theta}=(b,\bm{w})$, given that the bias term acting on the output can be thought as a weight multiplying a constant input.

The Fisher Information defined in equation 7 can be conveniently rewritten as $F(\bm{\theta})=B^{\top}B$ by defining the $2^{n}\times n$ matrix $B$ as $B^{\mu}_{i}=\sqrt{\nu(\mathbf{x}^{\mu})\cosh^{-2}(\bm{\theta}\cdot\mathbf{x}^{\mu})}x^{\mu}_{i}$, where we use subscripts and Roman letter for components indices (those ranging from 1 to $n$) whereas superscript and Greek letters for configurations indices (ranging from 1 to $2^{n}$). Moreover the matrix $B$ can be expressed as $B=DX$, namely as the product of the $2^{n}\times n$ matrix $X$ whose components are given by $x^{\mu}_{i}$ and the diagonal matrix $D$ whose elements are $D^{\mu\nu}=\sqrt{\nu(\mathbf{x}^{\mu})\cosh^{-2}(\bm{\theta}\cdot\mathbf{x}^{\mu})}\delta^{\mu\nu}$. These mathematical manipulations allow us to exploit Cauchy Binet formula and write the complexity in a more compact form. In fact, by Cauchy Binet formula, the determinant of the Fisher Information reads $\det F(\theta)=\det(B^{\top}B)=\sum_{\bm{\alpha}}{\det}^{2}(B^{\bm{\alpha}})$, where $B^{\bm{\alpha}}$ is now a square matrix obtained from $B$ by selecting $n$ configuration indices over the $2^{n}$ possible ones. The selection is labelled by the vector $\bm{\alpha}$ and the sum ranges over all ${2^{n}}\choose{n}$ possible $\bm{\alpha}$ vectors. Finally, since $B=DX$, the determinant of the Fisher Information can be written as $\det F(\theta)=\sum_{\bm{\alpha}}{\det}^{2}(X^{\bm{\alpha}})\prod_{i=1}^{n}\nu(\mathbf{x}^{\alpha_{i}})\cosh^{-2}(\bm{\theta}\cdot\mathbf{x}^{\alpha_{i}})$ and consequently the complexity takes the form

[TABLE]

From the above expression we can now find a mathematical approximation for the upper bound of the complexity. By employing the triangular inequality in equation 8, we move the sum over $\bm{\alpha}$ outside the square root so that $e^{C}\leq\sum_{\bm{\alpha}}\det(X^{\bm{\alpha}})\prod_{i=1}^{n}\sqrt{\nu(\bm{x}^{\alpha_{i}})}I(\mathbf{x}^{\alpha_{i}})$ where $I(\mathbf{x}^{\alpha_{i}})=\int d\bm{\theta}\cosh^{-1}(\bm{\theta}\cdot\mathbf{x}^{\alpha_{i}})$. The integral $I(\mathbf{x}^{\alpha_{i}})$ can be evaluated by making the change of variable $\phi_{i}=\bm{\theta}\cdot\mathbf{x}^{\alpha_{i}}$ and turns out to be equal to $\pi^{n}/\det(X^{\bm{\alpha}})$. Notice that if ${\det}(X^{\bm{\alpha}})=0$, there is no contribution to the summation in expression 8. Therefore the factors $\det(X^{\bm{\alpha}})$ cancel each other and the sum over $\bm{\alpha}$ is restricted to those selections $\bm{\alpha}\in\Omega$ such that $\det(X^{\bm{\alpha}})\neq 0$. The final result is the following simple inequality

[TABLE]

where the equal sign holds when only $n$ linearly independent configurations are present in the input data matrix (design matrix) as, in this case, there is only one possible selection ${\bm{\alpha}}$ in equation 8. Notice that equation 9 depends on both the dimensionality $n$ and the input distribution $\nu(\bm{x})$.

First, we consider the case in which the input distribution is localised on few configurations. In the limiting case in which the distribution is peaked around only one configuration $\xi$, i.e. $\nu(\mathbf{x}^{\mu})=\delta_{\mu\xi}$, where $\delta$ is the Kronecker delta function (most localised distribution), $e^{C}$ will be zero. As we have already observed in section 2.1 for models with $n=2$, in this extreme case the design matrix is not full rank. Consequently the parametrization is redundant: a model with only one parameter would have sufficed to explain the behaviour of the output variable in this particular case. Thus, in general the redundancy is caused by dependencies among the parameters induced by an “extremely localised” distribution of the input configurations. Those dependencies constrain the model to a lower dimensional manifold whose volume is of zero measure. Whenever this happens, the model represents a redundant description of the reality and one should look for a simpler explanation by reducing the number of parameters and/or combining inputs.

Moreover, it is rather intuitive and easy to see from equation 9 that the parametrization is not redundant if there are at least $n$ linearly independent configurations $\bm{x}^{\mu_{i}}$ with $i=1,...,n$ such that $\nu(\bm{x}^{\mu_{i}})\neq 0$ , that is if the design matrix is full rank. In the case of exactly $n$ linearly independent configurations there is only one selection $\bm{\alpha}$ such that $\det X^{\bm{\alpha}}\neq 0$ and equation 9 holds with the equal sign: the result is simply $e^{C}=\pi^{n}\prod_{i=1}^{n}\sqrt{\nu(\mathbf{x}^{\alpha_{i}})}$. In this case, given a sample of size $T$, the value of the complexity can vary between two values, $\pi^{n}\sqrt{(T-n+1)/T^{n}}\leq e^{C}\leq\pi^{n}/n^{n/2}$. The upper bound can be derived by maximizing $e^{C}$ while enforcing the normalisation of the frequencies through a Lagrange multiplier and is achieved when the $n$ linearly independent configurations are observed an equal number of times $\nu(\bm{x}^{\mu_{i}})=1/n$, $\forall i$. The lower bound can be verified numerically for small values of $n$ and $T$ and corresponds to the situation in which all configurations but one are observed only once, i.e. $\nu(\bm{x}^{\mu_{i}})=1/T$, $\forall i\neq j$, and $\nu(\bm{x}^{\mu_{j}})=(T-n+1)/T$. The difference between these two values tends to zero as $T$ approaches $n$ from above, while it grows large as $T\gg n$ and $n\gg 1$. In fact, their ratio scales as $(T/n)^{n/2}$.

Based on the fact that a more localised input distribution would lead to a redundant parametrization (zero value for the complexity) and given that equation 9 and numerical simulations (see Figure 1) suggests that a more spread distribution would generally lead to an equal or larger value for the complexity, we expect that the lower bound presented above, $e^{C}=\pi^{n}\sqrt{(T-n+1)/T^{n}}$, represents the minimum attainable (non zero) value of the complexity. In fact, when $n=2$, it matches the lower bound we found in section 2.1.

Interestingly, in this “most-localized scenario”, the complexity decreases instead of increasing as the dimensionality of the model increases (it goes at most as $e^{C}\sim\pi^{n}/n^{n/2}$). This means that, although adding a new parameter would still increase the dimensionality of the model (larger BIC penalization), it would also result in a more constrained model (decreasing of complexity).

In order to support these results and further investigate the dependence of the complexity on the input distribution, we resort to Monte Carlo simulations for estimating the integral $e^{C}=\int d\bm{\theta}\sqrt{\det F(\bm{\theta})}$ (details of the simulations are described in Supplemental Information, section 1).

As a first step, we assume an input distribution which is uniform and different from zero only in a limited number of configurations, $d$. The parameter $d$ tunes the localisation of the input distribution or, in other words, the strength of input correlations. We evaluated the complexity at varying $d$ over the entire range: from $d=1$ to $d=2^{n}$. The results are shown in Figure 1 for $n=2,3,4,5,6,8$.

As predicted, the complexity is zero if less than $n$ configurations have been observed at least once and, when the input distribution is uniformly distributed on exactly $n$ linearly independent configurations, the complexity jumps to $e^{C}=\pi^{n}/n^{n/2}$ (bottom black dashed line). For $d>n$, the value of the complexity depends on the degree of independence of the $d$ configurations, namely on the volumes spanned by each subset of $n$ configurations. Therefore for $d>n$ the complexity exhibits a high variability while generally increasing with $d$ from $e^{C}=\pi^{n}/n^{n/2}$ toward an upper bound which is achieved when the input distribution is uniform over all possible configurations (top black dashed line). This fact is in line with what has been observed in the case $n=2$ in section 2.1 where a uniform distribution leads to the model with the largest complexity.

Unfortunately, when the distribution is uniform, the inequality 9 provides a trivial upper bound, namely an upper bound larger than expected from simple arguments i.e. $e^{C}\leq\pi^{n}$, as argued in 2.1.

The upper bound of the complexity, obtained by evaluating the integral with Monte Carlo simulations assuming a uniform distribution on the inputs, is shown in Figure 2 a) at varying $n$, the number of inputs. Interestingly, the trend is consistent, within error bars, with the simple function $e^{C}=\pi^{n}/n$ (black dotted line) which has been also reported in Figure 1 (black dashed line) for comparison. As expected, the upper bound is found to be smaller than $\pi^{n}$.

Next, we assume that the input distribution follows an Ising model, namely $\nu(\bm{x})\propto\exp(-\sum_{i<j}\beta J_{i,j}x_{i}x_{j}-\sum_{i}\beta h_{i}x_{i})$, where $J_{i,j}=1/n$ $\forall i\neq j$, $J_{i,i}=0$ and $h_{i}=0.1$, $\forall i$333The small parameter $h_{i}$ has been introduced to break the symmetry between the two configurations with all inputs equal.. The localisation of the distribution is regulated by the parameter $\beta$: very small values of $\beta$ produce an almost uniform distribution whereas for larger values the distribution becomes more and more localised around the configuration with all inputs being one. As before, we used Monte Carlo simulations to estimate the integral for $n=4,5,6$ and $8$ for different values of $\beta$. The results are shown in Figure 2 b) with red lines. As expected, for small $\beta$ the curves approach the upper bound $e^{C}=\pi^{n}/n$ (top black dashed line) and when increasing $\beta$, the complexity decreases until some point around the limit $e^{C}=\pi^{n}/n^{n/2}$ (bottom black dashed line).

As a final remark, in Appendix A we prove that, due to the parity of the hyperbolic cosine function, the complexity reaches the upper bound already when the input distribution is uniform over the configurations of a subset of $n-1$ inputs. This is evident in Figure 1 for the cases $n=2,3$ and $4$: the complexity reached the upper bound already when all possible $2^{n-1}$ patterns relative to a subset of $n-1$ inputs were observed the same number of times. For larger values of $n$ this is not observed because none of the 10 different random paths in the simulation passes through the point in which the distribution is uniform over the $2^{n-1}$ patterns relative to a subset of $n-1$ inputs. In particular, as observed in Appendix A, this fact implies that the upper bound of the complexity for a model with $n$ inputs and without the bias term on the output variable is equal to that of a model with $n-1$ inputs and the bias.

2.3 The emergence of statistical constraints in regularized models

We now consider logistic models whose parameters $\theta_{k}$ are defined on a finite support, $|\theta_{k}|\leq\Theta$, $\forall k$, with $\Theta$ being any finite positive number, and we study the effect of this effective regularization on the model complexity. In particular, we focus on the case of uniformly distributed inputs, $\nu(\bm{x}^{\mu})=2^{-n}$, $\forall\mu$, and compare the results obtained in the last section with un-regularized models.

In Appendix B we derive a fully analytical expression for the complexity of such models in the limit of an indefinitely large input layer. In this regime, the complexity can be approximated by the following equation

[TABLE]

As it is clear from the above expression, if the parametric family is defined in a box with half length $\Theta$, the volume occupied by the parametric family in the space of distribution will start shrinking after a certain critical dimension in contrast with the steady increase observed for un-regularized models, namely $e^{C}=\pi^{n}/n$ . More precisely, beyond the critical dimension, $n_{c}=8\pi\Theta^{2}e^{2}$, increasing the number of inputs will make the model more constrained although still adding further degrees of freedom.

Intuitively this is due to the emergence of statistical constraints, namely the applicability of Central Limit arguments as shown in Appendix B, which induces effective dependencies among parameters. In fact, in the proof of equation 10 given in Appendix B, the assumptions of a large number of uniformly distributed inputs and parameters defined on a finite support allowed us to invoke Central Limit Theorem 444Notice that the requirement of the finite support is the crucial assumption in the derivation and makes the difference with the previous cases of un-regularized models. which implies that, for any choice of the parameter values, their linear combination with the inputs is constrained to follow a Gaussian distribution with zero mean. This is similar to what was observed in the previous section in the case of the most-localized input distribution where the complexity decreases with increasing the dimensionality of the model. In that case, the effective dependencies among parameters were induced by the correlations among inputs while, in this case, by statistical arguments.

Finally, the fact that the complexity starts decreasing after a certain critical dimension in regularized models is expected to be true for any non trivial input distribution given that in this section we proved it to be true for the uniform distribution which corresponds to the upper bound in the un-regularized case.

3 Model selection

In this section, we are going to apply the outcomes of the previous sections to the problem of model selection with logistic regression models and binary variables. As outlined at the beginning of section 2, we are interested in solving the following problem: given $T$ measurements of the input and output variables in a logistic model with binary inputs, what is the model $\mathcal{M}$, i.e. the set of non-zero parameters, that has most likely generated the output data given the inputs? Since each weight parameter multiplies a different input in a logistic model, this is also an input selection problem: which subset of the inputs is relevant for predicting the output? In the following, we will first derive a novel Model Selection criterion based on the posterior expansion and the results on the complexity obtained in section 2. Then we will apply this criterion to investigate the quality of model recovery and compare the performance with those achieved with other standard model selection recipes. As before and without any loss of generality, we will neglect the bias and consider only weight parameters, i.e. $b=0$ in equation 1 and $\bm{\theta}=\bm{w}$, since, if needed, the bias can always be introduced as a weight parameter multiplying an additional constant input.

3.1 A novel model selection criteria

Given $N$ potential predictors (inputs) for the output variable, there are $2^{N}$ possible models each of which trying to explain the output with a different subset of $n\leq N$ inputs. In other words, each model corresponds to a different subset of $n\leq N$ non-zero weight parameters. In a Bayesian Model Selection framework, we compare models according to their posterior probability. Assuming that all these models are a priori equally likely, this coincides with ranking models based on their likelihood $p(\hat{y}|\hat{\bm{x}},\mathcal{M})$. In section 2 we saw that $p(\hat{y}|\hat{\bm{x}},\mathcal{M})$ can be approximated as $p(\hat{y}|\hat{\bm{x}},\mathcal{M})=T\ell(\bm{\theta^{\ast}})+\log r+O(1/T)$ where the first term is the likelihood of the data calculated at the Maximum Likelihood estimates of the parameters and the second term is a penalization factor. The latter is given by $\log r=-\frac{n}{2}\log\left(\frac{T}{2\pi}\right)-C$ and consists of the sum of the BIC and the complexity term respectively. In the last section we found that the value of $C$ depends on both the localisation of the input distribution and the dimensionality of the model; it is upper bounded by $\log(\pi^{n}/n)$; finally, it decreases down to some value around $\log(\pi^{n}/n^{n/2})$ as the localisation of the input distribution increases. In the first case, the penalisation would be $\log r=-\frac{n}{2}\log\left(\frac{\pi T}{2}\right)+\log n$ which is dominated by a linear dependence on n and a logarithmic dependence on T, similar to the case of BIC. In the second case, the penalisation would be $\log r=-\frac{n}{2}\log\left(\frac{\pi T}{2n}\right)$ which is a weaker penalisation with respect to the previous case.

We now introduce the entropy $H_{n}(\bm{x})=-\sum_{\bm{x}}\nu(\bm{x})\log_{2}\nu(\bm{x})$ as a measure of the localisation of the input distribution: the value of the entropy approaches $n$ (the dimension of the model) when the distribution of the inputs approaches the uniform distribution; whereas it shrinks to some value around $\log_{2}n$ when the distribution is localised on approximatively $n$ linearly independent configurations.

Based on the results of the previous section and on the aforementioned measure of localisation, we introduce the following heuristic for calculating the penalization of a model with $n$ inputs whose distribution has entropy $H_{n}$:

[TABLE]

where $H_{N}$ is the entropy of the full design matrix with all $N$ predictors taken into account. We further motivate this formula in Supplemental Information (section 2), while, in the following, we provide the intuition behind it and discuss its properties and performance on model selection tasks.

Equation 11 allows us to approximately interpolate between the expected trends as derived in the limiting cases of a uniform and localised distribution. In fact, when $H_{n}\simeq n$, the largest contribution to the penalization exhibits the expected BIC-like trend, $\log r\sim-\frac{n}{2}\log\left(T\right)$, which represents the leading contribution in this regime given that in order to have $H_{n}\simeq n$, a sample size larger than the dimensionality of the configuration space is needed, i.e. $T\geq 2^{n}$, and given that $H_{N}\leq\log_{2}T$. Alternatively, one might think about this as the type of penalization which affects mostly models with $n\ll\log_{2}T$ (or sparsity $s=1-n/N\gg 1-\log_{2}T/N$).

On the other hand, when increasing the localisation of the distribution, the entropy growth with $n$ would be drastically reduced. This is typically the case for models with a large number of inputs, namely with $\log_{2}T\ll n\leq N$, where $H_{n}\simeq H_{N}$ (see also Figure 3 in Supplemental Information). In this case, the factor $T$ in the argument of the logarithm in equation 11 is counterbalanced by a large value of $n$ leading to a weaker penalization, similarly to what was obtained when substituting the lower bound of the complexity into the posterior expansion. Moreover, for very dense models $n\approx N$ and small datasets $T\approx N$, the leading term of the penalization tends to an AIC-like one of order $O(n)$.

3.2 Model selection tests

In general, we expect that results delivered by different model selection recipes will differ more among each other as the dimensionality of the models grows large. For example, the difference between the AIC and BIC penalization factors grows as $O(n\log T)$. This is an important issue for real world datasets for which $T=k\cdot N$ with $N\gg 1$ and $k\geq 1$. In our tests, we consider $N=50$ binary inputs ($N=100$ in Supplemental Information) and $k=5,50$ and $200$. As before, in section 2.2, we generate samples for the inputs according to an Ising model, $p(\bm{x})\propto\exp(-\sum_{i<j}\beta J_{i,j}x_{i}x_{j}-\sum_{i}\beta h_{i}x_{i})$, where $J_{i,j}=1/N$, $\forall i\neq j$, $J_{i,i}=0$ and $h_{i}=0.1$, $\forall i$, where the small parameters $h_{i}$ have been introduced to break the symmetry between the two configurations with all inputs being equal. The choice of an Ising distribution is motivated by the fact that it is easy to tune the localisation of the distribution by varying the parameter $\beta$.

Given the input data, we generate the output according to the conditional probability distribution $p(y|\bm{x})$ of equation 1 employing only a random subset $n\leq N$ of the predictors. This is achieved by assigning non zero values in the parameter vector only to the weights multiplying the selected predictors. The number of selected predictors, $n$, is decided according to the desired level of sparsity, $s=1-n/N$. Each non zero weight is drawn independently from the distribution defined as: $P(w)=1/2$, for $|w-1|\leq 1/2$ and $|w+1|\leq 1/2$; $P(w)=0$ otherwise. After that, the parameter vector is normalized by $\sqrt{n}$ to ensure that the effective field acting on the output is of $O(1)$ regardless of the number of non-zero parameters in the model.

Thus, after having defined the ground truth and generated $T$ samples from it, we aim at recovering its structure from the data by doing model comparison: we rank the models according to their maximum likelihood value and penalise them for their complexity by employing equation 11. We compared the performance with those obtained with some of the most common approaches, i.e BIC (Schwarz,, 1978), AIC (Akaike,, 1974), $\ell_{1}$ regularization (Tibshirani,, 1996; Koh et al.,, 2007) with a K-fold cross validation procedure for selecting the regularizer ($K=5$). Except for the $\ell_{1}$ method where the search in the space of models is embodied in the optimization, for the other methods we need to define the pool of candidate models. In fact, the set of all possible candidates grows exponentially with the number of predictors so that an extensive search among all possibilities is impracticable for large values of $N$. We therefore resort to approximate procedures to walk through this huge model space. We use a decimation procedure (LeCun et al.,, 1990; Hassibi and Stork,, 1993; Decelle and Ricci-Tersenghi,, 2014): we start from a full model with all predictors $N$ connected to the output and estimate its parameters; then the next model to include in the pool would employ all predictors as in the previous model except for the one whose corresponding interaction parameter had been inferred as having the least absolute value; at each step, we define new models by removing the parameter with the smallest absolute value until we reach the model with no active interactions at all, namely the model with the output independent from all candidate predictors. We have also tried a “forward” approach by starting from the null model and adding at each step the parameter that would cause the biggest jump in the likelihood. The reconstruction errors with both procedures are found to be almost indistinguishable in our simulations and, therefore, we report only the results with the decimation procedure. These errors are plotted in Figure 3 versus the sparsity of the ground truth and for different levels of the localisation of the input distribution $\beta$ (columns) and sample size $T$ (rows). The errors are the mean fraction of misclassified non-zero/zero parameters over $100$ independent realisations of the same experiments and the error bars are the corresponding standard deviations. In the figure, the sparsity level varies from [math] to $0.8$ meaning that the ground truth in the first case is a model with all weights being non-zero, whereas in the second case $80\%$ of them are set to zero.

Clearly, the reconstruction error decreases as the sample size $T$ increases for all methods, but more slowly for AIC. On the other hand, at increasing $\beta$, model retrieval becomes harder. Consistently, the worst reconstruction is observed for the largest value of $\beta$ and the smallest sample size ($\beta=1.5$ and $T=5N$ in Figure 3). In the latter case, the likelihood is quite flat so that a very sparse model is typically selected by any model selection criterion. In fact, the information on the activity of the output provided by a new predictor which is highly correlated with the ones already employed in the model is expected to be small because much of its variability is already encoded in the set of predictors in use. This is true even if this predictor is indeed present in the ground truth. On the other hand, when all the predictors are maximally different from each other (maximal information conveyed by each predictor) and the data length is large, the likelihood is expected to be quite sharp around the ground truth, allowing for a very good reconstruction. This is indeed the case of $\beta=0.01$ and $T=200N$ in Figure 3 where each model selection recipe achieved its minimal reconstruction error.

BIC and $\ell_{1}$ regularization employ a stronger penalization than AIC which explains why in the figure the latter methods are always better than AIC in recovering sparse models, but worse when it comes to unveiling dense models. The proposed penalization in equation 11 acts against this unbalance in the recovery errors between sparse and dense models. In fact, it tends to match the BIC trend for high sparsity and the AIC one for low sparsity, resulting in the best model selection recipe in our simulations: in Figure 3, for any level of sparsity, the reconstruction error with our method is the smallest or gets close to the smallest one achieved by the other methods whose performances instead strongly depend on the sparsity of the ground truth. In Supplemental Information (Figure 4), numerical experiments with $N=100$ further substantiates these results. This result is very important in applications given that normally we are not given any information on the expected number of predictors (sparsity) and therefore we might want to treat sparse and dense models on the same footing. As shown in this section, our method achieves this goal whereas standard alternatives typically introduce a bias toward some sparsity range.

4 Degenerate models

In this section we consider degenerate models, namely those models for which some components of the parameter vector ${\bm{\theta}}$ are constrained to be equal. Evaluating the complexity for degenerate models is particularly interesting for model selection given that usually model selection criteria, such as Bayesian Information Criteria and Akaike Information Criteria, can distinguish models only based on the number of parameters employed: for instance, a model with one parameter connecting one input to the output would be equivalent to a model connecting $n$ inputs through the same parameter. For logistic regression models, it is easier to see that any degenerate model, i.e. any model built from a non-degenerate model by explicitly constraining a subset of parameters to be equal, can be mapped into a non-degenerate model with categorical inputs 555From equation 1, the predictor of a degenerate parameter is the sum of all binary predictors connected to the output through the same parameter, i.e. a categorical predictor.. Therefore the additional complexity of a degenerate model with respect to a non degenerate one with the same number of parameters is related to the additional complexity conveyed by categorical inputs with respect to binary ones.

In the following we will consider a simple example of such models and discuss an application on a real dataset.

4.1 1-parameter degenerate models

The simplest example of such models is a logistic regression model with $n$ binary inputs $x_{i}$ with $i=1,...,n$ and the corresponding $n$ weights tied together through the same parameter $\theta$:

[TABLE]

where $b$ is a fixed parameter which works as a threshold. Given $T$ observations of the inputs and the output $X^{(t)}=\{y^{(t)},\bm{x}^{(t)}\}$, for $t=1,...,T$, the normalised log-likelihood of the model is given by

[TABLE]

where $X_{n}=\sum_{i=1}^{n}x_{i}$. As the form of the likelihood suggests, this model can be mapped into a model with only one categorical input $X_{n}\in\{-n,-n+2,...,n-2,n\}$, connecting to a binary output $y$ through the parameter $\theta$. Consequently the Fisher Information takes the form $F(\theta)=\sum_{X_{n}}\nu(X_{n})\,{(X_{n}-b)}^{2}\cosh^{-2}(\theta(X_{n}-b))$, where $\nu(X_{n})$ represents the distribution of the sum of random variables $x_{i}$, and the complexity is simply the integral of its square root. As we did in section 2.2 for non-degenerate models, we resort to triangular inequality to study the bounds of the complexity by varying the distribution $\nu(X_{n})$ and the dimension $n$. This procedure leads to the simple inequality:

[TABLE]

It can be easily shown that the upper bound in the right hand side of the previous equation is maximised when the distribution is uniform over all possible outcomes of the variable $X_{n}$ which means $\nu(X_{n})=1/(n+1)$. By substituting this form of the distribution in equation 14, we find that $e^{C}\leq\pi\sqrt{n+1}$. Thus the complexity grows at most as the square root of the number of inputs which is a much slower growth compared with the exponential one in non-degenerate models. Therefore, it is clear that the additional complexity gained by growing the alphabet of a categorical variable is generally very small compared to the additional complexity obtained by adding more independent parameters.

On the other hand, the minimum value for the right hand side of equation 14 will be reached when only one particular value, say $\bar{X}$, of the variable $X_{n}$ is observed, namely when $\nu(X_{n})=\delta_{X_{n},\bar{X}}$, where $\delta$ stands for the Kronecker delta function. In this case, the resulting complexity is exactly $e^{C}=\pi$, as we saw for a one-parameter model with a binary input in section 2.1.

An interesting case is when all inputs are uniformly distributed so that in a large sample ($T\gg 2^{n}$) we expect to observe all possible $2^{n}$ configurations the same number of times. In this case, the number of occurrences of a particular value of $X_{n}$ is given by $\binom{n}{(n+X_{n})/2}$ 666The number of plus ones in a pattern with $X_{n}=x$ is $(n+x)/2$; it follows that the combinatorial factor $\binom{n}{(n+x)/2}$ counts all possible patterns with the same number of plus ones, i.e. the same value of $X_{n}$.. It follows that the distribution of $X_{n}$ will be $\nu(X_{n})=\binom{n}{(n+X_{n})/2}2^{-n}$. By plugging this distribution into equation 14, we find $e^{C}\leq\pi\,2^{-n/2}\sum_{k=0}^{n}\sqrt{\binom{n}{k}}$. The right hand side of the inequality is very well approximated by the function $\pi(a_{1}\cdot n+a_{2})^{1/4}$ with $a_{1}=2\pi$ and $a_{2}=\pi/\sqrt{2}$ which gives a better intuition of the behaviour of the complexity in this particular case: it increases with $n$ but much slower than the upper bound; it is always above the lower bound identified by $\pi$ and, as $n$ goes to zero, it tends to the same value taken by the upper bound. Notice that both the upper bound and this particular solution do not tend to $e^{C}=\pi$ for $n=1$ because these results are obtained by employing triangular inequality.

We present these analytical results in Figure 4 (dotted lines) at varying the size of the input layer $n$. In the same figure, we compare these trends with simulations obtained by evaluating the one dimensional integral $e^{C}=\int d\theta\sqrt{F(\theta)}$ numerically (solid lines). In evaluating the integral we fixed $b=0$ and studied it at varying $n$. As it is clear from the figure, the complexity oscillates while very slowly increasing with $n$ in both cases of uniformly distributed random inputs (red line) and uniformly distributed sum of inputs (blue line). The oscillation is due to the fact that when $n$ is even, the categorical variable $X_{n}$ can take exactly the value of the threshold $b$, which has been set to zero. The complexity increases with the dimension $n$ more rapidly for models with a small number of predictors, whereas its value becomes almost independent of $n$ for large models. Moreover, it is interesting to notice that the complexity of a one parameter degenerate model for $n<100$, as shown in 4 (dotted line), is always smaller than the maximum value of the complexity for a two parameter non degenerate model which is $e^{C}=\pi^{2}/2$, as we have seen in section 2.1.

As a conclusion, given that a one parameter degenerate model with $n$ inputs can be mapped into a model with only one categorical input with an alphabet of $n+1$ symbols, the additional complexity conveyed by categorical inputs with respect to binary ones increases very slowly with the cardinality of the alphabet of the categorical variable, showing a larger effect for small values of $n$, and it was found to be lower than the increment in complexity which would follow from adding another non degenerate parameter to the model, even for large values of $n$.

4.2 The 13 keys to the White House

An interesting application is the question of predicting the U.S. presidential elections from the values of some predictors. A successful method for forecasting the verdict of US presidential elections is “The 13 keys to the White House” (Lichtman and Keilis-Borok,, 1981). The method was developed in the eighties and, since then, has been able to predict correctly presidential elections (Lichtman,, 2016). It is based on 13 binary questions (keys) and each time the answer to a key is false, the key is turned against the party in power. If six or more keys are false, the challenging party is predicted to win the elections otherwise the party holding the White House will be more likely to win. Since the predictions depend only on the sum of the keys, i.e. how many of them are false, the method can be interpreted as the deterministic counterpart of a 1-parameter degenerate model (see Supplemental Information, section 3).

The question we would like to ask from a Bayesian Model Selection point of view is: is there any “simpler” explanation of the outcomes of presidential elections which would perhaps involve less predictors? We compare all possible models corresponding to all possible ways of selecting a subset of predictors from the original set of the 13 keys and use that subset for forecasting the presidential elections (see Supplemental Information). This amounts at comparing $2^{13}$ models. Since all models, except the model employing no predictors at all, can be cast into degenerate models with only 1 parameter, they will all have the same BIC penalization factor. It follows that the largest penalization factor which differs among models is represented by the complexity.

In Figure 5 a) we plot the posterior probability $p({\cal M}|\hat{X})$ for all models ranked according to their maximum likelihood values. The figure clearly shows that the effect of taking into account the complexity on model selection is local, i.e. it concerns models which are very close in likelihood. In fact, the posterior seems to be monotonically increasing on a larger scale while it presents fluctuations on a finer scale (insets). Therefore, the complexity induces a local reordering of the models with respect to the likelihood ranking. We measure the range of such an effect by introducing the quantity $\Delta rank$. Specifically, given all models ranked in a non-decreasing order according to their likelihood, $\Delta rank$ is the range of models surrounding a given one within which the fluctuations of the complexity are approximately equal to the variation of the likelihood. In some sense $\Delta rank$ is a measure of the uncertainty induced by the complexity in the likelihood ranking. It varies across models with a mean value of $\overline{\Delta rank}\approx 26$, as shown in Figure 5 b).

Moreover, Figure 5 a) reveals that, interestingly, there are some models in our pool that are attached a posterior significatively larger than the others. There are 42 of such models and they are reported in Figure 6 a). These models fit the dataset exactly ($\ell\approx 0$). It follows that, given that the BIC term is the same overall, the complexity is the only term varying across models. In Figure 6 b) the models are ranked according to their posterior (or equivalently from the most to the least complex model, given that the likelihood and the BIC term are constant). As mentioned in the previous section, the complexity does not depend merely on the number of predictors but also slightly on the distribution of the corresponding categorical variable attached to them. In fact, the model with all 13 keys active turns out to be not the most complex model although still among the most complex ones. On the other hand, the model with the largest posterior overall, contains only 8 out of the 13 keys (Keys 1,2,3,4,6,7,8,13) where the first seven keys are among the most frequent keys over all 42 models. Interestingly, by looking at the set of predictors not included in the best model (Keys 5,9,10,11,12), we could see that: according to this model, it does not matter for the forecast of the next presidential election whether the economy is in recession during the electoral campaign (key 5), but what does matter is rather the economic growth during the whole previous term (key 6); military failures or successes are irrelevant as predictors (keys 10 and 11); finally, it does not matter whether the incumbent administration is untainted by a major scandal (key 9) or whether the incumbent-party candidate is charismatic or a national hero while it is relevant whether the challenging-party candidate is charismatic (key 13)777We refer to (Lichtman,, 2016) for a detailed explanation of the meaning of the keys.

As a conclusion, by comparing all possible degenerate models with different predictors through a Bayesian model selection approach, we found that the original method of the 13 keys to the White House is not attached a posterior probability larger than many other candidates. Instead, a model with only 8 predictors is the one that achieves the largest posterior. However it is worth noting that the differences in the values of the posterior among these 42 most likely models are very tiny and eventually the result of the model comparison would strongly depend on the choice of the prior over models. Our choice of the prior favours the least complex models among the best candidates, i.e. the one which account for the minimum number of distinguishable probability distributions.

5 Conclusions

In this paper we address the question of estimating the complexity of a logistic regression model with binary inputs. As shown in Balasubramanian, (1997), the complexity of a model can be related to the number of distinguishable distributions that a model can represent. In this perspective, it is defined as the integral of the square root of the determinant of the Fisher Information matrix. We evaluate analytically the integral in special cases and corroborate the analysis with Monte Carlo simulations. We find that the integral depends on the degree of correlation of the inputs: correlations among inputs induce effective dependencies on the parameters thus reducing the volume that the model occupies in the space of distributions. It follows that when the inputs are maximally uncorrelated, the model achieves its maximal value of the complexity, whereas when the inputs are highly correlated, it becomes highly constrained and it is characterized by a low value of the complexity. The upper bound of the complexity has been estimated with Monte Carlo methods. We find that it can be very well approximated by the simple relation $e^{C}=\pi^{n}/n$, where $n$ is the number of parameters. The lower bound instead is reached when the number of unique patterns observed in the input layer is equal to the number of parameters. In the latter case, we find an exact solution for the value of the complexity, $e^{C}=\pi^{n}\prod_{i=1}^{n}\sqrt{\nu(\bm{x}^{i})}$. The solution depends on the distribution $\nu(\bm{x}^{i})$ of these $n$ patterns and ranges from $e^{C}=\pi^{n}\sqrt{(T-n+1)/T^{n}}$ to $e^{C}=\pi^{n}/n^{n/2}$. Monte Carlo estimates of the integral at varying degrees of correlation (localisation of input distribution) further confirm these results. Finally when the number of unique patterns observed in the input layer is less than the number of parameters, the data matrix is not full rank and $e^{C}=0$: the model is constrained into a lower dimensional manifold in the space of distribution induced by the deterministic relation among parameters.

Moreover, we investigate the complexity of models whose parameters are defined on a finite support (regularized models). In the case of maximally uncorrelated inputs, we find that the complexity starts decreasing with the dimensionality of the model after a certain critical dimension as opposed to the steady exponential increase, i.e. $e^{C}=\pi^{n}/n$, found for models with parameters defined on the entire real axis. This is due to the gradual appearance of dependencies among parameters induced by emergent statistical constraints. Thus, for large values of $n$, increasing the number of parameters in regularized models will always increase the degrees of freedom of the model, but it will also make the model more constrained (less complex), at odds with the corresponding un-regularized case.

We apply these results on the model complexity for devising a novel model selection criterion that would take into account the correlations among inputs when deciding how much a model should be penalised in a Bayesian model selection framework. To this purpose, we introduce a heuristic for the penalization cost based on the entropy of the input distribution. Consistently, our proposal assigns a strong penalization, similar to a BIC-like one, to models with a large value of the complexity and a weaker one, more similar to an AIC-like penalization, to extremely constrained models. We then tested this ansatz by simulating data from logistic models with various levels of sparsity, sample size and input correlations, and then asking different model selection criteria to retrieve the correct model. Numerical tests clearly show that our proposal always achieves (or gets very close to) the smallest reconstruction error for all levels of sparsity, data size and input correlations whereas performances with the other tested methods strongly depend on the sparsity of the ground truth. As a consequence, our method qualifies as the best choice in our tests and would be especially useful in those cases in which a priori we have no information about the sparsity of the ground truth. As an example, this might be the case for a population of neurons simultaneously recorded from some region of the brain. In fact, the sparsity of the network of dependencies among neurons strongly varies according to if the recorded cells belong to a functionally connected module or many independent ones as well as according to how the behavioural task during the recording shapes the activity of the cells (Dunn et al.,, 2015). These information might often be not known to the scientist a priori.

Finally, we address the question of how the alphabet size of the inputs affect the value of model complexity. The problem is equivalent of studying the complexity of a model with binary inputs, but with degenerate parameters, namely when subsets of parameters are forced to be equal. In this respect, we study analytically and numerically the simple case of a model with only one parameter by varying its degree of degeneracy or, alternatively, the alphabet size of the corresponding effective input. We find that the additional complexity conveyed by categorical inputs with respect to binary ones increases very slowly with the alphabet size of the categorical variable, showing a slightly larger rise for small sizes, and it is always lower than the increment in complexity which would follow from adding another non degenerate parameter to the model.

Studying the complexity of degenerate models is also interesting because many commonly used model selection recipes are based on penalising models according to some cost function which often depends only on the number of parameters and sample size. A Bayesian Model selection approach would instead employ a penalization which might differ even for models with the same number of parameters. As an application, we study the dataset of the 13 keys to the White House (Lichtman and Keilis-Borok,, 1981; Lichtman,, 2016). In fact, the question of deciding if a model with less predictors would be better at explaining the variability in the dataset with respect to the full model with all 13 predictors, can actually be cast into a model selection problem with 1-parameter degenerate models. As one would expect, the contribution of the complexity is very small as compared with that coming from the likelihood. Yet, we show that it becomes relevant for deciding among the subset of models with the largest value for the likelihood.

It would be interesting to extent this study of the complexity to a more general class of models known as Generalized Linear Models (Nelder and Wedderburn,, 1972) which is widely employed in applications, for instance in neuroscience for studying the neural code of cells or populations of cells in the brain (Roudi et al.,, 2015).

Another interesting application is related to the fact that one of the most powerful approaches for inference and structure learning in Ising model is based on the so-called pseudo-likelihood approximation (Ravikumar et al.,, 2010; Aurell and Ekeberg,, 2012; Decelle and Ricci-Tersenghi,, 2014). This approximation basically consists of substituting the exact likelihood function, whose optimization is hard because of the partition function, with a sum of conditionally independent logistic regression functions with binary inputs and outputs. Within this approximation, based on the results of section 3, we would expect that our proposed criterion (equation 11) will deliver close to optimal performances for nodes with any level of degree distribution (sparsity) and, therefore, on a wide range of network topologies, including topologies where highly connected nodes coexist with almost isolated ones, as for instance in scale free and small world networks.

Finally, it must be noticed that equation 11 is a heuristic inspired by the truncated posterior expansion of equation 5 so it should inherit the limitations buried in the truncation. In particular, terms of order 1/T, which are neglected in equation 5, might be relevant in some regimes (Balasubramanian,, 1997). As shown also by our numerical tests, when the sample size approaches the dimensionality of the model or for large input correlations, model selection represents a hard task and performances are generally poor with any approach. It would be interesting to study higher order terms in the expansion and see whether the information contained in them could be used to improve model selection in this regime.

6 Acknowledgments

The authors are grateful to Prof. Benjamin Dunn for reading through the paper.

Appendix A Upper bound for models including the bias term

Here we prove that due to the parity of the hyperbolic cosine, the complexity of a logistic regression model with $n$ binary inputs, where $n-1$ of them are uniformly distributed, is equal to the complexity of the same model when all $n$ inputs are uniformly distributed, which achieves the upper bound, as we found in section 2.2.

To prove the statement, let’s consider the Fisher Information elements for a model with $n$ inputs connected through a $n$ dimensional vector of parameters $\bm{\theta}$

[TABLE]

Let’s perform the sum over an arbitrary variable $x_{k}$,

[TABLE]

where $k\neq i,j$, $\bm{x}_{/k}$ and $\bm{\theta}_{/k}$ are the $n-1$ dimensional vectors obtained by removing the $k$-th index and $\nu^{+}(\bm{x}_{/k})=\nu(x_{k}=1;\bm{x}_{/k})$ and $\nu^{-}(\bm{x}_{/k})=\nu(x_{k}=-1;\bm{x}_{/k})$. By changing $\bm{x}_{/k}$ with $-\bm{x}_{/k}$ in the second summation we obtain

[TABLE]

Notice that $\sum_{\bm{x}_{/k}}$ and $\sum_{-\bm{x}_{/k}}$ is a sum over the same elements but in a reversed order. It follows that we can write it as

[TABLE]

We have considered here $k\neq i,j$, but the latter equation is true for any elements of the Fisher Information matrix. In fact, it is not difficult to see that equation 18 holds in the case when $k=i=j$. For the cases $k=i\neq j$ and $k=j\neq i$, there will be a minus sign in front of the second summation in equation 16 which will become again a plus sign in equation 17 after transforming $\bm{x}_{/k}$ in $-\bm{x}_{/k}$. Now, if the distribution of the inputs satisfy the relation

[TABLE]

then the elements of the Fisher Information matrix can be written as if the k-th input was absent and the corresponding weight was the bias term acting on the output, namely

[TABLE]

In particular, equation 19 is true if the distribution factorize over the spin variable $k$, i.e. $\nu(\bm{x})=\nu(x_{k})\nu(\bm{x}_{/k})$ and if all possible patterns and anti-patterns of the remaining $n-1$ spin are observed with the same frequency, i.e. $\nu(\bm{x}_{/k})=\nu(-\bm{x}_{/k})$. One instance in which this is true is when $n-1$ out of the $n$ inputs are uniformly distributed, i.e. $\nu(\bm{x}_{/k})=2^{-n+1}$, regardless of $\nu(x_{k})$. It follows that, the upper bound for the complexity is reached already in the case in which the input distribution is flat over the configuration space of only $n-1$ inputs, as it has been observed in the simulations and discussed in the paper in section 2.2.

As a consequence, since the bias term can be thought as a weight parameter connecting a constant input (always equal to $1$ or $-1$), the complexity of a logistic model with $n-1$ inputs (weights) and the bias term reaches the same upper bound as the model with $n$ inputs (weights) and no bias. The upper bound is attained when the $n-1$ inputs are uniformly distributed.

Appendix B Proof of equation 10

We consider logistic models with parameters defined on a finite support $|\theta_{k}|\leq\Theta$, $\forall k$ with $\Theta$ being any finite positive number, in the case of uniformly distributed inputs, $\nu(\bm{x})=2^{-n}$, $\forall\bm{x}$. Recall that the Fisher Information matrix elements are given by $F_{i,j}(\bm{\theta})=\sum_{\bm{x}}\nu(\bm{x})\cosh^{-2}(\bm{\theta}\cdot\bm{x})x_{i}x_{j}$.

We examine first the off-diagonal terms. We carry out the summations over the variables $i$ and $j$ and define the variables $\hat{\bm{\theta}}$ and $\bm{\hat{x}}$ by removing the indices $i$ and $j$ from the variables $\bm{\theta}$ and $\bm{x}$ respectively,

[TABLE]

where we have first exploited the parity of the hyperbolic cosine and then rearranged the terms in the summations. We think of the average over observations as the average over a random process with the same limiting distribution, i.e. $p(\bm{x})=\nu(\bm{x})$. This will allow us to exploit central limit theorem and perform calculations. It follows that, since the distribution over $\bm{x}$ is uniform, the quantity $z\equiv\hat{\bm{\theta}}\cdot\bm{\hat{x}}=\sum_{k\neq i,j}\theta_{k}x_{k}$ is the sum of $n-2$ independent random variables, $z_{k}=\theta_{k}x_{k}$, distributed as $p(z_{k}=\pm\theta_{k})=1/2$. Therefore each random variable $z_{k}$ has zero mean and variance $\theta_{k}^{2}$. Asymptotically, i.e. as $n$ goes to infinity, $z$ becomes a Gaussian random variable with zero mean and variance $r^{2}=\sum_{k\neq i,j}\theta_{k}^{2}$. This is ensured by Lyapunov’s condition which is always satisfied when the parametric family is bounded (see for instance Billingsley,, 1995). Thus, in the last equation, we substitute the expectation over the configuration states $\hat{\bm{x}}$ with an expectation over $z$,

[TABLE]

where $p(z)=e^{-z^{2}/2r^{2}}/\sqrt{2\pi r^{2}}$ is the gaussian density and where we used the change of variable $z\rightarrow z+\theta_{i}$ for passing from the first to the second line. The integrand in the last equation is the difference between two functions centered at $z=\pm\theta_{j}$ and weighted according to a gaussian density centered at $z=\theta_{i}$. However the standard deviation of the Gaussian density grows with $n$ and eventually for $n\gg 4\Theta^{2}$, the two contributions will be weighted equally cancelling each other. More specifically, for $r\gg 2\Theta$, namely when $n\gg 4\Theta^{2}$, the density can be approximated as $p(z)\sim 1/\sqrt{2\pi r^{2}}+O(1/r^{3})$. The term in $O(1/r)$ has a zero contribution to the integral so the first non zero contribution comes from the $O(1/r^{3})$ term. Thus off-diagonal elements go to zero faster than $O(n^{-1/2})$ which ensures that these terms will not bring a relevant contribution to the determinant as $n$ goes to infinity.

Following the same argument, it is easy to obtain a simple expression also for the diagonal terms $F_{i,i}(\bm{\theta})=\sum_{\mu}\nu(\bm{x}^{\mu})\cosh^{-2}(\bm{\theta}\cdot\bm{x}^{\mu})$, $\forall i$. The latter will become $F_{i,i}(\bm{\theta})=\int_{-\infty}^{\infty}dz\,p(z)\cosh^{-2}(z)$ where $p(z)=e^{-z^{2}/2r^{2}}/\sqrt{2\pi r^{2}}$ and $r^{2}=\sum_{k}\theta_{k}^{2}$. The elements on the diagonal are all equal and their value depends only on $r$, namely $F_{i,i}(\bm{\theta})=\lambda(r)$ where

[TABLE]

The square root of the determinant of the Fisher Information is simply given by $\lambda^{n/2}(r)$ and, given that $\lambda(r)$ depends only on the euclidean distance $r$ in the space of parameters, it’s easier to perform the integral in spherical coordinates. Thus, $e^{C}=\int d\bm{\theta}\lambda^{n/2}(r(\bm{\theta}))$ becomes $e^{C}=\int d\Omega dr\,r^{n-1}\lambda^{n/2}(r)$ and by integrating over the solid angle $\Omega$, it becomes

[TABLE]

where $\Gamma(x)$ is the Euler gamma function. In the large $r$ limit, since $\lambda(r)\sim\sqrt{2/\pi r^{2}}+O(1/r^{3})$, the product $r^{n-1}\lambda(r)^{n/2}\sim({2/\pi})^{n/4}r^{n/2-1}$. Employing the latter approximation, the above integral becomes $e^{C}\sim 2(2\pi\Theta^{2}n)^{n/4}/\Gamma(n/2+1)$. Finally, using Stirling’s approximation for the gamma function, $\Gamma(n/2+1)\sim\sqrt{\pi n}(n/2e)^{n/2}$, the complexity in the large $n$ regime can be approximated by equation 10

[TABLE]