Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
49,711 characters · 25 sections · 33 citation commands
Structured Probit Estimation with Expectation Propagation
TODO: 0. Sort out issues with zero and non-zero coefficients. 1. Subsampling in the E step to speed up. 2. Apply our algorithm to the sign-based probit model as well.
Main contributions: EP for probit; increase computational efficiency; theory? Should having covariates help or hurt the estimation?
We have $n$ individuals, faced with $m$ choices. Individual $i$ has latent choice utility $Z_i = X_i\beta + \epsilon_i$, observed choice-covariate matrix $X_i \in \mathbb{R}^{m \times p}$ and decision $Y_i \in \{1, 2, \ldots, m\}$. $\beta \in \mathbb{R}^p$ is the unknown covariate coefficient and $\epsilon_i \sim \mathcal{N}_m(0, \Sigma)$ is the utility shock with choice covariance $\Sigma$. We are interested in the setting where the individual picks only one option, $j$--the one with greatest latent utility such that $Z_{ij} \ge Z_{ik}$, $k \in \{1, \ldots, m\}$. We assume sparsity in $\Sigma$, and impose this via a penalty term on the elements of $\Sigma$ in the likelihood. Sparsity structure is naturally motivated by applications. For example, categories of distinct goods result in a block diagonal structure of $\Sigma$. Moreover, an additional potential benefit of imposing sparsity is that it may help with the identification issue, since the model is scale and level invariant and hence only a subset of parameters of a general $\Sigma$ can be identified.
loaiza2021scalable uses MCMC based on a factor structure with trace constraint to achieve scalability and identification. They adopt a Bayesian approach and find priors that enforce the trace constraint (scale normalization), and use the reduced system with respect to a reference alternative for level normalization.
We also need to compare with the literature on simulated maximum likelihood approaches, which basically uses importance sampling to approximate the gradient of the likelihood, then use Newton/gradient based methods to maximize the likelihood directly, instead of using EM algorithm to maximize indirectly. For example, the GHK sampler uses the Cholesky factorization of the covariance matrix to simplify the task of sampling from multivariate truncated normal to sampling from one dimensional truncated normals, which is easy. However, the main shortcoming of this approach that makes it challenging to handle large choice sets is the fact that later samples of one-dimensional truncated normal depends on earlier samples.
We may also be able to use EP to estimate the gradient of the log likelihood, and perform direct MLE instead of EM. How would this compare with simulated likelihood?
We should also add back the penalty approaches (need to update the Newton method to reflect the added penalty).
One last thing to add is the normalization constraint that latent utilities sum to zero.
So it seems to me that there is a 2x2 combination table of methods that can be used to perform MLE for probit in discrete choice. In one dimension, we have sampling vs. EP to estimate moments and gradients of the likelihood; in the other, we have direct/gradient based optimization methods (such as Newton) vs. EM to perform the actual optimization part. They give rise to four methods. Right now we are focusing on one of these four, using EM algorithm for optimization, and then EP for estimating moments. The simulated maximum likelihood approach popular in the 90s in economics uses direct optimization methods with sampling. I'm wondering how the other two methods would perform. For example, what if we use EP to estimate the gradient of the likelihood (this should be possible) and use direct method? How does EP compare with GHK sampler at estimating moments in terms of speed? That'd be important for us to justify our approach.
Given observations $(X_i, Y_i)$, the covariance penalized log likelihood is
where $\mathcal{A}(Y_i) = \{z_i : z_{iY_{i}} \ge z_{ik}, \, k \in \{1, \ldots, m\}\}$. The complete penalized data log likelihood is
The E step requires computing $Q(\beta, \Sigma) = E_{Z|X,Y,\beta, \Sigma}\left[ \log L(\beta, \Sigma | Z, X, Y) \right]$, and hence computing
(ref) requires computing the first and second moments of the latent utilities $Z_i$ given $Y_i$ and current estimates of $\beta$ and $\Sigma$. Given $Y_i = j$, we know that $Z_{ij} \ge Z_{ik}\,, k \in \{1, \ldots, m\}$, so $Z_i | Y_i=j, X_i, \beta, \Sigma \sim \mathcal{TN}_{\mathcal{A}(Y_i)}(X_i\beta, \Sigma)$, a truncated multivariate normal (TMVN) distribution.
After the E step we have the lower bound of the log likelihood
where $S$ is the sample covariance
The M step requires optimizing (ref). This objective is of the form studied in bien_sparse_2011. We use their method to find a new sparse estimate of $\Sigma$. Other possible methods for optimizing (ref) include thresholding and banding bickel_covariance_2008.
Computing the truncated multivariate normal (TMVN) moments (ref) is expensive for large choice sets $m$. The expectation propagation (EP) approximation of cunningham_gaussian_2013 can provide a fast and relatively accurate estimate of the moments. The idea is to approximate the target TMVN density with an unrestricted multivariate normal density using EP minkaExpectationPropagationApproximate2001a. Once we have the approximate multivariate normal's mean and covariance, we use those as estimates of $E[Z_i|Y_i,X_i,\beta, \Sigma]$ and $Var[Z_i|Y_i,X_i,\beta, \Sigma]$ in the E step.
The method proposed in cunningham_gaussian_2013 applies an approximate Bayesian inference algorithm to computing multivariate normal density integrals.
In Bayesian inference we are interested in probabilistic queries involving the posterior density over latent variables $z$ given data $x$,
For example we may want to compute posterior expectations $\int z p(z|x)dz$. The challenge is that (ref) is often intractable because of the marginal likelihood $p(x)$. Markov chain Monte Carlo (MCMC) methods avoid this problem by drawing samples from $p(z|x)$ without needing knowledge of $p(x)$, which can be used for Monte Carlo estimates of posterior quantities of interest. MCMC methods will eventually sample from the true posterior after enough iterations. However they are computationally intensive.
Variational inference methods bleiVariationalInferenceReview2017 are a more efficient, albeit approximate, alternative to MCMC methods. The idea is to find a density $q^*(z)$ from a tractable family $\mathcal{Q}$ that best approximates $p(z|x)$ (Figure (ref)).
$q^*(z)$ is determined by minimizing a divergence $D$ between the true posterior $p(z|x)$ and approximation $q(z) \in \mathcal{Q}$, where $\mathcal{Q}$ is a family of distributions over $z$:
We will need to specify the approximating family $\mathcal{Q}$ and the choice of divergence $D$. The approximating family $\mathcal{Q}$ needs to trade off flexibility, which allows accurate approximation, with computational speed. The choice of divergence $D$ also has an impact on the approximation. The forward Kullback-Leibler (KL) divergence
is an inclusive divergence that promotes approximations $q^*$ that cover as much of the nonzero probability regions of $p(z|x)$ as possible, and also well approximates the normalizing constant of $p(z|x)$ minkaDivergenceMeasuresMessage2005.
However this objective is still intractable since it involves $p(z|x)$. Expectation propagation (EP) minkaExpectationPropagationApproximate2001a is a method that circumvents this issue by minimizing local forward KL divergences. It sets $\mathcal{Q}$ to be the family of multivariate normal distributions. This combination of choices leads to an approximation that captures the zeroth (normalizing constant), first, and second moments of the target posterior.
cunningham_gaussian_2013 use EP to approximate truncated multivariate normals with unconstrained multivariate normal distributions (Figure (ref)).
cunningham_gaussian_2013 derives a deterministic expectation propagation approximation for estimating truncated multivariate normal (TMVN) first and second moments under linear domain constraints. The idea is to approximate the target TMVN with an un-truncated multivariate normal (MVN) that has tractable moments.
Let the target $d$-dimensional TMVN density of interest be $p(x)$, where
and $m$ is the mean, $K$ is the covariance, $\zeta$ is the normalizing constant, $l,u \in \mathbb{R}^{c}$ are the lower and upper bounds, and $A \in \mathbb{R}^{c \times d}$ is the linear constraint matrix. We rewrite $p(x)$ as a product of factors:
where $p_0(x) = \frac{1}{\zeta} \exp \left(-\frac{1}{2}(x - m)^T K^{-1}(x-m) \right)$, and
where $a_k$ is the $k$th row of $A$, the $k$th linear constraint.
We want to approximate $p(x)$ with $q(x)$, where $q(x)$ has the same factorized form as $p(x)$:
and
$p_0(x)$ and the $\tilde{t}_k(x)$ are all normal densities, so the product of their densities will also be normal. Therefore $q(x)$ has a multivariate normal density,
Our goal is to find $\mu$ and $\Sigma$ such that $q(x)$ is a good approximation to $p(x)$.
cunningham_gaussian_2013 use expectation propagation minkaExpectationPropagationApproximate2001a to achieve this. The algorithm iterates over the factors $k = 1, \ldots, c$:
After iterating through all the factors this procedures repeats, until $\mu$ and $\Sigma$ converge. Upon termination, we use $\mu$ and $\Sigma$ as approximations of the mean and covariance of the TMVN $p(x)$.
EP does not have convergence guarantees. Empirically, cunningham_gaussian_2013 find that it is fast and accurate for multivariate normal probability estimation, especially for rectangular constraints where $A = I_d$. They found some failure cases for EP, for example when there are many redundant constraints.
The forward KL divergence objectives in EP ($D_{KL}(p(x)|| q(x)$) are inclusive, meaning the approximating distribution tries to cover all nonzero probability regions of the true target density. This is in contrast to the exclusive reverse KL divergence ($D_{KL}(q(x)|| p(x)$) minimization, meaning the approximating distribution seeks a mode of the target density. See minkaDivergenceMeasuresMessage2005.
The EP algorithm minimizes local forward KL divergences, instead of a global objective $D_{KL}(p(x)|| q(x))$. This is because minimization of the global objective is intractable.
The EP approximation for truncated multivariate normal has been used for the E-step in EM for probit estimation before, in mandt_sparse_2017. That work proposes a binary probit model with fixed latent covariance. Therefore the constraint matrix in that case is the identity matrix. This variant of the probit model has received most of the attention in the statistics literature, whereas the version of probit models that are used in economics to model choice behavior is distinct from the binary version.
In our setting each observation has an associated constraint that depends on the response $Y_i$. If $Y_i=j$, then the $j$th column is all ones, and the constraint matrix is \setcounter{MaxMatrixCols}{20}
with bounds
In our experiments we fix $\beta = 0$. We sample $n=2000$ observations. We set the maximum number of EM iterations to 300. The penalty parameter is fixed at $\lambda = 0.1\sqrt{\frac{\log m}{n}}$, and we set the thresholding value to $\tau = 3\sqrt{\frac{\log m}{n}}$. We set the number of choices $m=3$, with
We also consider estimating $\Sigma^{-1}$ using the graphical lasso friedman2008sparse in the M-step, and all other settings remain the same. In this case the true $\Sigma^{-1}$ is the same as the true $\Sigma$ in (ref).
Figure (ref) illustrates the good efficiency-accuracy trade off of EP. Here, we compare various TMVN sampling methods for Monte Carlo (MC) estimation to EP for TMVN moment estimation in the E step. For each MC method we draw 25 samples from $Z_i|Y_i,X_i,\beta, \Sigma$. EP recovers the true precision, while taking less than a quarter of the time that the fastest MC method uses.
\FloatBarrier
We compare several sparse covariance estimators for the M-step, and use the EP approximation for the E step for all of these methods. As seen in (ref), soft thresholding rothman_generalized_2009 and sparse covariance estimation bien_sparse_2011 can recover the true non-zero elements.
Reminders on a few issues related to identification.
As discussed in train2009discrete, $\Sigma$ is generally not identified due to the scale and level invariance of the probit choice model, and often the reduced $\Sigma$ is used. This involves the following two normalizations due to invariance of the probit model:
The standard identification strategy first subtracts the latent utility of (usually) the last alternative from the other alternatives for level invariance, then normalize the top diagonal entry (or the trace) to be 1 for scale invariance. The resulting covariance matrix $\tilde \Sigma$ for the reduced system is identifiable. In the next Section, we will focus on this identification strategy in the unpenalized setting, and replicate and compare with results from natarajan2000monte.
However, this identification strategy is not suitable if we want to impose sparsity structures on the unormalized covariance matrix $\Sigma$. We may therefore need to explore alternative identification strategies.
In order to preserve the sparsity structure, we first consider the identification strategy requiring the diagonal of $\Sigma$ to be 1.
Following Train's normalization, if we restrict the covariance matrix $\Sigma$ to have normalized diagonals, the following parameters are identifiable:
Unfortunately, we cannot obtain $\sigma$s from $\theta^{\ast}$s.
In general, models with $J$ alternatives will have $\frac{J^{2}-J}{2}-1$ identifiable parameters, as opposed to $\frac{J^{2}+J}{2}$ unnormalized parameters. The difference is $J+1$ parameters. This means that it is not sufficient to normalize the diagonal to 1. However, if we impose the additional constraint that an off-diagonal matrix is equal to 0, the model is identifiable. Sparsity-inducing penalties do not explicitly impose the constraint, but since they are approximations of the $\ell^{0}$ penalty, we should expect sparsity-inducing parameters to help with identification as well. The question is whether there is a way to impose this additional constraint without specifying the particular off-diagonal entry that is zero.
One identification one might think would work is requiring $\Sigma_{11}=1$, which takes care of scale. Then there are $J$ unidentified parameters. We may consider requiring that there are $J$ zero off-diagonal entries to take care of these. However, it is not the case that imposing such constraints will always yield an identifiable model. Consider the following covariance specification:
It can be shown that as long as $\rho_{1}+\rho_{2}=\rho$, the resulting models will be the same.
Another specification that guarantees identification is the one used in burgette2012trace, which assigns a default utility of 0 to one of the choices to handle level invariance. This means that if all other alternatives have negative latent utilities, then the detault alternative is chosen. Otherwise, the alternative with the highest positive latent utility is chosen. The usual $\Sigma_{11}=1$ requirement is used to handle scale invariance. This specification preserves the sparsity structure and ensures identification.
There is also an additional computational concern regarding imposing constraints on the covariance matrix, since the optimization problem in the M step requires solving a problem which is easy in $\Sigma^{-1}$, but generally hard in $\Sigma$. This implies that we need to choose a specification that guarantees identification, and at the same time also allows a straightforward constraint expression in terms of $\Sigma^{-1}$. If we are not requiring sparsity, then a specification in terms of $\Sigma^{-1}$ instead of $\Sigma$ may be more advantageous from an optimization point of view.
In our paper, we adopt the symmetric identification strategies proposed by burgette2012trace and burgette2021symmetric. More specifically, to address scale invariance, burgette2012trace propose to restrict the trace of the covariance matrix. To facilitate the optimization routine, we will adapt their approach and instead place the trace restriction on the precision matrix. Similarly, to address the level invariance, burgette2021symmetric propose to require that the latent utilities sum to zero. Compared to the standard identification strategy, which chooses a base alternative and subtracts its latent utility from those of the other alternatives, this strategy has the advantage that it retains the original covariance matrix structure without reducing to a sub-system. This is desirable because the original covariance matrix may possess specific structures that correspond to product correlations, such as clustering or block diagonal structures. The caveat is that now the “covariance matrix” is not full rank, and if there are $m$ alternatives, it is exactly of rank $m-1$. This rank constraint needs to be reflected in the optimization routine. However, rank constraints are generally not convex, and hence hard to incorporate into optimization softwares such as CVX. To address this problem, we propose to encode the constraint in a different way in the E step instead of the M step. More precisely, recall the target $d$-dimensional TMVN density of interest:
where, if alternative $j$ is chosen in the $i$-th observation, we have
with bounds
A key insight of burgette2021symmetric is that their identification constraint that latent utilities sum to zero corresponds exactly to transforming the unnormalized utilities with the matrix $T_s$ given by
It can be shown that with the transformation $Z_i \rightarrow T_sZ_i$, the resulting utilities satisfy the normalization constraint while respecting the original choice decisions. This suggests that when using EP to compute the first two moments in the E step of our algorithm, we can replace $A_i$ by $A_iT_s$. This allows us to encode the level invariance constraint directly in the E step. In the M step, we then only impose a scale invariance constraint, such as on the trace or top diagonal entry of the covariance or precision matrix.
Alternatively, we can impose the constraint that latent utilities sum to zero as an additional truncation constraint, i.e.,
where $l'=u'=\mathbf{0}$.
In addition, we may also consider replacing $\Sigma$ with $T_s\Sigma T_s^T$ in the M step
to enforce the identification constraint. However, $T_s\Sigma T_s^T$ is not invertible, so we need to think about how to implement it.
We consider the simulation setting of natarajan2000monte, where the latent covariance is compound symmetric with off diagonal values of 0.5, and the univariate choice specific coefficient is 1:
The data is generated according to
We consider choice set sizes $m \in \{3, 5, 10\}$. We generate $n = 500$ observations. For identifiability, we use the first choice category as the baseline and infer the covariance matrix $\tilde{\Sigma}$ for the difference in latent utilities, $\tilde{Z}_{ij} = Z_{ij} - Z_{i1},\; j \in \{2, \ldots, m\}$, and we set $\tilde{\Sigma}_{11}=1$. For this problem, the true $\tilde{\Sigma}_{11}$ is
For the M step of the EM algorithm we use the nonlinear constrained optimization routine in the R package NlcOptim, which implements the sequential quadratic programming (SQP) method described in nocedal_numerical_2006. We maintain identifiability by restricting the $1,1$ element of the covariance to 1 in the optimization process. We do not include a penalty term for this experiment, and the estimated covariance matrix is dense.
In the E step we compare the EP method and Gibbs sampling for Monte Carlo estimation, specifically the Gibbs sampler of li2015efficient.
Figure (ref) compares the parameter recovery performance of the EM algorithm with EP and Gibbs sampling approximations to the truncated multivariate normal moments. EP generally leads to better estimates of the parameters.
We have found that nonlinear optimization with respect to $\Sigma$ of the conditional likelihood in the M step is numerically unstable, preventing us from dealing with problems larger than 10 dimensions. Instead, we propose to solve a constrained optimization problem with respect to the precision matrix $\Omega=\Sigma^{-1}$ in the M step, which is convex. More precisely, we impose a trace constraint on $\Omega$ instead of $\Sigma$ to ensure scale normalization:
where the trace constraint $Tr(\Omega)=m$ can be replaced with other normalization constants. This approach differs from that in burgette2012trace, which imposes the constraint on the covariance matrix $\Sigma$, and proposes a prior that enforces this constraint, since they take a Bayesian approach. In contrast, we enforce the constraint directly in the optimization problem in the M step.
Problem (ref) is a trace constrained version of the usual maximum likelihood estimation of covariance matrix, where $S$ takes the role of the empirical covariance matrix. With $\lambda$ the optimal dual variable, the convex optimization problem (ref) has optimal conditions
which can be solved using bisection to find the root $y$ in a subinterval of $y\in(-\infty,\lambda_{1})$ where $\lambda_{1}$ is the smallest eigenvalue of $S$. Alternatively, if we find the eigenvalues $\lambda_{1}\leq\dots\leq\lambda_{m}$ of $S$, the problem becomes $y<\lambda_{1}$, $\sum_{i}\frac{1}{\lambda_{i}-y}=m$, which can be solved using Newton's method. Both methods are guaranteed to find the optimal $\Omega=(S-yI)^{-1}$, and we recommend using a hybrid approach that starts with a bisection procedure to locate a starting point to the right of the optimal $y$, and then proceed using Newton ye1992new. Instead of relying on optimization packages, we implemented a simple customized solver based on the Newton method, which is able to achieve significant speedup over off-the-shelf convex optimization packages such as CVX. Moreover, for large problems $(m\geq50)$, our customized solver also achieves much better solution accuracy for (ref).
Our simulation experiments in this setting have $n=500$ observations, $\beta=1$, and true covariance
The data is generated according to
To ensure identification, we reduce the system by subtracting a reference utility to address level invariance. For scale invariance, we compare two types of normalization approaches. The first is the trace normalization constraint on $\Omega$, while the other is to set the diagonal elements of $\Omega$ to be 1.
Figure (ref) shows that the estimates are converging to the true value in the $m=100$ case. Here there are a few entries that seem to be increasing, but note that we are only using sample size $n=500$. By increasing the sample size as number of alternatives increases, we observe convergence. More detailed experiments to come.
Another observation we had was that the choice of normalization/identification constraint is important for numerical performance. For example, suppose the ground truth $\Omega$ satisfies $Tr(\Omega)=m$ and $\Omega_{11}=1$, in practice the former constraint on the trace performs much better than the latter, because it constrains the entire diagonal. Similarly, if the ground-truth $\Omega$ has diagonal entries equal to $1$, then imposing this constraint results in better convergence.
However Figure (ref) reveals that the estimates diverge after further iterations. This is because we need to invert the precision matrix and provide the covariance to the EP estimation procedure. The covariance values explode, causing numerical problems. After about 18 iterations the procedure fails.
Figure (ref) shows an attempt to resolve this issue by adding a nugget term to the precision to guarantee its invertibility and allow computation of the covariance for all iterations. This does allow computation of the covariance, but the resulting estimates are far off.
Code:
Experiments:
Simulation setting:
The data is generated according to
Longer running simulations show when the algorithm currently converges--when the true coefficient is 0:
When the coefficient is not 0, the precision estimates diverge:
Initializing $\beta$ at different results does not change this behavior. For example when the true $\beta=0$ and we vary the initial value we still get good estimates:
\FloatBarrier Similarly when $\beta=1$ and we vary the initial value we still get poor estimates.
Thought process:
The centering: currently we are computing the mean and covariance of a truncated normal, $Z_i$, that is parameterized with mean $X\beta$. We instead want to compute the mean and covariance of a truncated normal $\epsilon_i$, that has mean 0. Then we can recover the mean of $Z_i$ by undoing the transformation.
We have a choice $Y$, covariate $X$, coefficient $\beta$, choice covariance $\Sigma$. If we observe the choice is the base case $Y = 1$, then the constraints are
If the choice is $Y = j$, $j\neq 1$, then the constraints are
where