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.
129,465 characters · 34 sections · 39 citation commands
Torch-Choice: A PyTorch Package for Large-Scale Choice Modeling with Python
Choice models are ubiquitous in economics, computer science, marketing, psychology, education and more fields. In consumer demand modeling, the researcher considers a setting where a consumer selects a product or brand from a set of alternatives, and the goal is to analyze how product characteristics affect consumer purchase decisions. These estimates can in turn be used to answer counterfactual questions about how purchases would change if, for example, prices or availability change. Choice models can be applied to aggregate market level demand data {berry1995}, cross-sectional individual choice data {berry2004}, and panel data where customer purchase events are tracked over time {chintagunta2001}. In transportation choice modeling, choice models are used to analyze intercity passenger travel mode choice {mcfadden1974measurement, wilson1990}, and in spatial modeling to model location choice {miyamoto2004}. In labor economics, choice models are used to model labor supply and childcare decisions {kornstad2007}, and family labor supply {van1995}. In education, Item Response Theory (IRT) and related models are choice models {baker2001basics, embretson2013}. In psychology, choice models are used to model decision-making under uncertainty {bonoma1979}, and psychiatric treatment diagnoses {fleiss1986}. In computer science, choice models are used to build recommendation systems {zhang2002} as well as for classification tasks such as spam filtering {chang2008}, image and character recognition {deng2012mnist}, and many other applications.
Multinomial logistic regression and extensions are the most commonly used choice models. More recently, there has been work using choice models and their extensions (matrix-factorization, multilayer neural networks for choice modeling, Bayesian embeddings based choice models) to model restaurant choice {athey2018}, grocery shopping {donnelly2021}, inter-dependencies among choices when putting together a shopping basket {ruiz2020}, the welfare effects of recommendations {donnelly2022}, and to model labor market transitions {vafa2022}. For each of the models, either models based on multinomial logistic regression are used, or they are important baselines for further model development, creating a strong need for effective choice modeling software.
Recent years have seen the proliferation of various tools and software packages that go a long way in making it easy to estimate the parameters of choice models. While each of these have their particular strengths, there are important limitations of these packages when they are used with the goal of building flexible, fast choice models which scale to massive data.
In the discussion that follows, when we refer to \code{R} packages, we mean \code{glm} and \code{mlogit}. When we refer to \code{Stata} packages, we mean the entire suite of choice models offered in the packaged software.\footnote{See \url{https://www.stata.com/manuals/cm.pdf}, \url{https://www.stata.com/help.cgi?xtmelogit} and \url{https://www.stata.com/manuals/cmnlogit.pdf}} When we refer to flexible functional forms, we mean model specifications with user- and item-specific latent coefficients (latent variables). When we refer to regularization, we mean the ability to incorporate penalization or “shrinkage” of regression coefficients {krogh1991simple}.
We propose \code{torch-choice}, a library for flexibly, fast choice modeling with PyTorch; it implements conditional and nested logit models, designed for both estimation and prediction. This package aims to address these limitations and fill the gap between tools used by econometricians/statisticians (R-style, end-to-end model training) and those computer scientists who are familiar with (PyTorch).
A limitation of \code{torch-choice} is that while it natively supports estimating user-specific coefficients in the setting of panel data, \code{torch-choice} estimates a constant coefficient (instead of distribution) for each user, and the same coefficient is applied across all sessions involving this user. It does not, however, learn the distribution over user-level parameters, and does not model that distribution. As such, it does not support random coefficients in cross-sectional data settings, where one can learn the parameters of distribution over user-level parameters in \code{R}, \code{Stata}, or \code{xlogit} in Python {arteaga2022xlogit}. Please see section (ref) for more details.
What are the advantages of using \code{torch-choice}?
The \code{torch-choice} library is publicly available on Github at: \href{https://github.com/gsbDBI/torch-choice}{\code{torch-choice} github repository}. Researchers can install our packages via pip following instructions at \href{https://pypi.org/project/torch-choice}{\code{torch-choice} on pip}. We have prepared tutorials on the documentation website: \href{https://gsbdbi.github.io/torch-choice}{\code{torch-choice} documentation website}, where researchers can find the full documentation of APIs. Jupyter notebook {Kluyver2016jupyter} versions of these tutorials are available on Github at: \href{https://github.com/gsbDBI/torch-choice/tree/main/tutorials}{\code{torch-choice tutorial notebooks}}. Code demonstrated in this paper can be found in this \href{https://github.com/gsbDBI/torch-choice/blob/main/tutorials/paper_demo.ipynb}{Jupyter Notebook}.
Section 2 outlines the choice problem \code{torch-choice} was built to model. Section 3 covers the data structures. Section 4 covers the models implemented in \code{torch-choice}. Section 5 shows performance benchmarks for \code{torch-choice} against another open source package, \code{mlogit} in \code{R}, and demonstrates its performance advantages. Finally, Section 6 concludes this paper.
We start by outlining the building blocks of the choice problem that will be processed with \code{torch-choice}. Models in \code{torch-choice} predict which item among all available items a user will choose given the context (where the context is referred to as a session). The paper refers to the basic unit of analysis as a choice record.
Figure (ref) illustrates the layout of the \code{ChoiceDataset} data structure. Each row of the central dataframe in Figure (ref) corresponds to a choice record.
Choice records are the fundamental unit of observation in choice modeling describing who (i.e., which user) chooses what (i.e., which item) under what context (i.e., in which session). Let $N$ denote the size of the dataset and $n \in \{1, 2, \dots, N\}$ denote the index for records in the dataset. Let $\mathcal{D}$ denote the dataset and $\square^{(n)}$, the super-script with parentheses, denote the $n$-th record in the dataset, where $\square$ is a placeholder for the object of interest. Therefore, $\mathcal{D}^{(n)} = (u^{(n)}, i^{(n)}, s^{(n)})$ denotes the item chosen, the user involved, and the corresponding session index of the $n$-th record. Similarly, to represent subsets of a dataset, we can use superscripts. For instance, $\mathcal{D}^{(\text{train})}$ denotes the training set and $\mathcal{D}^{(\text{test})}$ denotes the test set.
The set of $U$ users is indexed by $u \in \{1, 2, \dots, U\}$, while the set of $I$ items is indexed by $i \in \{1, 2, \dots, I\}$. The set of sessions $S$ is indexed by $s \in \{1, 2, \dots, S\}$. The set of items is partitioned into $C$ categories indexed by $c \in \{1,2,\dots,C\}$.\footnote{The data structure in \code{torch-choice} does not limit how many categories a user could choose each time. The notion of categories is fully general and can be any partition of the item set. However, the nested-logit model implemented in \code{torch-choice} allows for modeling only one category, with two levels of nesting. See the model section for more details. However, modeling multiple categories in a single model is not currently supported for nested-logit models implemented in \code{torch-choice}. This is left to future work.} Let $I_c$ denote the collection of items in category $c$, $I_c$'s are disjoint sets and $\bigcup_{c=1}^C I_c = \{1, 2, \dots I\}$. Note that when we learn over multiple categories simultaneously, parameters that are not item-specific are the same across categories {rudolph2017structured}. To learn per category parameters at the user and constant level, one needs to specify separate datasets and learn separate models over them. Let $\mathcal{C}(i) \in \{1, 2, \dots, C\}$ denote the category of item $i$, $I_{\mathcal{C}(i)}$ represents items in the same category as item $i$.\footnote{Since we will be using PyTorch to train our model, we represent user, item, session identities with consecutive integer values instead of the raw human-readable names of items (e.g., Dell 24-inch LCD monitor). Similarly, we encode user indices and session indices as well. Raw item names can be encoded to integers easily with \href{https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.LabelEncoder.html}{\code{sklearn.preprocessing.LabelEncoder}} or \href{https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.OrdinalEncoder.html}{\code{sklearn.preprocessing.OrdinalEncoder}}.}
The model assumes each user must choose exactly one item from each category (where one item may be defined to represent not selecting any items) and that choices are independent across categories and conditional on user preferences. The user $u$ chooses an item from each category to maximize their utility $\mathcal{V}_{uis}$ from item $i$ in the context of session $s$. We assume the utility can be decomposed into a deterministic utility $\mu_{uis}$ and an error term $\varepsilon_{uis}$ following the Gumbel distribution as in Equation ((ref)).\footnote{The Gumbel distribution is a continuous distribution with a single parameter $\beta$ and the probability density function is given by $f(x) = \frac{1}{\beta} \exp(-x/\beta - \exp(-x/\beta))$. It is easy to show that with Gumbel noise, the probability density of $\max_{i}\mathcal{V}_{uis}$ follows the softmax function of $\mu_{uis}$.}
The same user will make choices from different categories in the same session, so that the dataset will contain multiple purchasing records with the same indices for user and session $(u, s)$ but different item indices $i$.
The utility $\mu_{uis}$ depends on the characteristics of users and items and other factors (e.g., the utility from iced coffee drops during winter time). We introduce the notion of session indexed by $s \in \{1,2,\dots, S\}$ to capture factors that vary by time and context.
The flexibility of sessions enables easier data management (e.g., observables and item availability) as we will see later.
The session can be used to capture the fact that item availability may change over time. Our dataset admits a $S \times I$ binary matrix $A$ to indicate which items are available to be selected in each session.\footnote{The item availability acts as a constraint in users' choice problems. Computationally, we assign $-\infty$ utility to unavailable items to operationalize their exclusion in the choice expression.} We use $A(s)$ to denote the set of available items in session $s$. All items are available all the time by default, i.e., $A(s) = \{1, 2, \dots, I\}$ for all session $s$. Researchers should define the granularity of this matrix guided by the granularity with which item availability changes, along with the granularity of how users make decisions. For example, if the availability of items is the same in a day and a user faces at most one choice decision in a day, sessions may be defined at a day level; however, if users make more than one choice per day or if the availability of items is personalized to users, sessions will need to be more granular. Less granular session definitions encourage more efficient data management computationally. Note that in some applications, the prominence with which an item is displayed can also change over time; that can be captured as a time-varying item characteristic (i.e., an item-session observable).
The objective of our models is to predict utilities $\mu_{uis}$ as a function of user, item, and session characteristics called observable tensors.\footnote{Observables are also referred to as features in machine learning literature, independent variables in statistics, they are often denoted as $X$ with appropriate super/sub-script.}
The package supports all six types of observables with various dependencies on user $u$, item $i$, and session $s$. Observables are essentially high-dimensional arrays called tensors in machine learning literature. For example, an array $X$ of shape $U \times I \times 3$ can encapsulate three user-item-specific variables, indexing \code{X[u,i,:]} retrieves a vector in $\mathbb{R}^3$ for these observables corresponding to user \code{u} and item \code{i}. The following is a list of observable types supported and expected shapes of tensors/arrays encapsulating them.
Note on Computational Efficiency The user-session-item observable tensor can be extremely large and cause out-of-memory issues even at moderate levels of $U$, $I$, and $S$. For example, suppose we are modeling shopping data with customers as users and dates as sessions; there is a user-session-item-specific observable $x_{uis}$. Storing the complete $x_{uis}$ requires $U \times S \times I$ entries. Often, each customer only shops on a couple of dates, say three dates on average and $S \gg 3$. We can instead define sessions to be all existing pairs of (customer, date), and $X$ becomes a session-item-specific observable instead. Storing $X$ as a session-item-observable only requires $3 \times S \times I$ entries, a factor of $U/3$ reduction.
Since \code{torch-choice} leverages PyTorch and GPU accelerations, all observables need to be numerical, and indices of users/items/sessions should be consecutive integers starting from 0. Raw item names can be encoded easily with \href{https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.LabelEncoder.html}{\code{sklearn.preprocessing.LabelEncoder}} or \href{https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.OrdinalEncoder.html}{\code{OrdinalEncoder}}.
We assume each user selects at most one item from each item category each session by maximizing his/her utility $\mathcal{V}_{uis}$ within each category. When users are choosing over items, we can write utility $\mathcal{V}_{uis}$ that user $u$ derives from item $i$ in session $s$ as:
where $\varepsilon$ is an unobserved random error term. If we assume i.i.d. extreme value type I errors for $\varepsilon_{uis}$, which leads to the logistic probabilities of user $u$ choosing item $i$ among $I_{C(i)}$ (i.e., items in the same category as $i$) in session $s$, as shown by McFadden{mcfadden1973conditional}, and as often studied in Econometrics.
Given the distributional assumption on $\varepsilon_{uis}$, the probability that user $u$ chooses item $i$ in session $s$ is captured by the logistic function
where $A(s)$ is the set of items available in session $s$ (i.e., the choice set), and $I_{C(i)}$ is the set of items in the same category as item $i$. The probability in Equation (ref) is normalized over available items only; therefore, probabilities of choosing unavailable item are zero. The mean utility $\mu_{uis}$ for a user, item, and session is specified according to a functional form, and we consider different combinations of observable and latent variables, each of which corresponds to a distinct version of the model.
Equation (ref) describes the general form of the utility function in conditional logit models. Greek letters with tilde and prime (e.g., $\alpha, \tilde{\alpha}_i, \alpha_u'$) are coefficients to be estimated; capital letters denote observables; subscripts of them denote their dependencies on user $u$, item $i$, and session $s$.
For example, the specification in Equation (ref) allows for user observables ($W_u$), item observables ($X_i$), and item-session observable ($P_{is}$), and intercepts that vary with the item:
The Greek letters $\alpha_i$, $\beta_u$, $\gamma$, and $\delta_i$ denote coefficients to be estimated. These coefficients can be interpreted as the individual user fixed effect, impacts of item characteristics, user characteristics, and (session, item)-specific characteristics on the mean utility for user $u$ for item $i$. Researchers can add user-item observable interactions by constructing the corresponding item-session specific observable. The \code{torch-choice} package offers flexible ways to estimate the mean utility $\mu_{uis}$ beyond this example; the model section will cover these aspects in detail.
Beyond the basic and widely implemented multinomial logit model, the \code{torch-choice} library provides two widely used classes of choice models, the conditional logit model (CLM) and the nested logit model (NLM) to model $\mu_{uis}$. The standard multinomial logit model is nested in the conditional logit model, and so is treated as a special case. The nested logit model relaxes the i.i.d. assumption on $\varepsilon_{uis}$.
In each of the conditional logit and nested logit models, researchers can specify the functional form they want to use for observables. Each model takes user and session information $(u^{(n)}, s^{(n)})$ together with observables and outputs the predicted probability of purchasing each $\ell \in \{1, 2, \dots, I\}$, denoted as $\hat{P}(\ell \mid u^{(n)}, s^{(n)})$. In cases when not all items are available, the model sets the probability of unavailable items to zero and normalizes probabilities among available items. By default, model coefficients are estimated using the gradient descent algorithm to maximize log-likelihood in Equation ((ref)). Since \code{torch-choice} is built upon PyTorch, any other optimization algorithms supported by PyTorch can be used as well. Please refer to our online documentation for details regarding applying alternative optimization algorithms.
The models implemented in \code{torch-choice} offer flexible ways to estimate the deterministic part $\mu_{uis}$. Given the Gumbel distribution of $\varepsilon_{uis}$ and the independence of irrelevant items assumption, the predicted probability for choosing item $i$ among all items in the same category $I_{\mathcal{C}(i)}$ is a function of utilities $\{\mu_{ui's}\}_{i' \in I_{\mathcal{C}(i)}}$ from the model. The exact functional form differs by model, and we defer the discussion to the model section. The form of predicted probabilities for any given specificaton of $\mu$ is a normalized exponential function (“softmax”) function in Equation ((ref)).
This section covers the \code{ChoiceDataset} class in the \code{torch-choice} library.\footnote{Please see our online documentation \href{https://gsbdbi.github.io/torch-choice/data_management/}{Data Management Tutorial} and \href{https://gsbdbi.github.io/torch-choice/easy_data_management/}{Easy Data Wrapper Tutorial} for more details. We provide Jupyter notebooks offering hands on experience on our \href{https://github.com/gsbDBI/torch-choice/tree/main/tutorials}{Github repository}.}\footnote{Code examples in this section are for illustrating the usage only; please refer to online materials for executable notebooks.} \code{ChoiceDataset} is a specialized instantiation of PyTorch's \code{Dataset} class. It inherits support for all of the functionality of PyTorch's \code{Dataset} class, including subsetting, indexing, and random sampling. Most importantly, GPU memory is often much more scarce on modern machines, \code{ChoiceDataset} is manage information in a memory-efficient manner. Existing packages such as \code{xlogit} in Python arteaga2022xlogit and \code{mlogit} JSSv095i11 in R require the data to be in a long-format or wide-format, which store duplicate information. \code{ChoiceDataset} minimizes its memory footage via its specialized structure eliminating redundant information. Moreover, all information in \code{ChoiceDataset} is stored as PyTorch tensors,\footnote{Tensors are simply PyTorch's implementation of multi-dimensional arrays. Tensor operations are very similar to operations on \code{numpy} arrays; the \href{https://pytorch.org/tutorials/beginner/blitz/tensor_tutorial.html}{PyTorch tensor tutorial} overviews basics of tensors. This manuscript doesn't assume researchers' knowledge on PyTorch tensor operations. Two-dimensional tensors are simply matrices; we will use tensor and matrix interchangeably in this manuscript.} which allows researchers to leverage half-precision floats to further improve memory efficiency and seamless data transfer between CPU and GPU.
The \code{torch-choice} package provides two methods to build \code{ChoiceDataset}; the first method is to use the \code{EasyDataWrapper} class, which is a wrapper class that converts pandas \code{DataFrame} in long-format to \code{ChoiceDataset}.
The \code{EasyDataWrapper} only requires that the researcher (1) load dataframes to Python (for example, using \code{pandas} tools to load various types of data files including CSV, TSV, Stata database, and Excel spreadsheet) and (2) usage of \code{pandas} to pre-process the dataset.
We load the car-choice dataset to demonstrate \code{torch-choice}'s data wrapper utility. The car-choice dataset is a synthetic dataset on consumers' choices regarding the nationality of cars; the dataset was modified from Stata's tutorial dataset.
The dataset is in “long” format: multiple rows constitute a single choice record.
Table (ref) presents the first 30 rows of the car-choice dataset, in which each user (consumer) only made one choice. The first four rows with \code{record_id == 1} correspond to the first purchasing record. It means that consumer 1 was deciding on four types of cars (i.e., items) and chose American car (since the \code{purchase == 1} in that row of American car). Not all cars were available all the time. For example, there is no row in the dataset with \code{session_id == 4 & car == "Korean"}, which means only American, Japanese, and European cars were available in the fourth session.
The \code{EasyDataWrapper} requires a long-format main dataset with the following columns:
The next step is to identify observables to be incorporated in the model. In the car-choice dataset, we want to add (1) \code{gender} and \code{income} columns as user-specific observables; (2) \code{speed} as item-specific observable; (3) \code{discount} as session-specific observable; (4) \code{price} as (session, item)-specific observable.
The above car-choice dataframe already includes observables we want to add. The first approach to adding observables is to simply identify these columns while initializing the \code{EasyDatasetWrapper} object. This is accomplished by supplying a list of column names to each of \code{\{user, item, session, itemsession\}_observable_columns} keyword arguments. For example, we assign
to specify user-specific observables from the \code{gender} and \code{income} columns of the main dataframe. The code snippet below shows how to initialize the \code{EasyDatasetWrapper} object.
The \code{data_wrapper.summary()} method can be used to print out a summary of the dataset. The \code{ChoiceDataset} object can be accessed by the wrapper using \code{data_wrapper.choice_dataset}.
Researchers can also construct dataframes and use dataframes to supply different observables. This is useful in the case of longitudinal data or when users choose items from multiple categories, so that there are multiple choice records by the same user. Suppose there are $K$ choice records for each of the $U$ users. Using a single dataframe for all variables requires a lot of memory and repeated data, while using a separate dataframe mapping user index to user observables only requires $U$ entries for each observable. Similarly, storing item information or item-session information in long-format data leads to inefficient storage.
Dataframes for user-specific observables need to contain a column of user IDs (e.g., \code{consumer_id}), this column should have exactly the same name as the column containing user indices. Similarly, for item-specific observables, the dataframe should contain a column of item names (i.e., \code{car} in this car-choice example); session identifier column \code{session_id} needs to be present for session-specific observable dataframes, and both item and session identifiers need to be present for item-session-specific observable dataframes.
The following code snippet extracts dataframes for observables from the car-choice dataframe.
Then, researchers can initialize the \code{EasyDatasetWrapper} object with these dataframes via \code{\{user, item, session, itemsession\}_observable_data} keyword arguments.
The \code{EasyDataWrapper} also support the mixture of above methods: in the following example, we supply user, item, and session-specific observables (i.e., gender, income, speed, and discount) as dataframes but the item-session specific observable (i.e., price) using the corresponding column name.
As discussed above, storing data in a long format is not memory efficient due to repeated data. \code{torch-choice} enables storage of observables as separate tensors.
This section demonstrates an example of creating a dataset from observable and index tensors. We use a randomly generated dataset with $N=10,000$ records from $U=10$ users, $I=4$ items and $S=500$ sessions. We generate random observable tensors for user, item, session, and item-session observables from independent multivariate Gaussian distributions following Equation ((ref)).\footnote{We chose 128, 64, 10, and 12 as the number of observables, this choice is arbitrary and researchers can choose any number of observables.} For simplicity, we generate data where all items are available in all sessions.
The code snippet below builds a \code{ChoiceDataset} object from randomly generated $(u^{(n)}, i^{(n)}, s^{(n)})$ and observable tensors. This construction assumes that all items are available in all sessions, since the \code{item_availability} argument is not specified.
The \code{ChoiceDataset} is expecting the following arguments at its initialization.
As mentioned previously, all observable tensors should meet shape requirements in Table (ref).
Further, the data structure in \code{torch-choice} package utilizes variable name prefixes to keep track of observables' variation. Observable tensors need to follow naming conventions for \code{ChoiceDataset} to build correctly. While constructing the \code{ChoiceDataset} with observable tensor keyword arguments, observable tensors are named with prefixes indicating how they depend on the user (with \code{user_} prefix), item (with \code{item_} prefix), session (with \code{session_} prefix), and both of (item, session) (with \code{itemsession_} prefix) .
For example, suppose we have a tensor \code{X} for user income of shape $(U, 1)$. The income observable should be introduced as \code{ChoiceDataset(..., user_income=X, ...)}, initializing the dataset with \code{ChoiceDataset(..., income_of_user=X, ...)} using an incorrect keyword argument will introduce errors. Table (ref) summarizes the naming requirement for different types of observable tensors.
When multiple observable tensors of the same level of variation are needed, researchers can include as many keyword arguments as needed, as long as they follow the prefix requirement. For example, \code{ChoiceDataset(..., user_income=X, user_age=Z, user_education=E...)} includes three different user-specific observable tensors. Such a feature is particularly useful when you want, for example, an item-specific effect on users' income but a constant effect on users' education levels.
This section overviews the functionalities of the \code{ChoiceDataset} class. The command \code{print(dataset)} provides a quick overview of shapes of tensors included in the object as well as where the dataset is located (i.e., in CPU memory or GPU memory). The \code{dataset.summary()} method offers more detailed information about the dataset.
The \code{dataset.num_\{users, items, sessions\}} attribute can be used to obtain the number of users, items, and sessions, which are determined automatically from the \code{\{user, item, session\}_obs} tensors provided while initializing the dataset object. The \code{len(dataset)} reveals the number of records $N$ in the dataset.
The \code{ChoiceDataset} offers a \code{clone} method to make a copy of the dataset; researchers can modify cloned dataset safely without changing the original dataset.
\code{ChoiceDataset} objects can be moved between host memory (i.e., CPU memory) and device memory (i.e., GPU memory) using the \code{dataset.to()} method. The \code{dataset.to()} method moves all tensors in the \code{ChoiceDataset} to the target device. The following code runs only on a machine with a compatible GPU as well as a GPU-compatible version of PyTorch installed. The \code{dataset._check_device_consistency()} method checks if all tensors are on the same device. For example, an attempt to move the \code{item_index} tensor to CPU without moving other tensors will result in an error message.
The \code{ChoiceDataset.x_dict} attribute reconstructs the long-format representation of the dataset, which is useful for data exploration, visualization, and designing customized models upon \code{ChoiceDataset}. \code{dataset.x_dict} is a dictionary mapping observable names to tensors with shape \code{(N, num_items, *)}, where \code{N} is the length of dataset, \code{num_items} is the number of items, and \code{*} is the dimension of the corresponding observable. The helper function \code{print_dict_shape} prints the shapes of tensors in a dictionary.
\code{dataset[indices]} can be applied with \code{indices} as an integer-valued tensor or array to retrieve the corresponding rows of the dataset.\footnote{In Python, the square bracket indexing operation implicitly calls the \code{__getitem__} method of the object. For example, \code{dataset[indices]} is equivalent to \code{dataset.__getitem__(indices)}, we occasionally call the subset method as \code{__getitem__}.} The subset \code{ChoiceDataset} contains \code{item_index}, \code{user_index}, and \code{session_index} corresponding to $(u^{(n)}, i^{(n)}, s^{(n)})_{n \in \text{selected indices}}$ and full copies of observable tensors. The sub-setting method allows researchers to easily split dataset into training, validation, and testing subsets.
The example code block below randomly queries five records from the dataset.
It is worth noting that the subset method automatically creates a clone of the datasets so that any modification applied to the subset will not be reflected on the original dataset. The code block below demonstrates this feature.
The nested logit model requires two separate datasets, one for the categorical level information and the other for the item level information. The \code{JointDataset} class offers a simple way to combine multiple \code{ChoiceDataset} objects. The following code snippet demonstrates how to create a joint dataset from two separate datasets, in which one dataset is named as \code{item_level_dataset} and the other one is named as \code{cateogry_level_dataset}. Researchers can name these datasets arbitrarily as long as names are distinct.
The joint dataset structure allows easy and consistent indexing on multiple datasets. The indexing method of \code{JointDataset} will return a dictionary of subsets of each dataset encompassed in the joint dataset structure. Specifically, \code{joint_dataset[indices]} returns a dictionary \code{\{"item_level_dataset": item_level_dataset[indices], "category_level_dataset": category_level_dataset[indices]\}}. The chaining functionality is particularly helpful when the researcher wishes to experiment with the model on multiple datasets.
The \code{ChoiceDataset} is a subclass of PyTorch's dataset class and researchers with PyTorch expertise can use the PyTorch data loader to iterate over the dataset with customized batch size and shuffling options using PyTorch's \code{Sampler} and \code{DataLoader}.\footnote{This section covers advanced materials for researchers who want to build models upon \code{ChoiceDataset}. The \code{torch-choice} offers utilities for model estimation and researchers can safely skip this section.} This section is intended for researchers who wish to develop their model training pipeline; researchers who are interested in the standard training loop provided by \code{torch-choice} and \code{bemb} can skip this section.
For demonstration purposes, we turned off the shuffling option. We first create an index sampler and then use it to create a data loader.
The following code snippet iterates through the \code{dataloader} to get the mini-batches.
The \code{torch-choice} package implements two classes of models, namely the Conditional Logit Models (CLM) and Nested Logit Models (NLM). Researchers specify the utility as a function of coefficients and observables. There are multiple ways to handle user random coefficients in existing packages; for example in \code{Stata}, one can model a distribution $\mathcal{D}$ and treat user-specific coefficients as samples from this distribution, $\theta_u \sim \mathcal{D}$, and \code{Stata} then learns the parameters of the distribution $\mathcal{D}$. The \code{torch-choice} package estimates a real-valued coefficient for every user, $\bm{\theta}_u$, instead of a distribution of coefficients over all users. This is useful when learning user level coefficients in panel data settings, and an advantage of \code{torch-choice} is that in these settings, it directly reports per user coefficients.\footnote{In contrast, users need to write their own code to estimate these per user coefficients in \code{Stata}.} In contrast, since \code{torch-choice} does not support learning a distribution over user level parameters, it cannot be used to learn a random effects model with cross-sectional data, since user level parameters would not be well identified in this case with only one record per user. For example, suppose the researcher specifies a functional form $\mu_{uis} = \theta_u \times \text{Price}_i$, models in \code{torch-choice} will estimate $U$ scalar coefficients, one $\theta_u \in \mathbb{R}$ for each user $u$. The coefficient specific to user $u$ is estimated using all choice records from user $u$, and the coefficient is the same across all sessions in which user $u$ is involved. The researcher can use the model's summary method to obtain coefficients estimated for all users.
The CLM admits a linear form of the utility function. Equation ((ref)) demonstrates the general form of the utility function for the conditional logit model, in which each $\bm{\theta}_p$ is a latent coefficient of the model to be estimated, and $\bm{x}_p$ is an observable vector. The coefficient $\bm{\theta}_p$ can be fixed for all $(u, i, s)$, user-specific, or item-specific and the observable $\bm{x}_p$ can be user-specific, item-specific, session-specific, or (session, item)-specific. Please note that the researcher can set $\bm{x}_1 = 1 \in \mathbb{R}$ to include an intercept term (including user-fixed effect and item-fixed effect) in the utility function.
The CLM assumes that all items are in the same category, and the probability for user $u$ to choose item $i$ under session $s$ is given by the logistic function in Equation ((ref)).
One notable property of the CLM is the independence of irrelevant alternatives (IIA). The IIA property states that the relative probability for choosing item $i$ and $j$ is independent from any other item, which is evident from Equation ((ref)).
This section demonstrates the usage of CLM model with the Mode-Canada transportation choice dataset bhat1995heteroscedastic-modecanada The transportation choice dataset contains 3,880 travelers' choices of traveling method for the Montreal-Toronto corridor. The dataset includes the cost, frequency, in-vehicle and out-vehicle time (IVT and OVT) of four transportation methods (air, bus, car, and train), and individual income as observables.
We assign a separate session to each choice record for the transportation dataset. The cost, frequency, OVT, and IVT are session-item-specific observables. The income of travelers only depends on the session.\footnote{We can also include the income variable as an user-specific observable. However, having income as a session-specific observable is more flexible and allows for different income levels of the same user but in different sessions.} We model the utility as the following.
where $X^{itemsession: (cost, freq, ovt)}_{is} \in \mathbb{R}^3$ denotes the vector of cost, frequency, and in-vehicle time of transportation choice $i$ in session $s$. $X^{session:income}_s \in \mathbb{R}$ denotes the income of decision maker in session $s$, and $X_{is}^{itemsession:ivt} \in \mathbb{R}$ represents the in-vehicle time of choice $i$ in session $s$.
The \code{torch-choice} offers two methods for creating the model. Researchers can specify the model in Equation (ref) with a \code{formula} similar to that used in the R programming language, as follows:
More generally, the \code{formula} string specifies the functional form of $\mu_{uis}$ as additive terms \code{(observable|variation)} as the following.
The \code{observable_i} is the observable name (i.e., name of $\bm{x}_p$ in \code{ChoiceDataset}), The \code{variation_i} specifies how the coefficient (i.e., $\bm{\theta}_p$) depends on the user, item, and session. The \code{variation} takes the following possible values. The coefficient can be (1) \code{constant} for all users and items; (2) \code{user} fits a different coefficient for each user; (3) \code{item-full} setting fits a different coefficient for each item (i.e., the item with index 0); or (4) \code{item} configuration is similar to the \code{item-full} setting, but the coefficient to the first item will be set to zero. The researcher needs to pass in the \code{dataset} so that \code{model} can infer dimensions of $\bm{x}_p$'s.
Alternatively, the code snippet below shows how to create the \code{ConditionalLogitModel} class with Python dictionaries. The dictionary-based method is more suitable if researchers want to create multiple model configurations systematically.
The researcher needs to specify three keyword arguments while initializing the \code{ConditionalLogitModel} class. To fit a CLM on \code{ChoiceDataset}, there should be an attribute of the choice dataset for every observable $\bm{x}_p$. For each observable, the user can specify the coefficient variation (i.e., the coefficient type) and the number of parameters for the observable, which is the same as the dimension of the observable due to the inner product.
Our estimation procedure maximizes the log-likelihood in Equation (ref) by updating coefficients iteratively using algorithms such as stochastic gradient descent.
The \code{torch-choice} package currently supports optimization algorithms in PyTorch, which covers first-order methods such as Adam and quasi-Newton methods like LBFGS. Since \code{torch-choice} is designed for large-scale datasets and models, second-order methods are not currently supported due to computational constraints.
Researchers often want to regularize magnitudes of coefficients, especially when the number of parameters is large. The \code{torch-choice} implements \code{ConditionalLogitModel} with the support for the commonly used $L^1$ and $L^2$ regularization {krogh1991simple}. With $L^p$ regularization and regularization weight $\lambda > 0$, the objective function to be maximized changes from log-likelihood to
where $\bm{\theta}$ denotes the set of all coefficients of the model.\footnote{When regularization is used, the model is trained by minimizing the Cross-Entropy loss objective (along with an added regularization term), instead of maximizing an MLE (Maximum Likelihood Estimation) objective, which is what is done without regularization. This is achieved in practice in the package by maximizing the negative of the regularized Cross Entropy objective. The Cross-Entropy objective, without regularization, minimizes the negative likelihood for the choice models we are working with.}
Researchers can add regularization terms by simply specifying \code{regularization} $\in$ \{\code{"L1"}, \code{"L2"}\} and \code{regularization_weight} $\in \mathbb{R}_{+}$ in the model creation call. As usual, researchers may select the value of the regularization parameter using cross-validation. The following model initialization code adds $L^1$ regularization with $\lambda = 0.5$ to the model estimation.
The \code{torch-choice} package offers a model training pipeline. The following code snippet shows how to train a CLM with the \code{run()} function, in which \code{batch_size=-1} means the gradient descent is performed on the entire dataset.\footnote{If you don't know what optimizing with mini-batches means, you leave the \code{batch_size} argument to be -1, which is the default setting.} We use the LBFGS optimizer in this example since we are working on a small dataset with only 2,779 choice records and 13 coefficients to be estimated. We recommend using the Adam optimizer (the default optimizer) instead for larger datasets and more complicated models.
The \code{run(...)} function reports the model provided, the estimation progress, the final log-likelihood, and table summarizing the estimated coefficients. Due to the limited space, we only present the coefficient estimation table in this paper; readers can refer to the supplementary Jupyter Notebook for the complete output.
After the model estimation begins, \code{torch-choice} creates a TensorBoard log to illustrate the training progress; researchers can start a TensorBoard using the following command and navigate to \code{http://localhost:6006} in the browser to obtain training progress; Figure (ref) shows an example of the training progress curve, the x-axis is the number of epochs, and the y-axis is the log-likelihood. Please note that Figure (ref) is only an illustrative example, the actual training progress curve may vary depending on the dataset and model configuration.
The \code{torch-choice} package provides utilities to help researchers easily retrieve estimated coefficients from fitted models. Consider the following stylized model configured with the code block.
After the model estimation, the \code{get_coefficient()} method allows researchers to retrieve the coefficient values from the model using the general syntax \code{model.get_coefficient(COEFFICIENT_NAME)} with \code{COEFFICIENT_NAME} from the \code{formula} string. For example,
Researchers can refer to our online documentation \href{https://gsbdbi.github.io/torch-choice/conditional_logit_model_mode_canada/}{Tutorial: Conditional Logit Model on ModeCanada Dataset} for a complete Jupyter Notebook tutorial.
When items can naturally be clustered into nests, the nested logit model becomes a suitable option.\footnote{A set of nests is a partition of item sets. When the model requires the user to choose only one item across all partitions each time, these partitions are called nests. When such a requirement is removed, the model allows users to choose an item from every partition, these partitions are referred to as categories.} The nested logit model assumes a weaker version of independence of irrelevant alternatives (IIA): the IIA assumption holds for each pair of items belonging to the same nest but not for pairs from different nests.
We implement a two level nested logit model, where at the higher nesting level the user chooses one of many nests, and at the lower level, the user chooses an item within a nest. Thus, the user selects exactly one item across all nests. In the current version of the software, the nested logit is available only in a setting where there is a single category.\footnote{Our nested logit setup can be constrasted with modeling multiple categories using the multinomial logit in \code{torch-choice}, where we model choice over multiple categories simultaneously, and the user chooses one item from each category simultaneously. In the future we aim to support modeling choice over multiple categories where each category can include a two level nested logit.} The nested logit model decomposes the utility of choosing item $i$ into the item-specific values and inclusive values for each nest as in Equation ((ref)). This effectively allows correlation among the $\varepsilon_{uis}$ within the choices in a nest. Nests are indexed by $k \in \{1, 2, \dots, K\}$, $\mathcal{K}(i) \in \{1, 2, \dots, K\}$ denotes the nest of item $i$, and $I_k \subseteq \{1, 2, \dots, I\}$ denotes the set of items in nest $k$. The \code{NestedLogitModel} constructor requires a \code{nest_to_item} dictionary mapping $k \mapsto I_k$, where keys of \code{nest_to_item} are nest indices $k$'s and \code{nest_to_item[j]} is a list consisting of item indices in $I_k$. Researchers do not need to supply a \code{num_items} argument to the \code{NestedLogitModel} class, the number of items is automatically inferred from the \code{nest_to_item} dictionary.
This section briefly discusses the derivation of the log-likelihood formula \code{torch-choice} maximizes. Researchers can refer to Discrete Choice Methods with Simulation by Kenneth Train for a more detailed treatment {train2009discrete}.
Let $\mathcal{V}_{uis} = \mu_{uis} + \varepsilon_{uis}$ denote the utility user $u$ gains from item $i$ in the context of session $s$. The NLM allows $\varepsilon_{uis}$ of items from the same nest to be correlated, but $\varepsilon_{uis}$ of items from different categories should be independent. The $\lambda_k$ denotes the degree of independence among $\varepsilon_{uis}$'s of items in nest $k$. Equation (ref) illustrates the likelihood of item $i$ from the nested logit model.
Without loss of generality, the utility $\mathcal{V}_{uis}$ can be decomposed into a nest-level component $W_{u \mathcal{K}(i) s}$ and an item-level component $T_{uis}$.
The likelihood in Equation (ref) decomposes into a marginal probability of choosing nest $\mathcal{K}(i)$ and a conditional probability of choosing $i$ given nest $\mathcal{K}(i)$ is chosen.
The inclusive value of nest $k$ (for user $u$ in session $s$), denoted as $\mathcal{I}_{uks}$, is defined as $\log \sum_{j \in I_k} \exp(T_{ujs}/\lambda_k)$, which is the expected utility from choosing the best alternative from nest $k$ given Gumbel error term in utility.
The softmax form of marginal and conditional probabilities in Equation (ref) suggests a way to model $\mathcal{V}$: estimate both $W_{u \mathcal{K}(i) s}$ and $T_{uis}$ using linear functional forms as in the conditional logit model with nest observables and item observables, respectively.
where $\bm{\theta}_p$'s are coefficients of the nest-level model that can be (1) fixed for all items and users, (2) user-specific, and (3) item-specific (i.e., nest-specific since we are talking about the nest level model). $\bm{x}_p$'s are nest level observables. The nest-level model in NLM is similar to a CLM for choices over categories as items. The only difference between specifying the nest-level model and a CLM is that an item-specific coefficient or observable in the NLM nest-level model is indeed nest-specific. For example, by specifying \code{nest_formula="(1|item)"} for the nest-level model, the model specification in Equation (ref) is $W_{uks} = \lambda_c$. $\bm{\gamma}_p$'s denote coefficients of item-level model and $\bm{z}_p$'s are item observables. Specification of the item-level model is exactly the same as in the CLM.
The configuration of NLM requires specifying both the item-level CLM and nest-level CLM. Similar to the CLM, researchers can configure models using Python dictionaries or R-like formulas. For example,
specify the following model
It is worth mentioning that everything tagged with \code{item} in the \code{nest_formula}, such as item-specific intercepts/coefficients (i.e., \code{(1|item)} and \code{(user_obs1|item)}) and item-specific observables (i.e., \code{(item_obs1|constant)}) are indeed nest-specific with $k$ sub-script in Equation (ref). For example, by specifying nest-level formula \code{"(1|item)+(user_obs1|item)"}, it is interpreted as nest-specific fixed effect together with nest-specific random effect from user observable 1. Equivalently, $W_{uks} = \tau_k + \delta^\top_k \bm{x}^{(\text{user obs 1})}_u$.
Similar to specifying CLMs, researchers can use \code{\{nest, item\}_coef_variation_dict} and \code{\{nest, item\}_num_param_dict} to specify the functional form of $W_{uks}$ and $T_{uis}$ respectively.
Due to the decomposition in Equation (ref), the log-likelihood for user $u$ to choose item $i$ in session $s$ can be written as Equation (ref), and \code{torch-choice} jointly optimizes all coefficients to maximize Equation (ref).
\code{torch-choice} allows for empty nest-level models from an empty dictionary or empty formula string. In this case, $W_{uks} = \varepsilon_{uks}$ since the inclusive value term $\lambda_k I_{uks}$ will be used to model the choice over categories. In this special case, the model's log-likelihood function reduces to Equation ((ref)).
However, by specifying an empty item-level model ($T_{uis} = \varepsilon_{uis}$), the nested logit model reduces to a conditional logit model of choices over categories. Hence, researchers should never use the \code{NestedLogitModel} class with an empty item-level model.
The $\lambda_k$'s capture the degree of independence among unobserved utilities $\varepsilon$ of items in nest $k$ {train1978goods}. \code{torch-choice} estimates $\lambda_k$'s jointly with other coefficients. By default, these paramaters differ over categories, indicating different levels of correlations among unobserved utilities of items in different categories. Researchers can specify \code{shared_lambda=True} to constrain $\lambda_k = \lambda$ for all categories, indicating the same level of correlations in all categories. Whether to add such a constraint is an empirical question depending on specific datasets; researchers can use hypothesis testing to determine the right configuration {train2009discrete}. While working with large datasets, setting the same $\lambda$ for all categories helps reduce the number of parameters when the number of categories is potentially very large.
The following code snippets demonstrate the construction of nested logit models in \code{torch-choice}.
Similar to the CLMs, researchers can add $L^1$ or $L^2$ regularizations to NLMs by specifying \code{regularization} $\in$ \{\code{"L1"}, \code{"L2"}\}and \code{regularization_weight} $\in \mathbb{R}_{+}$ keywords. For example, the following code creates a model that penalizes $\|\bm{\theta}\|_2$ during estimation.
The \code{NestedLogitModel} admits a \code{JointDataset} structure, which consists of separate \code{ChoiceDataset} objects holding observables for categories and items. Naming conventions for \code{ChoiceDataset} for item observables are exactly the same as before: researchers can specify item observables \code{item_<obs_name>}, user observables \code{user_<obs_name>}, session observables \code{session_<obs_name>}, and (item, session)-specific observables \code{item_session_<obs_name>}. For the dataset holding nest information, researchers need to add the \code{item_index} array $(i^{(n)})_{n=1}^N$, and all nest-specific observables should be named as \code{item_<obs_name>} because categories are considered as items in the nest-level model. Similarly, researchers can add (nest, session)-specific observables via the \code{item_session_<obs_name>} keyword. Lastly, the researcher can link two \code{ChoiceDataset} objects using the \code{JointDataset} class as in the following code block.
For NLM estimation, \code{torch-choice} maximizes the likelihood defined in Equation (ref) by updating all coefficients iteratively using optimization algorithms like stochastic gradient descent. To estimate NLM, researchers need to provide a joint dataset encompassing both nest and item-level datasets. The \code{run()} helper function works for NLM and joint datasets.
Similar to the conditional logit model, researchers can specify regularization on coefficients by passing \code{regularization} and \code{regularization_weight} into the \code{run()} method.
Researchers can refer to our online documentation here - \href{https://gsbdbi.github.io/torch-choice/nested_logit_model_house_cooling/}{Random Utility Model (RUM) Part II: Nested Logit Model} for complete tutorials on applying various NLM to the Heating and Cooling System Choice in Newly Built Houses in California dataset {train2009discrete}.
The nested logit model has a similar interface for coefficient extraction to the conditional logit model demonstrated above. Consider a nested logit model consisting of an item-level model in Equation (ref) and a nest-level model incorporating user-fixed effect, category-fixed effect (i.e., specified by \code{(1|item)} in the \code{nest_formula}), and the user-specific coefficient on a 64-dimensional nest-specific observable (i.e., specified by \code{(item_obs|user)} in the \code{nest_formula}).
Researchers need to retrieve the coefficients of the nested logit model using the \code{get_coefficient()} method with an additional \code{level} argument.
For example, researchers can use the following code snippet to retrieve the coefficient of the user-fixed effect in the nest-level model, which is a vector with \code{num_users} elements.
With \code{level="item"}, the researcher can obtain the coefficient of user-specific fixed effect in the item level model (i.e., $\gamma_u$'s in Equation (ref)), which is a vector with length \code{num_users}.
Such API generalizes to all other coefficients listed above, such as \code{itemsession_obs[item]} and \code{itemsession_obs[user]}. One exception is the coefficients for inclusive values (often denoted as $\lambda$). Researchers can retrieve the coefficient of the inclusive value by using \code{get_coefficient("lambda")} without specifying the \code{level} argument. In fact, \code{get_coefficient} will disregard any \code{level} argument if the coefficient name is \code{lambda}. The returned value is a scalar if \code{shared_lambda=True}, and a 1D array of length \code{num_nests} if \code{shared_lambda=False}.
Leveraging parallel processors in modern graphical processing units (GPU), the \code{torch-choice} package is particularly efficient dealing with a large amount of parameters. For example, estimating the item-specific slope $\beta_i^{(2)}$ in Equation (ref) becomes computationally expensive when we have a large number of items. This section examines the computational efficiency of \code{torch-choice} as the dataset grows larger in different dimensions.\footnote{Through this section, \code{torch-choice} uses full-batch Adam optimization with a learning rate of 0.01 to maximize the sample likelihood of the dataset. The optimization stops at an epoch when the sample likelihood fails to improve for more than 50 epochs. The optimization runs for up to 5,000 epochs to ensure full algorithm convergence. Experiments related to \code{torch-choice} are run on an Nvidia RTX3090 GPU. R benchmarks were conducted on a server with 16 cores and 128GiB of memory.}
For the analysis in this section, we create a synthetic seed dataset of consumers' choices; then, we subsample datasets of various sizes to measure the runtime performance of packages.
\paragraph{User Features}
Let $\mathbf{x}_u \in \mathbb{R}^{d}$ denote the $d$-dimensional feature vector for user $u \in \{1,2,\ldots,U\}$, where we simulate $U = 500$ users and $d = 30$ features. User features $\bm{x}_u$ are generated from a mixture of $K_u = 5$ Gaussian clusters:
where $\sigma_u = 0.5$ is the within-cluster standard deviation. Cluster centers are sampled from a Gaussian distribution $\textbf{c}_k^u \sim \mathcal{N}(\mathbf{0}, \delta_u^2\mathbf{I})$ with between cluster distance $\delta_u = 3$. Cluster assignments $\pi_k^u$ are determined by uniform allocation with random perturbation so that the $\{\pi_1^u, \dots, \pi_{K_u}^u\}$ vector is one-hot for each user $u$. We selected $K_u = 5$ to represent distinct market segments with clear separation ($\delta_u = 3$) while maintaining heterogeneity within cluster ($\sigma_u = 0.5$). This parameterization yields well-formed clusters.
\paragraph{Item Features} Item features are generated analogously to these user features. Let $\mathbf{z}_i \in \mathbb{R}^{d}$ denote the feature vector for item $i \in \{1,2,\ldots,I\}$, we simulated $I = 500$ items with $d = 30$ item-specific features.
where $\sigma_i = 0.3$ is the within-cluster standard deviation. Cluster centers $\mu_k^i \sim N(\mathbf{0}, \delta_i^2 \mathbf{I})$ with $\delta_i = 2.5$. The cluster assignments $\pi_k^i$ are determined by uniform allocation with random perturbation and each item belongs to a single cluster. The smaller within-cluster variance ($\sigma_i = 0.3$) and between-cluster variance ($\delta_i = 2.5$) reflect that product categories are typically more tightly defined than consumer segments. The 30-dimensional feature space was chosen to capture the complex preference structure while remaining computationally feasible.
\paragraph{User Price Sensitivity Coefficients} In our simulation, each user $u$ has a price sensitivity coefficient $\beta_u$ generated as:
where $\mu_\beta = -1.5$ and $\sigma_\beta = 0.5$.
\paragraph{Item Prices} For each item $i$ in session $s \in \{1,2,\ldots,S\}$, where $S = 1,000$, prices $p_{s,i}$ are generated from the following steps. We first generate a base-price $p_i^\text{base} = \max(\eta_i, 1.0)$ with $\eta_i \sim \mathcal{N}(10, 3^2)$ for each item $i$. This base price is constant across sessions. In each session $s$, the price of item $i$, $p_{s, i}$, is constructed by adding symmetric noise with much smaller variance $\varepsilon_{s,i} \sim \mathcal{N}(0, (0.6)^2)$ to the base-price: $p_{i, s} = p_i^\text{base} + \varepsilon_{s,i}$. Finally, we take $p_{s, i} = \max(p_{s, i}, 1)$ to make sure all prices are positive.
\paragraph{Item Availabilities} We assume a 95% availability of items in our simulation. Specifically, item availability $A_{s,i} \in \{0,1\}$ in session $s$ is:
with the constraint that $\sum_i A_{s,i} \geq 1$ for all $s$ (i.e., the item is at least available in one session).
\paragraph{Choice Generation Process} We generate $3,000,000$ choice records, each associating a randomly selected (i.e., sampling from uniform distributions) user $u^{(n)}$, session $s^{(n)}$. The corresponding chosen item $i^{(n)}$ is generated using a hybrid choice model combining exploration and exploitation to mimic real-world consumer choices. The utility of item $i$ for user $u^{(n)}$ in session $s^{(n)}$ is:
where $\varepsilon_{u_n,i,s_n} \sim \mathcal{N}(0, 0.5^2)$ is a random utility shock.
With probability 0.15, the user makes a purely exploratory choice by randomly picking an available item:
Otherwise, with a probability of 0.85, we further create two nested cases. With probability 0.35 ($0.85 \times 0.35 = 0.2975$ total probability), the item chosen is sampled from a softmax with temperature $T = 2.0$:
With probability 0.65 ($0.85 \times 0.35 = 0.5525$ total probability), choice follows the maximum utility rule:
The exploration rate of 15% and mixed decision rules capture the dual-process nature of consumer decision-making, from utility maximization to more stochastic choice patterns, ensuring the dataset exhibits realistic choice behavior across various market conditions.
The seed dataset includes 3,000,000 choice records from 500 users and 500 items; each user ($u$) and item ($i$) is associated with a 30-dimensional vector of characteristics generated randomly, denoted as $X_u \in \mathbb{R}^{30}$ and $Z_i \in \mathbb{R}^{30}$ respectively. There is no session-specific observable. We consider the models in Equations (ref) to (ref) to cover model specifications with and without random effects.
As a convention, the coefficient of the first item $\bm{\alpha}_0$ is set to zero.
Experiments in this section demonstrate the scalability of \code{torch-choice} when researchers expand (1) the number of choice records (i.e., size of dataset), (2) the dimension of user/item covariates, and (3) the cardinality of the item set. For each of these, and for each specification above, two sets of experiments are conducted: (1) we compare the performance of \code{torch-choice} and an open source implementation, \code{mlogit} implemented in R, on datasets with small-to-medium sizes;\footnote{The \code{torch-choice} and R comparison only focuses on small-scale datasets because the size of long-format CSV file grows fast as the dataset expands. For example, 200,000 records of 500 items result in 100 million rows of data with 60 columns of user/item characteristics, which takes more than 100GiB of disk space and memory, which leads to out-of-memory issues.} (2) we then examine the scalability of \code{torch-choice} when we expand to even larger datasets. We run each of these benchmarks five times and report the average and standard deviation of run times across five trials.
The first set of experiments focuses on the performance of \code{torch-choice} as the number of records ($N$) grows while the dimensions of user/item observables and the item set remain the same.
Scalability on Small Datasets. We sampled $30$ items and only used the first 10 dimensions in user and item characteristics through this experiment. The top panels in Figure (ref) show the time cost of \code{torch-choice} as the number of records $N$ grows from 3,000 to 100,000. The vertical axis shows the ratio between the time taken at each level of $N$ and the baseline time cost at $N=3,000$. Both \code{torch-choice} and \code{mlogit} scale linearly as $N$ grows, but the runtime of the R-implementation increases by 30x to 40x as the sample size rises from 3,000 to 100,000 while \code{torch-choice} only increases by 20x.
Scalability on Larger Datasets. We then turn to the complete dataset with $500$ items and $30$-dimensional user/item characteristics. The bottom panel in Figure (ref) illustrates the time taken relative to the baseline case ($N=3,000$) as $N$ increases from $3,000$ to $100,000$, which leads to a 30x increase in run time. The runtime of model $\mu_{uis} = \bm{\beta}^\top \mathbf{z}_i$ becomes vastly different from the other two because models with random effects now have 500x more parameters than model $\mu_{uis} = \bm{\beta}^\top \mathbf{z}_i$. The linear cost of \code{torch-choice} is preserved on the large-scale dataset, which provides evidence of scalability.
Absolute Time Costs. Figure (ref) depicts the time cost in seconds of \code{mlogit} and \code{torch-choice} on different datasets. The Figure illustrates that in these datasets, \code{torch-choice} provided up to a 2x reduction in run time.
The number of parameters plays an important role in computational efficiency; as the number of parameters grows, the time cost of computing the gradient of the likelihood function also grows. In this section, we examine the performance of \code{torch-choice} as the number of covariates grows. We investigate the models in Equations (ref) to (ref) using only the first $P \in \{3, 5, 10, 15, 20, 30\}$ dimensions of user features $\mathbf{x}_u$ and item features $\mathbf{z}_i$.
Scalability on Small Datasets. The first part examines different packages' time costs at each level of $P$ normalized by the baseline runtime ($P=3$) on a dataset with $N=10,000$ records and $I=50$ items. The top panels of Figure (ref) compare computational efficiencies of R and \code{torch-choice} as the number of covariates grows. For the two models with item-specific coefficients $\alpha_i$, the time taken is more than 20x higher when $P=30$ compared to the baseline case while using the R-implementation; in contrast, the time cost was less than 3x higher when $P=30$ with \code{torch-choice}. The time ratio for the $\beta^\top \mathbf{z}_i$ model using \code{torch-choice} is higher because the model converges extremely fast in the baseline setting, which can be verified in Figure (ref).
Scalability on Larger Datasets. The bottom panel of Figure (ref) illustrates the run time with different $P$'s on the large-scale dataset with 3 million records, $500$ users/items. Model $ \mu_{uis} = \bm{\alpha}_i^\top X_u + \bm{\beta}^\top Z_i$ with both random and common coefficients scaled worse than linear scaling, the run time increased by up to 4x when $P=30$. The scaling factor in time cost as $P$ expands depends on the number of items, which we will explore later, partly because the total number of parameters is $\mathcal{O}(P\times I)$. To summarize, \code{torch-choice} scales well when the number of items is moderate. Researchers should use caution while specifying item-specific random effects on high dimensional covariates when many items are present.
Absolute Time Costs. Finally, Figure (ref) compares time cost in seconds of R and \code{torch-choice} in various scenarios, in which \code{torch-choice} delivers up to 30x speed-up compared to \code{mlogit}.
The number of parameters grows as the item set expands as well when item-specific coefficients are specified. In this section, we examine the performance of \code{torch-choice} as the number of items grows.
Scalability on Small Datasets. We first compare performances, in terms of time cost relative to baseline case ($I=10$), of \code{torch-choice} and R-implementation as $I$ increases from $10$ to $200$ on a dataset with $N=10,000$ records and 5-dimensional user/item characteristics.\footnote{We choose a smaller $N$ in this experiment because we have a limited number of observations when $I=10$ and we wish to keep $N$ constant when $I$ changes.}\footnote{We stop at $I=200$ because \code{mlogit} encountered an out-of-memory error with $I=250$ on a server with 128GiB memory.} We investigate the models in Equations (ref), (ref), and (ref). The top panels in Figure (ref) show the time cost of \code{torch-choice} and \code{mlogit} as the number of items $I$ grows. Since the complexity of model $\mu_{uis} = \bm{\beta}^\top \mathbf{z}_i$ does not depend on the number of items, its runtime remains relatively constant as the number of items grows. For the other two models, \code{torch-choice} exhibits a mild runtime increase as $I$ grows, and the costs of estimating these models are only 3x-4x higher when the number of items increases by 20x from 10 to 200. In contrast, the curvatures in the panel for R implementation suggest an exponential growth trend in the time cost of \code{mlogit}; the computational cost required explodes by more than 400x-500x when the number of items rises by 20x. These observations again highlight the computational efficiency of parallel computing capability in PyTorch and \code{torch-choice} while handling a large number of parameters.
Scalability on Larger Datasets. The bottom panel in Figure (ref) demonstrates the runtime of \code{torch-choice} as $I$ increases up to 500 with $N=30,000$ and all 30 user/item features. As we expected, the cost of estimating model $\mu_{uis} = \bm{\beta}^\top \mathbf{z}_i$ grows slowly as we expand the item set. The computational costs of both model $\mu_{uis} = \bm{\alpha}_i^\top \mathbf{x}_u$ and model $\mu_{uis} = \bm{\alpha}_i^\top \mathbf{X}_u + \bm{\beta}^\top \mathbf{z}_i$ rise around 4x as we expand the item set. This provides further evidence that the computational cost using \code{torch-choice} follows a sub-linear trend with the expansion of the item set, demonstrating \code{torch-choice}'s scalability.
Absolute Time Costs. Lastly, Figure (ref) demonstrates the run time of \code{torch-choice} and \code{mlogit} in seconds through experiments considered in this section, showing superior performance of \code{torch-choice}. Especially, we observe up to 60x speed-up from the \code{torch-choice} compared to \code{mlogit} for models with item-specific coefficients.
Choice modeling problems are ubiquitous in multiple disciplines, including Economics, Psychology, Education and Computer Science. As a result, access to effective choice modeling software is crucial for both researchers and practitioners. The \code{torch-choice} package aims to greatly improve this modeling landscape by providing a flexible and scalable choice modeling implementation.
In this paper, we proposed \code{torch-choice}, a PyTorch-based choice modeling package. The \code{torch-choice} package offers a memory-efficient data structure, \code{ChoiceDataset} for storing choice datasets. Towards this end, \code{torch-choice} also provides functionalities allowing researchers to build the \code{ChoiceDataset}, from databases in various formats. For the advanced user, the \code{ChoiceDataset} object is completely flexible and supports all functionalities in PyTorch's dataset data structure.
The \code{torch-choice} package implements two widely-used classes of models, the \code{ConditionalLogitModel} and the \code{NestedLogitModel} using \code{PyTorch}. It allows researchers to fully specify flexible functional forms and changing item availability, which overcomes some of the difficulties when modeling with existing libraries. Models in \code{torch-choice} support regularization during model estimation, making it usable for engineering applications as well, a feature not provided in standard choice modeling libraries used by econometricians and statisticians. Such a feature vastly increases the scope of choice modeling options as compared to standard libraries used by engineers in terms of flexibility and scale.
The package has the advantage of scaling to massive datasets while being computationally efficient by using GPUs for estimation. We examined the scalability of \code{mlogit} in R and \code{torch-choice} while (1) increasing the number of observations, (2) increasing the number of covariates, and (3) expanding the item set with various model specifications. Experiment results show that \code{torch-choice} scales significantly better than R on small-to-medium-scale datasets. Then we test \code{torch-choice} on large-scale datasets that R cannot handle, on which \code{torch-choice} demonstrates its ability to scale up with computational efficiency.
Moving forward, we aim to maintain the reference implementation of \code{torch-choice} on the Python Package Index (PyPI) via Github, developing relatively few new features. We foresee the following important extensions to our models. First, \code{torch-choice} currently does not allow researchers to specify Bayesian priors over coefficients and learn posteriors instead of point estimates. With Bayesian coefficients, in addition to using item/user characteristics as observables, researchers can incorporate these characteristics into the coefficient prior and let them influence the coefficient posterior directly. Second, currently our nested logit model only supports two levels of nesting. Another direction is to implement multi-level nested logit models (i.e., nests within a nest), which are particularly useful when users make sequential decisions. In addition, allowing for modeling multiple categories with nested logit will be an important extension. We are actively developing these extensions and aim to release them soon.
\quad
We want to express our sincere appreciation and gratitude to the following researchers who have contributed to the completion of this paper. We thank Charles P\'ebereau, Keshav Agrawal, Aaron Kaye, Rob Kuan, Emil Palikot, Rachel Zhou, Rob Donnelly, and Undral Byambadalai for their valuable comments and feedback that help improve the quality of our package and paper.