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.
91,554 characters · 22 sections · 60 citation commands
On the Estimation of Multinomial Logit and Nested Logit Models: A Conic Optimization Approach
\newcolumntype{C}{>{\arraybackslash}p{4em}}
\spacing{1.5}
\affil[1]{\it School of Computing and Information Systems, Singapore Management University} \affil[2]{\it SLSCM and CADA, Faculty of Data Science and Artificial Intelligence, College of Technology, National Economics University, Hanoi, Vietnam} \affil[*]{\it Corresponding author, [email removed]}
{\bf Keywords:} Discrete Choice Models; Multinomial Logit; Nested Logit; Tree Nested Logit; Exponential Cone Programming; Parameter Estimation
Discrete choice models are fundamental tools for analyzing individual decision-making in domains such as transportation, marketing, industrial organization, and operations research. Since the seminal contributions of McFadden1974 and BenAkivaLerman1985, the multinomial logit (MNL), nested logit (NL), and tree--nested logit (TNL) models have become standard due to their behavioral interpretability, closed-form choice probabilities, and analytical tractability Train2009. Estimation is typically carried out via maximum likelihood, solved with gradient-based routines such as Newton, quasi-Newton, or sequential quadratic programming methods BerndtEtAl1974,NocedalWright2006. While these methods are effective in small- to medium-scale applications, they can be sensitive to step-size selection and initialization, and their numerical performance often deteriorates in large-scale settings with high-dimensional covariates, large choice sets, or complex nesting structures BorschSupanHajivassiliou1993,Train2009,HessDaly2014.
This paper revisits estimation of MNL, NL, and TNL through the lens of convex conic optimization. We show that the MLE for the MNL model admits an exact reformulation as an exponential-cone program (ECP) BoydVandenberghe2004,Nesterov1994interior,MOSEKcookbook. For NL and TNL, when the dissimilarity (scale) parameters are fixed, the MLE objective is convex and can likewise be written as an ECP. Building on these observations, we propose a two-stage procedure that alternates between (i) an outer update of the nesting scale parameters and (ii) an inner ECP that updates the utility coefficients.
Casting estimation as an ECP brings algorithmic and practical benefits. First, it enables the use of primal–dual interior-point methods with iteration counts that grow only logarithmically with the desired accuracy Nesterov1994interior,BoydVandenberghe2004. Second, off-the-shelf solvers implementing these methods MOSEKcookbook can exploit the block structure induced by observations and alternatives, delivering robust performance without delicate hand-tuning of step sizes or line-search heuristics. Third, the conic viewpoint yields unified formulations across MNL, NL, and TNL, simplifying implementation and facilitating large-scale estimation.
We evaluate the proposed approach on synthetic datasets spanning a wide range of dimensions, numbers of nests, and choice-set sizes. Across all regimes, the ECP reformulations consistently deliver better MLE solutions and substantially lower wall-clock times than strong gradient-based baselines, with the largest gains appearing on medium- and large-scale instances. The results underscore that exponential-cone programming is a practical, scalable, and numerically stable alternative to standard gradient-based MLE for a broad class of discrete-choice models.
\paragraph{Contributions.} This paper makes the following key contributions. First, we provide an exact ECP reformulation of the maximum-likelihood estimator for the MNL model. Second, we develop convex ECP reformulations for NL and TNL when the dissimilarity parameters are fixed, and propose a two-stage procedure that alternates between updating the scale parameters and solving the convex inner problem for the utility coefficients. Third, we present an extensive empirical study showing that our approach delivers superior scalability, robustness, and solution quality compared to standard gradient-based estimators. To the best of our knowledge, this is the first work to formulate and solve the estimation of several widely used discrete choice models using exponential-cone programming, achieving state-of-the-art performance compared to standard gradient-based methods. Moreover, building on these findings, we are the first to establish polynomial-time guarantees for estimating the MNL model, as well as the NL and TNL models with fixed scale parameters.
\paragraph{Organization.} The remainder of the paper is organized as follows. Section (ref) presents a literature review. Section (ref) introduces the MNL formulation and its ECP reformulation. Sections (ref) and (ref) extend the approach to NL and TNL under fixed scale parameters and describes the two-stage estimation scheme. Section (ref) reports numerical results, and Section (ref) concludes with a discussion of implications and future research directions. The appendix contains proofs omitted from the main paper, along with detailed numerical analyses.
\paragraph{Notation:} Boldface characters represent matrices (or vectors), and $a_i$ denotes the $i$-th element of vector $\textbf{a}$. We use $[m]$, for any $m\in \mathbb{N}$, to denote the set $\{1,\ldots,m\}$.
Discrete-choice models have a long history as empirical workhorses in transportation, marketing, and operations. The MNL model occupies a central place because of its interpretability and analytical tractability under i.i.d.\ Type-I extreme value errors McFadden1974,BenAkivaLerman1985,Train2009. The nested logit and tree–nested logit models generalize the framework to accommodate richer substitution patterns when alternatives share unobserved components, using generating functions from the Generalized Extreme Value (GEV) family Williams1977,DalyZachary1978,Train2009. Across these models, estimation is typically framed as maximum likelihood, with algorithms built around smooth optimization of the log-likelihood and its derivatives.
For the MNL model, the log-likelihood is strictly concave in the utility coefficients, which makes maximum-likelihood estimation well-posed and enables the use of Newton, Fisher scoring/BHHH, and quasi-Newton methods such as L-BFGS BerndtEtAl1974,LiuNocedal1989,Train2009. Practical implementations emphasize efficient evaluation of utilities and choice probabilities across large numbers of observations and alternatives, as well as robust variance estimation Train2009. When the universe of alternatives is very large, sampling-of-alternatives with correction terms becomes important to reduce computational burden while preserving consistency BenAkivaLerman1985. Regularization has also become common, with $\ell_2$ penalties stabilizing estimates and $\ell_1$ penalties enabling high-dimensional specifications HastieTibshiraniFriedman2009. Despite these advances, performance can deteriorate when covariates are poorly scaled, when choice sets are massive, or when gradients are ill-conditioned, making step-size selection, line searches, and preconditioning critical in practice NocedalWright2006.
It is well known that the MNL model satisfies the independence from irrelevant alternatives (IIA) property, which states that the relative odds of choosing between two alternatives are unaffected by the presence or characteristics of other options. This assumption is often unrealistic in practice. The NL model relaxes IIA by allowing alternatives within a nest to share unobserved components Williams1977,DalyZachary1978. From an estimation perspective, the likelihood is concave in the utility coefficients conditional on the nest dissimilarity (scale) parameters, but the joint problem over both sets of parameters is not generally concave Train2009,DaganzoKusnic1993_convex_nested. This has led to two common strategies: fixing the scale parameters based on prior evidence and estimating the utilities by standard MLE, or alternating between updates of the scales and conditional concave updates of the utilities Train2009. To ensure identification, the NL model requires normalizing the overall scale and restricting the dissimilarity parameters to lie within valid bounds. In practice, estimation routines impose these conditions using box or nonlinear constraints Train2009,NoceWrig06. The use of analytical “inclusive value” (logsum) formulas helps keep computations manageable, but the extra constraints and nonconvexity make estimation more challenging and less straightforward than in the MNL model BenAkivaLerman1985,Train2009.
The TNL extends the NL model to multi-level correlation structures, attaching scale parameters to internal nodes along each branch of a choice hierarchy Train2009,Daly1987. This added flexibility improves behavioral realism in settings with deep category structures, but it also complicates estimation: the likelihood remains convex in the utility coefficients for fixed scales, while the full problem (utilities plus all scale parameters) is generally nonconvex and subject to additional admissibility restrictions along each path of the tree Train2009,DaganzoKusnic1993_convex_nested. As the depth and breadth of the nesting structure grow, conditioning worsens and gradient-based methods can become sensitive to initialization and step-size rules. In applied work, it is therefore common either to fix a subset of scales or to employ alternating or two-step procedures that separate scale updates from utility estimation to preserve conditional concavity Train2009,Daly1987.
Beyond the MNL and NL models, two important generalizations are the cross-nested logit (CNL) Bier06 and the mixed logit (MMNL) models McFaddenTrain2000. The CNL model extends the standard nesting structure by allowing an alternative to belong to multiple nests with fractional allocation parameters, thereby capturing flexible substitution patterns across overlapping groups VovshaBekhor2002. The MMNL model, also known as the random parameters logit, introduces random taste heterogeneity by allowing coefficients in the utility function to vary across individuals according to a specified distribution McFaddenTrain2000,Train2009. While both models are behaviorally attractive and widely applied, their estimation is significantly more challenging. In the CNL case, the log-likelihood is generally non-concave because the allocation parameters interact nonlinearly with the inclusive values, complicating identification and leading to multiple local optima. Similarly, in the MMNL model, the likelihood involves high-dimensional integrals over the random coefficients that are approximated through simulation, resulting in a simulated log-likelihood that is also non-concave in the parameters. These non-convexities mean that standard maximum-likelihood solvers may converge slowly or become trapped in local maxima, and that robust and scalable estimation of CNL and MMNL remains substantially more delicate than for MNL or NL.
The estimation of these discrete choice models has also been implemented in several widely used software tools and libraries. In transportation research, packages such as BIOGEME Bierlaire2003,Bierlaire2020 and Alogit Daly1999 are standard, while in economics and marketing, implementations are available in econometric toolboxes such as Stata, NLOGIT/LIMDEP Greene2016, and more recently in open-source environments such as R (mlogit, Croissant2019) and \texttt{Python} (\texttt{PyLogit}, Brathwaite2018). Despite differences in interfaces and supported model classes, the underlying estimation routines across these platforms remain largely gradient-based, typically relying on Newton-type or quasi-Newton algorithms to solve the maximum-likelihood problem.
Our work explores a novel estimation approach for the aforementioned discrete choice models by leveraging convex conic optimization. In this context, it worth mentioning that conic optimization provides a unified framework for solving a broad class of convex programs by expressing them as linear optimization problems over convex cones BenTalNemirovski2001,BoydVandenberghe2004. A key advantage of this formulation is that it enables the use of primal--dual interior-point methods with strong polynomial-time complexity guarantees NesterovNemirovskii1994. Within this framework, the exponential cone has emerged as a powerful modeling tool, as it allows nonlinear functions such as the exponential, relative entropy, and the log-sum-exp to be represented exactly in a conic form Chares2009. Moreover, casting problems as exponential-cone programs makes them amenable to highly robust and efficient conic solvers such as MOSEK MOSEKcookbook, which can handle large-scale instances without the delicate step-size tuning required by first-order methods. Exponential-cone programming underpins applications ranging from geometric programming BoydGP2007 and logistic regression McCullaghNelder1989 to entropic optimal transport Cuturi2013 and distributionally robust optimization NamkoongDuchi2017, highlighting its flexibility and wide-ranging impact. While conic optimization, and in particular exponential-cone programming, has been applied in several domains, to the best of our knowledge, this work is the first to employ these techniques for the estimation of discrete choice models.
The multinomial logit (MNL) model is perhaps the most widely used specification in discrete choice analysis McFadden1974,Train2009. It assumes that the utility of each alternative consists of a systematic component, linear in observed attributes, plus an unobserved error term following the Type-I extreme value distribution. This assumption yields closed-form choice probabilities, simple behavioral interpretations of parameters, and tractable likelihood functions. Despite its limitations (e.g., the independence of irrelevant alternatives property - IIA), the MNL model remains the workhorse of empirical studies across transportation, marketing, and operations.
We first present the standard maximum likelihood estimation of the MNL model. Let $m$ be the number of alternatives, and denote by $[m] = \{1,\ldots,m\}$ the universal set of alternatives. For each individual $i$, the deterministic component of the utility associated with alternative $j \in [m]$ is $v_{ij} = \pmb{\beta}^\top \textbf{a}_{ij},$ where $\pmb{\beta}$ is the parameter vector to be estimated, and $\textbf{a}_{ij}$ is the vector of observed attributes of alternative $j$ for individual $i$. The full random utility is defined as $u_{ij} = v_{ij} + \mu \epsilon_{ij},$ where $\epsilon_{ij}$ are i.i.d.\ disturbances following the Type I extreme value (Gumbel) distribution. For notational simplicity, we set $\mu = 1$ without loss of generality, since any $\mu \neq 1$ can be absorbed into the coefficient vector by defining $\tilde{\pmb{\beta}} = \pmb{\beta} / \mu$.
Under this distributional assumption, the choice probabilities take the familiar logit form \[ P_{ij}(S_i) = \frac{\exp(\pmb{\beta}^\top \textbf{a}_{ij})}{\sum_{k \in S_i}\exp(\pmb{\beta}^\top \textbf{a}_{ik})}, \] where $S_i \subseteq [m]$ denotes the choice set available to individual $i$.
Given $N$ independent observations $(j_n, S_n)$, where $S_n \subseteq [m]$ is the choice set of individual $n$ and $j_n \in S_n$ is the chosen alternative, the maximum likelihood estimator (MLE) solves
Equivalently, the log-likelihood can be expressed as
which is a concave function problem involving log-sum-exp terms BoydVandenberghe2004.
A typical approach to estimation is to compute the gradient of the log-likelihood, \[ \nabla_{\pmb{\beta}} {\mathcal{L}}^{\textsc{MNL}}(\pmb{\beta}) = \sum_{n=1}^N \left( \textbf{a}_{nj_n} - \sum_{j \in S_n} P_{ij}(S_n)\, \textbf{a}_{nj} \right), \] and then apply a gradient-based routine such as Newton’s method, quasi-Newton methods (e.g., BFGS or L-BFGS), or other iterative schemes NoceWrig06. While effective for small- to medium-scale problems, such methods require careful step-size selection and can become computationally demanding when both $N$ and $m$ are large.
We now reformulate the MLE problem under the MNL model as an exponential cone program. The key observation is that the log-sum-exp terms in the likelihood can be expressed exactly using the exponential cone, a convex set that captures exponential and logarithmic relationships in a conic form. The exponential cone is typically defined as
This cone provides a convex representation of nonlinear expressions involving exponentials, logarithms, and relative entropy, and has become a standard tool in modern convex optimization Chares2009,MOSEKcookbook.
Figure (ref) illustrates the boundary of the exponential cone. The curved surface corresponds to points satisfying $x_1 = x_2 e^{x_3/x_2}$ for $x_2>0$, while the flat region at $x_2=0$ corresponds to the closure of the cone for $x_3 \leq 0$. Together, these surfaces delineate the convex region ${\mathcal{K}}_{\exp}$ that forms the foundation of our reformulation.
In the context of MLE for the MNL model, the exponential cone provides an exact convex representation for exponential and log-sum-exp functions. In particular, problem (ref) can equivalently be written as:
Moreover, the equality constraints can be relaxed to the inequality: \[ t_n \;\geq\; \log \Big(\sum_{j\in S_n}\exp(\pmb{\beta}^\top \textbf{a}_{nj})\Big), \] since in the maximization problem the optimizer will always drive $t_n$ to the smallest feasible value, so equality holds at the optimum. Rearranging, this inequality is equivalent to \[ \sum_{j\in S_n} \exp\!\big(\pmb{\beta}^\top \textbf{a}_{nj} - t_n\big) \;\leq\; 1, \qquad \forall n \in [N]. \] Introducing auxiliary variables $z_{nj} = \exp(\pmb{\beta}^\top \textbf{a}_{nj} - t_n)$, the constraints $\sum_{j\in S_n} \exp\!\big(\pmb{\beta}^\top \textbf{a}_{nj} - t_n\big) \;\leq\; 1 $ can be equivalently written as
The nonlinear inequalities $z_{nj} \geq \exp(\pmb{\beta}^\top \textbf{a}_{nj} - t_n)$ can in turn be expressed using exponential cone constraints: \[ (z_{nj},\, 1,\, \pmb{\beta}^\top \textbf{a}_{nj} - t_n) \;\in\; {\mathcal{K}}_{\exp}, \qquad \forall n \in [N], \; j \in S_n. \] Thus, problem (ref) admits the ECP formulation:
Figure (ref) illustrates the key transformation of the log-sum-exp constraint into its equivalent conic representation. The ECP formulation in (ref) is equivalent to the classical MLE but is expressed directly in the language of convex conic optimization. As a result, it can be solved efficiently by modern conic solvers (e.g., MOSEK) using interior-point methods, which provide robustness, polynomial-time complexity guarantees MOSEKcookbook,Nesterov1994interior.
As noted earlier, a key advantage of the ECP reformulation is that it enables the MLE to be solved using primal--dual interior-point methods, which enjoy polynomial-time complexity guarantees that are difficult to achieve with standard gradient-based approaches. In this section, we analyze the computational complexity of solving the exponential cone program in (ref). We begin by examining the size of the program, in terms of its variables and constraints, and then provide a formal statement on the complexity of attaining an $\varepsilon$-optimal solution.
The exponential-cone reformulation in (ref) introduces a total of $p+N+Z$ decision variables, where $p$ denotes the number of utility coefficients (the dimension of $\pmb{\beta}$), $N$ is the number of individuals, and $Z=\sum_{n=1}^N |S_n|$ represents the total number of alternative occurrences across all observations. These variables correspond respectively to the parameter vector $\pmb{\beta}$, the auxiliary scalars $t_n$, and the exponential-cone variables $z_{nj}$. The constraints consist of $N$ linear inequalities of the form $\sum_{j\in S_n} z_{nj}\leq 1$ together with $Z$ exponential-cone membership constraints $(z_{nj},1,\pmb{\beta}^\top \textbf{a}_{nj}-t_n)\in{\mathcal{K}}_{\exp}$. When expressed in the standard form used by conic solvers, each linear inequality is represented as an equality with an additional nonnegative slack, yielding in total $p+Z+2N$ variables, $N$ linear equalities, and $Z$ exponential-cone blocks of dimension three. In terms of scaling, if $\bar m=\tfrac{1}{N}\sum_n |S_n|$ denotes the average choice-set size, then $Z=N\bar m$, so the problem size grows linearly with both $N$ and $\bar m$, and in the full-availability worst case one has $Z\leq Nm$.
The following proposition summarizes the overall complexity of solving (ref) via interior-point methods.
The proof (provided in the appendix) follows directly from the standard complexity analysis of interior-point methods for self-concordant barriers NesterovNemirovskii1994,BenTalNemirovski2001. In this framework, the iteration complexity is governed by the barrier parameter $\nu$, which in our case scales as $\nu = \Theta(Z)$ because each exponential-cone block contributes a constant barrier parameter. Thus, a path-following method requires $\mathcal{O}(\sqrt{Z}\log(1/\varepsilon))$ Newton steps to obtain an $\varepsilon$-optimal solution. The cost per iteration is dominated by assembling and factoring the Schur complement in the global variables $(\pmb{\beta},\textbf{t})$, which requires $\mathcal{O}(Zp^2+(p+N)^3)$ arithmetic operations. Combining these bounds yields the stated overall complexity.
The complexity stated in Proposition (ref) implies that the conic reformulation of MNL estimation enjoys the same strong polynomial-time guarantees as other problems in convex optimization. In particular, the number of interior-point iterations grows only with the square root of the total number of exponential-cone blocks $Z$, and only logarithmically with the target accuracy $\varepsilon$. This means that the bulk of the computational burden lies in the per-iteration cost of assembling and solving the Schur complement system, which scales as $\mathcal{O}(Zp^2+(p+N)^3)$. From a practical perspective, this structure is highly favorable: the $\mathcal{O}(Zp^2)$ assembly term is embarrassingly parallel across individuals and alternatives, while the $\mathcal{O}((p+N)^3)$ factorization step depends only on the number of parameters and individuals, not directly on the choice set size. Hence, the conic approach is particularly well suited for large-scale applications where $N$ and the average choice set size are large, and provides robustness and scalability advantages over traditional gradient-based maximum likelihood routines, which lack such complexity guarantees and may suffer from convergence issues in ill-conditioned settings.
The NL model is one of the most widely used extensions of the multinomial logit because of its ability to relax the restrictive IIA property and capture correlations among groups of similar alternatives Williams1977,DalyZachary1978,Train2009. Its popularity stems from both its behavioral appeal and its analytical tractability, making it a standard tool in transportation research for mode and route choice analysis, in marketing for product differentiation and market share prediction, and in economics and operations for demand estimation and policy evaluation. Estimation of the NL model typically proceeds by maximum likelihood, and it is well known that when the dissimilarity (scale) parameters are fixed, the log-likelihood is concave in the utility coefficients DaganzoKusnic1993_convex_nested. However, joint estimation of the utility and dissimilarity parameters leads to a non-convex problem, raising concerns of local optima and numerical instability BenAkivaLerman1985,Train2009. As a result, robust and scalable estimation of the NL model remains more challenging than that of the simpler MNL, despite its broader applicability and flexibility.
In the NL model, alternatives are partitioned into nests, which form disjoint subsets of the overall choice set. The choice probability of each alternative can be decomposed into two components: first, the probability of choosing the nest that contains the alternative, and second, the conditional probability of choosing the alternative within that nest. Figure (ref) illustrates a simple nested structure with two nests, each containing two alternatives. For instance, the choice probability of alternative $j_1$ can be expressed as \[ P(j_1) \;=\; P(\text{Nest A}) \cdot P(j_1 \mid \text{Nest A}), \] where $P(\text{Nest A})$ is the probability of selecting Nest A at the upper level, and $P(j_1 \mid \text{Nest A})$ is the conditional logit probability of choosing $j_1$ among the alternatives within Nest A. Each nest is associated with a dissimilarity (or scale) parameter $\lambda \in (0,1]$, which governs the degree of correlation among the alternatives inside the nest: values of $\lambda$ closer to one indicate weaker correlation (approaching the standard MNL), while smaller values of $\lambda$ allow for stronger correlations within the nest.
We first present the MLE of the NL model. In the NL model, the full set of alternatives is partitioned into \(L\) disjoint subsets (or nests) \({\mathcal{N}}_1, \ldots, {\mathcal{N}}_L \subset [m]\). Each nest \(l \in [L]\) is associated with a dissimilarity (scale) parameter \(\lambda_l \in (0,1]\), which measures the degree of independence in unobserved utility among the alternatives within that nest.
A central feature of the NL model is that the probability of choosing an alternative can be decomposed into two components: (i) the probability of selecting the nest that contains the alternative, and (ii) the conditional probability of choosing the alternative within that nest. Formally, for an observed choice $j_n \in S_n$ made by individual $n$, we can write \[ P_n(j_n \mid S_n) \;=\; P_n(l_n \mid S_n)\,\cdot\,P_n(j_n \mid l_n, S_n), \] where $l_n$ denotes the nest containing $j_n$. The probability of choosing nest $l$ is given by \[ P_n(l \mid S_n) \;=\; \frac{W_{nl}^{\lambda_l}}{\sum_{l' \in [L]} W_{nl'}^{\lambda_{l'}}},~~\text{where } W_{nl} \;=\; \sum_{j \in {\mathcal{N}}_l \cap S_n} \exp\!\left(\frac{\pmb{\beta}^\top \textbf{a}_{nj}}{\lambda_l}\right) \] is the so-called inclusive value (or log-sum term) for nest $l$. Conditional on nest $l$, the probability of choosing alternative $j \in {\mathcal{N}}_l$ is \[ P_n(j \mid l, S_n) \;=\; \frac{\exp\!\left(\pmb{\beta}^\top \textbf{a}_{nj}/\lambda_l\right)}{W_{nl}}. \] Combining the two terms, the overall choice probability of $j_n$ can be written as \[ P_n(j_n \mid S_n) \;=\; \frac{W_{n l_n}^{\lambda_{l_n}-1}\, \exp\!\left(\pmb{\beta}^\top \textbf{a}_{n j_n}/\lambda_{l_n}\right)} {\sum_{l' \in [L]} W_{nl'}^{\lambda_{l'}}}. \] The NL model reduces to the standard MNL model when all scale parameters satisfy $\lambda_l = 1$ for every $l \in [L]$, in which case nests collapse to single-level choice sets.
Given the above choice probabilities, the MLE of the NL model can be formulated as
where $P_n(j_n \mid S_n)$ is the nested logit probability of observing choice $j_n$ from the offered set $S_n$ for individual $n$. Expanding the probability expression, the log-likelihood can be written equivalently as
The NL log-likelihood has a more complicated structure than that of the MNL model. Moreover, it can be shown that if the scale parameters ${\pmb{\lambda}}$ are fixed, then ${\mathcal{L}}^{\textsc{\tiny NL}}(\pmb{\beta},{\pmb{\lambda}})$ is concave in $\pmb{\beta}$ DaganzoKusnic1993_convex_nested. Standard approaches therefore typically rely on iterative optimization methods that compute the gradient (and sometimes the Hessian) of the log-likelihood with respect to $\pmb{\beta}$ and update the parameters using Newton or quasi-Newton schemes such as BFGS or L-BFGS. When both $\pmb{\beta}$ and ${\pmb{\lambda}}$ are estimated simultaneously, the problem becomes non-convex, and specialized algorithms or two-step procedures are often employed to improve stability and convergence.
We now present an exponential cone programming (ECP) reformulation of the NL maximum likelihood problem in (ref). The key idea is to replace the nested log-sum-exp and log-sum terms that appear in the likelihood with equivalent exponential-cone constraints, thereby obtaining a tractable conic representation of the estimation problem. To this end, we introduce auxiliary variables that explicitly represent the inclusive values within each nest and the top-level aggregation across nests. Specifically, for individual $n$ and nest $l \in [L]$, define
where $z_{nl}$ represents the log of the inclusive value within nest $l$, and $y_n$ represents the top-level log-sum across all nests. With these variables, the NL log-likelihood can be rewritten as the following constrained optimization problem:
Here we can see that the equality constraints (ref)--(ref) can be safely relaxed to inequalities
without changing the optimal solution. This is because $\lambda_l \leq 1$, so in the maximization problem the optimizer will always push $z_{nl}$ and $y_n$ to their smallest feasible values, ensuring the inequalities are tight at optimality. Next, we rewrite these inequalities in exponential-cone form. Rearranging gives:
By introducing additional auxiliary variables:
the problem can be written as the following exponential cone program:
This formulation replaces the log-sum-exp and log-sum expressions in the NL likelihood with linear and exponential-cone constraints, thereby making the entire estimation problem amenable to modern conic solvers. It preserves the convexity structure when the scale parameters $\{\lambda_l\}$ are fixed and provides a principled way to incorporate NL estimation into the exponential cone programming framework.
Similar to the analysis in the previous section, we now discuss the computational complexity of solving (ref) using interior-point algorithms. For ease of notation, similar to the case of MNL model, let \[ Z \;:=\; \sum_{n=1}^N |S_n|, \qquad L_n \;:=\; \big|\{\,l\in[L] : {\mathcal{N}}_l \cap S_n \neq \emptyset\,\}\big|, \qquad \Lambda \;:=\; \sum_{n=1}^N L_n. \] Here, $Z$ counts the total number of alternative appearances across all observations, while $\Lambda$ counts the total number of active nests across all individuals (i.e., nests that contain at least one available alternative).
In (ref), the main decision variables are as follows: $\pmb{\beta} \in \mathbb{R}^p$ ($p$ variables), $y_n$ ($N$ variables), $z_{nl}$ ($\Lambda$ variables), $k_{njl}$ (one for each $(n,j)$ with $j \in {\mathcal{N}}_l \cap S_n$, contributing a total of $Z$ variables since nests are disjoint), and $h_{nl}$ ($\Lambda$ variables). Hence, the model has in total $p \;+\; N \;+\; Z \;+\; 2\Lambda$ decision variables. The linear constraints consist of: (i) one inequality $\sum_{j \in {\mathcal{N}}_l \cap S_n} k_{njl} \leq 1$ for each active pair $(n,l)$, and (ii) one inequality $\sum_{l \in L_n} h_{nl} \leq 1$ for each $n$, for a total of $\Lambda + N$ linear inequalities (each of which introduces a nonnegative slack variable in solver-ready form).
The exponential-cone constraints comprise: (i) one block for each $k_{njl}$ ($Z$ blocks), and (ii) one block for each $h_{nl}$ ($\Lambda$ blocks), so that there are in total $Z \;+\; \Lambda$ exponential-cone blocks. In the worst case of full availability ($S_n = [m]$ for all $n$) and all nests active for every individual, we have $Z = N m$ and $\Lambda = N L$.
Given this problem size, the following proposition states the polynomial-time complexity of solving (ref) with a path-following interior-point method Nesterov1994interior.
The proof of Proposition (ref) can be found in the appendix. The complexity analysis for the NL model highlights both the benefits and the challenges of the exponential cone reformulation. On the one hand, once the scale parameters are fixed, the estimation problem is convex and the inner step can be solved to global optimality with polynomial-time guarantees using interior-point methods. On the other hand, the overall computational cost grows with the total number of observations, the size of the choice sets, and the number of active nests, reflecting the richer structure of the NL compared to the MNL. This scaling underscores why estimation of NL models is substantially more demanding than MNL, and also explains the practical importance of efficient solvers and structural exploitation (e.g., parallel assembly or sparse linear algebra) to handle large-scale applications.
We discuss the estimation of the TNL model. The TNL model extends the standard NL by allowing a hierarchical structure of nests organized as a tree Train2009, Daly1987. This generalization is particularly useful in applications where alternatives share unobserved components at multiple levels of aggregation, for example in transportation where travelers first choose a travel mode, then a service type within that mode, and finally a specific route, or in marketing where products can be classified by brand, subcategory, and item. The tree structure provides greater flexibility than the two-level NL model, as it can capture correlations across alternatives at different depths of the hierarchy. However, this added flexibility comes at the cost of increased complexity in estimation. In particular, the likelihood involves nested inclusive values at multiple levels of the tree, and while it remains concave in the utility coefficients conditional on fixed scale parameters at each node, the joint estimation of both utilities and scale parameters is non-convex. As a result, estimation of TNL models is computationally more demanding and numerically less stable than for standard NL, often requiring multi-stage procedures to achieve convergence. This complexity has limited the widespread use of TNL relative to simpler logit models, despite its strong behavioral appeal in contexts with naturally hierarchical choice structures. In this section, we show that, similar to the NL case, the hierarchical estimation structure of the TNL model can also be reformulated as an exponential cone program, thereby enabling efficient solution with modern conic optimization methods.
In the TNL model, a choice is represented as a path through a tree structure: starting at the root node, the decision-maker selects a sequence of intermediate nodes until reaching a leaf node, which corresponds to the chosen alternative. The tree-nested logit thus provides greater flexibility than the standard nested logit, as many sets of alternatives can be naturally organized hierarchically.
Figure (ref) illustrates an example of a three-level tree-nested logit structure. For example, in the context of transportation mode choice, the hierarchy may take the following form:
This hierarchical organization captures correlations at multiple levels of aggregation, allowing for more realistic substitution patterns between alternatives.
We now turn to the mathematical formulation of the TNL model. Let the tree have $T$ levels, with the root node at level 1 and internal nodes at levels $1$ through $T-1$ corresponding to nests. The leaf nodes at level $T$ correspond to the set of available alternatives. Let $\mathcal{N}$ be the set of all nodes and $\mathcal{I}$ be the set of all internal nodes. For any internal node $k \in \mathcal{I}$, let $C(k)$ denote the set of child nodes of $k$. Moreover, let $\mathcal{S}$ be the set of all leaf nodes, which correspond to the alternatives in the choice set. Let $r$ denote the root node of the tree structure.
Each internal node $k \in \mathcal{N}$ in the tree is associated with a scale (or dissimilarity) parameter $\lambda_k > 0$. For leaf nodes $k \in \mathcal{S}$, corresponding to the actual alternatives, we normalize $\lambda_k = 1$. The choice probability generating function (CPGF) Fosgerau2013 is then defined recursively as
where $v_k$ is the deterministic utility of alternative $k$, $\mathcal{I}$ is the set of internal nodes, and $C(k)$ denotes the set of child nodes of $k$. To ensure random utility maximization (RUM) consistency, we require $\lambda_s \geq \lambda_k$ for every parent--child pair $(k,s)$ Train2009. Since $\lambda_s=1$ for all $s \in \mathcal{S}$, this condition implies $\lambda_k \leq 1$ for all $k \in \mathcal{N}$.
The TNL choice probability can be decomposed into conditional probabilities of selecting a child node at each level of the tree. Specifically, for any internal node $k \in \mathcal{I}$ and child $s \in C(k)$, the probability of selecting $s$ given $k$ is \[ P(s \mid k) \;=\; \frac{V_s^{\lambda_s/\lambda_k}}{V_k}. \] Thus, the probability of reaching a particular alternative $i \in \mathcal{S}$ can be expressed as the product of conditional probabilities along the unique path from the root $r$ to $i$. If this path is denoted $\{k_1 = r, k_2, \ldots, k_T = i\}$, then \[ P(i \mid \mathcal{S}) \;=\; \prod_{t=1}^{T-1} P(k_{t+1} \mid k_t) = \prod_{t=1}^{T-1} \frac{V_{k_{t+1}}^{\lambda_{k_{t+1}}/\lambda_{k_t}}}{V_{k_t}}. \] We now describe MLE of the TNL model. Suppose we observe $N$ individuals. For each individual $n$, the data consists of a pair $(j_n, S_n)$, where $S_n \subseteq [m]$ is the offered choice set and $j_n \in S_n$ is the chosen alternative. Let $\{k^n_1 = r, k^n_2, \ldots, k^n_T = j_n\}$ denote the unique path from the root to $j_n$. The probability of observing choice $j_n$ is
where the values $V^n_k$ are computed recursively as \[ V^n_k =
\] The log-likelihood for parameters $(\boldsymbol{\beta}, \boldsymbol{\lambda})$ is therefore
This formulation highlights how the hierarchical structure of the TNL model decomposes the likelihood into contributions from each level of the decision tree.
We now reformulate the estimation of the TNL model under fixed scale parameters $\{\lambda_k\}_{k\in{\mathcal{N}}}$ as an ECP. Recall that the choice probability for observation $n$ and chosen alternative $j_n$ can be expressed as
where $\{k^n_1,\ldots,k^n_T=j_n\}$ is the unique path from the root to the chosen leaf node $j_n$. Note that $k^n_T$ is a leaf node, hence $\lambda_{k^n_T}=1$ and \[ V^n_{k^n_T} = \exp(\pmb{\beta}^\top \textbf{a}_{n j_n}). \] The log-likelihood can thus be written as
where $V^n_r$ denotes the value function at the root for individual $n$. Recall that the value functions $V^n_k$ are defined recursively as
Because $\lambda_s \geq \lambda_k$ for all parent--child pairs $(k,s)$, the coefficients of $\ln V^n_k$ in the log-likelihood are non-positive. Thus, in maximizing the log-likelihood, the optimizer will push $V^n_k$ to be as small as possible, forcing tightness at the optimum. Therefore, the equalities in (ref)--(ref) can be relaxed to inequalities:
We now let $z^n_k = \ln V^n_k$ for all $k \in {\mathcal{N}}$. The problem can then be reformulated as
Constraint (ref) can be rewritten as \[ 1 \;\geq\; \sum_{s\in C(k)} \exp\!\left(\tfrac{\lambda_s}{\lambda_k} z^n_s - z^n_k\right), \qquad k \in {\mathcal{I}},\; n\in[N]. \] We introduce auxiliary variables \[ y^n_{ks} = \exp\!\left(\tfrac{\lambda_s}{\lambda_k} z^n_s - z^n_k\right), \qquad \forall k\in {\mathcal{I}}, \; s \in C(k), \; n\in [N], \] so that the problem becomes the following exponential cone program:
This ECP formulation explicitly replaces the nested log-sum-exp constraints of the TNL likelihood with linear inequalities and exponential cone constraints, making the problem amenable to modern interior-point conic solvers (when all the scale parameters $\lambda_k$ are fixed).
We now discuss the computational complexity of solving (ref) (with fixed scale parameters) using interior-point methods. Analogous to the analysis for the MNL and NL models, we begin by defining the aggregate counts: \[ Z \;:=\; \sum_{n=1}^N |S_n|, \qquad A_n \;:=\; \{\,k\in{\mathcal{I}} : \text{the subtree of }k\text{ contains some } s\in S_n\,\}, \qquad \Gamma \;:=\; \sum_{n=1}^N |A_n|. \] Thus, $Z$ is the total number of leaf appearances across all observations and $\Gamma$ is the total number of active internal nodes across individuals (i.e., internal nodes whose subtrees intersect the offered set). Let $E_n$ be the set of active parent--child edges in the minimal subtree that spans $S_n$ and the root; then $|E_n| = |A_n|+|S_n|-1$, and the total number of active edges across observations is \[ E \;:=\; \sum_{n=1}^N |E_n| \;=\; \Gamma + Z - N . \] In the TNL-ECP (ref), the decision variables are: $\pmb{\beta}\in\mathbb{R}^p$ ($p$ vars), node logs $\{z_k^n\}$ for all active nodes ($\Gamma+Z$ vars), and edge auxiliaries $\{y^n_{ks}\}$ for all active edges ($E$ vars). Hence the model has \[ p \;+\; (\Gamma+Z) \;+\; E . \] variables. The linear inequalities are one constraint $1\ge \sum_{s\in C(k)} y^n_{ks}$ for each active internal node $k\in A_n$ (total $\Gamma$ constraints, adding $\Gamma$ nonnegative slacks in solver form) and the affine leaf bounds $z_k^n\ge \pmb{\beta}^\top \textbf{a}_{nk}$ for each $(n,k)$ with $k\in S_n$ (total $Z$ affine rows). The exponential-cone blocks consist of one block per active edge $(k,s)\in E_n$, so there are $E$ three-dimensional ${\mathcal{K}}_{\exp}$ blocks in total. In the worst case of full availability ($S_n={\mathcal{S}}$ for all $n$) on a tree with $m:=|{\mathcal{S}}|$ leaves and $q:=|{\mathcal{I}}|$ internal nodes, we have $\Gamma=Nq$ and $E=N(q+m-1)=\Theta(Nm)$.
The following proposition states the computational complexity of solving (ref) using a path-following interior-point algorithm NesterovNemirovskii1994.
The proof can be found in the appendix. The complexity formulation for the TNL model underscores the additional computational burden introduced by the hierarchical tree structure. Compared to the NL model, the per-iteration cost of the interior-point method now depends not only on the number of observations and choice set sizes, but also on the number of active internal nodes and edges in the decision tree. This reflects the richer substitution patterns captured by TNL, but also implies that estimation can become considerably more expensive for deep or wide trees. At the same time, the analysis shows that the problem remains polynomial-time solvable when scale parameters are fixed.
In this section, we present a general framework for estimating the MNL, NL, and TNL models, building on their ECP reformulations. When the scale parameters are fixed, as in the MNL model and in the NL and TNL models under fixed $\pmb{\lambda}$, the estimation problem is convex and can be solved efficiently to global optimality using off-the-shelf conic solvers such as MOSEK MOSEKcookbook.
In the more general case where the scale parameters $\{\pmb{\lambda}\}$ must be estimated jointly with the utility coefficients $\pmb{\beta}$, the problem becomes non-convex due to the nonlinear dependence of the likelihood on $\pmb{\lambda}$. To address this, we adopt a two-stepstep procedure that alternates between updating the scale parameters and optimizing the utility parameters:
For the outer update step, the gradient with respect to the parameter $\pmb{\lambda}$ can be obtained by differentiating the log-likelihood function while holding $\pmb{\beta}$ fixed at its optimal value. Specifically, let $ \pmb{\beta}^*(\pmb{\lambda}) = \text{argmax}_{\pmb{\beta}} \; {\mathcal{L}}^{\textsc{\tiny NL}}(\pmb{\beta}, \pmb{\lambda}),$ and define $ {\mathcal{L}}^{\textsc{\tiny NL}*}(\pmb{\lambda}) = {\mathcal{L}}^{\textsc{\tiny NL}}(\pmb{\beta}^*(\pmb{\lambda}), \pmb{\lambda}).$ According to the envelope theorem milgrom2002envelope, the gradient of ${\mathcal{L}}^{\textsc{\tiny NL}*}(\pmb{\lambda})$ with respect to $\pmb{\lambda}$ can then be expressed as
A similar gradient formulation can be derived for the TNL model. These gradients admit a closed-form expression and can be supplied to a nonlinear constrained optimization solver to efficiently solve the outer problem.
Figure (ref) illustrates the pseudo-code of our two-step procedure. This alternating two-step procedure decouples the non-convex estimation problem into a sequence of convex inner problems and constrained nonlinear outer updates. The inner conic programs leverage the computational advantages of exponential cone reformulations, while the outer updates adjust the scale parameters under their structural constraints. In practice, we iterate between the outer and inner steps until convergence of the log-likelihood or until successive updates fall below a prescribed tolerance. This framework provides a flexible and scalable estimation method applicable to the MNL, NL, and TNL models, encompassing both fixed and estimated scale parameter settings.
We conduct experiments to evaluate the performance of our conic optimization approach in comparison with standard gradient-based methods. The experiments are based on generated model structures and datasets of varying sizes, with the objective of assessing how the ECP reformulation performs, particularly in large-scale instances.
To examine the performance of the ECP reformulations, we randomly generate three sets of instances corresponding to the MNL, NL, and TNL models.
\paragraph{MNL instances.} We generate $36$ sets of MNL instances with parameter configurations $\text{dim}(\pmb{\beta}) \in \{5,10,20,50\}$, $(N,m) \in \{(500,50), (1000,100), (2000,200)\}$, and $|S_n| \in \{0.2m, 0.5m, 0.8m\}$. Each set consists of $5$ instances. For each instance, every element of the vector $\mathbf{a}_{ij}$ ($i \in N, j \in [m]$) is drawn independently from $U[0,3]$. The choice set $S_n$ and the chosen alternative $j_n$ are randomly selected from $[m]$, subject to the condition that $j_n \in S_n$.
\paragraph{NL instances.} For the NL formulation, we generate $72$ sets of instances using the same specifications for $\dim(\pmb{\beta})$, $(N,m)$, $|S_n|$, and $\mathbf{a}_{ij}$ as in the MNL case. In addition, the number of nests is set to $L \in \{2,5\}$.
\paragraph{TNL instances.} For the TNL model, we also generate $72$ sets of instances. Here we consider a tree structure with depth $T=4$, where the number of child nodes at levels $1$ and $2$ is chosen from $\{2,3\}$.
For cases where $\pmb{\lambda}$ is fixed, we draw the upper bound $u_{\lambda}$ from $U[0.8,0.9]$, the lower bound $l_{\lambda}$ from $U[0.1,0.2]$, and generated each $\lambda$ uniformly from the interval $[l_{\lambda}, u_{\lambda}]$.
We compare our approach with the nonlinear optimization solvers implemented in SciPy.optimize---a state-of-the-art library that provides gradient-based methods for minimizing (or maximizing) continuous nonlinear objective functions, possibly subject to constraints 2020SciPy-NMeth. The package includes solvers for general nonlinear problems and supports both local and global optimization algorithms.
To select the most suitable solver from the optimize module, we conducted a preliminary experiment on several instances of the MNL, NL, and TNL datasets with five attributes. The results (reported in the Appendix) indicate that L-BFGS-B and SLSQP achieve the lowest runtime, number of iterations (nit), and number of function evaluations (nfev) among the solvers available in SciPy.optimize. Consequently, \texttt{L-BFGS-B} is selected as the baseline solver for estimating the MNL, NL, and TNL models with fixed ${\pmb{\lambda}}$, while \texttt{SLSQP} is employed for estimating the NL and TNL models (jointly estimating $\pmb{\beta}$ and ${\pmb{\lambda}}$) due to its ability to handle constraints on $\pmb{\lambda}$.
All methods in our experiments are implemented in Python, using MOSEK 11.0 and SciPy 1.15.2, with a solving time limit of one hour per run. The experiments are conducted on a PC equipped with an Intel(R) Core(TM) i7-9700 CPU @ 3.00GHz, 16 GB of RAM, and the Windows 11 operating system.
We present comparative results between our ECP-based method and gradient-based baselines across various estimation instances of the MNL, NL, and TNL models, considering both the fixed $\pmb{\lambda}$ case and the joint estimation of $(\pmb{\beta}, \pmb{\lambda})$.
Table (ref) compares the performance of L-BFGS-B and the proposed ECP solver across the MNL, NL, and TNL models (with fixed ${\pmb{\lambda}}$). The instances are grouped by numbers of attributes, $\{5,10,20,50\}$ (denoted by “\#Att”), and dataset sizes $(N,m) \in \{(500,50), (1000,100), (2000,200)\}$ (denoted by $\{S, M, L\}$). The detailed results of each set in 36 sets are visualized in the Appendix (ref). In the table, the columns “\#Opt” report the number of optimal solutions found by each method for a given set of instances, while the columns “AveTime(s)” present the average solving time in seconds. The best results in each row are highlighted in bold.
The results demonstrate that ECP consistently outperforms L-BFGS-B in both solution reliability and computational efficiency. Specifically, for the {MNL instances,} across all $180$ instances, ECP achieves a perfect optimization rate of $180/180$, compared to $168/180$ for L-BFGS-B. For the {NL instances,} ECP solves all $360$ instances successfully, whereas L-BFGS-B solves only $177$. For the {TNL instance,} ECP again achieves a perfect success rate of $360/360$, while L-BFGS-B solves $256$ instances. In terms of runtime, ECP is significantly faster. In the NL model with $50$ attributes and large dataset size ($L$), ECP requires on average only $26.95$ seconds, compared to $3372.16$ seconds for L-BFGS-B. In the TNL model under the same configuration, ECP completes in $6.39$ seconds, while L-BFGS-B requires $2820.66$ seconds.
Overall, these results highlight the superior scalability and robustness of ECP, establishing it as the preferred method for efficiently solving large-scale and nested logit models.
Figures (ref) illustrate the average solution times of the optimization methods for estimating MNL, NL, and TNL instances (with fixed $\pmb{\beta}$). The results show that ECP consistently outperforms L-BFGS-B, particularly on medium and large instances, where the runtime of L-BFGS-B increases sharply.
Table (ref) presents the comparative performance of three estimation approaches---L-BFGS-B, Mixed L-BFGS-B, and ECP+L-BFGS-B---for estimating the NL model across varying sizes of the parameter vector $\pmb{\beta}$ (column “\#Att”), the number of nests (column “\#Nest”), and dataset sizes (see Appendix (ref) for more details of each set in 72 sets). For the L-BFGS-B baseline, we directly apply the solver to estimate the MLE. The Mixed L-BFGS-B baseline adapts our two-stage procedure (Section (ref)), where both stages are solved using L-BFGS-B. Finally, our proposed method, \texttt{ECP+L-BFGS-B}, employs \texttt{L-BFGS-B} in the outer step and the ECP solver in the inner step to solve the ECP reformulations. Because global optimality cannot be guaranteed for the joint estimation problems, the columns “\#Opt” are replaced by “\#Best,” which report the number of best solutions found. In addition, the columns “AveGap(%)” present the average objective gaps of \texttt{Mixed L-BFGS-B} and \texttt{ECP+L-BFGS-B} relative to \texttt{L-BFGS-B}.
The results show that our approach ECP+L-BFGS-B approach consistently achieves the highest number of best solutions ($342/360$), far surpassing Mixed L-BFGS-B ($262/360$) and L-BFGS-B alone ($136/360$). Moreover, ECP+L-BFGS-B not only attains superior solution quality but also offers faster computation times in nearly all settings. For example, with $50$ attributes, $5$ nests, and a large dataset, ECP+L-BFGS-B finds all $15$ best solutions in an average of $1205.79$ seconds, compared to $3600.00$ seconds for Mixed L-BFGS-B (which finds only one best solution) and $3254.23$ seconds for \texttt{L-BFGS-B} (which finds only two).
In terms of average gaps, ECP+L-BFGS-B again demonstrates a clear advantage, particularly in high-dimensional settings. For instance, with $50$ attributes and large datasets, it achieves average gaps of $70.61\%$ and $56.00\%$, significantly outperforming Mixed L-BFGS-B, which fails to yield meaningful gaps (reported as $-\infty$).
Overall, these results demonstrate that our approach provides a powerful and scalable solution for the NL estimation problem, excelling in both efficiency and reliability.
Table (ref) summarizes the performance of three optimization strategies on the TNL dataset across various problem sizes and tree structures. We consider the following methods:
The results demonstrate that ECP+SLSQP clearly outperforms the other methods, achieving the highest number of best solutions ($340/360$), with the largest average gaps and competitive runtimes. While SLSQP alone is efficient on small problems, it fails to maintain solution quality as the problem size increases. L-BFGS-B+SLSQP improves accuracy but often fails to converge or exceeds time limits on large instances, leading to infeasible or divergent solutions (e.g., infinite or negative gaps). In contrast, ECP+SLSQP remains stable and efficient even in the most challenging settings, demonstrating superior scalability, robustness, and solution quality.
Figure (ref) presents the results for joint estimation of the NL and TNL models, comparing the SciPy solvers with the exponential cone reformulation approaches. In both cases, ECP-based methods achieve significantly lower runtimes and exhibit superior scalability, clearly outperforming the SciPy solvers as model complexity increases.
Our experiments comprehensively evaluate the proposed exponential cone reformulation approach against standard gradient-based solvers from SciPy across the MNL, NL, and TNL models. The results consistently demonstrate the superiority of ECP in terms of solution quality, reliability, and scalability. For the MNL model, ECP achieves a perfect optimization rate across all instances, while L-BFGS-B (best gradient-based method for the problem settings) fails on a subset of cases. The performance gap becomes even more pronounced for the more complex NL and TNL models: ECP-based methods reliably solve all instances, whereas the gradient-based solvers frequently fail to converge, return infeasible solutions, or hit time limits. In terms of runtime, ECP is not only significantly faster on medium and large-scale instances, but also scales more gracefully with increasing numbers of attributes, nests, and dataset sizes. Furthermore, for the joint estimation problems, our approaches such as ECP+L-BFGS-B and ECP+SLSQP leverage the strengths of the ECP solver, offering the best balance between accuracy and efficiency in joint estimation problems. Overall, the results highlight that ECP-based reformulations provide a robust and scalable framework for large-scale logit model estimation, outperforming state-of-the-art gradient-based solvers.
We revisited parameter estimation for MNL, NL, and TNL models through convex conic optimization, showing that their maximum-likelihood problems can be reformulated as ECPs. Building on this insight, we analyze the polynomial-time computational complexity of solving the equivalent ECP reformulations for the MNL model, as well as for the NL and TNL models with fixed scale parameters. we proposed a two-stage procedure that alternates between updating scale parameters and solving ECPs for utility coefficients. Experiments on synthetic datasets demonstrate that ECP-based methods consistently outperform gradient-based solvers in terms of likelihood values, robustness, and runtime, especially on large-scale and high-dimensional instances. These results establish exponential-cone programming as a practical, scalable, and reliable alternative for discrete-choice model estimation.
A natural extension is the development of conic optimization techniques for the global estimation of nested logit, cross-nested logit model, or generalized network-based choice models Train2009,Bier06, Mai16. Unlike the fixed-parameter setting studied here, these formulations involve jointly optimizing over both utility coefficients and dissimilarity parameters, leading to highly non-convex landscapes. Additional challenges include ensuring identifiability of parameters, handling cross-nesting structures where alternatives belong to multiple nests, and designing scalable algorithms that can cope with the resulting high-dimensional optimization problems. Addressing these challenges would mark a significant step toward robust and efficient estimation of the most general forms of discrete-choice models.