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.
109,088 characters · 15 sections · 65 citation commands
Equilibrium-Constrained Estimation of Recursive Logit Choice Models
\newcolumntype{C}{>{\arraybackslash}p{4em}}
\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:} {Recursive Logit, Dynamic Discrete Choice, Maximum Likelihood Estimation, Conic Optimization, Transportation Networks}
The recursive logit (RL) model, introduced by FosgFrejKarl13, has emerged as a powerful framework for modeling sequential decision-making processes on networks. In transportation research, the RL model has been widely applied to analyze route choice behavior, capturing how travelers make link-by-link decisions that ultimately determine complete paths. Beyond route choice, RL-type formulations have proven useful in activity-based travel demand modeling, where individuals plan sequences of activities throughout the day zimmermann2018mixedRL, as well as in multiple discrete choice settings, where selections of composite alternatives can be represented as paths on suitably defined directed acyclic graphs tran2024network. By decomposing complex combinatorial choice problems into recursive structures, the RL model provides both theoretical elegance and practical tractability, making it a cornerstone in the field of behavioral modeling and prediction.
Despite these advantages, the practical application of RL models critically depends on estimation methods. Existing approaches overwhelmingly rely on the nested fixed-point (NFXP) algorithm Rust87,Rust94, in which model parameters are estimated by repeatedly solving a dynamic program to evaluate the likelihood function. In practice, this requires computing the value function as the fixed point of a Bellman operator through iterative methods such as value iteration. However, this procedure is computationally demanding and often numerically unstable. When the Bellman recursion is ill-conditioned, value iteration can fail to converge, diverging instead to infinity or producing undefined values. As a result, the computed log-likelihood may become invalid, leading to failed estimation runs or unreliable parameter estimates FosgFrejKarl13,MaiFrejinger2022. These limitations represent a serious bottleneck for large-scale applications, where estimation speed and numerical robustness are critical.
To address these challenges, we propose a novel reformulation of the RL maximum likelihood estimation problem. Our key idea is to represent the likelihood estimation task as an optimization problem with equilibrium constraints, in which both the structural parameters and the expected downstream utilities (i.e., value functions) are treated as joint decision variables. This formulation allows us to transform the problem into a convex conic program, specifically an exponential-cone optimization model, that can be solved efficiently by modern commercial solvers such as MOSEK mosek2023. By avoiding repeated fixed-point computations and instead embedding equilibrium conditions directly into the optimization problem, our approach provides a robust and scalable alternative to classical NFXP.
\paragraph{Contributions.} This paper makes several methodological and practical contributions to the literature on RL model estimation:
Overall, this work is the first to show that the estimation the RL model can be cast as an exact convex optimization problem, thereby replacing decades-old reliance on unstable fixed-point iterations with a theoretically principled and computationally efficient alternative. This shift has both methodological significance—connecting dynamic discrete choice estimation with conic optimization—and practical impact, enabling stable large-scale estimation in domains where RL models have so far been under-utilized.
\paragraph{Paper Outline:} The remainder of this paper is organized as follows. Section (ref) reviews relevant literature on route choice modeling, recursive logit estimation methods, and related optimization-based approaches. Section (ref) introduces the RL framework and revisits its estimation using the classical NFXP algorithm. In Section (ref), we present the proposed equilibrium-constrained reformulation and discuss its theoretical properties, including convexity and equivalence to the standard RL formulation. Section (ref) reports numerical experiments on both synthetic and real transportation networks, highlighting the computational efficiency and stability of the proposed approach compared to NFXP. Finally, Section (ref) concludes the paper and discusses potential directions for future research. Appendix extends the framework to the Nested Recursive Logit model MaiFosFre15.
Route choice models play a central role in transportation research and travel behavior analysis, with applications ranging from network design and congestion management to demand forecasting and intelligent transportation systems BenABier99,Trai03,Domencich1975. The literature on route choice modeling can broadly be categorized into two main approaches: path-based and link-based recursive formulations. Path-based models Prat09 rely on sampling feasible paths between each origin--destination pair, which makes the resulting parameter estimates sensitive to the path sampling procedure. Even when correction terms are introduced to ensure consistent estimation, prediction remains computationally demanding, particularly in large-scale networks where the number of feasible paths grows exponentially.
In contrast, link-based RL models Fosgerau2013,MaiFosFre15,mai2021RL_STD,oyama2017discounted build upon the dynamic discrete choice framework of Rust87 and are mathematically equivalent to discrete choice models defined over the set of all feasible paths. These models offer several advantages: they can be consistently estimated without path enumeration, and they enable efficient prediction via dynamic programming. The first RL model was introduced by Fosgerau2013 and has since been extended in multiple directions, including modeling correlated utilities across overlapping paths MaiFosFre15,Mai_RNMEV,mai2018decomposition, accounting for dynamic network conditions de2020RL_dynamic, incorporating stochastic time-dependent link costs mai2021RL_STD, capturing discounted behavior oyama2017discounted, or modeling multiple-discrete choice behaviors tran2024network. RL models have been successfully applied in various transportation and network design contexts, including traffic management BailComi08,Melo12, network pricing zimmermann2021strategic, and network interdiction mai2024stackelberg. A comprehensive overview of these developments and related estimation techniques is provided in zimmermann2020tutorial.
The RL model belongs to the class of dynamic discrete choice models, where estimation is typically performed using the NFXP algorithm Rust87. In this framework, the inner loop solves the Bellman equation to obtain the value function, while the outer loop maximizes the likelihood with respect to model parameters. Despite its conceptual simplicity, several studies have noted that NFXP is highly sensitive to initialization, scaling, and network topology, which often lead to numerical instability or infeasible value function estimates in practice Fosgerau2013,mai2018decomposition,MaiFrejinger2022. Several efforts have been made to improve the NFXP procedure for RL model estimation. For instance, mai2018decomposition proposed a decomposition-based method that can reduce estimation time by up to a factor of thirty. However, this approach cannot accommodate individual-specific attributes, which limits its applicability in general route choice settings. More recently, MaiFrejinger2022 introduced a least-squares-based technique to accelerate the computation of the value function within the NFXP framework. Despite its computational efficiency, this method suffers from the same numerical instability as standard NFXP—specifically, the value function may fail to converge during the parameter search process, leading to incomplete or invalid likelihood evaluations. These challenges motivate the search for alternative estimation frameworks that preserve the theoretical rigor of the RL model while improving numerical robustness.
Our work contributes to this line of research by connecting route choice estimation with the field of conic optimization. Specifically, we propose an equilibrium-constrained optimization formulation that casts the RL estimation problem as a convex conic program. Compared to standard gradient-based optimization approaches, conic formulations offer several advantages, including guaranteed global optimality under convexity, well-established interior-point solution methods, and improved numerical stability in ill-conditioned problems BoydVandenberghe2004,nesterov2006conic,mosek2023. This approach bridges dynamic discrete choice modeling and modern convex optimization, providing a new computational perspective for stable and scalable estimation of RL models.
It is worth emphasizing that our proposed equilibrium-constrained formulation framework is conceptually related to the Mathematical Programming with Equilibrium Constraints (MPEC) approaches developed for structural and dynamic discrete choice estimation SuJudd2012,IskhakovEtAl2016,qeconomics2023. Despite this connection, our method differs fundamentally in both formulation and scope. MPEC methods jointly estimate model parameters and value functions by imposing the Bellman equation as an equilibrium constraint. However, these formulations typically preserve the nonlinear and nonconvex structure of the original problem, resulting in large-scale constrained nonlinear programs that are computationally demanding and often difficult to solve to global optimality. In contrast, our convex formulation leverages the specific analytical structure of the RL model to derive a conic representation of the equilibrium conditions. This convexification guarantees global optimality, enhances numerical stability, and enables the use of efficient interior-point solvers such as MOSEK. Moreover, as a side note, prior studies have shown that MPEC approaches are generally less computationally efficient than traditional NFXP algorithms, particularly in large-scale settings where the number of equilibrium constraints increases rapidly with the state space IskhakovEtAl2016. In contrast, as shown extensively in our experiments, our proposed method achieves both numerical robustness and computational tractability, offering a practical and theoretically grounded alternative to the NFXP.
The recursive logit (RL) model, introduced by FosgFrejKarl13, is a dynamic discrete choice framework designed to represent sequential decision-making on a transportation network. Unlike classical multinomial logit formulations that treat entire paths as atomic alternatives, the RL model decomposes path choice into a sequence of link- or node-level decisions, thereby avoiding the combinatorial explosion of path enumeration.
Let ${\mathcal{S}}$ denote the set of states, where each state $s \in {\mathcal{S}}$ corresponds to a link or node in the network. For each state $s$, let $A(s) \subseteq {\mathcal{S}}$ represent the set of directly reachable successor states. For a trip with destination $d$, the destination is modeled as an absorbing state, so the effective state space is $\widetilde{{\mathcal{S}}} = {\mathcal{S}} \cup \{d\}$. At each step, a traveler located in state $s$ chooses a successor state $s' \in A(s)$. The instantaneous random utility associated with this transition is given by
where $v(s'|s)$ denotes the deterministic component of utility, typically expressed as a linear function of observable attributes (e.g., travel time, cost, or congestion) with parameters to be estimated, $\boldsymbol{\beta}$. Hence, we can write this component as $v(s'|s;~\boldsymbol{\beta})$. The term $\epsilon(s')$ represents an i.i.d.\ extreme value type I disturbance, and $\mu > 0$ is a scale parameter. At each state $s \in {\mathcal{S}}$, the individual chooses the next state by maximizing the sum of two components: (i) the instantaneous random utility of moving from $s$ to $s' \in A(s)$, and (ii) the expected maximum utility of continuing from $s'$ to the destination $d$. Formally, the next-state decision is characterized by
where $V(s')$ denotes the expected maximum utility of reaching $d$ starting from state $s'$. By the dynamic programming principle, the value function $V(s)$ can be expressed as
where the expectation is taken with respect to the random shocks $\{\epsilon(s')\}_{s' \in A(s)}$.
In the RL model, the error terms $\epsilon(s')$ are assumed to be i.i.d.\ extreme value type I (Gumbel). Under this distributional assumption, the maximum operator in (ref) admits a closed-form representation, and the resulting conditional choice probabilities follow a multinomial logit structure. Consequently, the value function satisfies the log-sum-exp recursion:
Equation (ref) illustrates how the RL model inherits the tractability of the multinomial logit: the continuation value at each state is obtained as a smooth “log-sum-exp” operator of the utilities of feasible successor states. This recursive structure eliminates the need to enumerate paths explicitly and ensures that choice probabilities are consistent with random utility maximization.
Moreover, conditional on being at state $s$, the probability of choosing successor $s' \in A(s)$ is
This expression illustrates the recursive nature of the model: the attractiveness of an action depends not only on its immediate attributes but also on the continuation value $V(s')$.
For model estimation, consider an observed path $\sigma = \{s_0, s_1, \ldots, s_T\},$ which represents a sequence of states from an origin $s_0$ to the destination $d$. The probability of observing this path under parameter vector $\pmb{\beta}$ is the product of the conditional choice probabilities at each step:
Using the recursive logit structure, this probability can be expressed in closed form as
where $v(\sigma) = \sum_{t=0}^{T-1} v(s_{t+1} \mid s_t)$ is the total deterministic utility accumulated along the path $\sigma$, and $V(s_0)$ is the value function at the origin state. Given a dataset of $N$ observed paths $\{\sigma_n\}_{n=1}^N$, the log-likelihood function for parameter estimation is
where $s_0^n$ denotes the origin of the $n$-th observed path, and $V^n(s)$ represents the expected maximum utility of reaching the destination starting from state $s$, corresponding to observation $n$. The model parameters $\pmb{\beta}$ can then be estimated by maximum likelihood estimation (MLE), which consists of maximizing the log-likelihood function $\mathcal{L}(\pmb{\beta})$ with respect to $\pmb{\beta}$.
\paragraph{Estimation via NFXP.} A central difficulty in maximizing the log-likelihood function in (ref) lies in the fact that the value functions $\{V(s)\}$ are not given in closed form but are instead implicitly defined through the fixed-point system (ref). For a given parameter vector $\pmb{\beta}$, solving this system may not always succeed: depending on the parameter values, the fixed-point equations can be ill-conditioned, which may lead to numerical instability or even the absence of a valid solution. This lack of robustness makes the computation of $V(s)$ itself a non-trivial task.
The standard estimation strategy is the NFXP algorithm Rust87. In this approach, the estimation problem is organized into two loops. In the outer loop, a candidate parameter vector $\pmb{\beta}$ is proposed, and the corresponding likelihood $\mathcal{L}(\pmb{\beta})$ must be evaluated. In the inner loop, for the given $\pmb{\beta}$, the system (ref) is solved to recover the value functions $\{V(s)\}$, typically through iterative methods such as value iteration. Once $V(s)$ is obtained, the outer loop updates $\pmb{\beta}$ using a numerical optimization routine, for example quasi-Newton methods.
While this nested structure guarantees statistical consistency, it comes with a significant computational burden. Each evaluation of the likelihood requires repeatedly solving a high-dimensional dynamic program defined over the entire network. When this process must be performed many times across successive updates of $\pmb{\beta}$, the computational cost grows prohibitively large. This combination of potential instability in computing value functions and the repeated fixed-point evaluations forms the main challenge of applying the RL model in practice. These difficulties motivate the search for alternative formulations, such as the equilibrium-constrained optimization approach developed in the next section.
We now reformulate the maximum likelihood estimation problem as an optimization problem with equilibrium constraints. The key idea is to treat both the model parameters $\pmb{\beta}$ and the value functions $\{V^n(s)\}$ as decision variables to be jointly estimated. For notational convenience, we denote the value function $V^n(s)$ by $V^n_s$, moving the state index to a subscript. The destination $d$ also depends on the observed path $n$, but the path index is omitted for clarity. This representation allows us to handle $\pmb{\beta}$ and $V$ more naturally as vectors of optimization variables in the reformulated problem.
We write the MLE as the following equilibrium-constrained optimization problem:
In this formulation, the equilibrium conditions defining the value functions appear explicitly as equality constraints. The goal is therefore to simultaneously optimize over the structural parameters $\pmb{\beta}$ and the value functions $\{V^n_s\}$. However, problem (ref) is highly impractical in its current form: the constraints are nonlinear equalities involving exponential terms, which makes the optimization problem non-convex and computationally challenging. In the next section, we describe our approach to reformulate (ref) into a tractable convex optimization problem. This convex reformulation allows the problem to be efficiently solved using modern commercial solvers.
We now describe our approach to reformulate the non-convex equilibrium-constrained problem (ref) into a convex optimization problem. The central idea is to relax the nonlinear equality constraints into inequalities that preserve the equilibrium structure while ensuring tractability. To make this reformulation valid, we impose two mild assumptions.
Assumption (ref) guarantees that all states are reachable from the observed origins. This condition is natural in most transportation or choice network applications, where the network is typically strongly connected or can be restricted to its reachable component.
Assumption (ref) ensures that the estimation problem is well-posed, i.e., the optimization problem has a feasible solution that maximizes the likelihood. This is a standard regularity condition in maximum likelihood estimation: without it, the optimization could be ill-defined (for instance, if the likelihood is unbounded or the feasible region is empty). In practice, this assumption is usually satisfied under mild conditions on the utility specification and the observation data.
The following theorem establishes a key result that enables us to overcome the nonlinearity of (ref): it shows that the equality constraints in (ref) can be safely relaxed to inequalities without changing the optimal solution. This relaxation not only preserves the theoretical integrity of the model but also paves the way for a convex reformulation that can be efficiently solved using modern optimization techniques.
For notational convenience, let us first define the Bellman operator $\mathcal{T}_{\boldsymbol{\beta}}$ acting on a vector of values $\boldsymbol{V}^n$ as \[ (\mathcal{T}_{\boldsymbol{\beta}}[\boldsymbol{V}^n])_s =
\] To prove the theorem, we first introduce the following lemma, which shows that the Bellman operator $\mathcal{T}_{\boldsymbol{\beta}}[\boldsymbol{V}^n]$ is monotonic in $\boldsymbol{V}^n$. This property is critical for relaxing the equality constraint in the theorem.
We now return to the proof of Theorem (ref). The general idea is to show that the equality constraint in the original formulation can be relaxed to an inequality constraint without altering the optimal solution. This is achieved by leveraging the monotonicity of the Bellman operator $\mathcal{T}_{\boldsymbol{\beta}}$, established in Lemma (ref), which ensures that iterative applications of $\mathcal{T}_{\boldsymbol{\beta}}$ preserve the order of value functions. Consequently, any feasible solution satisfying the relaxed inequality will converge to the same fixed point as that of the equality-constrained formulation.
To further elucidate the implications of Theorem (ref), we now provide a series of remarks that analyze specific scenarios in which the equivalence between the equality- and inequality-constrained formulations either holds, breaks down, or requires additional structural conditions. These discussions aim to underscore the critical role of Assumption (ref) (network conectivity) in ensuring the validity of the relaxation and illustrate the pathological cases that may arise when this assumption is violated, thereby deepening our understanding of the theoretical framework.
We now present our approach to further reformulate the inequality-constrained formulation in (ref) into a convex optimization problem that can be solved efficiently by modern commercial solvers. The key step is to introduce auxiliary variables \[ Q^n_{ss'} = v(s'|s;\pmb{\beta}) + V^n_{s'}, \qquad \forall n\in[N], ~ s\in {\mathcal{S}}, ~ s'\in A(s). \] These variables capture the one-step transition utilities, combining both the deterministic part of the utility and the continuation value. With this definition, the value function recursion can be written in the compact form \[ V^n_s = \log\!\left(\sum_{s'\in A(s)} \exp(Q^n_{ss'})\right). \] This representation makes it possible to directly embed the equilibrium constraints into the optimization problem. Substituting the expression above into (ref) leads to the following formulation:
It is important to emphasize that (ref) is a convex optimization problem. The concavity of the objective function follows from the fact that $-\log \sum \exp$ is concave, while $v(\sigma_n \mid \pmb{\beta})$ is linear in $\pmb{\beta}$. On the other hand, the constraints consist of inequalities of the form \[ Q^n_{ss'} \;\geq\; v(s'|s;\pmb{\beta}) + \log\!\left(\sum_{t\in A(s')} \exp(Q^n_{s't})\right), \] whose right-hand side is convex in $(\pmb{\beta},\textbf{Q})$. This guarantees that the feasible set is convex. Together, these two properties imply that the problem is convex and therefore admits a globally optimal solution. This convex formulation can in principle be solved by generic gradient-based convex optimization solvers such as projected gradient descent or quasi-Newton methods (e.g., L-BFGS). These methods guarantee convergence to the global optimum under mild conditions. However, they may suffer from slow convergence, particularly in large-scale networks where the log-sum-exp terms create ill-conditioned subproblems. Moreover, gradient-based approaches typically require careful step-size tuning and may be sensitive to scaling of the problem data. These limitations motivate the use of more specialized reformulations, in particular casting the problem as an exponential cone program, as discussed in the following.
\paragraph{Exponential Cone Reformulation} An additional advantage of the reformulation in (ref) is that the log-sum-exp terms can be represented exactly using exponential cone constraints. Recall that the exponential cone is defined as MOSEKcookbook,DahlAndersenExpCone2022: \[ \mathcal{K}_{\exp} = \Big\{(x,y,z)\in \mathbb{R}^3 : y>0, \; y \exp(x/y) \leq z\Big\} \;\cup\; \Big\{(x,0,z) : x \leq 0, \; z \geq 0\Big\}. \]
The convex formulation in (ref) involves multiple log-sum-exp expressions, each of which can be equivalently modeled using exponential cone constraints. To make this explicit, we introduce additional epigraph variables
which capture the log-sum-exp values at each state. With these variables, problem (ref) can be equivalently written as
Here we note that the equality definition in (ref) can be safely replaced by the inequality constraints in (ref). The reason is that the objective function decreases monotonically in $u^n_{s^n_0}$, so at optimality there is no incentive to make $u^n_s$ larger than the log-sum-exp value. In other words, any slack in (ref) would strictly worsen the objective, and hence all inequalities must bind at the optimum. This guarantees that the relaxation is exact and that (ref) is equivalent to (ref).
The constraints in (ref) are of the generic form \[ x \;\geq\; \log\!\left(\sum_{i=1}^m \exp(w_i)\right). \] It is standard in conic optimization that such log-sum-exp constraints can be reformulated exactly using exponential cones. Specifically, we have the equivalence \[ x \;\ge\; \log\!\Big(\sum_{i=1}^m e^{w_i}\Big) \quad\Longleftrightarrow\quad
\] Indeed, the condition $(w_i-x,\,1,\,r_i)\in \mathcal{K}_{\exp}$ implies that $r_i \ge e^{w_i-x}$. Summing over $i$ yields \[ \sum_{i=1}^m r_i \;\ge\; \sum_{i=1}^m e^{w_i-x} \;=\; e^{-x} \sum_{i=1}^m e^{w_i}. \] Thus the constraint $\sum_i r_i \le 1$ enforces $e^{-x} \sum_{i=1}^m e^{w_i} \;\le\; 1,$ which is equivalent to $x \ge \log \sum_i e^{w_i}$. This provides an exact exponential-cone representation of the log-sum-exp inequality. Applying this transformation to (ref), we obtain the following exponential-cone constraints using auxiliary variables $r^n_{ss'}$:
Collecting everything, the exponential-cone reformulation of (ref) is:
Problem (ref) is a convex conic program: the objective is linear in $(\pmb{\beta},\textbf{Q},\textbf{u})$, the coupling constraints are affine, and the only nonlinearities arise through memberships in the exponential cone $\mathcal{K}_{\exp}$, which is convex. By construction of the epigraph reformulation, (ref) is equivalent to (ref), and hence also equivalent to (ref). Therefore, any global optimum returned by a modern exponential-cone solver (e.g., MOSEK) is a globally optimal estimator for the RL MLE problem.
The main advantage of solving the exponential-cone program (ref), compared to generic gradient-based convex optimization methods, is its numerical robustness and guaranteed polynomial-time complexity. Interior-point solvers nesterov1994interior,mosek2023 for conic programs can exploit the exact conic structure of the log-sum-exp constraints, avoid issues of ill-conditioning common in gradient descent, and typically converge in far fewer iterations. In addition, solvers such as MOSEK MOSEKcookbook handle scaling automatically and return certificates of optimality, whereas first-order methods often require delicate step-size tuning and may converge slowly on large-scale networks.
\paragraph{Computational Complexity of the Exponential-Cone Reformulation: } As mentioned earlier, one of the main advantages of the exponential-cone reformulation proposed in (ref) is that it can be efficiently solved using a modern interior-point method with polynomial-time complexity. In the following, we now delve deeper into this complexity analysis.
Interior-point methods nesterov1994interior are the state of the art for solving exponential-cone programs. For a generic conic optimization problem with $n$ decision variables and $m$ exponential-cone constraints, the iteration complexity is $\mathcal{O}\!\left(\sqrt{m}\,\log\frac{1}{\epsilon}\right),$ where $\epsilon$ denotes the required accuracy. Each iteration requires forming and solving a Newton system whose dimension is proportional to $n+m$. The cost of solving such a system via direct methods is typically $\mathcal{O}((n+m)^3)$, though structure-exploiting solvers can often reduce this to nearly quadratic complexity. Hence, the overall worst-case complexity of an interior-point solver for the exponential-cone reformulation is $\mathcal{O}\!\Big((n+m)^3\sqrt{m}\,\log(1/\epsilon)\Big).$
For problem (ref), the problem size scales with the network and the data as follows. Each observation $n$ introduces a set of variables $Q^n_{ss'}$ for every arc $(s,s')$, one variable $u^n_s$ for each state $s$, and auxiliary epigraph variables $r^n_{ss'}$ for each arc. Thus, the total number of decision variables is $n = \mathcal{O}\big(N(|{\mathcal{A}}|+|{\mathcal{S}}|)\big),$ where $N$ is the number of observed paths, $|{\mathcal{A}}|$ is the number of arcs in the network, and $|{\mathcal{S}}|$ is the number of states. The number of exponential-cone constraints is also proportional to $N|{\mathcal{A}}|$, since each arc contributes one cone constraint. Therefore, the solver complexity in terms of the problem parameters is on the order of $\mathcal{O}\!\left( (N(|{\mathcal{A}}|+|{\mathcal{S}}|))^3 \sqrt{N|{\mathcal{A}}|}\,\log(1/\epsilon)\right).$ It is worth emphasizing that gradient-based methods such as projected gradient descent or stochastic gradient descent have a per-iteration cost of $\mathcal{O}(N|{\mathcal{A}}|)$, since each iteration requires traversing all observed arcs to compute gradients of the log-sum-exp terms. Their iteration complexity to achieve an $\epsilon$-accurate solution is typically $\mathcal{O}(1/\epsilon)$, leading to an overall complexity of $\mathcal{O}\!\left(\frac{N|{\mathcal{A}}|}{\epsilon}\right)$ Bubeck2015. In practice, this can result in a very large number of iterations, especially when high accuracy is required or when the problem is ill-conditioned. By contrast, interior-point methods are usually shown to converge within a modest number of iterations—often, regardless of the target accuracy $\epsilon$, but at the expense of a much higher per-iteration computational cost BoydVandenberghe2004.
In summary, the exponential-cone reformulation yields a convex conic program with polynomial-time solvability and strong global optimality guarantees. The theoretical worst-case complexity is cubic in the number of variables and constraints, which can be demanding for very large-scale networks. Nevertheless, for moderate-sized networks, interior-point solvers such as MOSEK are highly effective and often far more reliable than first-order methods, especially when accuracy and robustness are critical.
Compared to the classical NFXP algorithm, solving the exponential-cone formulation offers several important advantages. The exp-cone approach is considerably more robust with respect to initialization, since feasibility is directly enforced by convex constraints, whereas NFXP requires carefully chosen initial parameters $\pmb{\beta}$ to ensure that the Bellman equations admit a solution throughout the estimation process MaiFrejinger2022. Moreover, NFXP is computationally demanding because it requires solving multiple systems of linear equations both to evaluate the value function $V$ and to compute its Jacobian, with the number of such systems growing proportionally to the dimension of $\pmb{\beta}$. In contrast, the computational burden of the exp-cone method depends primarily on the size of the network rather than the number of parameters, making it more scalable in that respect.
While the exponential-cone reformulation is robust and guarantees global optimality, the program (ref) can become prohibitively large in real-world transportation networks. In such settings, the solver must process a vast number of cone constraints, often resulting in heavy memory consumption and long computation times. To address these challenges, it is useful to apply additional pre-processing steps, such as network reduction or approximate formulations, in order to improve scalability. In the following, we introduce a network trimming procedure that helps reduce the problem size while maintaining the statistical validity of the estimation.
The idea is that, given a specific origin–destination pair $(o,d)$, not all states in the network are equally relevant. In particular, states that are far from both the origin and the destination are very unlikely to appear in feasible paths, and thus contribute negligibly to the likelihood. Removing such states from the optimization problem can reduce the size of the exponential-cone program significantly without affecting estimation outcomes. To operationalize this idea, we use the concept of network flow under a reference parameter vector $\pmb{\beta}^0$. We first choose $\pmb{\beta}^0$ such that the deterministic utilities $v(s'|s;\pmb{\beta}^0)$ are sufficiently negative (e.g., by scaling travel times or costs), which guarantees that the Bellman equation $\textbf{V} = {\mathcal{T}}_{\pmb{\beta}^0}[\textbf{V}]$ admits a well-defined solution. Based on this solution, we can compute approximate choice probabilities at each state $P_{\pmb{\beta}^0}(s'|s)$. Specifically, to compute flow values from the orginal state $o$, let $P_{\pmb{\beta}^0}$ denote the $|{\mathcal{S}}|\times|{\mathcal{S}}|$ transition matrix with entries \[ \big(P_{\pmb{\beta}^0}\big)_{s,s'} =
\] Let $F \in \mathbb{R}^{|{\mathcal{S}}|}$ be the vector of flow values, and let $e_o$ be the unit vector corresponding to the origin state $o$. Injecting one unit of flow at $o$ and propagating it through the network yields the fixed-point system $F = e_o + P_{\pmb{\beta}^0}^\top F.$ Rearranging gives the linear system $(I - P_{\pmb{\beta}^0}^\top) F = e_o.$ Provided the spectral radius of $P_{\pmb{\beta}^0}$ is strictly less than one, this system admits a unique finite solution BailComi08. Each entry $F(s)$ then represents the expected flow through state $s$, starting from origin $o$ under the reference policy $P_{\pmb{\beta}^0}$.
For the trimming criterion, states with flow values below a prescribed threshold $\varepsilon > 0$ are deemed irrelevant and can be removed from the network, along with their associated arcs: ${\mathcal{S}}_{\text{trim}} = \{ s \in {\mathcal{S}} : F(s) \geq \varepsilon \}$. The optimization problem is then restricted to the trimmed network $({\mathcal{S}}_{\text{trim}}, {\mathcal{A}}_{\text{trim}})$, where ${\mathcal{A}}_{\text{trim}}$ denotes the set of arcs between states in ${\mathcal{S}}_{\text{trim}}$. This results in a reduced exponential-cone program with significantly fewer variables and constraints. It is worth noting that, beyond reducing the size of the network, the trimming procedure preserves the structural property required in our convex reformulation. In particular, the trimmed network $({\mathcal{S}}_{\text{trim}},A_{\text{trim}})$ automatically satisfies Assumption (ref) (network connectivity with respect to the origin), even if the original network does not satisfy this assumption. We formalize this in the proposition below.
This trimming procedure has two important properties: (i) it preserves all high-probability states and arcs, ensuring that the resulting estimator remains accurate; and (ii) it substantially reduces the size of the convex reformulation in practice, particularly for large networks with many peripheral states. In applications, the choice of the threshold $\varepsilon$ provides a trade-off between computational tractability and approximation fidelity. Our later experimental results will demonstrate that combining the proposed network trimming procedure with the convex reformulation enables efficient estimation of the RL model on a large-scale, real-world network.
To evaluate the performance of the proposed convex reformulation for RL model estimation, we focus on two fundamental applications of the recursive logit framework: route choice modeling FosgFrejKarl13 and multiple discrete choice modeling tran2024network. These two domains feature networks with distinct structural characteristics, providing complementary test cases for assessing the effectiveness and generality of our approach. The proposed ECP method is compared against the standard NFXP algorithm Rust87, with the MOSEK conic solver mosek2023 used to handle the exponential-cone reformulation. Although several alternative methods have been proposed as substitutes for NFXP, they are not directly suitable for inclusion in our experiments. For instance, the value decomposition method of mai2018decomposition can substantially accelerate RL model estimation but cannot accommodate general route choice settings with individual-specific attributes. Similarly, the MPEC-based approaches developed for dynamic discrete choice models SuJudd2012,SchjerningNotesMPEC have been shown to be less computationally efficient than the standard NFXP, particularly in large-scale problems IskhakovEtAl2016. Consequently, our comparisons focus on the classical NFXP benchmark to ensure both fairness and interpretability in evaluating the advantages of the proposed convex reformulation.
All experiments are conducted on a Google Cloud TPU v3-8 server equipped with a 96-core Intel Xeon 2.00 GHz CPU. In the following, we present our experimental results for each of the two application domains.
Multiple discrete choice (MDC) problems Bhat05,Bhat08, tran2024network are central to a wide range of economic and behavioral applications, as individuals or households often make simultaneous selections of multiple alternatives rather than a single option. Examples include consumers purchasing a basket of goods during a shopping trip, families choosing several entertainment or leisure activities, or commuters selecting a combination of transport modes and transfer options. The analysis of MDC behavior is therefore critical in domains such as marketing, urban planning, and transportation policy, as it enables richer modeling of preferences, substitution patterns, and complementarities among alternatives.
Traditional approaches to MDC within the RUM framework, dating back to McFa81, typically model each combination of alternatives as a distinct composite choice. While theoretically sound, this formulation quickly becomes computationally intractable due to the combinatorial explosion in the number of possible choice sets. For instance, when $m$ elemental alternatives are available, the number of potential composite alternatives grows exponentially as $2^m$, rendering classical RUM estimation infeasible even for moderately sized $m$.
To address this scalability barrier, tran2024network recently proposed the Logit-based Multiple Discrete Choice (LMDC) model, a RUM-based framework that exploits a network-based representation of MDC problems. In LMDC, composite alternatives are represented as feasible paths in a DAG. Each node or link in the DAG corresponds to a decision regarding the inclusion or exclusion of an elemental alternative, and complete paths encode valid multi-item choice bundles. This construction transforms the original combinatorial problem into a RL route choice model on the DAG. Importantly, the recursive structure enables efficient computation of choice probabilities using RL techniques developed in the context of route choice modeling. As a result, the LMDC model preserves the theoretical rigor of the RUM framework while maintaining tractability and scalability in high-dimensional multiple discrete choice settings.
There are two DAG representations proposed for the LMDC problem, as illustrated in Figure (ref). Each topology has its own properties. BiC-DAG has significantly fewer arcs compared to MuC-DAG, which reduces the graph density but also puts more pressure on individual arcs and causes computational instability.
In this experiment, we compare the performance of the proposed equilibrium-constrained reformulation, solved using the MOSEK conic optimizer, against the classical NFXP method originally employed in tran2024network and prior studies on RL models. Recall that the NFXP approach estimates the parameters $\pmb{\beta}$ by repeatedly solving a system of linear equations via value iteration in an inner loop, while updating $\pmb{\beta}$ in an outer loop using gradient-based optimization. By contrast, our method reformulates the MLE problem as a convex exponential-cone program, which can be solved in a single stage using interior-point algorithms fujisawa1997interior. Specifically, we consider the following two estimation methods for comparison:
This comparison underscores a fundamental methodological trade-off. The NFXP algorithm relies on repeated evaluations of an inner fixed-point system, and in practice it may suffer from numerical instability and the absence of guaranteed convergence. By contrast, our conic optimization approach exploits a convex reformulation that ensures global convergence within polynomial time. In the following, we present comparative results on both synthetic and real-world datasets to illustrate these trade-offs in practice.
The synthetic dataset, together with all settings and instances, is taken from the original paper tran2024network. Problem instances vary by the number of elemental alternatives $m \in \{5,10,20,30,50\}$ and by the bounds $[L,U]$, which represent the minimum and maximum number of alternatives an individual may simultaneously consider. The lower bound is fixed as $L=0$ for all instances, allowing individuals to select nothing, while the upper bound $U$ is specified for each size: $U=3$ for $m=5$, $U \in \{3,5\}$ for $m=10$, $U \in \{5,10\}$ for $m=20$, $U \in \{10,20\}$ for $m=30$, and $U \in \{20,30\}$ for $m=50$. The elemental alternative values are identical to those in the original experiment. The “ground-truth” model is the LMDC framework with parameter vector $\pmb{\beta}^0 = [-0.5, -0.02, -0.1]$. For each instance, two sets of observations are generated: one estimation set of size $1000$ for parameter estimation, and one prediction set of size $300$ to evaluate the average log-likelihood achieved by the estimated parameters. For comparison, log-likelihood values closer to zero indicate better predictive performance.
It is worth noting that the DAG networks used in this application are cycle-free. Consequently, both the equality-constrained formulation (ref) and the exponential-cone reformulation (ref) always admit a unique optimal solution. In particular, solving the value functions in the NFXP framework is always successful under these settings.
\paragraph{Estimation time comparison.} Since both the NFXP and ECP formulations are theoretically equivalent, and the underlying networks are acyclic, the estimation process always converges successfully. Consequently, both methods yield identical parameter estimates. We therefore focus our comparison exclusively on estimation times. Table (ref) and Figure (ref) reports the average and standard deviation of the estimation times obtained by both approaches over 40 independent runs. The results clearly highlight the computational advantages of the proposed ECP method. Although the NFXP approach using the MuC-DAG representation is slightly slower than the BiC-DAG variant---as also observed in tran2024network---even the faster NFXP+BiC configuration remains significantly slower than our ECP approach using BiC-DAG. The superior efficiency of ECP in this setting can be attributed to the low-density topology of the BiC-DAG representation, which substantially reduces the size of the optimization problem. Conversely, when the high-density MuC-DAG is employed, the running time of ECP increases notably due to the larger number of arcs and decision states.
\paragraph{Scalability with respect to the number of model parameters.} We compare the scalability of the proposed ECP and the classical NFXP method in terms of the dimensionality of the parameter vector $\boldsymbol{\beta}$. As illustrated in Figure (ref), the estimation time of NFXP increases almost linearly as the number of parameters grows, whereas ECP remains largely unaffected. This contrast arises from their computational structures: NFXP must repeatedly solve a sequence of linear systems for both value function evaluation and Jacobian computation, with the number of such systems scaling proportionally to $\dim(\boldsymbol{\beta})$. In contrast, ECP incorporates the parameters as additional decision variables within a convex conic program whose size depends primarily on the network topology rather than on $\boldsymbol{\beta}$. Consequently, increases in parameter dimensionality introduce negligible computational overhead for ECP, as the underlying constraint structure and sparsity pattern remain fixed. Empirical results consistently confirm this behavior across all tested datasets: NFXP exhibits a steady increase in runtime with larger parameter sizes, occasionally showing minor fluctuations due to multithreading optimizations in low-level numerical libraries, while ECP maintains stable and consistent computation times. Overall, these findings demonstrate that ECP achieves superior scalability with respect to $\boldsymbol{\beta}$, making it well-suited for high-dimensional estimation problems.
We further evaluate the practical performance of the proposed ECP against the classical NFXP approach using two real-world datasets---the Jewelry-Store and Book-Crossing datasets---provided in tran2024network. Both experiments are conducted under identical experimental settings. The Jewelry-Store dataset contains purchase records from an online jewelry retailer collected over a three-year period. It includes 436 elemental alternatives (products), each characterized by four attributes: Quantity, Price, Gender, and a constant term equal to one across all products. The estimation set consists of 34,291 observations. Table (ref) summarizes the estimation results. The estimated parameters obtained by NFXP and ECP are nearly identical, confirming the theoretical equivalence established in earlier sections. The slight difference in log-likelihood values is attributed to minor numerical rounding effects inherent in the value iteration computations. In contrast, the computational efficiency difference is substantial: ECP achieves convergence in only 28 seconds, compared to 333 seconds for NFXP, representing more than an order-of-magnitude improvement in runtime.
The second dataset, Book-Crossing, contains four weeks of book purchase histories and user ratings from Amazon. This dataset comprises 100 best-selling books, each described by four attributes: Popularity, Rating, Lifespan, and a constant term. The sample includes 16,989 observations. As shown in Table (ref), both methods yield identical parameter estimates and log-likelihood values, demonstrating the robustness of the ECP approach. However, ECP again achieves a remarkable computational advantage, completing estimation in approximately 7 seconds compared to more than 70 seconds required by NFXP. This performance gap stems primarily from ECP’s convex reformulation, which eliminates the repeated nested evaluations of the Bellman equations inherent in NFXP. Overall, these real-world experiments confirm that ECP not only matches the statistical accuracy of the classical NFXP framework but also provides substantial gains in computational efficiency, particularly in large-scale settings with many product alternatives.
The RL model was originally proposed for route choice problems in transportation networks Fosgerau2013. It aims to model how travelers select paths between origin-destination pairs based on factors such as travel time or distance. Traditional multinomial logit models treat entire routes as discrete alternatives; however, this formulation quickly becomes intractable as the number of feasible paths grows combinatorially with network size. The RL model overcomes this limitation by framing route choice as a sequence of link-level decisions, enabling the computation of choice probabilities without explicitly enumerating all possible paths. Despite this theoretical advantage, the estimation of RL models remains computationally demanding. The NFXP algorithm is the conventional approach for parameter estimation in the route choice context, but as discussed earlier, it frequently fails to recover valid value functions $\mathbf{V}$, particularly when the magnitudes of the parameters $\boldsymbol{\beta}$ are small or when the network exhibits ill-conditioned structures MaiFrejinger2022.
In this experiment, synthetic traffic networks are constructed using random geometric graph methods - a topology commonly used in graph algorithm benchmarks penrose2003random. Two types of transport networks are considered: DAGs and undirected graphs. DAG-based networks have predetermined origin and destination nodes. Undirected graph-based networks are undirected versions of DAG instances (by converting all arcs into undirected edges). Four variants of transport networks are spawned (20, 30, 40 and 50 nodes), with 5 instances for each variant. In the RL model, each link represents a transition between a pair of connected arcs and is characterized by four attributes. The first, $TT(a)$, denotes the travel time on link or route $a$, which is proportional to the length of arc $a$. The second attribute, $LT(a \mid k)$, is a left-turn indicator that equals 1 if the turning angle from node $k$ to arc $a$ lies between $30^{\circ}$ and $150^{\circ}$. Similarly, $RT(a \mid k)$ is a right-turn indicator that equals 1 when the turning angle lies between $-150^{\circ}$ and $-30^{\circ}$. Finally, $UT(a \mid k)$ is a U-turn indicator that equals 1 if the absolute turning angle exceeds $150^{\circ}$ in either direction. These attributes collectively capture key geometric and behavioral features of route choice, allowing the RL model to reflect travelers’ preferences for smoother trajectories and to penalize sharp or reversing turns.
The following specification for the link utility function is employed:
To generate synthetic observations, we specify a ground-truth parameter vector $\widehat{\boldsymbol{\beta}}$ in which all coefficients are negative, reflecting travelers’ general preference to avoid delays and difficult maneuvers. Among these attributes, travel time exerts the strongest negative influence, while turning penalties decrease in magnitude from U-turns to left turns to right turns. Accordingly, the following parameter values are adopted across all experiments: $\hat{\beta}_{TT} = -4,~ \hat{\beta}_{LT} = -0.1,~ \hat{\beta}_{RT} = -0.05,~ \hat{\beta}_{UT} = -0.3.$ A total of 3,000 observations are generated for parameter estimation and in-sample prediction. Further details on the data generation process for each experimental scenario are provided in the corresponding subsections.
\paragraph{Estimation time comparison.} Table (ref) presents the estimation results obtained on DAG-based synthetic networks. Across all test instances, the NFXP and ECP approaches produce identical log-likelihood values, confirming their theoretical equivalence and the numerical consistency of the proposed equilibrium-constrained formulation. This alignment demonstrates that the ECP method recovers the same optimal solution as the conventional NFXP algorithm. However, the runtime comparison clearly highlights the computational advantage of ECP. For every tested network size, ECP achieves convergence in a fraction of the time required by NFXP—often more than an order of magnitude faster. This improvement stems from the convex reformulation of the estimation problem, which eliminates the need for nested fixed-point iterations inherent to NFXP. Consequently, ECP not only preserves the statistical accuracy of the original recursive logit estimation but also significantly enhances its scalability and numerical stability, particularly for larger and more complex network structures.
We now compare the performance of the NFXP and the proposed ECP approaches under undirected cyclic network settings, which are known to be substantially more challenging for RL model estimation MaiFrejinger2022. In such networks, both methods occasionally fail to produce valid value functions and link choice probabilities due to the presence of cycles, which complicates the convergence of the Bellman operator. To overcome this issue when generating synthetic observations, we temporarily convert each undirected graph into a DAGs. Specifically, the conversion process limits the maximum path length to the number of nodes in the original network. For each undirected graph, a corresponding multi-layer DAG is constructed following these rules: each layer contains one duplicate of all nodes (without intra-layer edges), and a node at layer $k$ can only connect to nodes at layer $k+1$ if the two nodes are connected in the original undirected graph. The origin node is positioned in the first layer and the destination node in the final layer.
We evaluate the two methods using two performance metrics: (i) the successful rate, defined as the percentage of runs (out of 40) where the RL model is successfully estimated without numerical failure, and (ii) the estimation time, computed as the average runtime over successful runs. A method is deemed unsuccessful if it terminates prematurely due to numerical instability (e.g., the solver is forced to stop) or produces invalid log-likelihood values such as NaN or positive values. The comparative results are summarized in Table (ref).
As shown in the table, the proposed ECP method demonstrates a clear advantage over NFXP in terms of both stability and computational efficiency. Across all network sizes, the success rate of ECP remains close to 100%, whereas NFXP exhibits highly variable and often poor stability, with success rates dropping as low as 2.5% in some instances. This instability arises because NFXP relies on iterative fixed-point computations of the value function, which are highly sensitive to initialization and scaling in cyclic network structures. In contrast, ECP formulates estimation as a single convex optimization problem, avoiding the nested iterations that typically lead to divergence in NFXP. Regarding computational performance, ECP also achieves competitive or faster runtimes in most cases despite solving a larger conic program. For smaller networks, both methods exhibit comparable speeds; however, as network size grows, NFXP’s runtime increases more erratically due to repeated failures and restarts, while ECP maintains consistent convergence behavior. For example, with 50-node networks, ECP achieves nearly full stability while completing estimation in practical time, whereas NFXP succeeds in less than half the runs and displays large runtime fluctuations. Overall, these results highlight that ECP provides a significantly more robust and reliable estimation framework for RL models in complex cyclic network environments, preserving accuracy while offering superior numerical stability and scalability.
\paragraph{Scalability with respect to the number of parameters.} Similar to the MDC experiments discussed earlier, we evaluate the scalability of both the NFXP and ECP approaches with respect to the dimensionality of the parameter vector $\pmb{\beta}$. The experiments are conducted on both DAGs and undirected cyclic networks. For each combination of graph size and attribute dimensionality, five instances are generated, and each instance produces 3,000 observations using the same data generation procedures described previously for the DAG and undirected cases. To ensure statistical reliability, each configuration is repeated 40 times, and the reported estimation times correspond to the mean runtime over all successful runs.
For DAG instances, the number of nodes ranges from 20 to 50, and the number of attributes per link varies from 10 to 100. The results, shown in Figure (ref), exhibit a consistent pattern with those observed in the LMDC framework: ECP’s estimation time remains largely unaffected by the increase in parameter dimensionality, while NFXP shows an approximately linear growth in runtime.
For undirected graphs, where the estimation process is considerably more computationally demanding, we restrict the analysis to 20- and 30-node networks. The corresponding results are presented in Figure (ref). Only runs where both NFXP and ECP successfully converged are included, since NFXP frequently encountered numerical instabilities, as demonstrated in earlier experiments. Unlike the DAG case, the results on undirected networks reveal a distinct behavioral pattern. While NFXP continues to scale linearly with the number of parameters, ECP exhibits irregular and less predictable timing behavior, failing to scale smoothly. This degradation in performance is primarily attributed to the underlying cyclic topology: the equilibrium constraints in ECP are more sensitive to feedback loops, which exacerbate the computational burden in solving the conic program. Consequently, ECP’s runtime grows sharply as network complexity increases, even though it still achieves higher stability and convergence reliability compared to NFXP.
In summary, ECP maintains near-constant scalability with respect to parameter dimensionality in DAG settings but suffers in cyclic graphs due to the structural complexity of undirected networks. Conversely, NFXP demonstrates better time scalability in cyclic settings but remains less stable overall. These findings collectively highlight the trade-off between stability and computational efficiency across different network topologies.
We evaluate the performance of the proposed method using a real-world dataset collected from the Borlänge transportation network in Sweden. This dataset has been widely employed in previous studies on route choice modeling FosgFrejKarl13,MaiFosFre15,mai2018decomposition. The underlying network consists of $3{,}077$ nodes and $7{,}459$ links. Since traffic conditions in this network are generally uncongested, link travel times can be treated as static and deterministic. The sample includes $1{,}832$ observed trips, each represented by a simple path containing at least five links. In total, the dataset covers $466$ destinations, $1{,}420$ distinct origin--destination pairs, and more than $37{,}000$ link-level choices. The instantaneous link utility function is specified as:
where $TT(a)$ denotes the travel time on link $a$, $LT(a \mid k)$ is a left-turn indicator that equals one if the transition from link $k$ to link $a$ involves a left turn, $LC$ is a constant attribute equal to one for all links, and $UT(a \mid k)$ is a U-turn indicator for transitions that reverse direction. The topology of the Borlänge transportation network derived from this dataset is illustrated in Figure (ref).
For the real-world Borlänge network, the large number of links results in an extremely high-dimensional ECP formulation. Consequently, solving the full-scale conic program directly can be computationally demanding, requiring substantial memory resources and often leading to numerical instability. To mitigate these issues, we apply the network trimming procedure introduced in Section (ref) to reduce the network size before full estimation. The main objective of this preprocessing step is to obtain a smaller yet representative subnetwork in which the ECP can be efficiently solved. By estimating the RL model on the trimmed network, we aim to recover a set of parameter estimates that closely approximate those of the full-scale model. These estimates serve as high-quality initial values for the NFXP algorithm, which can then be applied to the complete network. Starting from parameters near the true optimum significantly enhances the numerical stability of NFXP and reduces the likelihood of divergence in the value function iteration. In essence, this two-stage estimation strategy leverages the robustness and convexity properties of ECP to provide a reliable initialization for NFXP, thereby combining the strengths of both methods: ECP ensures stable convergence in reduced networks, while NFXP refines the estimates to achieve final accuracy on the full network.
We apply the network trimming procedure described in Section (ref), where the threshold parameter $\epsilon$ is chosen such that approximately $90\%$ of the links are removed from the original network. This produces one subnetwork for each origin--destination pair, with a size roughly one-tenth that of the full network. We then estimate the RL model on each subnetwork using both the NFXP and ECP methods. To ensure a fair comparison, we deliberately exclude stabilization heuristics commonly employed in prior RL studies (e.g., manually replacing invalid values with finite approximations during value function computation). The estimation results obtained from both methods on the trimmed networks are summarized in Table (ref). In terms of computation time, NFXP demonstrates a clear advantage, requiring only $46$ seconds compared to $363$ seconds for the convex ECP approach. However, NFXP suffers from numerical instability, as indicated by the positive log-likelihood values, while the ECP estimates remain numerically stable and yield parameter values consistent with those reported in the original study by Fosgerau2013. These findings are consistent with our earlier observations from synthetic experiments on undirected graphs, further confirming the robustness of ECP under challenging numerical conditions.
Table (ref) summarizes the estimation results of the RL model on the full Borlänge network using three different estimation strategies, namely: (i) direct estimation using NFXP on the full network, (ii) NFXP applied to the full network with initialization from parameters estimated on trimmed subnetworks using NFXP, and (iii) NFXP applied to the full network with initialization from parameters estimated on trimmed subnetworks using the ECP approach. For the first approach, the initial parameter vector is set to $-1.5$ for all coefficients $\beta$, following the setup used in FosgFrejKarl13. As shown in the table, direct estimation via NFXP on the full network fails to converge, as indicated by the missing parameter values. This outcome is expected given the extremely large number of links and the inherent instability of the nested fixed-point approach when applied to high-dimensional, cyclic networks. The fixed-point iterations required by NFXP often diverge or return infeasible value functions, leading to numerical failure.
In contrast, the two-stage approaches—where model estimation is first performed on trimmed subnetworks—yield stable results. When NFXP is initialized with parameters estimated from trimmed networks (second row), the algorithm successfully converges on the full network, producing reasonable parameter estimates and a normalized log-likelihood of $-4.37$. The signs and magnitudes of the estimated coefficients are consistent with behavioral expectations: all parameters are negative, reflecting aversion to travel time, turning movements, and U-turns, with travel time ($\hat{\beta}_{TT}$) and U-turns ($\hat{\beta}_{UT}$) having the strongest influence.
The hybrid approach, where ECP is first applied on trimmed subnetworks before refinement with NFXP on the full network, achieves the best balance between numerical stability and model fit. The resulting coefficients remain behaviorally plausible and closer in magnitude to those reported in previous literature Fosgerau2013, while the corresponding log-likelihood ($-3.44$) indicates a superior fit compared to the purely NFXP-based initialization. This improvement demonstrates that ECP provides more reliable and well-conditioned initial estimates, enabling NFXP to converge more efficiently to a high-quality local optimum.
Overall, these results reinforce the advantages of using ECP as a preprocessing step for large-scale RL model estimation. Direct NFXP estimation remains impractical for complex real-world networks, but the combination of ECP-based initialization and subsequent NFXP refinement offers a computationally feasible and statistically robust alternative for achieving accurate model parameters.
Here we note that the estimation of the NRL model is not included in the present experiments. This omission is primarily due to methodological and computational considerations. First, the proposed ECP framework is not directly applicable to NRL under the general setting, as the introduction of nested structures breaks the global convexity of the underlying value function formulation. In particular, the equilibrium constraints that make ECP tractable for standard RL rely on a single-level, additive utility structure, which no longer holds in the presence of hierarchical or correlated error components inherent to NRL. Second, the classical NFXP estimation approach for NRL suffers from the same numerical and convergence issues as those encountered in the standard RL model, such as instability in the value function computation and divergence during the inner fixed-point iteration. In practice, researchers often use the simpler RL model—estimated via NFXP—to obtain a reasonable set of initial parameter values before proceeding with NRL estimation on the same network MaiFrejinger2022,MaiFosFre15. Therefore, including NRL results here would not yield additional insights beyond those already demonstrated for RL, as both methods share similar computational challenges while the ECP framework, by design, cannot be extended to NRL in a straightforward manner.
This paper proposed an equilibrium-constrained optimization (ECP) framework for estimating RL models, offering a convex and numerically stable alternative to the classical NFXP approach. Theoretical analysis established the equivalence between ECP and NFXP, while experiments on synthetic and real transportation networks demonstrated that ECP achieves similar log-likelihood performance with improved numerical stability and scalability. Empirical results revealed that ECP performs particularly well on directed acyclic networks, where it maintains stable convergence and near-constant scalability as parameter dimensionality grows. In contrast, NFXP remains faster in small-scale problems but exhibits instability and divergence in complex or cyclic networks. When used together, ECP provides reliable initialization that significantly enhances NFXP convergence on large-scale networks.
Future research will focus on extending the ECP formulation to handle the NRL model and other hierarchical structures Mai_RNMEV, as well as developing decomposition-based or distributed algorithms to further enhance scalability. In addition, applying the ECP framework beyond transportation—such as to facility location and assortment optimization problems mai2020multicut,Talluri2004revenue—represents a promising direction for expanding its applicability.