A Benders decomposition for maximum simulated likelihood
estimation of advanced discrete choice models
T. Haeringa,∗ , C. Bongiovannia and M. Bierlairea
a École Polytechnique Fédérale de Lausanne, Transport and mobility laboratory,

Lausanne, Switzerland
tom.haering@epfl.ch, claudia.bongiovanni@epfl.ch, michel.bierlaire@epfl.ch∗ Corresponding author
Extended abstract submitted for presentation at the 11th Triennial Symposium on
Transportation Analysis conference (TRISTAN XI)
June 19-25, 2022, Mauritius IslandApril 11, 2022

Keywords: maximum likelihood estimation, discrete choice, simulation, mixed integer linear programming, benders decomposition

1INTRODUCTION

Maximum likelihood estimation (MLE) is a broadly used method to estimate the parametersof a previously specified distribution, given observed data. It finds its use in many areas ofphysics (e.g. Hauschild & Jentschel (2001)), machine learning (e.g. Goodfellow et al. (2016))
and discrete choice modelling (e.g. Bierlaire (2003)). Advanced discrete choice specifications likelatent class or probit models are challenging to estimate, as the choice probabilities resulting fromsuch models are highly non-convex and do typically not have a closed-form expression. For thisreason, optimization approaches have relied on simulation techniques, i.e. maximum simulated
likelihood estimation (MSLE), see Train (2009). A general approach for MSLE has been proposed
in Fernández Antolín (2018), where the problem is formulated as a mixed integer linear program(MILP). This allows to define the model in terms of its error components, enabling the approach
to be flexibly applied to any advanced discrete choice model. It is furthermore independent fromthe complexity of the error term distributions, the only requirement being that its possible totake draws. If the number of draws is sufficiently large, the MILP formulation is guaranteed toconvergence to a globally optimal solution. However, since the complexity of the MILP scales
exponentially with the number of draws, the approach can currently only be applied to solving
small-scale instances (i.e., with few individuals and alternatives).
In this work, we extend the range of applicability of the MILP approach in Fernández Antolín(2018) by means of a Benders decomposition, which speeds-up the MILP solution process for
the MLSE drastically and enables to scale-up the tackled instances. Our designed Benders decomposition exploits total unimodularity to keep the master problem linear, thus eliminating thebottleneck in computational time usually associated with Benders decomposition. The proposed
approach is benchmarked against the full MILP and PandasBiogeme. Preliminary results on asimple logit model show that the Benders decomposition approach solves instances up to 60xtimes faster than the MILP, while retaining high quality solutions.
TRISTAN XI Symposium

Original abstract submittal

T. Haering, C. Bongiovanni and M. Bierlaire

22METHODOLOGY

In this section, we formally introduce an MILP formulation for the MSLE problem, based on thework in Fernández Antolín (2018), and a problem-specific Benders decomposition approach.

2.1

MILP formulation

Consider a set of n = {1, . . . , N } individuals
choosing exactly one alternative among a set
β,ω,s,z,U,Hniof i = {1, . . . , I} alternatives. Such choice iss.t.
depicted by a binary decision variable yin . AsX
ωinr = 1(µnr ) sume that each individual n selects the alteri
native i corresponding to the maximal utilX
Hnr =Uinr ωinr
(ζnr ) ity Uin , i.e. yin = 1 ⇔ Uin = maxj Ujn .
iThe utility function is defined in constraints
Hnr ≥ Uinr
(αinr ) (κinr ) and depends on k parameters β whichX
are to be estimated. The simulated choices for
sin =
ωinr
(θin )
each scenario r = {1, . . . , R} are captured byr
the binary variables ωinr , together with conzin ≤ Lr − Kr sin(ξinr )straints (µnr ), limiting the choice to a sinX
Uinr =βk xink + ϵinr (κinr ) gle alternative. Constraints (ζnr ) and (αinr )kguarantee that the chosen alternative correωinr ∈ {0, 1}
sponds to the one with the highest utility. Theobjective is to maximize
the simulated logβ, s, z, U, H ∈ R
Q Q
likelihood, given by ln( n i P̂n (i)yin ), whereP̂n (i) represents the estimator for the probFormulation 1 – MSLE as an MILP
ability of individualn choosing alternative i
1 P
and is given by R r ωinr . Constraints (θin )and (ξinr ) model a piece-wise linearization of the log-transformation, utilizing constants Lr =(1 + r) ln(r) − r ln(1 + r) and Kr = ln(r) − ln(1 + r) for the intercepts and slopes.
max

