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.
116,401 characters · 22 sections · 69 citation commands
Amortized Inference for Correlated Discrete Choice Models via Equivariant Neural Networks
Discrete choice models are a widely used tool in management science, economics, marketing, and other fields for understanding how individuals and organizations make decisions among finite sets of alternatives mcfadden1974conditional. These models provide a structural framework that delivers interpretable parameters---such as willingness-to-pay and demand elasticities---and that also enables predictions about the effects of hypothetical interventions, including pricing changes, product introductions, and policy modifications.
The dominant discrete choice model in applied work is the multinomial logit (MNL), which assumes that the utility consumer $i$ obtains from alternative $j$ takes the form $U_j = v_j + \epsilon_j$ where $v_j$ is the deterministic component of utility, typically written as a linear function $\boldsymbol{x}_j^{\top} \boldsymbol{\beta}$ of alternative $j$'s attributes $\boldsymbol{x}_j$, and $\epsilon_j$ is a Type I extreme value (Gumbel) error that is independent and identically distributed (iid) across alternatives $j = 1, \ldots, K$. This setup yields simple closed-form choice probabilities: the probability of choosing alternative $j$ is given by the “softmax” function $\exp(v_j) / \sum_{k=1}^{K} \exp(v_k)$. The resulting computational convenience has made MNL the default choice in countless applications.
However, this convenience comes at a well-known cost. The iid Gumbel error assumption implies the restrictive independence of irrelevant alternatives (IIA) property: the odds ratio between any two alternatives is unaffected by the presence or attributes of other alternatives, which can lead to unrealistic substitution patterns. A popular generalization of MNL is the mixed logit model (mixed-MNL), in which the utility weight $\boldsymbol{\beta}$ is allowed to be heterogeneous across consumers. Estimation is straightforward via simulation methods train2009discrete. The mixed-MNL relaxes IIA at the aggregate level, but IIA still holds at the level of consumer $i$'s individual choices given the corresponding random utility weight, $\boldsymbol{\beta}_i$.
The multinomial probit (MNP) is a well-known alternative to MNL that relaxes IIA and allows for flexible substitution patterns by assuming $\epsilon_1,\ldots,\epsilon_{K}$ come from a multivariate normal distribution with cross-alternative correlations. Despite this advantage, MNP has seen relatively limited adoption, largely due to its computational demands. MNP choice probabilities have no closed form and require evaluating multivariate normal rectangle probabilities. Estimation relies on more sophisticated simulation methods than needed for mixed-MNL: the GHK simulator in a classical setting geweke1989bayesian, hajivassiliou1990simulation, keane1994solution or MCMC in a Bayesian setting McCullochRossi1994.
We propose a fundamentally different approach: Rather than simulating choice probabilities anew for each likelihood evaluation, we train a neural network emulator to directly approximate the choice probability function. This strategy---known as amortized inference---shifts the computational burden from inference time to a one-time training phase. Once trained, the emulator provides rapid deterministic approximations to choice probabilities via simple function calls. The amortized inference framework has proven highly successful in multiple scientific domains for approximating computationally intensive simulation models lueckmann2019likelihood,cranmer2020frontier. Here, we adapt it to discrete choice models, extending beyond MNL and MNP to general correlated error distributions.
In the econometrics literature, Norets2012 employs a related strategy, approximating the expected value function in dynamic discrete choice models via a neural network. He applies this strategy in the context of a specific parametric model, so changes in the dynamic model require modifying and retraining the neural network. In contrast, our emulator operates on the generic choice problem given only the deterministic utilities and the parameters of the error distribution, enabling modifications to the parametric model (e.g., the functional form of the deterministic utilities) without altering or retraining the emulator.
Our methodology relies on several new contributions. First, we develop a neural network architecture specifically designed for choice models. It respects their fundamental invariance properties, which we formalize as symmetries (group actions) of the choice-probability map.\footnote{If these properties are not embedded in the architecture, the NN must learn them during training. As a continuum of models generate equivalent choice probabilities, this would slow learning considerably.} Choice probabilities are invariant to location shifts (adding a constant to all utilities) and scale transformations (multiplying utilities by a positive constant), and equivariant with respect to permutations of choice alternatives. Our architecture incorporates a preprocessing transformation that reduces the size of the feature space by enforcing location and scale invariance. The processed features are then passed to a per-alternative encoder module based on the DeepSet architecture zaheer2017deep. The output is then concatenated across alternatives and processed through equivariant layers that impose a sum-to-one constraint on the output probabilities while respecting the permutation equivariance property.
Second, we establish the theoretical foundations of our architecture. We prove that it can universally approximate choice probabilities on compact subsets of the parameter space, outside a measure-zero exceptional set. This result connects our architecture to the theory of orbit separation under group actions, extending the universal approximation results of blumsmith2025machine on symmetric matrices to the joint space of utility vectors and centered covariance matrices. blumsmith2025machine employ Galois theory to prove generic orbit separation, while we provide a direct proof based on an invariant reconstruction argument.
Third, we establish the statistical properties of maximum likelihood estimators (MLEs) formed using an emulator approximation of the true discrete choice probabilities. We show that if the emulator approximates the true log likelihood sufficiently well---specifically, if the average approximation error is $o_p(n^{-1})$---then the emulator-based estimator inherits the consistency and asymptotic normality of the exact MLE. When this condition is not met, we show that valid inference can still be obtained via sandwich standard errors, treating the emulator as a working model in the quasi-maximum likelihood framework.
By careful design, our preprocessing transformations and architecture enforce smoothness of choice probabilities with respect to the inputs: the deterministic utilities and corresponding scale (covariance) matrix. We use Sobolev training czarnecki2017sobolev to match emulator gradients to those of the log choice probabilities. Together, these design choices enable reliable use of automatic differentiation for downstream model-fitting and inference tasks.
With a pretrained emulator, generalizing from logit to probit (or any other error distribution) requires only the replacement of closed-form softmax probabilities with emulator evaluations of choice probabilities when forming the likelihood. In the probit case, given a fixed computational budget, we show via simulations that an ML estimator using our amortized inference procedure matches or exceeds the performance of one using the GHK algorithm in terms of estimation error and coverage rates. Furthermore, our approach can easily handle models other than MNP where an efficient simulation algorithm like GHK is unavailable.
For concreteness, much of our exposition focuses on MNP models. However, the emulator approach and supporting theory are largely agnostic to the assumed parametric form of the errors. Generalizing to other error distributions (e.g., correlated Gumbel) is straightforward, requiring only modest changes to the training data generation and emulator inputs.
The remainder of this paper is organized as follows. Section (ref) summarizes related literature on choice modeling. Section (ref) describes the MNP model and our inferential goals. Section (ref) presents the emulator architecture and training procedure. Section (ref) establishes the theoretical properties of the architecture and emulator-based estimators. Section (ref) presents simulation results, and Section (ref) concludes.
The literature on applying machine learning methods to discrete choice can be divided into two broad streams. The first stream maintains the classic random utility model (RUM) of choice, in which choice probabilities are generated by a population of rational consumer types with different (but stable) preference orderings over the universe of choice objects, as explained by BlochMarschak1960 and McfaddenRichter1990.\footnote{Suppose there are $\bar{K}$ objects in the universe and a consumer is presented with a choice set (or “assortment”) that contains $K \le \bar{K}$ elements. A key implication of RUM is that the utility a consumer derives from product $k$ is invariant to the choice set, ruling out context or assortment effects.} MNL, MNP, and mixed-MNL are all members of the RUM class provided one maintains the RUM assumption that the utility of an option depends only on its own characteristics. Key papers in this strand start from this basic MNL structure and use neural networks to generalize the functional form of utility: Bentz2000, Sifringer2020, WangMoZhao2020, han2022, and singh2023choice. Two important recent papers extend this work: Aouad2025 develop an architecture that implements the mixed logit model with a flexible distribution of taste heterogeneity (RUMnet), and Bagheri2025 develop another architecture that generalizes the Gumbel error assumption (RUM-NN). Both papers use softmax-smoothed sample averages to approximate choice probabilities within the loss function, resulting in increased computation time relative to pure logit-based models.
Aouad2025 build on the generic approximation property of mixed-MNL models shown in mcfadden2000mixed. This result relies on a flexible basis expansion of the utility function and a potentially nonparametric mixing distribution. In practice, these model attributes are unknown to analysts. As emphasized in both Aouad2025 and mcfadden2000mixed, this approximation property is not unique to Gumbel errors. Rather, it holds more generally for a larger class of hierarchical discrete choice models, including mixed-MNP models. mcfadden2000mixed use this approximation property to justify adoption of mixed-MNL as a computationally convenient alternative to other models lacking closed-form choice probabilities, such as MNP. But this justification of mixed-MNL becomes less compelling given a practical, general-purpose alternative like the emulators we propose here.
The second stream dispenses with the RUM structure, often motivated by a desire to relax restrictive MNL assumptions like IIA and to allow more flexible substitution patterns. Some papers in this stream view discrete choice as a general classification problem that is amenable to machine learning methods. This is exemplified by WangRoss2018, Lheritier2019, Rosenfeld2020, ChenMisic2022, and Chen2025. Others maintain an MNL structure at the top level---that is, choice probabilities are determined by a vector of alternative-specific utilities that enter a softmax function---but a neural net is used to construct the alternative-specific utility functions in flexible ways that deviate from RUM assumptions. For instance, the utility of alternative $j$ is allowed to depend on attributes of other alternatives to generate context effects. This is exemplified by WangMoZhao2021, Wong2021, Cai2022, Pfannschmidt2022, and Berbeglia2025.
Between these two streams, a fundamental tension emerges: While RUM models are favored for their interpretability and grounding in economic theory, the dominant RUM models in empirical work (linear MNL and mixed-MNL) enforce constraints on deterministic utilities and substitution patterns that may be too restrictive in some applications; in particular, they assume IIA at the level of individual choices. The first stream aims to address this challenge by adding additional flexibility to RUM models. However, many of these proposals maintain the assumption of independent logit errors, constraining substitution patterns and potentially biasing utility estimates. The exceptions---namely, RUMnet and RUM-NN---address substitution, but they do so at the expense of interpretability and/or computational efficiency due to their (a) (potentially) nonlinear utility functions and (b) expensive sample-average approximations of the likelihood function. The second stream abandons the RUM structure altogether, making it difficult or impossible to obtain reliable inferences for many economically meaningful quantities, such as consumer welfare, willingness-to-pay measures, demand elasticities, and substitution effects. In summary, this tension results in a tradeoff between modeling flexibility and interpretation, ranging from rigid interpretable models like MNL on one hand to flexible atheoretic classification methods on the other.
MNP models constitute a notable exception to the above tradeoff as they allow for flexible substitution patterns via interpretable covariance relationships. Moreover, if additional flexibility is desired, they can be extended to allow the deterministic utilities to follow the functional form of a neural network Hruschka2007. Conversely, if researchers desire a more parsimonious or interpretable model, analysts can constrain the covariance structure via penalization methods as in Jiang2025 or via factor structures as we illustrate in Section (ref). Despite these virtues, empirical applications of MNP are relatively sparse in the literature due, in large part, to the difficulty of evaluating MNP choice probabilities.
More broadly, MNP models are just one member of a larger class of RUM models featuring correlated error terms. Analogous to MNP, other members can be generated by assuming errors of the form $\boldsymbol{\epsilon} = \boldsymbol{\mathbf{\Sigma}}^{1/2} \boldsymbol{\epsilon}^*$, where $\boldsymbol{\mathbf{\Sigma}}^{1/2}$ is a square-root scale matrix and $\boldsymbol{\epsilon}^*$ is a vector of exchangeable errors. For example, $\epsilon_j^* \overset{iid}{\sim} \text{Gumbel}(0, 1)$ results in a correlated Gumbel distribution with scale matrix $\boldsymbol{\mathbf{\Sigma}}$. As another example, we could assume $\boldsymbol{\epsilon}^* \sim \text{Multivariate-}t(\boldsymbol{0}, \boldsymbol{\mathbf{I}}, \nu)$ to capture heavy tail behavior via the degrees-of-freedom parameter $\nu > 0$ ($\boldsymbol{\mathbf{I}}$ denotes an identity matrix). Bagheri2025 consider correlated error distributions formed in this manner, resorting to expensive simulation methods to approximate the choice probabilities because an efficient algorithm (GHK) is available only in the Gaussian case.
The primary contribution of this work is an amortized inference framework that produces accurate, reusable, and computationally efficient estimates of choice probabilities for RUM models with general error distributions, including those featuring nontrivial correlation structures. The framework presents a potential resolution to the flexibility--interpretability tradeoff highlighted above, enabling practical estimation of flexible choice models without sacrificing the economically meaningful insights and parsimony of RUM models. The framework is supported by strong theoretical justification, including a universal approximation guarantee and asymptotic inference results under mild approximation conditions (see Section (ref)). Moreover, the framework is largely complementary to recent advances in machine learning methods, including the two streams discussed above. In many cases, these methods could be enhanced with emulator-based likelihood evaluations, resulting in flexible models with interpretable substitution patterns and manageable computational demands.
This section summarizes the problem setup, including the invariance properties of discrete choice models and our inferential goals.
We consider a decision-maker choosing among $K$ mutually exclusive alternatives. The decision-maker assigns a latent utility $U_j$ to each alternative $j \in \{1, \ldots, K\}$ and selects the alternative with the highest utility:
The latent utilities decompose into deterministic and stochastic components:
where $v_j$ is the deterministic (or systematic) utility that depends on observable characteristics, and $\epsilon_j$ is a centered random error capturing unobserved factors. As explained in Section (ref), the deterministic utility typically takes a linear form $v_j = \boldsymbol{x}_j^\top \boldsymbol{\beta}$, where $\boldsymbol{x}_j$ is a vector of alternative-specific attributes and $\boldsymbol{\beta}$ is a parameter vector to be estimated.
Different distributional assumptions on the error vector $\boldsymbol{\epsilon} = (\epsilon_1, \ldots, \epsilon_K)^\top$ yield different choice models. The MNL model assumes that the $\epsilon_j$ are independent and identically distributed with the $\text{Gumbel}(0, 1)$ distribution. In contrast, the MNP model assumes that $\boldsymbol{\epsilon}$ follows a multivariate normal distribution with mean zero and covariance matrix $\boldsymbol{\mathbf{\Sigma}}$. Below, we consider correlated models in which $\boldsymbol{\epsilon} = \boldsymbol{\mathbf{\Sigma}}^{1/2} \boldsymbol{\epsilon}^*$, where $\boldsymbol{\mathbf{\Sigma}}^{1/2}$ is a square-root scale matrix as described in Section (ref). We impose the following assumption on $\boldsymbol{\epsilon}^*$.
This setup and assumption include the MNP as a leading example with $f_{\boldsymbol{\epsilon}^*}$ equal to the PDF of the $\mathcal{N}(\boldsymbol{0}_{K}, \boldsymbol{\mathbf{I}}_{K})$ distribution. However, the setup is general enough to include many other possibilities, such as the correlated Gumbel and multivariate-$t$ distributions described in Section (ref). We defer treatment of additional regularity conditions until Section (ref) when we develop the theoretical properties of our framework in the context of an assumed parametric model with parameter vector $\boldsymbol{\theta} \in \Theta \subset \mathbb{R}^p$.
We denote the choice probability for alternative $j$ as
where $\boldsymbol{v} = (v_1, \ldots, v_K)^\top$. In practice, we may also choose to include additional parameters governing the distribution of $\boldsymbol{\epsilon}^*$, such as the degrees-of-freedom parameter, $\nu$, described in Section (ref), but these parameters are suppressed in the notation for simplicity. The choice probability can be expressed as the following integral:
where $\mathcal{R}_j$ is the region of the error space where alternative $j$ is chosen and $\big[\boldsymbol{\mathbf{\Sigma}}^{1/2}\big]_{j\mathpalette\bigcdot@{.5}}$ denotes the $j$th row of $\boldsymbol{\mathbf{\Sigma}}^{1/2}$ so that $\big[\boldsymbol{\mathbf{\Sigma}}^{1/2}\big]_{j\mathpalette\bigcdot@{.5}} \boldsymbol{\epsilon}^* = \epsilon_j$. In general, this integral is not analytically tractable. In MNP models, in particular, it does not have a closed-form solution for $K \geq 3$, necessitating numerical methods, such as the GHK algorithm, which simulates these probabilities via recursive conditioning. For each evaluation, GHK requires $R$ simulation draws, with the approximation error decreasing at rate $O(R^{-1/2})$.
Choice probabilities satisfy three fundamental invariance properties.
Location invariance. Adding a constant to all utilities does not change choice probabilities:
where $\boldsymbol{1}_K$ is the $K$-vector of ones and $C \in \mathbb{R}$ may be random. This follows because choices depend only on utility differences.
Scale invariance. Scaling all utilities and the scale matrix preserves choice probabilities:
This follows because the scaled model is equivalent to the original model with $\boldsymbol{v}$ and $\boldsymbol{\epsilon}$ replaced by $\alpha \boldsymbol{v}$ and $\alpha \boldsymbol{\epsilon}$, respectively.
Permutation equivariance. Relabeling alternatives permutes the choice probabilities correspondingly. For any permutation $\pi$ of $\{1, \ldots, K\}$ with associated permutation matrix $\boldsymbol{\mathbf{P}}_\pi$:
A naive attempt at constructing an emulator may not respect these properties. Translating, rescaling, or permuting the alternatives could produce different emulator predictions, resulting in slow learning and poor generalization because the emulator would need to learn these properties manually. In Section (ref), we propose a solution that automatically satisfies these properties based on two design components: a preprocessing step that fixes the location and scale, and an equivariant architecture that ensures permutation equivariance without sacrificing expressivity.
The above invariance properties imply that the choice model is not identified without normalization. We assume throughout that the model has been appropriately normalized so that the parameter vector $\boldsymbol{\theta}$ is identified. This is typically achieved in practice by setting the utility of one alternative to zero deterministically and fixing one of the nonzero diagonal elements in $\boldsymbol{\mathbf{\Sigma}}$ (see keane1992). Similarly, we assume that the data exhibit sufficient variation to identify the remaining parameters, which typically requires variation in alternative-specific covariates, such as prices.
Given observations $\{(y_i, \boldsymbol{\mathbf{X}}_i)\}_{i=1}^n$, where $y_i \in \{1, \ldots, K\}$ is the observed choice and $\boldsymbol{\mathbf{X}}_i$ represents the covariates, we seek to estimate the parameter vector, $\boldsymbol{\theta}$. This parameter vector includes both the coefficients affecting deterministic utilities and additional parameters governing the scale matrix.
The (averaged) log-likelihood function is
where $\boldsymbol{v}_i(\boldsymbol{\mathbf{X}}_i, \boldsymbol{\theta})$ and $\boldsymbol{\mathbf{\Sigma}}_i(\boldsymbol{\mathbf{X}}_i, \boldsymbol{\theta})$ denote the utility vector and scale matrix for observation $i$ as functions of the covariates and parameters. Maximum likelihood estimation requires maximizing $\ell_n(\boldsymbol{\theta})$ over $\Theta$, which in turn requires repeated evaluation of $\ell_n(\boldsymbol{\theta})$ and its gradient $\nabla_{\boldsymbol{\theta}} \log P_j(\boldsymbol{v}, \boldsymbol{\mathbf{\Sigma}})$. Similarly, Bayesian methods such as Hamiltonian Monte Carlo require repeated evaluations of both $\ell_n(\boldsymbol{\theta})$ and $\nabla_{\boldsymbol{\theta}} \log P_j(\boldsymbol{v}, \boldsymbol{\mathbf{\Sigma}})$, often thousands of times throughout the estimation process.
Our goal is to construct a neural network emulator $\hat{P}_j(\boldsymbol{v}, \boldsymbol{\mathbf{\Sigma}}; \boldsymbol{\phi})$ that provides fast, accurate approximations of the true choice probabilities, along with analytic gradients via automatic differentiation. The emulator is trained once on simulated data spanning the relevant input space, after which it can be used for rapid inference on many datasets and models.
This section outlines the emulator preprocessing transformations, neural network architecture, training procedure, and inference process.
The emulator must satisfy four desiderata:
We achieve these goals through a combination of preprocessing transformations and a carefully designed neural network architecture.
We address location and scale invariance through preprocessing transformations that project the inputs onto a canonical subspace.
Centering. Define the centering matrix
which we employ to project the vector of utilities, $\boldsymbol{U} \in \mathbb{R}^{K}$, onto the subspace orthogonal to $\boldsymbol{1}_K$, resulting in transformed utilities $\tilde{\boldsymbol{U}} = \boldsymbol{\mathbf{M}} \left(\boldsymbol{v} + \boldsymbol{\epsilon}\right) = \boldsymbol{v} + \boldsymbol{\epsilon} - \boldsymbol{1}_{K} (\bar{v} + \bar{\epsilon})$, where $\bar{v} = \frac{1}{K} \sum_{j=1}^{K} v_j$ and $\bar{\epsilon} = \frac{1}{K} \sum_{j=1}^{K} \epsilon_j$. These transformed utilities sum to zero within each realization: $\sum_{j=1}^{K} \tilde{U}_j = 0$. In effect, this transformation results in modified emulator inputs:
Scaling. We normalize by the trace of the centered scale matrix, producing transformed utilities $\boldsymbol{U}^* = \sqrt{K / \mathrm{tr}\big(\tilde{\boldsymbol{\mathbf{\boldsymbol{\mathbf{\Sigma}}}}}\big)} \tilde{\boldsymbol{U}}$. The final transformed parameters are then
We can omit the case $\mathrm{tr}(\tilde{\boldsymbol{\mathbf{\Sigma}}}) = 0$ as this results in known, deterministic choices. The complete transformation produces normalized inputs $(\boldsymbol{v}^*, \boldsymbol{\mathbf{\Sigma}}^*)$ satisfying:
This normalization reduces the size of the emulator's input space, accelerating learning and improving generalization. Importantly, this transformation commutes with permutation matrices: for any permutation matrix $\boldsymbol{\mathbf{P}}_\pi$, applying the transformation to $(\boldsymbol{\mathbf{P}}_\pi \boldsymbol{v}, \boldsymbol{\mathbf{P}}_\pi \boldsymbol{\mathbf{\Sigma}} \boldsymbol{\mathbf{P}}_\pi^\top)$ yields $(\boldsymbol{\mathbf{P}}_\pi \boldsymbol{v}^*, \boldsymbol{\mathbf{P}}_\pi \boldsymbol{\mathbf{\Sigma}}^* \boldsymbol{\mathbf{P}}_\pi^\top)$. This ensures that the preprocessing preserves permutation equivariance.
After preprocessing, we construct a neural network that maps $(\boldsymbol{v}^*, \boldsymbol{\mathbf{\Sigma}}^*)$ to choice probabilities while maintaining permutation equivariance. The architecture consists of three components: a per-alternative encoder, permutation-equivariant layers, and an output layer. Below, we interpret $\boldsymbol{\mathbf{\Sigma}}$ as a covariance matrix to simplify the exposition, referring to its diagonal and off-diagonal elements as variances and covariances, respectively. This interpretation is correct for MNP models. For other models, it is only illustrative as these parameters may not exactly coincide with the variances and covariances.
Per-alternative encoder. For each alternative $j \in \{1, \ldots, K\}$, we construct a representation $\boldsymbol{z}_j \in \mathbb{R}^{d_z}$ that captures how alternative $j$ relates to all other alternatives. This representation is built from two complementary DeepSet networks---a diagonal DeepSet and an off-diagonal DeepSet---whose outputs are combined with alternative $j$'s own features. The diagonal DeepSet processes pairwise relationships between $j$ and each other alternative, while the off-diagonal DeepSet summarizes the covariance structure among alternatives other than $j$. In the descriptions below, the features labeled as “base inputs” are sufficient for universal approximation of choice probabilities (see Section (ref)). We include additional features to improve the expressivity of the emulator, as we find this improves its predictions.
Diagonal DeepSet. For each pair $(j, k)$ with $k \neq j$, we construct a feature vector $\boldsymbol{d}_{jk}$ containing differentiable functions of the utilities and covariances associated with alternatives $j$ and $k$. The base inputs are the utilities $v_j^*$ and $v_k^*$, the variances $\Sigma_{jj}^*$ and $\Sigma_{kk}^*$, and the covariance $\Sigma_{jk}^*$. From these, we derive additional feature transformations that facilitate learning, including standard deviations $\sigma_j = \sqrt{\Sigma_{jj}^*}$ and $\sigma_k = \sqrt{\Sigma_{kk}^*}$, the correlation $\rho_{jk} = \Sigma_{jk}^*/(\sigma_j \sigma_k)$, and the standardized utility difference (z-score):
The diagonal DeepSet processes these features using the architecture of zaheer2017deep:
where $\phi_{\text{diag}}$ is a multi-layer perceptron (MLP) applied to each pair, the sum aggregates over all $k \neq j$, and $\rho_{\text{diag}}$ is another MLP that produces a learned nonlinear representation of the sum. We use the notation $\phi_{\text{diag}}$ and $\rho_{\text{diag}}$ following zaheer2017deep; the MLPs $\rho_{\text{diag}}$ and $\rho_{\text{off}}$ (defined below) can be distinguished from the correlation $\rho_{jk}$ by their subscripts.
singh2023choice also employ a variant of the DeepSet architecture. However, their method assumes iid errors and operates on problem-specific features, which requires their neural network to be modified and retrained to accommodate new data structures. In contrast, our architecture relies on generic RUM model inputs, such as deterministic utilities.
Off-diagonal DeepSet. For each pair $(k, l)$ with $k < l$ and $k, l \neq j$, we construct a feature vector $\boldsymbol{o}_{kl}$ containing differentiable, symmetric functions comparing alternatives $k$ and $l$; we require symmetry because the pairs do not possess a natural ordering. The base input is $\Sigma_{kl}^*$, and we include other derived features, such as the correlation $\rho_{kl}$, squared utility difference $(v_k^* - v_l^*)^2$, the squared z-score $z_{kl}^2$, and the variance sum $\Sigma_{kk}^* + \Sigma_{ll}^*$.
The off-diagonal DeepSet has an analogous structure to that of the diagonal DeepSet:
Combining network. The outputs of both DeepSets are concatenated with additional “pass-through features” $\boldsymbol{s}_j$, such as $v_j^*$, $\sigma_j$, and summary statistics of alternative $j$'s covariances with other alternatives. As with the additional features above, the pass-through features are not required for universal approximation, but we find that they improve performance empirically. The concatenated result is then processed by the combining MLP, $\zeta$, to produce the per-alternative representation:
Permutaition-equivariant layers. The per-alternative representations are stacked into a matrix $\boldsymbol{\mathbf{Z}} \in \mathbb{R}^{K \times d_z}$ and processed through linear permutation-equivariant layers zaheer2017deep, allowing information exchange across alternatives. These layers take the form
where $\boldsymbol{\mathbf{A}}, \boldsymbol{\mathbf{B}} \in \mathbb{R}^{d_{\text{in}} \times d_{\text{out}}}$ and $\boldsymbol{c} \in \mathbb{R}^{d_{\text{out}}}$ are learnable parameters, and $\sigma$ is a smooth activation function. The first term applies a per-alternative transformation, the second aggregates information across alternatives, and the third adds a shared bias (intercept).
Output layer. The final layer produces one logit per alternative:
where $\boldsymbol{\mathbf{Z}}^{(L)}$ is the output of the last equivariant layer and $\boldsymbol{a} \in \mathbb{R}^{d_L}$ so that the output $\textbf{logit}$ is a vector in $\mathbb{R}^{K}$. We omit the other terms because they produce a constant shift that does not affect the output probabilities.
Choice probabilities are obtained by applying the softmax function to the output layer:
For any vector of valid choice probabilities $(P_1, \ldots, P_{K})$ on the probability simplex $\Delta^{K-1}$, there exists a corresponding vector of logits yielding $(P_1, \ldots, P_{K})$ via softmax, so this parameterization entails no loss of generality. Figure (ref) illustrates the overall architecture.
Activation functions. Throughout the architecture, we use smooth (infinitely differentiable) activation functions. This ensures that the emulator output is a smooth function of the inputs, enabling accurate gradient computation for downstream model fitting. We found that the Swish activation function $\text{Swish}(x) = x / (1 + e^{-x})$ works well empirically ramachandran2018searching, though other smooth activation functions could also be applied, such as $\tanh$ or $\operatorname{softplus}$ dugas2000softplus.
Computational complexity. A naive implementation of the per-alternative encoder would require $O(K^3)$ forward-pass operations due to the off-diagonal DeepSet computation. However, by precomputing the sum over all off-diagonal pairs and using subtraction to obtain the sum excluding pairs involving alternative $j$, we achieve $O(K^2)$ complexity. Because $\boldsymbol{\mathbf{\Sigma}}$ has $O(K^2)$ unique elements influencing the choice probabilities, this scaling is optimal for general $\boldsymbol{\mathbf{\Sigma}}$.
Network size. The architecture can be compact. In our simulations, we use networks with 8--64 hidden units and 0--2 hidden layers per component. In general, larger values of $K$ require larger networks for accurate approximation.
Relation to blumsmith2025machine. Our proposed architecture can be seen as an extension of the DS-CI architecture of blumsmith2025machine with three key differences. First, whereas blumsmith2025machine require a “linking feature” combining diagonal and off-diagonal elements, our centering procedure makes this feature redundant, simplifying the set of required inputs. Second, our problem includes an additional vector input ($\boldsymbol{v}$) not considered by blumsmith2025machine; this additional complexity requires us to couple deterministic utility values with the diagonal variance elements. Third, the equivariant nature of our problem requires us to encode each alternative separately, prompting the separation into a single distinguished alternative ($j$) and rival alternatives ($k \neq j$) in the per-alternative encoder.
We train the emulator using simulated data spanning a diverse range of utility and covariance configurations, employing Sobolev training to ensure accurate approximation of both choice probabilities and their gradients.
Data generation. To generate each training example for our simulations in Section (ref) involving the MNP model, we follow the procedure below:
Steps 3 and 4 ensure diversity in the relationship between deterministic utilities and error covariances. When $\omega = 1$, the utilities are perfectly aligned with the covariance structure; when $\omega = 0$, they are independent. Intermediate values produce partial correlation.
This procedure can be generalized to other error distributions by replacing the MNP model in the final step with another model with a different error distribution. Similarly, the other distributions listed and the number of choices ($10^6$) can be modified based on the problem domain and computational resources.
Loss function. We employ Sobolev training czarnecki2017sobolev, which augments the standard likelihood loss with a gradient-matching penalty. The total loss is:
where $\lambda_{\text{grad}} > 0$ controls the weight of the gradient-matching term.
The cross-entropy component is the multinomial negative log-likelihood (up to constants):
where $\hat{P}_j^{\text{sim}} = n_j / \sum_{k=1}^{K} n_k$ is the simulated choice frequency for alternative $j$.
For the gradient-matching component, we follow the stochastic Sobolev training approach of czarnecki2017sobolev, which avoids computing full Jacobian matrices by instead matching directional derivatives along random directions. Let $\boldsymbol{d} = (\boldsymbol{d}_v, \operatorname{vech}(\boldsymbol{\mathbf{D}}_\Sigma))$ denote a random direction in the joint space of utilities and covariances, where $\boldsymbol{d}_v \in \mathbb{R}^K$, $\boldsymbol{\mathbf{D}}_\Sigma \in \mathbb{R}^{K \times K}$ is symmetric, and $\operatorname{vech}(\cdot)$ denotes the half-vectorization (upper triangular elements including the diagonal). The gradient-matching loss is:
where $\nabla_{\boldsymbol{d}}$ denotes the directional derivative along $\boldsymbol{d}$. Weighting by the simulated choice frequencies $\hat{P}_j^{\text{sim}}$ focuses learning on alternatives that are frequently chosen, whose gradients most influence parameter estimation in downstream inference tasks.
Constraint-respecting directions. The canonical inputs $(\boldsymbol{v}^*, \boldsymbol{\mathbf{\Sigma}}^*)$ lie on a constrained manifold:
Random directions must respect these constraints to remain in the tangent space of $\mathcal{X}_K$. We construct valid directions as follows:
Computing directional derivatives. We compute the emulator's directional derivatives via finite differences:
where $\epsilon > 0$ is a small step size (we use $\epsilon = 10^{-5}$). This approach is computationally efficient, requiring only two additional forward passes per direction regardless of the input dimension. Alternatively, one could compute the directional derivative via automatic differentiation as in czarnecki2017sobolev.
Target directional derivatives are computed from pre-stored Jacobians. During data generation, we compute and store the full Jacobians $\nabla_{\boldsymbol{\theta}} \log P_j^{\text{target}}$ with respect to $\boldsymbol{v}^*$ and $\boldsymbol{\mathbf{\Sigma}}^*$. During training, the target directional derivative is obtained via the inner product:
where $\langle \cdot, \cdot \rangle_F$ denotes the Frobenius inner product.
Target gradient computation. Computing target gradients for Sobolev training requires differentiating through the choice simulation process. Since discrete choices are inherently non-differentiable, we employ a smooth relaxation using a temperature-scaled softmax. For each Monte Carlo draw $r$ with utility vector $\boldsymbol{U}^{(r)} = \boldsymbol{v}^* + \boldsymbol{\epsilon}^{(r)}$, we approximate the hard choice indicator with:
where $\tau > 0$ is a temperature parameter. As $\tau \to \infty$, this converges to the hard argmax; for finite $\tau$, it provides a differentiable approximation. The target log-probability is then:
and the target gradients $\partial \log P_j^{\text{target}} / \partial \boldsymbol{v}^*$ and $\partial \log P_j^{\text{target}} / \partial \boldsymbol{\mathbf{\Sigma}}^*$ are obtained via automatic differentiation through this expression. This approach is closely related to the “Gumbel--softmax trick” frequently used for gradient estimation with neural networks jang2017categorical, maddison2017concrete. We use $\tau = 3$ in our experiments, which provides a balance between gradient accuracy (low bias) and computational stability.
Training details. Our training procedure makes use of two replay buffers, a common technique in deep learning Mnih2015. The first replay buffer stores precomputed training examples including inputs $(\boldsymbol{v}^*, \boldsymbol{\mathbf{\Sigma}}^*)$, simulated choice frequencies, and target Jacobians. The second stores random directions for the gradient-matching loss. At each training step, we sample random pairs from the replay buffers, compute the loss, and update the network parameters via stochastic gradient descent. The replay buffer accelerates training by allowing reuse of expensive high-precision simulations and gradient computations. We set the gradient penalty weight $\lambda_{\text{grad}} = 10^{-8}$ in our experiments; empirically, we find that this value provides sufficiently well-behaved derivatives.
Multi-$K$ training. Remarkably, none of the weight vectors and matrices has a dimension depending on $K$. The dimension $K$ enters the neural network only through the input data dimension; thus, we can design a single emulator capable of processing multiple values of $K$ simultaneously. The training procedures described above can then be adapted to iterate over a range of $K$ values, say $K \in \{3, 4, 5\}$, and the resulting emulator can then learn to approximate choice probabilities for a varying number of choices (see Section (ref)). Under this configuration, we recommend including $K$ as an element of $\boldsymbol{s}_j$.
Once trained, the emulator enables rapid inference. Given a new dataset $\{(y_i, \boldsymbol{\mathbf{X}}_i)\}_{i=1}^n$ and candidate parameter vector $\boldsymbol{\theta}$:
As the preprocessing transformation and neural network are both smooth functions, the gradients of the emulator log likelihood with respect to $\boldsymbol{\theta}$ can be computed exactly via automatic differentiation. This enables the use of gradient-based optimizers for maximum likelihood estimation, as well as gradient-based sampling methods such as Hamiltonian Monte Carlo for Bayesian inference. The smoothness of the emulator also facilitates computation of standard errors, which can be estimated via outer products of gradients (see Section (ref)).
This section establishes the theoretical foundations for the emulator architecture and the statistical properties of emulator-based discrete choice estimators.
We first show that the per-alternative encoder architecture can universally approximate choice probabilities. The key insight is that the multisets processed by the encoder generically determine the input up to the relevant symmetry group. For each alternative $j$, define the invariant function $g_j : \mathcal{X}_K \to (\text{Multiset}, \text{Multiset})$ by:
The multiset $\mathcal{T}_j$ consists of $K - 1$ triples, each coupling an alternative $k$'s utility and variance with its covariance with alternative $j$. The multiset $\mathcal{O}_j$ consists of the off-diagonal covariances among alternatives other than $j$. These are precisely the base inputs processed by the per-alternative encoder's diagonal and off-diagonal DeepSets.
The proof, given in Appendix (ref), proceeds in three steps:
Step 1: Genericity conditions. We define two genericity conditions on $(\boldsymbol{v}^*, \boldsymbol{\mathbf{\Sigma}}^*)$: (a) the triples $(v_k^*, \Sigma_{kk}^*, \Sigma_{jk}^*)$ for $k \neq j$ are pairwise distinct as elements of $\mathbb{R}^3$; and (b) for each $k \neq j$, there exists a unique subset of off-diagonal covariances that sums to $-\Sigma_{kk}^* - \Sigma_{jk}^*$. Condition (b) follows from the row-sum constraint $\boldsymbol{\mathbf{\Sigma}}^* \boldsymbol{1}_K = \boldsymbol{0}$: the off-diagonal entries in row $k$ (excluding $\Sigma_{ki}$) sum to $-\Sigma_{kk}^* - \Sigma_{jk}^*$.
Step 2: Recovery procedure. When both conditions hold, we can recover $(\boldsymbol{v}^*, \boldsymbol{\mathbf{\Sigma}}^*)$ from $g_j$ up to relabeling of alternatives other than $j$. First, order the triples in $\mathcal{T}_j$ lexicographically to assign labels to alternatives. Then, for each alternative $k$, identify its row of off-diagonal entries as the unique subset in $\mathcal{O}_j$ with the required sum. Finally, recover shared entries $\Sigma_{kl}^*$ as the singleton intersection of the subsets for rows $k$ and $l$.
Step 3: Measure-zero exceptional set. The set $\mathcal{B}_j$ where either condition fails is a finite union of zero sets of nontrivial analytic functions on the affine subspace $\mathcal{X}_K$. By the principle that proper analytic subvarieties have measure zero Mityagin2020, $\mathcal{B}_j$ has measure zero.
This separation result, combined with the Stone--Weierstrass theorem, yields universal approximation:
The exceptional set $\mathcal{B}$ has measure zero, so it does not affect approximation in practice: for any probability distribution over $\mathcal{X}_K$ that is absolutely continuous with respect to Lebesgue measure, inputs lie in $\mathcal{B}$ with probability zero.
Corollary (ref) guarantees that the per-alternative encoder alone is sufficient for universal approximation. The inclusion of the final equivariant layers ensures that the estimated probabilities sum to one and improves the expressivity of the neural network.
For large choice sets, we may opt for a simpler architecture that applies the per-alternative encoder to only the chosen alternative with no equivariant layers, resulting in a single output, $P_j$. This modification maintains the universal approximation property; however, it relaxes the sum-to-one constraint in return for computational efficiency.
Corollary (ref) could also be extended to hold simultaneously for a finite set of $K$ values under a suitable extension of the space $\mathcal{X}_K$ and its exceptional sets. Section (ref) presents numerical results showing that the proposed architecture can provide accurate approximations for multiple values of $K$ simultaneously.
We now establish that emulator-based estimators inherit the asymptotic properties of exact maximum likelihood estimators under appropriate conditions.
Let $\{(y_i, \boldsymbol{\mathbf{X}}_i)\}_{i=1}^n$ be iid observations from the true model with parameter $\boldsymbol{\theta}_0 \in \Theta$. Define the true and emulator log-likelihood contributions as
with corresponding score functions
The true and emulator (averaged) log likelihoods are $\ell_n(\boldsymbol{\theta}) = \frac{1}{n} \sum_{i=1}^n m(y_i, \boldsymbol{\mathbf{X}}_i; \boldsymbol{\theta})$ and $\hat{\ell}_n(\boldsymbol{\theta}) = \frac{1}{n} \sum_{i=1}^n \hat{m}_n(y_i, \boldsymbol{\mathbf{X}}_i; \boldsymbol{\theta})$. Define the MLE and emulator-based estimator as
We require that the emulator log likelihood converges to the true log likelihood.
The $n$ referenced in Assumption (ref) is the sample size of the observed data. Consequently, this assumption involves an asymptotic regime in which the emulator becomes increasingly accurate as $n \to \infty$, which is achievable by progressively (a) increasing the number of simulated training examples and (b) increasing the neural network complexity (i.e., the number of layers and hidden units). Appendix (ref) provides more primitive sufficient conditions for Assumption (ref) with references to closely related work on convergence rates of neural networks and their usage in econometric theory. Notably, for models with analytic choice probabilities, we expect this assumption to require only $O(n^{\alpha + \delta})$ training examples for any $\delta > 0$.
We also require sufficient regularity of the true model.
Under the stronger condition $\alpha \geq 1$, the emulator-based estimator is asymptotically normal.
The proof, given in Appendix (ref), shows that $\hat{\boldsymbol{\theta}}_n$ is asymptotically equivalent to $\tilde{\boldsymbol{\theta}}_n$.
To estimate the Fisher information, define the true and emulator outer-product-of-scores estimators:
Assumption (ref) is justified by Sobolev training, which penalizes discrepancies between emulator and target gradients. The penalization is required because training on log-likelihood values alone does not guarantee convergence of gradients. As with Assumption (ref), see Appendix (ref) for more primitive sufficient conditions and additional discussion.
The emulator score $\hat{\psi}_n(y_i, \boldsymbol{\mathbf{X}}_i; \boldsymbol{\theta})$ is computable via automatic differentiation through the smooth neural network. Standard errors for $\hat{\boldsymbol{\theta}}_n$ are obtained as the square roots of the diagonal elements of $\hat{\boldsymbol{\mathbf{J}}}_n(\hat{\boldsymbol{\theta}}_n)^{-1}/n$.
An alternative estimator based on the Hessian, $\tilde{\boldsymbol{\mathbf{J}}}_n = -\frac{1}{n} \sum_{i=1}^n \nabla_{\boldsymbol{\theta}}^2 \hat{m}_n(y_i, \boldsymbol{\mathbf{X}}_i; \hat{\boldsymbol{\theta}}_n)$, requires second-derivative accuracy not directly targeted by first-order Sobolev training. The score-based estimator $\hat{\boldsymbol{\mathbf{J}}}_n(\hat{\boldsymbol{\theta}}_n)$ is therefore preferred.
When the emulator does not perfectly approximate the true likelihood, we can still obtain valid inference by treating the emulator as a working model, following the quasi-maximum likelihood framework white1982maximum. In this section, we treat the neural network as fixed with $n$, dropping the $n$ subscript from the observation log likelihood, $\hat{m}(y, \boldsymbol{\mathbf{X}}; \boldsymbol{\theta})$, and score function, $\hat{\psi}(y, \boldsymbol{\mathbf{X}}; \boldsymbol{\theta})$. Define the pseudo-true parameter as the population maximizer of the emulator log likelihood:
When the emulator closely approximates the true likelihood, $\boldsymbol{\theta}_0^{\dagger} \approx \boldsymbol{\theta}_0$. The asymptotic distribution depends on two matrices:
Under correct specification (i.e., a perfect emulator), the information matrix equality gives $\boldsymbol{\mathbf{A}}(\boldsymbol{\theta}_0) = \boldsymbol{\mathbf{B}}(\boldsymbol{\theta}_0) = \boldsymbol{\mathbf{J}}(\boldsymbol{\theta}_0)$. Under misspecification, these matrices differ, necessitating the sandwich covariance form. Under regularity conditions stated in Appendix (ref), we obtain the following results.
The proofs follow standard M-estimator arguments; see Appendix (ref). The sandwich covariance $\boldsymbol{\mathbf{V}} = \boldsymbol{\mathbf{A}}(\boldsymbol{\theta}_0^{\dagger})^{-1} \boldsymbol{\mathbf{B}}(\boldsymbol{\theta}_0^{\dagger}) \boldsymbol{\mathbf{A}}(\boldsymbol{\theta}_0^{\dagger})^{-\top}$ can be consistently estimated by $\hat{\boldsymbol{\mathbf{V}}}_n = \hat{\boldsymbol{\mathbf{A}}}_n^{-1} \hat{\boldsymbol{\mathbf{B}}}_n \hat{\boldsymbol{\mathbf{A}}}_n^{-\top}$, where
Both quantities are computable via automatic differentiation. Note that $\hat{\boldsymbol{\mathbf{B}}}_n$ corresponds to $\hat{\boldsymbol{\mathbf{J}}}_n(\hat{\boldsymbol{\theta}}_n)$ from the previous section under an asymptotic regime with a fixed emulator.
Theorem (ref) provides valid inference for $\boldsymbol{\theta}_0^{\dagger}$, not necessarily $\boldsymbol{\theta}_0$. If we employ an asymptotic regime where the emulator varies with $n$, then the pseudo-true parameter also varies with $n$, so we denote it as $\boldsymbol{\theta}_n^{\dagger}$. When $\alpha \in [0, 1)$ in Assumption (ref), we have $\boldsymbol{\theta}_n^{\dagger} \xrightarrow{p} \boldsymbol{\theta}_0$ as $n \to \infty$, but this convergence may be too slow relative to $\sqrt{n}$ sampling variability. The sandwich standard errors remain valid for $\boldsymbol{\theta}_n^{\dagger}$, but confidence intervals should not be interpreted as targeting $\boldsymbol{\theta}_0$. Even if $\boldsymbol{\mathbf{A}}(\boldsymbol{\theta}_0^{\dagger})$ and $\boldsymbol{\mathbf{B}}(\boldsymbol{\theta}_0^{\dagger})$ both approximate $\boldsymbol{\mathbf{J}}(\boldsymbol{\theta}_0)$, the centering at $\boldsymbol{\theta}_n^{\dagger} \neq \boldsymbol{\theta}_0$ introduces asymptotic bias. Valid inference for $\boldsymbol{\theta}_0$ requires $\alpha \geq 1$ (Theorem (ref)).
Similarly, Proposition (ref) provides conditions under which $\hat{\boldsymbol{\mathbf{B}}}_n \xrightarrow{p} \boldsymbol{\mathbf{J}}(\boldsymbol{\theta}_0)$. However, our first-order Sobolev training procedure provides no guarantee that the emulator Hessian will approximate the Hessian of the true model; thus, $\hat{\boldsymbol{\mathbf{A}}}_n$ need not converge to any quantity from the true model. However, Theorem (ref) requires only that $\boldsymbol{\mathbf{A}}(\boldsymbol{\theta}_0^{\dagger})$ be nonsingular, not that it match the true model. When the emulator additionally provides accurate second derivatives, the sandwich covariance matrix reduces to $\boldsymbol{\mathbf{J}}(\boldsymbol{\theta}_0)^{-1}$.
We evaluate the performance of emulator-based MNP estimators across a range of sample sizes, numbers of alternatives, and model specifications. We compare their statistical performance and computational requirements to estimators based on the GHK simulator.
Model specification. We consider MNP models with $K \in \{3, 5, 10\}$ alternatives and $p = 2$ covariates. Alternative $K$ serves as the reference with utility normalized to zero, representing an outside option. For each non-reference alternative $j \in \{1, \ldots, K - 1\}$, the deterministic utility is $v_j = \boldsymbol{x}_j^\top \boldsymbol{\beta}$, where $\boldsymbol{x}_j \in \mathbb{R}^p$ is a vector of alternative-specific attributes and $\boldsymbol{\beta} = (0, 1)^\top$ is the true coefficient vector. We consider two specifications for $\boldsymbol{\mathbf{\Sigma}}$, denoting them as Dense and Factor.
The Dense specification assumes a dense structure for $\boldsymbol{\mathbf{\Sigma}}_{1:(K-1),1:(K-1)}$, the upper $(K-1) \times (K-1)$ submatrix of $\boldsymbol{\mathbf{\Sigma}}$. For identification, we set $\Sigma_{11} = 1$ and $\boldsymbol{\mathbf{\Sigma}}_{\mathpalette\bigcdot@{.5} K} = \boldsymbol{\mathbf{\Sigma}}_{K \mathpalette\bigcdot@{.5}}^{\top} = \boldsymbol{0}_{K}$; the latter equation sets the variance and covariances of the zero-utility option to zero.
The Factor specification assumes a one-factor structure for $\boldsymbol{\mathbf{\Sigma}}_{1:(K-1),1:(K-1)}$ of the form $\boldsymbol{\mathbf{\Sigma}}_{1:(K-1),1:(K-1)} = \boldsymbol{\mathbf{\Psi}} + \boldsymbol{\gamma} \boldsymbol{\gamma}^{\top}$, where $\boldsymbol{\gamma} \in \mathbb{R}^{K-1}$ and $\boldsymbol{\mathbf{\Psi}} \in \mathbb{R}^{(K - 1) \times (K - 1)}$ is a diagonal matrix of uniquenesses. For identification, we impose an anchored normalization that fixes the first uniqueness and first factor loading:
As with the Dense specification, we set $\boldsymbol{\mathbf{\Sigma}}_{\mathpalette\bigcdot@{.5} K} = \boldsymbol{\mathbf{\Sigma}}_{K \mathpalette\bigcdot@{.5}}^{\top} = \boldsymbol{0}_{K}$. These constraints imply $\boldsymbol{\mathbf{\Sigma}}_{11} = 1$, fixing the scale of the model, and $\boldsymbol{\mathbf{\Sigma}}_{1j} = \gamma_j$ for $j \geq 2$, so that the factor loadings directly encode covariances with the first alternative. Given these constraints, we parameterize the model in terms of $K - 2$ free uniquenesses $(\Psi_{22}, \ldots, \Psi_{K-1,K-1})$ and $K - 2$ free loadings $(\gamma_2, \ldots, \gamma_{K-1})$. This anchored normalization is designed to eliminate flat directions in the likelihood surface.
Table (ref) shows the number of nonredundant parameters by type for the Dense and Factor specifications. With the anchored normalization, the Factor specification is exactly identified for $K = 3$ and overidentified for $K \geq 4$. For $K=10$, the Factor specification results in fewer than half as many parameters as the Dense specification (18 vs.\ 46).
Data generation. For each replication, we generate covariates $\boldsymbol{x}_{ij} \stackrel{\text{iid}}{\sim} \mathcal{N}(\boldsymbol{0}, \boldsymbol{\mathbf{I}}_p)$ for observations $i = 1, \ldots, n$ and alternatives $j = 1, \ldots, K - 1$, with $\boldsymbol{x}_{iK} = \boldsymbol{0}$ for the reference alternative.
For the Dense specification, we draw the covariance matrix according to a $\text{Wishart}(\boldsymbol{\mathbf{I}}_{K-1}, K + 10)$ distribution and then normalize it so that $\Sigma_{11} = 1$, resulting in a probability density function (PDF) given by
where $s = (K + 10)(K-1)/2$ and $\Gamma_{K-1}(\cdot)$ denotes the $(K-1)$-dimensional multivariate gamma function (see Appendix (ref)).
For the Factor specification, we sample the free uniquenesses on the log scale, drawing $\log(\Psi_{jj}) \overset{iid}{\sim} \mathcal{N}(0, 0.1^2)$ for $j = 2, \ldots, K-1$, and exponentiate to obtain $\Psi_{jj} > 0$. For the factor loadings, we draw $\gamma_j \overset{iid}{\sim} \mathcal{N}(0, 0.2^2)$ for $j = 2, \ldots, K-1$. The fixed values $\Psi_{11} = 0$ and $\gamma_1 = 1$ complete the parameterization.
Given the covariates and true parameters, we compute deterministic utilities and simulate choices $y_i \in \{1, \ldots, K\}$ according to one of two data-generating processes (DGPs): Probit or Emulator. The Probit RUM generates choices from the exact MNP model, drawing errors $\boldsymbol{\epsilon}_i \sim \mathcal{N}(\boldsymbol{0}, \boldsymbol{\mathbf{\Sigma}})$ and selecting the utility-maximizing alternative. The Emulator RUM generates choices by sampling directly from the categorical distribution defined by the emulator's predicted choice probabilities. When coupled with the Emulator RUM, the emulator-based estimation method is correctly specified, enabling a comparison of emulator-based estimators under both correct specification (Emulator RUM) and approximate inference (Probit RUM).
Estimation. To ensure positive definiteness of $\boldsymbol{\mathbf{\Sigma}}$, we estimate the Dense specification on the scale of the Cholesky factor, with the diagonal elements optimized on the log scale to enforce positivity. Similarly, we estimate the Factor specification with the free uniquenesses on the log scale, optimizing $\log(\Psi_{jj})$ for $j = 2, \ldots, K-1$.
We form estimators using two methods, Emulator and GHK, corresponding to whether we approximate the log-likelihood function using a pretrained emulator or the GHK simulator. We fit both methods to data generated by the Probit RUM; we also fit the Emulator method to data generated by the Emulator RUM to analyze its performance under correct specification. We apply the GHK method with $R \in \{10, 50, 250\}$ draws, representing low, moderate, and high computational effort.
Both methods, Emulator and GHK, perform regularized maximum likelihood estimation. The regularization consists of two components: a barrier penalty and a pseudo-prior. The barrier penalty is applied to $\log(\Psi_{jj})$ for the Factor specification and to the log-diagonal Cholesky parameters for the Dense specification. For a latent parameter $\eta$, the barrier penalty is $\lambda_{\text{reg}} \cdot \log(\eta + 8)/n$. The pseudo-prior terms take the form $\lambda_{\text{reg}} \cdot \log\left\{p(\boldsymbol{\theta})\right\} /n$, where $p(\boldsymbol{\theta})$ is the density function of $\boldsymbol{\theta}$ under the data-generating process. We add these penalties together to form a combined penalty of the form $\lambda_{\text{reg}} \cdot r(\boldsymbol{\theta})/n$, where $r(\boldsymbol{\theta})$ sums the barrier and log-prior penalties and $\lambda_{\text{reg}}$ is a regularization hyperparameter.
We then estimate $\boldsymbol{\theta}$ as the maximizer of the sum of the (approximate) log likelihood and $\lambda_{\text{reg}} \cdot r(\boldsymbol{\theta})/n$. We take $\lambda_{\text{reg}} = 0.1$ so that the penalty is small and asymptotically negligible. In practice, we find that these penalties lead to more reliable convergence of the model-fitting algorithm and help address weak identification of the covariance parameters that arises when $K$ is large and $n$ is small. We account for these penalties in the standard error computation by estimating $\mathrm{Cov}(\hat{\boldsymbol{\theta}}_n)$ as $\big\{n \cdot \hat{\boldsymbol{\mathbf{J}}}_n(\hat{\boldsymbol{\theta}}_n) - \lambda_{\text{reg}} \cdot \nabla_{\boldsymbol{\theta}}^2\, r(\hat{\boldsymbol{\theta}}_n)\big\}^{-1}$. Again, the impact of the penalization is asymptotically negligible because the term $n \cdot \hat{\boldsymbol{\mathbf{J}}}_n(\hat{\boldsymbol{\theta}}_n)$ dominates asymptotically.
For all methods, we maximize the (approximate) regularized log likelihood using L-BFGS optimization liu1989. In practice, we found it necessary to limit the minimum step size (we used $10^{-6}$) and allow small decreases in the objective function to enable the optimizer to escape local minima.
Experimental design. We consider sample sizes $n \in \{1{,}000, 10{,}000, 100{,}000\}$ and conduct 100 replications for each configuration. We compute performance metrics aggregated across all parameter types. We compute Monte Carlo standard errors for the coverage and average timings using scaled empirical standard deviations. For the remaining metrics, we apply bootstrap resampling (across repetitions) with 1{,}000 bootstrap samples to estimate the Monte Carlo standard errors.
Computational environment. All simulations were implemented in PyTorch paszke2019pytorch and executed on a personal laptop with 48 GB RAM and an Apple M4 Pro processor. We fit the emulators using the Apple Metal Performance Shaders (MPS). To obtain comparable estimation timings, however, we performed the rest of the simulation on CPU.
We fit the emulators for $400{,}000$ training episodes using $400{,}000$ pregenerated training examples, each consisting of $10^6$ choices. Details on the optimizer, learning rate schedule, and network architecture are provided in Appendix (ref). This model-fitting procedure completed in 2.1, 4.6, and 16.6 hours for $K=3, 5, $ and $10$, respectively, representing a modest level of computational effort.
Figure (ref) displays the training and evaluation loss for the emulators over time. At each value of $K \in \{3, 5, 10\}$, the training loss decreases consistently over time even at the end of training, indicating that further training would likely lead to additional improvements. The evaluation loss closely tracks the training loss, implying that the emulator generalizes well with minimal overfitting to the simulated training probabilities. We compute the evaluation loss every 200 episodes; these periodic calculations produce significantly less visible variation in the evaluation loss compared to the training loss.
Table (ref) displays performance metrics for $K=3$ comparing the performance of the emulator and GHK estimation methods. The performance metrics are the root mean squared error (RMSE), root mean squared bias (RMS Bias), ratio of estimated to empirical standard errors (SE Ratio), 95% confidence interval coverage rates (Coverage), and elapsed model-fitting time in seconds (Time (s)). The emulator methods perform well across all sample sizes, matching the performance of GHK(50) and GHK(250) in terms of statistical performance but at a lower computational cost. GHK(10) significantly undercovers at $n=100{,}000$, exhibiting substantially higher RMS Bias than the other methods. This poor performance is likely an artifact of the GHK simulator providing biased estimates on the log scale due to Jensen's inequality; for sufficiently large $n$, this bias is eventually nonnegligible relative to $1/\sqrt{n}$ sampling error.
Table (ref) displays similar results for $K=5$. The emulator is competitive with GHK(50) and GHK(250) at all sample sizes. Compared to GHK(50), the compute time for the Emulator method is slightly lower for $n=1{,}000$, approximately equal for $n = 10{,}000$, and slightly higher for $n = 100{,}000$.
Table (ref) displays performance results for $K=10$. The emulator performs comparably to GHK(250) with much lower computational requirements. At $n=1{,}000$, the emulator fits the model in less than 3 seconds on average, outpacing even GHK(10). The computation time of the emulator falls between that of GHK(10) and GHK(50) at $n=10{,}000$ and between GHK(50) and GHK(250) at $n=100{,}000$. Compared to Tables (ref) and (ref), the standard error estimates are less reliable for the smaller sample sizes, presumably due to weak identification; the problem is amplified when we omit the penalty on $\boldsymbol{\mathbf{\Sigma}}$. The standard errors improve substantially with the sample size (with or without penalization), reaching SE Ratios of 0.92--1.07 for the Emulator and GHK(250) methods at $n=100{,}000$.
Table (ref) displays analogous results for the Factor specification with $K=3$. Performance metrics are averaged over the nonredundant elements of $\boldsymbol{\beta}$, $\boldsymbol{\mathbf{\Psi}}$, and $\boldsymbol{\gamma}$. The Emulator's performance mirrors that of GHK(250) but at a much lower computational cost.
Tables (ref) and (ref) display analogous results for 5 and 10 alternatives, respectively. In both cases, the Emulator performs comparably to GHK(50) and GHK(250) despite not being trained specifically for factor models. The Emulator is faster than GHK(50) for $n \in \{1{,}000, 10{,}000\}$. At $n = 100{,}000$, the Emulator's computation time falls between that of GHK(50) and GHK(250).
Across all scenarios, the Emulator method performs competitively with or significantly better than GHK. We expect that the relative performance of the emulator would increase with (a) additional training and (b) implementation on specialized deep-learning hardware, such as GPUs. The latter is likely to be especially advantageous for the emulator because its computations are trivially parallelizable and highly optimized on modern GPUs.
We conclude this section by providing numerical results from training a neural network emulator to approximate choice probabilities for $K \in \{3, 4, 5\}$ simultaneously. We use the same neural network as we used in the previous section for $K=5$. As above, we train for $400{,}000$ episodes, but we loop over training examples with $K=3$, $4$, and $5$ alternatives within each episode. The training process completes in 13.0 hours.
Figure (ref) shows the training and evaluation loss steadily decreasing throughout the training procedure, approaching levels similar to those of Figure (ref). Figure (ref) plots the estimated choice probabilities against estimation errors for $K \in \{3, 4, 5\}$ with 400 simulated examples each. The small errors indicate that the model accurately approximates the choice probabilities across all three values of $K$ simultaneously.
This paper proposes an amortized inference approach for discrete choice models. The approach uses a novel neural network architecture to approximate choice probabilities under general error distributions, including those featuring nontrivial cross-alternative correlations, such as MNP models. This generality enables flexible modeling of empirical choice behavior, including complex substitution patterns, without losing the theoretical grounding or interpretability of RUM models.
Because the emulator operates on the generic choice problem, it is agnostic to the assumed parametric form of the deterministic utilities or the accompanying scale (covariance) matrix. Given a pretrained emulator, this generality enables drop-in adoption without the need for custom architecture changes or additional training. In particular, generalizing from logit to probit requires only the replacement of closed-form softmax probabilities with emulator evaluations. Other correlated error distributions are similarly straightforward.
Our proposed architecture and training process feature a preprocessing transformation, DeepSet modules zaheer2017deep, equivariant layers, smooth activation functions, and Sobolev training czarnecki2017sobolev. Together, these design choices produce smooth approximations of choice probabilities that respect their invariance properties and accelerate learning via parameter sharing. We show that the architecture possesses a universal approximation property, extending recent theoretical results for symmetric matrices blumsmith2025machine to the joint space of utility vectors and centered scale (covariance) matrices.
In Monte Carlo simulations, we find that emulator-based estimators perform as well or better than GHK-based estimators in terms of both statistical performance and computational requirements when estimating MNP models. These findings support our theoretical results showing that emulator-based estimators inherit the properties of exact MLEs under mild approximation conditions. We focus on MNP models in the simulation study to enable comparison to the GHK simulator, which is directly applicable only to Gaussian errors. The computational benefits of the amortized inference approach are amplified by the fact that emulator evaluations are trivially parallelizable and highly optimized on modern deep learning hardware, such as GPUs.
Several limitations merit discussion. Practitioners must ensure sufficient overlap between the distribution of inputs in their application and the training distribution. Additionally, our asymptotic theory requires the emulator approximation error to vanish at an appropriate rate, meaning applications with very large samples may require more extensively trained emulators. Lastly, models with correlated errors require identifying restrictions not needed for independent logit-based models. Although our approach enables simple and efficient computation of choice probabilities for these models, the data requirements and modeling decisions to ensure proper identification remain a separate challenge.
The framework also allows numerous extensions. Our emulator architecture and training procedure could be specialized to specific settings, such as factor-structured scale matrices or mixture distributions. We expect these could be handled by reconfiguring the modules in our architecture and, in some cases, generalizing the preprocessing transformation. In particular, the two preceding examples would likely require an additional DeepSet operation over factors or mixture components. Our emulators could also be used in tandem with other modern choice modeling architectures, such as those described in Section (ref). Previous work has applied neural-network-assisted estimation approaches to other econometric models with intractable likelihood functions, such as dynamic discrete choice Norets2012 and consumer search wei2025, in the context of specific parametric models. The amortized inference approach presented here could likely be generalized to many of these settings, enabling efficient amortized inference without modifying or retraining the neural network.
The authors acknowledge the use of AI tools to assist in the development of this manuscript. All content has been manually reviewed for accuracy.