2.2

XX

yin zin

Benders decomposition approach

Combinatorial optimization problems that are characterized by complicating integer decisionvariables are typically tackled by a Benders decomposition approach (Benders (1962)). As theintegral master problem needs to be solved repeatedly, Benders is notorious for its slow convergence. In our case we can use an elegant trick to avoid this issue: by identifying the continuous
estimation parameters β as the complicating variables and fixing them in the subproblem, theutilities of all the alternatives become fixed as well. Thus the problem of choosing the highest
utility alternative simplifies to a knapsack problem, which is totally unimodular, which allows todrop the integrality constraints on the choice variables. Formulations 2 and 3 give the respective
definitions of the primal and dual of the subproblem, while Formulation 4 describes the masterproblem.
As the linearization of constraint (ζnr ) using a big-M approach no longer works when integrality constraints are relaxed, the formulation in the primal subproblem needs to be modified.
The product ηinrk = ωinr βk is modeled directly using constraints (πinr ), (λinrk ) and (φβnrk ). This
formulation is equivalent to Formulation 1. It is important to mention that, in order to preserve
the inegrality of the primal, information about β fixed had to be kept in the matrix, which implies
it also being contained in the matrix of the dual (constraints (χinr )). This means the feasible
region of the dual subproblem is not constant over iterations, which might distort the Bendercuts. Lastly, both the primal and the dual models are fully decomposable on the individuals n,
as individuals select alternatives independently from each other and across scenarios.
TRISTAN XI Symposium

Original abstract submittal

T. Haering, C. Bongiovanni and M. Bierlaire

2.2.1

3Subproblem

minβ,ω,χ,η,s,z,H−XXyin zinniXωinrXmaxµ,α,ζ,µ,λ,φβ ,θ,ξs.t.
=1(µnr )iXβk xink − Hnr≤−ϵinrX≤X(αinr )kHnr −

ηinrk xinkX

ωinr ϵinr (ζnr )ϵinr αinr +X

πinrnrinr
inr
XX
fixed+βk λinrk +Lr ξinr
inr
inrk

s.t.
µnr − ζnr ϵinr + πinr − θinX
πinr +
βkfixed λinrk≤0(ωinr )≤0(χinr )kiikµnr −−X

αinr + ζnr=0(Hnr )χinr + ωinr

=1(πinr )iηinrk + βkfixed χinr
X
βk −
ηinrk

=βkfixed(λinrk )=0(ηinrk )=0(φβnrk )−ζnr xink + λinrk − φβnrkXθin +Kr ξinr≤0(sin )

ξinr=−yin(zin )φβnrk=0(βk )µ, π, λ, θ, φβ∈Rα, ζ, ξ∈R≤0risin −Xωinr=0(θin )zin + Kr sin≤

Lr(ξinr )ω, χ, s∈R≥0β, η, z, H

∈

RXrr

X
inr

αinr xink +X
nrFormulation 2 – Primal subproblem
Formulation 3 – Dual subproblem
2.2.2Master problemmin LL,βs.t.
L≥L∗ +L≥LbestL, β∈RXXnϕ∗nk (βk − βkfixed )kFormulation 4 – Master problem

3

The master problem reduces to finding optimal values for the estimation parameters β.
For each β fixed , after solving the dual subprob(1) lem, a Benders cut of the same type as constraint (1) is added. Each optimal objective(2)
value serves as a new lower bound on the objective, enforced in constraint (2). The parameters ϕ∗nk of the Benders
P cuts in our case aredetermined by ϕ∗nk = ir λinrk .

PRELIMINARY RESULTS

To test our approach, we conduct experiments on a binary logit model. A mode choice problem
between two alternatives, public transport (pt) and car, is considered. The systematic utilitiesof the alternatives are:Vcar = βtime · traveltimecar
Vpt = βtime · traveltimept
The dataset is extracted from revealed preference data on mode choice collected in 1987 for theNetherlands Railways, consisting of 228 respondents (CASE, 2017). Experiments are performed
using GUROBI 9.5.0 (Gurobi Optimization, LLC, 2021) on a 2.6 GHz 6-Core Intel Core i7processor with 16 GB of RAM, with a three hour time limit per instance. Our proposed Benders
approach is benchmarked against PandasBiogeme (Bierlaire, 2020) and the full MILP, comparing
the objective values, the parameter values, and the runtimes. Biogeme’s objective function
TRISTAN XI Symposium

Original abstract submittal

T. Haering, C. Bongiovanni and M. Bierlaire4is the log-likelihood (LL), which is approximated by the simulated log-likelihood (sLL), theMILP objective. For comparison, the LL is also evaluated using the estimated parameters fromdecomposition and full MILP. We take random subsets of individuals from the population to getinstances that are manageable for the MILP. The results shown in Table 1 highlight the following:
1. On average, the decomposition solves the problem over 30 to 60 times faster, 2. comparing theoptimal solution values for the full MILP and our decomposition reveals differences, and 3. thesedifferences are small for the objective function and more significant for the estimated parameters.
Although a Benders decomposition is an exact approach, our formulation contains mathematicalN202020505050100100100
200200200R50100200
50100200
50100200
50100200

LL-Bio-12.3029-12.3029-12.3029-30.2645-30.2645-30.2645-64.8827-64.8827-64.8827-122.6885-122.6885-122.6885

LL-D-12.4934-12.4109-12.4633-30.683-30.481-30.2825-65.3962-66.0306-64.9245-122.6895-122.7386-122.7213

LL-M-12.4443-12.3948-12.3778-30.3258-30.3258-30.3254-64.8978-64.8828-64.8933-122.7352-122.92-123.3417

sLL-D-12.658-12.2581-12.6478-31.0302-31.0396-30.6917-65.8014-67.419-66.0181-124.0274-124.2428-124.1058

sLL-M-12.6074-12.2116-12.2834-30.8476-30.4611-30.5655-65.2044-65.7837-65.6991-123.5507-124.0001-124.7073

B-Bio-1.5579-1.5579-1.5579-1.4095-1.4095-1.4095-0.9481-0.9481-0.9481-1.31-1.31
-1.31

B-D
-0.9704-1.1097-2.1604-0.9354-1.7825-1.3072-0.6117-0.4513-0.8502-1.322-1.3929-1.377B-M-1.0481-1.1432-1.1823-1.2226-1.2226-1.2232-0.889-0.9431-0.8987-1.3901-1.4898-1.0213T-D
10.06139.902316.9385
29.7795
65.0063
98.2058
28.7805
274.1626
295.7408
120.5793
327.2528
1262.7548

T-M
64.9419403.6936
1117.0644
286.6792
1558.6037
5375.6553
2820.2287
4346.0671
10800+
1476.1851
10800+10800+

Table 1 – Comparing our decomposition method with the full MILP and PandasBiogeme(N = population size, R = number of draws, LL = log-likelihood , sLL = simulated log-likelihood,B = βtime , T = time (sec.), Bio = Biogeme, D = decomposition, M = MILP)
aspects that may currently prevent the convergence to the real global optimum. As mentioned in
the methodology, a possible explanation for the deviations is the fact that information about themaster variables is maintained in the matrix of the dual. Other explanations include numericalissues, stemming for example from the linearization of the logarithm or the way certain solvers
handle specific constraints.

4FUTURE WORKThis is ongoing work. We are currently investigating the cause of the deviations between theoptimal values from the decomposition and the full MILP, as well as testing extensions andperformance on various advanced discrete choice models. Additional results obtained by thetime of the conference will of course be included in the presentation.

References
Benders, Jacques F. 1962. Partitioning procedures for solving mixed-variables programming problems.
Numerische mathematik, 4(1), 238–252.
Bierlaire, Michel. 2003. BIOGEME: A free package for the estimation of discrete choice models. In: Swisstransport research conference.
Bierlaire, Michel. 2020. A short introduction to PandasBiogeme. A short introduction to PandasBiogeme.
CASE, NETHERLANDS MODE CHOICE. 2017. Data collection.
Fernández Antolín, Anna. 2018. Dealing with Correlations in Discrete Choice Models. Tech. rept. EPFL.
Goodfellow, Ian, Bengio, Yoshua, & Courville, Aaron. 2016. Machine learning basics. Deep learning,1(7), 98–164.
Gurobi Optimization, LLC. 2021. Gurobi Optimizer Reference Manual.
Hauschild, T, & Jentschel, M. 2001. Comparison of maximum likelihood estimation and chi-squarestatistics applied to counting experiments. Nuclear Instruments and Methods in Physics Research
Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 457(1-2), 384–401.
Train, Kenneth E. 2009. Discrete choice methods with simulation. Cambridge university press.

TRISTAN XI Symposium

Original abstract submittal