The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
86,708 characters
Iterative Distributed Multinomial Regression
\title{Iterative Distributed Multinomial Regression\thanks{We thank Jean-Marie Dufour, Hiro Kasahara, Artem Prokhorov, Dacheng
Xiu, Jin Yan, and participants of the seminar at the University of
British Columbia, the 6th International Conference on Econometrics
and Statistics, and 2024 Econometric Society North American Summer
Meeting for helpful discussions.}}
\author{Yanqin Fan\thanks{Department of Economics, the University of Washington; email: [email removed]}
\and Yigit Okar\thanks{Department of Economics, the University of Washington; email: [email removed]}
\and Xuetao Shi\thanks{School of Economics, the University of Sydney; email: [email removed]}}
\maketitle
\begin{abstract}
This article introduces an iterative distributed computing estimator for the multinomial logistic regression model with large choice sets. Compared to the maximum likelihood estimator, the proposed iterative distributed estimator achieves significantly faster computation and, when initialized with a consistent estimator, attains asymptotic efficiency under a weak dominance condition. Additionally, we propose a parametric bootstrap inference procedure based on the iterative distributed estimator and establish its consistency. Extensive simulation studies validate the effectiveness of the proposed methods and highlight the computational efficiency of the iterative distributed estimator.
\textbf{Keywords}: Distributed computing, Iterative methods, Maximum
likelihood, Multinomial logistic regression.
\textbf{JEL Codes}: C13, C25, C61, C63
\end{abstract}
\newpage{}
\section{Introduction}
\subsection{Motivation and Main Contributions}
Discrete choice models, including logit and multinomial-logit (MNL)
models, are widely used in applied social science research. With the
growing availability of diverse data types and the integration of
econometric models with textual and image data, researchers often encounter
applications of MNL models with a massive number of choices; see the
applications discussed in Section \ref{sec:literature}. In these cases, even
if the number of parameters for each choice is small, the total number
of parameters, which increases linearly with the number of choices,
will be large due to large choice sets. As a result, maximum likelihood
estimation (MLE) can become computationally intractable due to the high cost
of solving an optimization problem for a large number of parameters.\footnote{MLE of the MNL model refers to the conditional MLE given the covariate
and total counts.}
Researchers in diverse disciplines have explored various approaches and proposed numerous
methods to numerically solve the MLE; see the numerical algorithms discussed in Section \ref{sec:literature}. However, there is a notable lack of theoretical results ensuring the consistency
or asymptotic efficiency of the estimators derived from these methods. This gap in the literature between numerical computation and statistical inference motivates the present study.
Our proposed estimator utilizes the multinomial-Poisson (MP) transformation,
which reformulates the multinomial likelihood into a Poisson likelihood
by incorporating individual fixed effects into the MNL model. When all the covariates are \textit{categorical}, \citet{baker1994multinomial} establishes that the multinomial likelihood and the Poisson likelihood produce identical estimates of the parameters in the MNL model and advocates the computational advantage of maximizing the Poisson likelihood. However, the MP transformation has also been employed in MNL models with continuous or mixed discrete and continuous covariates. For instance, \citet{gentzkow2019measuring} note that they “... approximate
the likelihood of {[}their{]} multinomial logit model with the likelihood
of a Poisson model...”.
As the first contribution of
this paper, we establish an equivalence result for all types of covariates, \textit{continuous}, \textit{discrete}, or \textit{mixed}, justifying the application of the MP transformation in these cases.
To accomplish this, we re-interpret the Poisson likelihood from the MP transformation as a conditional quasi-log-likelihood function given the covariates. Maximizing this function provides a
quasi-maximum likelihood estimator (QMLE) for the MNL model. The equivalence result is established by demonstrating that the resulting QMLE is identical to the MLE of the MNL model.
While the QMLE is computationally more efficient than the MLE, it remains costly when applied to MNL models with large choice sets. To address
this computational challenge, \citet{taddy2015distributed} exploits
an important feature of the quasi-log-likelihood function: for any
given fixed effects, the function is additively separable in the parameters
across different choices of the MNL model. \citet{taddy2015distributed}
proposes to estimate parameters for each choice separately at a specific
value of the fixed effects and calls the resulting estimator \textit{distributed
computing} estimator. As noted in \citet{taddy2015distributed}, however,
his distributed computing estimator is inconsistent except in a few
very special cases.
To regain consistency and asymptotic efficiency, we adopt the idea
of iterative backfitting algorithms studied in \citet{pastorello2003iterative},
\citet{dominitz2005some}, and \citet{fan2015maximization} to compute the
QMLE of the parameters in the MNL model and fixed effects iteratively.\footnote{Similar iterative algorithms have also been developed for dynamic
discrete games of incomplete information where directly computing
the MLE using the nested fixed-point algorithm proves computationally
infeasible, see \citet{aguirregabiria2002swapping,aguirregabiria2007sequential}
for a nested pseudo-likelihood (NPL) algorithm. \citet{kasahara2012sequential}
further analyze the conditions necessary for the convergence of the
NPL algorithm and derive its convergence rate. }
During each iteration, we first solve for the parameters of interest
through \textit{distributed computing}, given approximate values
of the individual fixed effects. Then, we update estimates of the individual
fixed effects based on their expressions derived from maximizing the quasi-log-likelihood
function using the previous estimates of the model parameters. We
call our estimator \textit{iterative distributed computing} (IDC)
estimator. This is the second and main contribution of the paper.
The IDC algorithm is fast because of distributed computing,
even when the choice set is large.
We consider three IDC estimators based on three different initial
estimators: a consistent estimator based on pairwise binomial logistic
regression, \citet{taddy2015distributed}'s distributed computing
estimator, and a maximum likelihood estimator assuming that the distribution
of total counts is Poisson. The latter two estimators are inconsistent
in general. All three initial estimators are fast to compute because
they allow for distributed computing.
As the third contribution, we establish theoretical results on the
consistency and asymptotic efficiency of all three IDC estimators.
When the initial estimator is consistent, we show that the IDC estimator
with any finite number of iterations is always consistent and is asymptotically
efficient under an information dominance condition when the number
of iterations diverges with respect to the sample size $n$ at $\log\left(n\right)$
rate. With inconsistent initial estimators, the IDC estimators are
consistent and asymptotically efficient under a stronger contraction
mapping condition when the number of iterations diverges with respect
to $n$ at a polynomial rate.
When the number of choices is large, conducting inference via plug-in
estimation of the variance matrix becomes infeasible. This is because
the Fisher information matrix has a very large dimension, causing
the computation of the inverse of its estimator to be both time-consuming
and unreliable. The fourth contribution of the paper is that we propose
a parametric bootstrap inference procedure and show its consistency.
Because the IDC estimator is fast to compute, our inference procedure
is computationally feasible.
Lastly, we conduct extensive simulations to study the finite sample
performance of our estimator and inference procedure. We are particularly
interested in the computational time of the IDC estimator and its
accuracy compared to the maximum likelihood estimator. The
simulation results show that the IDC estimator is very fast to compute,
with a running time approximately linear in the number of
choices. Compared to the maximum likelihood estimator, our estimator
has a very similar mean squared error in all the different model settings
but is much faster to compute when the number of choices is large.
We also study the finite sample behavior of the proposed bootstrap
inference procedure. The results suggest that the procedure achieves
the correct size and is consistent.
\subsection{Related Literature\label{sec:literature}}
\paragraph{Applications}
The proposed IDC estimator can be applied to study various economics
and computer science topics such as text analysis, dimensionality
reduction, spatial choice models, image classification (\citet{russakovsky2015imagenet}),
and video recommendation (\citet{davidson2010youtube}).
\emph{Text analysis}: The integration of text data into econometric
models is increasingly prominent in economics. For example, \citet{baker2006investor}
analyze investor sentiment’s effect on stock returns, while \citet{chen2021sentiment}
explore how hedge funds capitalize on sentiment changes. Modeling
text data often involves treating word counts as a multinomial distribution,
as \citet{taddy2015distributed} demonstrates using Yelp reviews to
predict outcomes based on user and business attributes. \citet{gentzkow2019measuring}
use the distributed computing estimator in \citet{taddy2015distributed} for the multinomial regression
to measure polarization in Congressional speeches. \citet{NBERwTextSelection}
extend this approach, using Hurdle Distributed Multinomial Regression
to backcast, nowcast, and forecast macroeconomic variables from newspaper
text.
\emph{Dimensionality Reduction}: Our estimator aids in dimensionality
reduction for inverse multinomial regression, as \citet{taddy2013multinomial}
discusses. Instead of inferring sentiment from text, \citet{taddy2013multinomial}’s
approach estimates word distribution given sentiment. He introduces
a score based on word frequencies and regression parameters, which is useful
in forward-regression models.
\emph{Spatial Choice Models}: High-dimensional choices also appear
in spatial models. \citet{bucholztaxi2019} models taxi drivers’ location
choices with a dynamic spatial search, reducing dimensionality via
discretization. Similarly, \citet{pellegrinifotheringham2002} apply
hierarchical discrete choice models to immigration, while \citet{bettman1979memory}
addresses brand choices in limited-option settings, proposing hierarchical
selection for high-dimensional cases.
\paragraph{Numerical Algorithms}
To address the computational difficulty of solving the MLE of the
MNL model, researchers have proposed several numerical methods to
find approximate solutions. \citet{bohning1988monotonicity} and \citet{bohning1992multinomial}
propose replacing the Hessian matrix in the Newton-Raphson iteration
with its easy-to-compute global lower bound and show that the approximate
solution converges with the number of iterations. Because the convergence
rate depends crucially on the difference between the Hessian matrix
and its lower bound, the algorithm can be slow to run for certain
model parameters. Additionally, based on our simulation exercise,
if the choice probabilities vary significantly across different choices
with some being close to zero, the algorithm becomes unstable. In
comparison, our IDC algorithm is stable in all the simulation settings.
\citet{boyd2011distributed} introduce an alternating direction method
of multipliers, which reformulates the original optimization problem
by introducing redundant linear constraints. \citet{gopal2013distributed}
propose a log concavity method, which replaces the log partition function
of the multinomial logit with a parallelizable upper bound. \citet{recht2011hogwild},
\citet{raman2016ds}, and \citet{fagan2018unbiased} study a stochastic
gradient descent method, which uses random training samples to calculate
the gradient at each iteration. Although these methods can be computationally
efficient, to the best of the authors' knowledge, no consistency or
asymptotic efficiency result has been shown in these works.
Penalization methods have also been introduced to the MNL regression
and some of the numerical methods discussed above are adopted in solving
the penalized MNL regression, see e.g., \citet{friedman2010regularization},
\citet{simon2013blockwise}, and \citet{nibbering2022multiclass},
The proposed IDC procedure in this paper can be combined with the
aforementioned algorithms to further improve the performance of penalized
MNL regressions.
\paragraph{Organization of the rest of this paper}
The remainder of this paper is organized as follows. In Section \ref{sec:Multinomial-Logistic-Regression}
we present a comprehensive overview of the multinomial logistic regression
model and the MP transformation. Section \ref{sec:Iterative-Distributed-Computing-} introduces
our iterative distributed computing estimator along with some initial
values. In Section \ref{sec:Asymptotic-Theory} we provide the asymptotic
theory of the iterative distributed computing estimator. Section \ref{sec:Monte-Carlo-Simulation}
contains the simulation results. Finally, with Section \ref{sec:Conclusion}
we conclude. Appendix \ref{sec:Notations-and-Equalities} collects
the notations and equations used in the paper. All the technical proofs
are provided in Appendix \ref{subsec:Proofs}. The codes for implementing
the estimation and inference procedures are available \href{https://github.com/yigitokar/IDMR}{here}.
\paragraph{Notations}
Throughout the paper, we use index $i\in\left\{ 1,\ldots,n\right\} $
for individual, $j\in\left\{ 1,\ldots,p\right\} $ for covariate,
and $k\in\left\{ 1,\ldots,d\right\} $ for unique choice. Boldfaced
symbols such as $\boldsymbol{C}$ and $\boldsymbol{V}$ are used to
denote vectors; while elements of the vectors are denoted by plain
symbol such as $C_{k}$ and $V_{j}$. Denote $\sim$ as ``equality
up to a constant'', such that $f\left(\theta\right)\sim g\left(\theta\right)$
is equivalent to $f\left(\theta\right)=g\left(\theta\right)+h$, where
$h$ is a constant relative to $\theta$.
\section{Multinomial Logistic Regression \label{sec:Multinomial-Logistic-Regression}}
Let $\boldsymbol{C}_{i}\in\mathbb{R}^{d}$ denote the random vector
of counts on $d$ different choices for individual $i=1,...,n$, summing
up to $M_{i}=\sum_{k=1}^{d}C_{ik}$. We use the random vector $\boldsymbol{V}_{i}\in\mathbb{R}^{p}$
to denote the covariate vector that includes a constant.
Consider a \textit{correctly specified} multinomial-logit (MNL) model.
The conditional probability mass function is given by the following:
\begin{equation}
\Pr\left(\boldsymbol{C}_{i}\mid\boldsymbol{V}_{i},M_{i}\right)=\mathrm{MNL}\left(\boldsymbol{C}_{i};\boldsymbol{\eta}_{i}^{\ast},M_{i}\right)\equiv\frac{M_{i}!}{C_{i1}!\cdots C_{id}!}\left(\frac{e^{\eta_{i1}^{\ast}}}{\Lambda_{i}^{\ast}}\right)^{C_{i1}}\ldots\left(\frac{e^{\eta_{id}^{\ast}}}{\Lambda_{i}^{\ast}}\right)^{C_{id}},\label{eq: MNL model}
\end{equation}
where for $k=1,\ldots d$, we let $\eta_{ik}^{\ast}\equiv\boldsymbol{V}_{i}^{\prime}\boldsymbol{\theta}_{k}^{\ast}$
with unknown parameters $\boldsymbol{\theta}_{k}^{\ast}\equiv\left(\theta_{k1}^{\ast},\ldots,\theta_{kp}^{\ast}\right)^{\prime}$
and $\Lambda_{i}^{\ast}\equiv\sum_{k=1}^{d}e^{\eta_{ik}^{\ast}}$.
For the identification, we set $\boldsymbol{\theta}_{d}^{\ast}=\mathbf{0}$.
Let $\boldsymbol{\theta}^{\ast}\equiv\left(\boldsymbol{\theta}_{1}^{\ast\prime},\ldots,\boldsymbol{\theta}_{d}^{\ast\prime}\right)^{\prime}\in\Theta$
denote the parameter vector of interest. Throughout the paper, we
use the superscript $\ast$ to indicate the true value of the unknown
parameter. Denote the parameter space of $\boldsymbol{\theta}_{k}^{\ast}$
for $k=1,\ldots,d$ as $\varTheta_{k}$. We have that $\varTheta_{d}=\left\{ \mathbf{0}\right\} $
and $\Theta\equiv\prod_{k=1}^{d}\varTheta_{k}$.
In this paper, we focus on the case where $d$ is large (but fixed)
such that directly solving for the maximum likelihood estimator is
computationally costly. Applications include text corpora, where $\boldsymbol{C}_{i}$
represents the counts of $d$ different words/phrases in a text of
$M_{i}$ words; browser logs, where $\boldsymbol{C}_{i}$ indicates
the number of times a website among $d$ total websites is visited
by an individual; and location choices, where among $M_{i}$ number
of locations traveled by the driver, $\boldsymbol{C}_{i}$ contains
the number of times each location, among $d$ different ones, is visited.
\subsection{Maximum Likelihood Estimation (MLE) \label{subsec:Maximum-Likelihood-Estimator}}
Given a random sample of size $n$, let $\boldsymbol{\eta}_{i}\equiv\left(\eta_{i1},\ldots,\eta_{id}\right)^{\prime}\equiv\left(\boldsymbol{V}_{i}^{\prime}\boldsymbol{\theta}_{1},\ldots,\boldsymbol{V}_{i}^{\prime}\boldsymbol{\theta}_{d}\right)$
and $\Lambda_{i}=\sum_{k=1}^{d}e^{\eta_{ik}}$ for $i=1,...,n$. Ignoring
terms that are independent of the parameter $\boldsymbol{\theta}$,
the conditional log-likelihood function given the covariate $\boldsymbol{V}$
and total count $M$ takes the following form:
\begin{eqnarray}
l_{C\mid V,M}\left(\boldsymbol{\theta}\right) & \equiv & \sum_{i=1}^{n}\log\Pr\left(\boldsymbol{C}_{i}\mid\boldsymbol{V}_{i},M_{i}\right)\nonumber \\
& \sim & \sum_{i=1}^{n}\log\left[\left(\frac{e^{\eta_{i1}}}{\Lambda_{i}}\right)^{C_{i1}}\ldots\left(\frac{e^{\eta_{id}}}{\Lambda_{i}}\right)^{C_{id}}\right]\nonumber \\
& = & \sum_{i=1}^{n}\left\{ C_{i1}\left[\log\left(e^{\eta_{i1}}\right)-\log\left(\sum_{k=1}^{d}e^{\eta_{ik}}\right)\right]+\cdots+C_{id}\left[\log\left(e^{\eta_{id}}\right)-\log\left(\sum_{k=1}^{d}e^{\eta_{ik}}\right)\right]\right\} \nonumber \\
& = & \sum_{i=1}^{n}\left[C_{i1}\eta_{i1}+\cdots+C_{id}\eta_{id}-\left(C_{i1}+\cdots+C_{id}\right)\log\left(\sum_{k=1}^{d}e^{\eta_{ik}}\right)\right]\nonumber \\
& = & \sum_{i=1}^{n}\left[\boldsymbol{C}_{i}^{\prime}\boldsymbol{\eta}_{i}-M_{i}\log\left(\sum_{k=1}^{d}e^{\eta_{ik}}\right)\right].\label{eq:Conditional log likelihood}
\end{eqnarray}
Let $L_{C\mid V,M}\left(\boldsymbol{\theta}\right)$ denote the probability
limit of $\frac{1}{n}l_{C\mid V,M}\left(\boldsymbol{\theta}\right)$.
Denote $B\left(\boldsymbol{\theta},\varepsilon\right)$ as an open
ball in $\Theta$ centered at $\boldsymbol{\theta}$ with radius $\varepsilon$.
We make the following assumption throughout the paper.
\begin{assumption} \label{Assmp:: Maintained assumption} (i) The
true value $\boldsymbol{\theta}^{\ast}\in\Theta$ satisfies that $\sup_{\boldsymbol{\theta}\notin B\left(\boldsymbol{\theta}^{\ast},\varepsilon\right)}L_{C\mid V,M}\left(\boldsymbol{\theta}\right)<L_{C\mid V,M}\left(\boldsymbol{\theta}^{\ast}\right)$
for any $\varepsilon>0$. (ii) $\boldsymbol{\theta}^{\ast}$ is in
the interior of $\Theta$. \end{assumption}
Assumption \ref{Assmp:: Maintained assumption} (i) implies that $\boldsymbol{\theta}^{\ast}$
is identified as $\boldsymbol{\theta}^{\ast}=\arg\max_{\boldsymbol{\theta}\in\Theta}L_{C\mid V,M}\left(\boldsymbol{\theta}\right)$.
Define the following objective function:
\begin{equation}
Q_{n}^{\ast}\left(\boldsymbol{\theta}\right)\equiv-l_{C\mid V,M}\left(\boldsymbol{\theta}\right).\label{eq:Objective Function - original}
\end{equation}
Based on (\ref{eq:Objective Function - original}), the conditional
maximum likelihood estimator of $\boldsymbol{\theta}^{\ast}$ is:\footnote{The definition of $\widetilde{\boldsymbol{\theta}}$ implicitly assumes
that the solution to the minimization problem is unique. This can
be shown to hold with probability approaching one by the identification
of the model. See \citet{mcfadden1973conditional}. The same result
holds for all the estimators defined in the paper. We ignore such
mathematical subtlety for the remainder of the paper to simplify the
discussion. }
\begin{equation}
\widetilde{\boldsymbol{\theta}}=\arg\min_{\boldsymbol{\theta}\in\Theta}Q_{n}^{\ast}\left(\boldsymbol{\theta}\right).\label{eq:MLE}
\end{equation}
Solving the above optimization problem analytically is impossible.
In addition, due to the potentially large dimension $d$, numerical
algorithms such as the Newton-Raphson method are difficult to implement
either because they usually involve computing the inverse of the Hessian
matrix, which is of dimension $pd\times pd$, during each iteration.
In this paper, we propose an estimator that is both computationally
attractive and asymptotically efficient.
\subsection{Multinomial-Poisson Transformation}
In this section, we present the multinomial-Poisson (MP) transformation,
based on which we develop our estimator. We reinterpret the Poisson likelihood as a quasi-likelihood conditional on the covariates.
Let $\boldsymbol{1}_{d}\equiv\left(1,\ldots,1\right)^{\prime}\in\mathbb{R}^{d}$.
With a slight abuse of notation, define
\begin{align*}
l_{C\mid V,M}\left(\boldsymbol{\theta},\boldsymbol{\mu}\right) & \equiv l_{C\mid V,M}\left(\boldsymbol{\theta}\right)+\sum_{i=1}^{n}\left[\mu_{i}\boldsymbol{C}_{i}^{\prime}\boldsymbol{1}_{d}-M_{i}\log\left(e^{\mu_{i}}\right)\right]\\
& =\sum_{i=1}^{n}\left[\boldsymbol{C}_{i}^{\prime}\left(\boldsymbol{\eta}_{i}+\mu_{i}\boldsymbol{1}_{d}\right)-M_{i}\log\left(\sum_{k=1}^{d}e^{\eta_{ik}+\mu_{i}}\right)\right],
\end{align*}
where $\boldsymbol{\mu}\equiv\left(\mu_{1},\ldots,\mu_{n}\right)\in\mathbb{R}^{n}$.
The following lemma shows that the two functions $l_{C\mid V,M}\left(\boldsymbol{\theta}\right)$
and $l_{C\mid V,M}\left(\boldsymbol{\theta},\boldsymbol{\mu}\right)$
are the same for any $\boldsymbol{\theta}\in\Theta$ and $\boldsymbol{\mu}\in\mathbb{R}^{n}$.
In other words, argument $\boldsymbol{\mu}$ in $l_{C\mid V,M}\left(\boldsymbol{\theta},\boldsymbol{\mu}\right)$
does not affect the value of the function. The proof of the lemma
is straightforward by realizing that $\boldsymbol{C}_{i}^{\prime}\boldsymbol{1}_{d}=M_{i}$
by definition.
\begin{lemma} \label{lem:: Adding mu} $l_{C\mid V,M}\left(\boldsymbol{\theta}\right)=l_{C\mid V,M}\left(\boldsymbol{\theta},\boldsymbol{\mu}\right)$
for any $\boldsymbol{\theta}\in\Theta$ and $\boldsymbol{\mu}\in\mathbb{R}^{n}$.
\end{lemma}
Define the following two functions:
\begin{align}
f\left(\boldsymbol{\theta},\boldsymbol{\mu}\right) & \equiv\sum_{i=1}^{n}\left[M_{i}\log\left(\sum_{k=1}^{d}e^{\eta_{ik}+\mu_{i}}\right)-\sum_{k=1}^{d}e^{\eta_{ik}+\mu_{i}}\right]\textrm{ and }\label{eq:log likelihood l_M|V}\\
ql_{C\mid V}\left(\boldsymbol{\theta},\boldsymbol{\mu}\right) & \equiv l_{C\mid V,M}\left(\boldsymbol{\theta},\boldsymbol{\mu}\right)+f\left(\boldsymbol{\theta},\boldsymbol{\mu}\right)\nonumber \\
& =\sum_{i=1}^{n}\sum_{k=1}^{d}\left(C_{ik}\left(\eta_{ik}+\mu_{i}\right)-e^{\left(\eta_{ik}+\mu_{i}\right)}\right).\label{eq:Unconditional Likelihood with mu}
\end{align}
It is not difficult to see that $f\left(\boldsymbol{\theta},\boldsymbol{\mu}\right)$
takes the form of a log-likelihood function of $n$ conditional Poisson
distributions with means $\sum_{k=1}^{d}e^{\left(\eta_{ik}+\mu_{i}\right)}$,
$i=1,\ldots,n$ (after ignoring terms that are independent of $\boldsymbol{\theta}$
and $\boldsymbol{\mu}$). This in turn renders $ql_{C\mid V}$ a conditional
quasi-log-likelihood function of which $C_{ik}$ given $V_{i}$ is
drawn independently from a Poisson distribution with mean $e^{\left(\eta_{ik}+\mu_{i}\right)}$,
$k=1,...,d$. This property underlies the naming of the MP transformation.
Based on the conditional quasi-log-likelihood function $ql_{C\mid V}\left(\boldsymbol{\theta},\boldsymbol{\mu}\right)$,
we can compute a conditional quasi MLE (QMLE) of $\boldsymbol{\theta}^{\ast}$:
\begin{align}
\left(\widehat{\boldsymbol{\theta}},\widehat{\boldsymbol{\mu}}\right) & =\arg\min_{\boldsymbol{\theta}\in\Theta,\boldsymbol{\mu}\in\mathbb{R}^{n}}Q_{n}\left(\boldsymbol{\theta},\boldsymbol{\mu}\right),\label{eq:Definition of Distributed Multinomial Estimator}\\
\text{where }Q_{n}\left(\boldsymbol{\theta},\boldsymbol{\mu}\right) & \equiv-ql_{C|V}\left(\boldsymbol{\theta},\boldsymbol{\mu}\right).\nonumber
\end{align}
\citet{baker1994multinomial} shows that $\widehat{\boldsymbol{\theta}}=\widetilde{\boldsymbol{\theta}}$
if the covariate vector $\boldsymbol{V}$ contains only categorical
random variables. The following lemma demonstrates that $\widehat{\boldsymbol{\theta}}=\widetilde{\boldsymbol{\theta}}$
holds irrespective of the type of the covariate vector, thereby generalizing
the result of \citet{baker1994multinomial}.
\begin{lemma} \label{lem:: Equivalence conditional and unconditional}
It holds that $\widehat{\boldsymbol{\theta}}=\widetilde{\boldsymbol{\theta}}$.
\end{lemma}
It is important to note that Lemma \ref{lem:: Equivalence conditional and unconditional}
does not depend on the assumption that $f\left(\boldsymbol{\theta},\boldsymbol{\mu}\right)$
is the correct log-likelihood function of $M_{i}$, or equivalently that
$\Pr\left(M_{i}\mid\boldsymbol{V}_{i}\right)=\mathrm{Po}\left(\sum_{k=1}^{d}e^{\left(\eta_{ik}+\mu_{i}\right)}\right)$
or $\Pr\left(C_{ik}\mid\boldsymbol{V}_{i}\right)=\mathrm{Po}\left(e^{\left(\eta_{ik}+\mu_{i}\right)}\right)$,
where $\mathrm{Po}\left(\cdot\right)$ denotes the Poisson distribution.
No assumption about the conditional distribution of $M_{i}$ given $\boldsymbol{V}_{i}$
is required for any of the results in the paper to hold. As we show
in the following sections, the introduction of $f\left(\boldsymbol{\theta},\boldsymbol{\mu}\right)$
is merely a trick to achieve distributed computing.
\section{Iterative Distributed Computing Estimator \label{sec:Iterative-Distributed-Computing-}}
Lemma \ref{lem:: Equivalence conditional and unconditional} shows
that, instead of minimizing $Q_{n}^{\ast}\left(\boldsymbol{\theta}\right)$,
we can minimize $Q_{n}\left(\boldsymbol{\theta},\boldsymbol{\mu}\right)$
to obtain the QMLE of $\boldsymbol{\theta}^{\ast}$. However,
computing $\widehat{\boldsymbol{\theta}}$ remains computationally intensive, as solving (\ref{eq:Definition of Distributed Multinomial Estimator})
is impractical when $d$ is large. Nevertheless, the additive structure
of $Q_{n}\left(\boldsymbol{\theta},\boldsymbol{\mu}\right)$ enables the problem to be solved distributively.
\subsection{Distributed Computing Estimator in \citet{taddy2015distributed}}
As noted in \citet{taddy2015distributed}, although minimizing $Q_{n}\left(\boldsymbol{\theta},\boldsymbol{\mu}\right)$
with respect to $\boldsymbol{\theta}$ and $\boldsymbol{\mu}$ jointly
is computationally infeasible, given any value of $\boldsymbol{\mu}$,
solving $\arg\min_{\boldsymbol{\theta}}Q_{n}\left(\boldsymbol{\theta},\boldsymbol{\mu}\right)$
is much easier because the optimization can be done separately for
each $\boldsymbol{\theta}_{k}$ and be computed across machines. To
see this, we rewrite $Q_{n}\left(\boldsymbol{\theta},\boldsymbol{\mu}\right)$
as:
\begin{align}
Q_{n}\left(\boldsymbol{\theta},\boldsymbol{\mu}\right) & =\sum_{i=1}^{n}\sum_{k=1}^{d}\left(e^{\left(\eta_{ik}+\mu_{i}\right)}-C_{ik}\left(\eta_{ik}+\mu_{i}\right)\right).\nonumber \\
& =\sum_{k=1}^{d}\sum_{i=1}^{n}\left(e^{\left(\boldsymbol{V}_{i}^{\prime}\boldsymbol{\theta}_{k}+\mu_{i}\right)}-C_{ik}\left(V_{i}^{\prime}\boldsymbol{\theta}_{k}+\mu_{i}\right)\right)\nonumber \\
& \equiv Q_{1n}\left(\boldsymbol{\theta}_{1},\boldsymbol{\mu}\right)+\cdots+Q_{dn}\left(\boldsymbol{\theta}_{d},\boldsymbol{\mu}\right),\label{eq:Qn equals Q1n until Qdn}
\end{align}
where for $k=1,\ldots,d$, $Q_{kn}\left(\boldsymbol{\theta}_{k},\boldsymbol{\mu}\right)\equiv\sum_{i=1}^{n}\left(e^{\left(V_{i}^{\prime}\boldsymbol{\theta}_{k}+\mu_{i}\right)}-C_{ik}\left(V_{i}^{\prime}\boldsymbol{\theta}_{k}+\mu_{i}\right)\right)$.
In consequence, it holds that for any $\boldsymbol{\mu}$,
\begin{equation}
\arg\min_{\boldsymbol{\theta}\in\Theta}Q_{n}\left(\boldsymbol{\theta},\boldsymbol{\mu}\right)=\left[\arg\min_{\boldsymbol{\theta}_{1}\in\varTheta_{1}}Q_{1n}\left(\boldsymbol{\theta}_{1},\boldsymbol{\mu}\right),\ldots,\arg\min_{\boldsymbol{\theta}_{d}\in\varTheta_{d}}Q_{dn}\left(\boldsymbol{\theta}_{d},\boldsymbol{\mu}\right)\right]^{\prime}.\label{eq:Distributed Computing}
\end{equation}
Based on (\ref{eq:Distributed Computing}), solving $\arg\min_{\boldsymbol{\theta}\in\Theta}Q_{n}\left(\boldsymbol{\theta},\boldsymbol{\mu}\right)$
for any given $\boldsymbol{\mu}$ is equivalent to solving $d$ optimizations:
$\arg\min_{\boldsymbol{\theta}_{k}\in\varTheta_{k}}Q_{kn}\left(\boldsymbol{\theta}_{k},\boldsymbol{\mu}\right)$
for each $k=1,\ldots,d$, where each optimization is a Poisson regression.\footnote{Since $\varTheta_{d}=\left\{ \mathbf{0}\right\} $, solving for $\arg\min_{\boldsymbol{\theta}_{d}\in\varTheta_{d}}Q_{dn}\left(\boldsymbol{\theta}_{d},\boldsymbol{\mu}\right)$
is trivial.} Since $\boldsymbol{\theta}_{k}$ has only $p$ dimensions, $\arg\min_{\boldsymbol{\theta}_{k}\in\varTheta_{k}}Q_{kn}\left(\boldsymbol{\theta}_{k},\boldsymbol{\mu}\right)$
is easy to compute. In addition, the optimizations for $k=1,\ldots,d$
can be computed across machines allowing for distributed computing.
By Equation (\ref{eq:Definition of Distributed Multinomial Estimator}),
it is not difficult to see that $\widehat{\boldsymbol{\theta}}=\arg\min_{\boldsymbol{\theta}\in\Theta}Q_{n}\left(\boldsymbol{\theta},\widehat{\boldsymbol{\mu}}\right)$.
Because $\widehat{\boldsymbol{\theta}}$ is equivalent to the MLE
$\widetilde{\boldsymbol{\theta}}$ by Lemma \ref{lem:: Equivalence conditional and unconditional},
it has the desired properties such as being both consistent and asymptotically
efficient. As a result, we would hope to obtain $\widehat{\boldsymbol{\mu}}$
first and then compute $\widehat{\boldsymbol{\theta}}$ by distributed
computing. However, the value of $\widehat{\boldsymbol{\mu}}$ depends
on $\widehat{\boldsymbol{\theta}}$, which itself is difficult to
calculate. On the other hand, given any value of $\boldsymbol{\theta}$,
solving $\arg\min_{\boldsymbol{\mu}}Q_{n}\left(\boldsymbol{\theta},\boldsymbol{\mu}\right)$
is also fast, and the solution even has a closed form. Denote the
solution to $\arg\min_{\boldsymbol{\mu}}Q_{n}\left(\boldsymbol{\theta},\boldsymbol{\mu}\right)$
as $\overline{\boldsymbol{\mu}}_{n}\left(\boldsymbol{\theta}\right)$.
A simple calculation would show that
\begin{equation}
\overline{\boldsymbol{\mu}}_{n}\left(\boldsymbol{\theta}\right)=\left(\log\left(\frac{M_{1}}{\sum_{k=1}^{d}e^{\eta_{1k}}}\right),\ldots,\log\left(\frac{M_{n}}{\sum_{k=1}^{d}e^{\eta_{nk}}}\right)\right)^{\prime}.\label{eq:mu bar}
\end{equation}
Let $\widehat{\boldsymbol{\mu}}_{T}=\left(\log\left(M_{1}\right),\ldots,\log\left(M_{n}\right)\right)^{\prime}$.
Instead of solving for $\widehat{\boldsymbol{\mu}}$ using (\ref{eq:mu bar}),
\citet{taddy2015distributed} proposes an estimator $\widehat{\boldsymbol{\theta}}_{T}=\arg\min_{\boldsymbol{\theta}\in\Theta}Q_{n}\left(\boldsymbol{\theta},\widehat{\boldsymbol{\mu}}_{T}\right)$
and calls it the distributed computing estimator. Such an estimator
is fast to compute but fails to be consistent except in the special
cases discussed in \citet{taddy2015distributed}.
\subsection{Iterative Distributed Computing Estimator \label{sec_sub:Iterative-Distributed-Computing} }
We propose an iterative distributed computing (IDC) estimator, such
that during each iteration we solve (\ref{eq:Distributed Computing})
with $\boldsymbol{\mu}$ updated from the previous step estimate of
$\boldsymbol{\theta}$ via (\ref{eq:mu bar}). Our IDC estimator is
defined by the following steps.
\textbf{Step $\mathbf{0}$. }Compute an initial estimator of $\boldsymbol{\theta}^{\ast}$,
denoted as $\widehat{\boldsymbol{\theta}}^{\left(0\right)}$.
\textbf{Step $\mathbf{1},\ldots,\mathbf{S}$. }For step $s$, where
$s=1,\ldots,S$, we first update $\boldsymbol{\mu}$ using estimator
$\widehat{\boldsymbol{\theta}}^{\left(s-1\right)}$ from the previous
step via $\overline{\boldsymbol{\mu}}_{n}\left(\cdot\right)$. Then
we update $\boldsymbol{\theta}$ given the value of $\boldsymbol{\mu}$:
\begin{align}
& \widehat{\boldsymbol{\theta}}^{\left(s\right)}=\overline{\boldsymbol{\theta}}_{n}\left(\widehat{\boldsymbol{\theta}}^{\left(s-1\right)}\right)\nonumber \\
\equiv & \left[\arg\min_{\boldsymbol{\theta}_{1}\in\varTheta_{1}}Q_{1n}\left(\boldsymbol{\theta}_{1},\overline{\boldsymbol{\mu}}_{n}\left(\widehat{\boldsymbol{\theta}}^{\left(s-1\right)}\right)\right),\ldots,\arg\min_{\boldsymbol{\theta}_{d}\in\varTheta_{d}}Q_{dn}\left(\boldsymbol{\theta}_{d},\overline{\boldsymbol{\mu}}_{n}\left(\widehat{\boldsymbol{\theta}}^{\left(s-1\right)}\right)\right)\right]^{\prime}.\label{eq::iterative procedure}
\end{align}
The iterative estimator with $S$ iterations is defined as $\widehat{\boldsymbol{\theta}}^{I}\equiv\widehat{\boldsymbol{\theta}}^{\left(S\right)}$.
For any $\boldsymbol{\theta}$, the value of $\overline{\boldsymbol{\mu}}_{n}\left(\boldsymbol{\theta}\right)$
can be directly computed from (\ref{eq:mu bar}). In each step, we
compute $\arg\min_{\boldsymbol{\theta}_{k}}Q_{kn}\left(\boldsymbol{\theta}_{k},\overline{\boldsymbol{\mu}}_{n}\left(\widehat{\boldsymbol{\theta}}^{\left(s-1\right)}\right)\right)$
for $k=1,\ldots,d$ on $d$ parallel computers. This amounts to running
$d$ Poisson regressions with $p$ parameters, and the computational
burden for each step is low. The algorithm is described in Algorithm
\ref{alg::IDC algorithm}.
\begin{algorithm}[!h]
\DontPrintSemicolon
\KwInput{$S$} \KwOutput{$\widehat{\boldsymbol{\theta}}^{\left(S\right)}$}
Compute an initial estimator $\widehat{\boldsymbol{\theta}}^{\left(0\right)}$
\tcc{Start of Step 1}
Compute $\overline{\boldsymbol{\mu}}_{n}\left(\widehat{\boldsymbol{\theta}}^{\left(0\right)}\right)$
Solve for
\[
\widehat{\boldsymbol{\theta}}^{\left(1\right)}=\left[\arg\min_{\boldsymbol{\theta}_{1}\in\varTheta_{1}}Q_{1n}\left(\boldsymbol{\theta}_{1},\overline{\boldsymbol{\mu}}_{n}\left(\widehat{\boldsymbol{\theta}}^{\left(0\right)}\right)\right),\ldots,\arg\min_{\boldsymbol{\theta}_{d}\in\varTheta_{d}}Q_{dn}\left(\boldsymbol{\theta}_{d},\overline{\boldsymbol{\mu}}_{n}\left(\widehat{\boldsymbol{\theta}}^{\left(0\right)}\right)\right)\right]^{\prime}
\]
\tcc{End of Step 1. The output is $\widehat{\boldsymbol{\theta}}^{\left(1\right)}$}
\tcc{Start of Step 2}
Compute $\overline{\boldsymbol{\mu}}_{n}\left(\widehat{\boldsymbol{\theta}}^{\left(1\right)}\right)$
Solve for
\[
\widehat{\boldsymbol{\theta}}^{\left(2\right)}=\left[\arg\min_{\boldsymbol{\theta}_{1}\in\varTheta_{1}}Q_{1n}\left(\boldsymbol{\theta}_{1},\overline{\boldsymbol{\mu}}_{n}\left(\widehat{\boldsymbol{\theta}}^{\left(1\right)}\right)\right),\ldots,\arg\min_{\boldsymbol{\theta}_{d}\in\varTheta_{d}}Q_{dn}\left(\boldsymbol{\theta}_{d},\overline{\boldsymbol{\mu}}_{n}\left(\widehat{\boldsymbol{\theta}}^{\left(1\right)}\right)\right)\right]^{\prime}
\]
\tcc{End of Step 2. The output is $\widehat{\boldsymbol{\theta}}^{\left(2\right)}$}
...... Continue until Step $S$. The output of Step $S$ is $\widehat{\boldsymbol{\theta}}^{\left(S\right)}$
\caption{the iterative distributed computing procedure}
\label{alg::IDC algorithm}
\end{algorithm}
Unlike many existing algorithms that numerically solve for the MLE, such as gradient descent or stochastic
gradient descent, the IDC estimator does not inherently involve any tuning parameter.
This is advantageous because the performance of the classical (stochastic)
gradient descent is generally sensitive to the learning rate.
\subsection{Initial Estimators \label{subsec:Initial-Estimator}}
Similar to all iterative optimization procedures, the initial estimator
plays a critical role. In finite samples, a good initial guess of
$\boldsymbol{\theta}^{\ast}$ can improve the performance of the IDC
estimator. Asymptotically, a consistent initial estimator can lead
to consistent and asymptotically efficient iterative estimators under
weaker assumptions than an inconsistent initial estimator. In this
section, we propose three initial estimators: a consistent initial
estimator of $\boldsymbol{\theta}^{\ast}$ based on binomial MLE,
\citet{taddy2015distributed}'s estimator, and the MLE based on the
Poisson assumption of $M_{i}$. The latter two are inconsistent without
any assumption on the distribution of $M_{i}$.
\paragraph{A Consistent Initial Estimator }
Let $N_{ik}\equiv C_{ik}+C_{id}$. The following lemma results from
the MNL model defined in (\ref{eq: MNL model}).
\begin{lemma} \label{lem:: Ck Cd distribution} For any $k=1,\ldots,d-1$,
\begin{equation}
\Pr\left(C_{ik},C_{id}\mid\boldsymbol{V}_{i},N_{ik}\right)=\frac{N_{ik}!}{C_{ik}!C_{id}!}\left(\frac{e^{\eta_{ik}^{\ast}}}{e^{\eta_{ik}^{\ast}}+1}\right)^{C_{ik}}\left(\frac{1}{e^{\eta_{ik}^{\ast}}+1}\right)^{C_{id}}.\label{eq: Independent MNL}
\end{equation}
\end{lemma}
Lemma \ref{lem:: Ck Cd distribution} shows that we can consistently
estimate $\boldsymbol{\theta}_{k}^{\ast}$ based on a binomial logistic
regression with the log-likelihood function given by:
\begin{align*}
l_{C_{k},C_{d}\mid V,N_{k}}\left(\boldsymbol{\theta}_{k}\right) & \equiv\sum_{i=1}^{n}\log\Pr\left(C_{ik},C_{id}\mid\boldsymbol{V}_{i},N_{ik}\right)\\
& \sim\sum_{i=1}^{n}\left[C_{ik}\eta_{ik}-\left(C_{ik}+C_{id}\right)\log\left(e^{\eta_{ik}}+1\right)\right]\\
& =\sum_{i=1}^{n}\left[C_{ik}\boldsymbol{V}_{i}^{\prime}\boldsymbol{\theta}_{k}-\left(C_{ik}+C_{id}\right)\log\left(e^{\boldsymbol{V}_{i}^{\prime}\boldsymbol{\theta}_{k}}+1\right)\right].
\end{align*}
Let $\check{\boldsymbol{\theta}}_{k}=\arg\min_{\boldsymbol{\theta}_{k}\in\varTheta_{k}}-l_{C_{k},C_{d}\mid V,N_{k}}\left(\boldsymbol{\theta}_{k}\right)$
for $k=1,\ldots,d-1$ and $\check{\boldsymbol{\theta}}=\left(\boldsymbol{\check{\theta}}_{1}^{\prime},\ldots,\boldsymbol{\check{\theta}}_{d-1}^{\prime},\boldsymbol{\check{\theta}}_{d}^{\prime}\right)^{\prime}$
with $\boldsymbol{\check{\theta}}_{d}^{\prime}=\mathbf{0}$. The consistency
of $\check{\boldsymbol{\theta}}$ follows from standard arguments
in the maximum likelihood estimation.
Compared to $\widetilde{\boldsymbol{\theta}}$, the conditional probability
used in constructing the above binomial logistic log-likelihood function
does not use all the available information. Therefore, $\check{\boldsymbol{\theta}}$
is less efficient than $\widetilde{\boldsymbol{\theta}}$. However,
each component of $\check{\boldsymbol{\theta}}$, $\check{\boldsymbol{\theta}}_{k}$,
can be calculated independently, allowing for parallel computing.
The substantially short running time of $\check{\boldsymbol{\theta}}$
makes it a great candidate for the initial value $\widehat{\boldsymbol{\theta}}^{\left(0\right)}$.
\paragraph{Inconsistent Initial Estimators }
Even though \citet{taddy2015distributed}'s estimator fails to be
consistent in general, it could serve as a candidate for the initial
value in our Algorithm \ref{alg::IDC algorithm}. Another option is
to replace $\widehat{\boldsymbol{\mu}}_{T}$ with a zero vector to
obtain another estimator denoted as $\widehat{\boldsymbol{\theta}}_{P}=\arg\min_{\boldsymbol{\theta}\in\Theta}Q_{n}\left(\boldsymbol{\theta},\mathbf{0}\right)$.
It can also be computed across machines for each $k=1,...,d$. Like
\citet{taddy2015distributed}'s estimator, $\widehat{\boldsymbol{\theta}}_{P}$
could also serve as a candidate for the initial value in our algorithm.
Moreover, under an extra condition that $M_{i}$ follows a Poisson
distribution, $\widehat{\boldsymbol{\theta}}_{P}$ is the maximum
likelihood estimator of $\boldsymbol{\theta}^{\ast}$.
\begin{lemma} \label{lem:: consistency under stronger assumption}
If $\Pr\left(M_{i}\mid V_{i}\right)=\mathrm{Po}\left(\sum_{k=1}^{d}e^{\eta_{ik}^{\ast}}\right)$,
then $\widehat{\boldsymbol{\theta}}_{P}$ is the maximum likelihood
estimator of $\boldsymbol{\theta}^{\ast}$ based on the conditional
probability $\Pr\left(\boldsymbol{C}_{i}\mid V_{i}\right)$. \end{lemma}
Unlike $\check{\boldsymbol{\theta}}$, neither $\widehat{\boldsymbol{\theta}}_{T}$
nor $\widehat{\boldsymbol{\theta}}_{P}$ is a consistent estimator
of $\boldsymbol{\theta}^{\ast}$ without any additional assumption.
\subsection{Constrained Iterative Distributed Computing \label{sec:Extension}}
In some applications, researchers may have prior knowledge of some
linear equality constraints among parameters. Taking the constraints into consideration during the estimation would further improve asymptotic efficiency. In this section, we discuss how to modify our IDC estimator
introduced in Section \ref{sec_sub:Iterative-Distributed-Computing}
to incorporate equality constraints. The initial estimators introduced
in Section \ref{subsec:Initial-Estimator}, although they do not account
for the equality constraints, can be utilized to obtain the initial
$\widehat{\boldsymbol{\theta}}^{\left(0\right)}$.
We consider two different types of constraints: constraints on parameters
for the same choice and different choices. The procedures for
different types of constraints differ in the optimization problems
during each iteration. For each type, we use an example to illustrate
our procedure.
For the first type, the constraint is on components of individual
$\boldsymbol{\theta}_{k}^{\ast}$. Take the constraint $\theta_{k1}^{\ast}=\theta_{k2}^{\ast}$
for $k=1,\ldots d$ as an example. When computing $\widehat{\boldsymbol{\theta}}^{\left(s\right)}$
in Step $s$, we solve the constrained optimization problem:
\[
\arg\min_{\boldsymbol{\theta}_{k}\in\varTheta_{k},\theta_{k1}=\theta_{k2}}Q_{kn}\left(\boldsymbol{\theta}_{k},\overline{\boldsymbol{\mu}}_{n}\left(\widehat{\boldsymbol{\theta}}^{\left(s-1\right)}\right)\right)
\]
for each $k$. Because the constraint is on each $\boldsymbol{\theta}_{k}^{\ast}$,
the original distributed computing scheme remains.
The second type of constraints involves components of $\boldsymbol{\theta}_{k}^{\ast}$
across different choices. For example, researchers may impose restrictions
like $\theta_{11}^{\ast}=\cdots=\theta_{q1}^{\ast}$, where $q<d$.
For $k=1,\ldots,q$, let $\boldsymbol{\theta}_{k,-1}$ be the subvector
of $\boldsymbol{\theta}_{k}$ that excludes its first element. From
Steps $1$ to $S$, we first update $\boldsymbol{\mu}$ using estimator
$\widehat{\boldsymbol{\theta}}^{\left(s-1\right)}$ from the previous
step. Then, we compute $\widehat{\boldsymbol{\theta}}^{\left(s\right)}$
from the following optimization problems:
\begin{align}
& \left(\widehat{\boldsymbol{\theta}}_{1,-1}^{\left(s\right)},\ldots,\widehat{\boldsymbol{\theta}}_{q,-1}^{\left(s\right)},\widehat{\boldsymbol{\theta}}_{q+1}^{\left(s\right)},\ldots,\widehat{\boldsymbol{\theta}}_{d}^{\left(s\right)}\right)\nonumber \\
= & \biggl[\arg\min_{\boldsymbol{\theta}_{1,-1}}Q_{1n}\left(\widehat{\theta}_{11}^{\left(s-1\right)},\boldsymbol{\theta}_{1,-1},\overline{\boldsymbol{\mu}}_{n}\left(\widehat{\boldsymbol{\theta}}^{\left(s-1\right)}\right)\right),\ldots,\arg\min_{\boldsymbol{\theta}_{q,-1}}Q_{qn}\left(\widehat{\theta}_{q1}^{\left(s-1\right)},\boldsymbol{\theta}_{q,-1},\overline{\boldsymbol{\mu}}_{n}\left(\widehat{\boldsymbol{\theta}}^{\left(s-1\right)}\right)\right)\nonumber \\
& \arg\min_{\boldsymbol{\theta}_{q+1}}Q_{\left(q+1\right)n}\left(\boldsymbol{\theta}_{q+1},\overline{\boldsymbol{\mu}}_{n}\left(\widehat{\boldsymbol{\theta}}^{\left(s-1\right)}\right)\right),\ldots,\arg\min_{\boldsymbol{\theta}_{d}}Q_{dn}\left(\boldsymbol{\theta}_{d},\overline{\boldsymbol{\mu}}_{n}\left(\widehat{\boldsymbol{\theta}}^{\left(s-1\right)}\right)\right)\biggr],\label{eq: theta -1}\\
\widehat{\boldsymbol{\theta}}_{11}^{\left(s\right)} & =\arg\min_{\theta_{11}}\left[Q_{1n}\left(\theta_{11},\widehat{\boldsymbol{\theta}}_{1,-1}^{\left(s\right)},\overline{\boldsymbol{\mu}}_{n}\left(\widehat{\boldsymbol{\theta}}^{\left(s-1\right)}\right)\right)+\cdots+Q_{qn}\left(\theta_{11},\widehat{\boldsymbol{\theta}}_{q,-1}^{\left(s\right)},\overline{\boldsymbol{\mu}}_{n}\left(\widehat{\boldsymbol{\theta}}^{\left(s-1\right)}\right)\right)\right],\label{eq:theta 11}\\
\mbox{ and } & \left(\widehat{\boldsymbol{\theta}}_{21}^{\left(s\right)},\ldots,\widehat{\boldsymbol{\theta}}_{q1}^{\left(s\right)}\right)=\left(\widehat{\boldsymbol{\theta}}_{11}^{\left(s\right)},\ldots,\widehat{\boldsymbol{\theta}}_{11}^{\left(s\right)}\right).\nonumber
\end{align}
Optimization problems (\ref{eq: theta -1}) can be solved using parallel
computers. And (\ref{eq:theta 11}) is an optimization problem with
only one argument.\footnote{Let a general linear equality constraint of the second type be written
as $R\boldsymbol{\theta}^{\ast}=r$, where $R$ and $r$ are known
with dimensions $l_{R}\times pd$ and $l_{R}\times1$ respectively.
The matrix $R$ is assumed to have full row rank so that there is
no redundant constraint. We can always rearrange and decompose $R$
as $\left[R_{c},\mathbf{0}\right]$, where $R_{c}$ has dimension
$l_{R}\times q$ and has no zero column. The number of arguments in
the optimization problem (\ref{eq:theta 11}) for $R\boldsymbol{\theta}^{\ast}=r$
is $q-l_{R}\geq0$. The case where $q=l_{R}$ corresponds all $q$
number of elements in $\boldsymbol{\theta}^{\ast}$ having prespecified
values. The optimization problem (\ref{eq:theta 11}) is no longer
needed in this case.} In consequence, each step incurs a low computational burden. The
IDC estimator with such a constraint is also fast to compute.
The aforementioned two procedures can be generalized to accommodate
any linear equality constraint. In particular, the two procedures
can be combined in a straightforward way when constraints contain
both types.
\section{Asymptotic Theory\label{sec:Asymptotic-Theory}}
In this section, we establish the consistency and asymptotic normality
of our IDC estimator introduced in Section \ref{sec_sub:Iterative-Distributed-Computing}.\footnote{Asymptotic properties of the constrained estimators discussed in Section
\ref{sec:Extension} can be established in the same way. Moreover,
the procedures can be used to construct a likelihood ratio test for
testing the null hypothesis on linear equality constraint among parameters.} Technical proofs are collected in Appendix \ref{subsec:Proofs}.
We first impose the following two assumptions.
\begin{assumption}\label{Assm:: consistency 0} $\left\{ \left(\boldsymbol{C}_{i},\boldsymbol{V}_{i},M_{i}\right)\right\} _{i=1}^{n}$
are random samples of $\left(\boldsymbol{C},\boldsymbol{V},M\right)$.
\end{assumption}
\begin{assumption}\label{Assm:: consistency 1} (i) $\Theta$ is
compact and convex. (ii) $\mathbb{E}\left[e^{\boldsymbol{V}^{\prime}\boldsymbol{\theta}}\right]<\infty$
for all $\boldsymbol{\theta}\in\Theta$. (iii) $\mathbb{E}\left[M^{2}\right]<\infty$.\end{assumption}
Assumption \ref{Assm:: consistency 1} is standard. Assumption \ref{Assm:: consistency 1}
(ii) requires that the moment generating function of $\boldsymbol{V}$
exists within $\Theta$. Almost all commonly seen distributions satisfy
Assumption \ref{Assm:: consistency 1} (ii) and (iii).
\subsection{Consistency of the IDC Estimator \label{subsec:Consistency }}
In order to analyze the asymptotic properties of the IDC estimator,
we need to study the way in which the iterative estimator updates
itself in each step. In Step $s$ of the iteration, we first compute
$\overline{\boldsymbol{\mu}}_{n}\left(\widehat{\boldsymbol{\theta}}^{\left(s-1\right)}\right)$
using $\widehat{\boldsymbol{\theta}}^{\left(s-1\right)}$ from the
previous step and then calculate $\arg\min_{\boldsymbol{\theta}\in\Theta}Q_{n}\left(\boldsymbol{\theta},\overline{\boldsymbol{\mu}}_{n}\left(\widehat{\boldsymbol{\theta}}^{\left(s-1\right)}\right)\right)$.
To explicitly distinguish the argument in $\overline{\boldsymbol{\mu}}_{n}\left(\cdot\right)$
from the argument in $Q_{n}\left(\cdot,\boldsymbol{\mu}\right)$ for
any given $\boldsymbol{\mu}$, we introduce $\boldsymbol{\vartheta}$
and use it as the argument in $\overline{\boldsymbol{\mu}}_{n}\left(\cdot\right)$.
Define $Q_{n}^{\dagger}\left(\boldsymbol{\theta},\boldsymbol{\vartheta}\right)\equiv Q_{n}\left(\boldsymbol{\theta},\overline{\boldsymbol{\mu}}_{n}\left(\boldsymbol{\vartheta}\right)\right)$.
We have that
\[
Q_{n}^{\dagger}\left(\boldsymbol{\theta},\boldsymbol{\vartheta}\right)\equiv\sum_{i=1}^{n}\sum_{k=1}^{d}\left(\frac{M_{i}e^{\boldsymbol{V}_{i}^{\prime}\boldsymbol{\theta}_{k}}}{\sum_{k=1}^{d}e^{\boldsymbol{V}_{i}^{\prime}\boldsymbol{\vartheta}_{k}}}-C_{ik}\boldsymbol{V}_{i}^{\prime}\boldsymbol{\theta}_{k}-C_{ik}\log\left(\frac{M_{i}}{\sum_{k=1}^{d}e^{\boldsymbol{V}_{i}^{\prime}\boldsymbol{\vartheta}_{k}}}\right)\right).
\]
Further define function $Q^{\dagger}\left(\boldsymbol{\theta},\boldsymbol{\vartheta}\right)$
as the probability limit of $\frac{1}{n}Q_{n}^{\dagger}\left(\boldsymbol{\theta},\boldsymbol{\vartheta}\right)$:
\begin{align}
Q^{\dagger}\left(\boldsymbol{\theta},\boldsymbol{\vartheta}\right) & \equiv\sum_{k=1}^{d}\mathbb{E}\left[\left(\frac{Me^{\boldsymbol{V}^{\prime}\boldsymbol{\theta}_{k}}}{\sum_{k=1}^{d}e^{\boldsymbol{V}^{\prime}\boldsymbol{\vartheta}_{k}}}-C_{k}\boldsymbol{V}^{\prime}\boldsymbol{\theta}_{k}-C_{k}\log\left(\frac{M}{\sum_{k=1}^{d}e^{\boldsymbol{V}^{\prime}\boldsymbol{\vartheta}_{k}}}\right)\right)\right].\label{eq:Q infty Q_k}
\end{align}
In the following lemma, we provide some properties of $Q^{\dagger}\left(\boldsymbol{\theta},\boldsymbol{\vartheta}\right)$,
which are crucial for obtaining the consistency and asymptotic normality
of our IDC estimator.
\begin{lemma} \label{lem:: properties of Q star} Under Assumptions
\ref{Assm:: consistency 0} and \ref{Assm:: consistency 1}, the following
results hold.
(i) $\sup_{\boldsymbol{\theta},\boldsymbol{\vartheta}\in\Theta}\left|Q_{n}^{\dagger}\left(\boldsymbol{\theta},\boldsymbol{\vartheta}\right)-Q^{\dagger}\left(\boldsymbol{\theta},\boldsymbol{\vartheta}\right)\right|\stackrel{p}{\rightarrow}0$.
(ii) For any given $\boldsymbol{\vartheta}$, $Q^{\dagger}\left(\boldsymbol{\theta},\boldsymbol{\vartheta}\right)$
has a unique minimizer denoted as $\overline{\boldsymbol{\theta}}\left(\boldsymbol{\vartheta}\right)$.
(iii) $\overline{\boldsymbol{\theta}}\left(\cdot\right)$ is continuous
on $\Theta$.
(iv) $\overline{\boldsymbol{\theta}}\left(\boldsymbol{\theta}^{\ast}\right)=\boldsymbol{\theta}^{\ast}$,
i.e., the true value $\boldsymbol{\theta}^{\ast}$ is a fixed point
of the mapping $\overline{\boldsymbol{\theta}}:\Theta\rightarrow\Theta$.
(v) $\boldsymbol{\theta}^{\ast}$ is the unique fixed point of $\overline{\boldsymbol{\theta}}\left(\cdot\right)$.
\end{lemma}
Essentially, $\overline{\boldsymbol{\theta}}\left(\boldsymbol{\vartheta}\right)$
summarizes the operation in each step with $\boldsymbol{\vartheta}$
being the input and $\overline{\boldsymbol{\theta}}\left(\boldsymbol{\vartheta}\right)$
being the output when the sample size goes to infinity. By part (ii)
of Lemma \ref{lem:: properties of Q star}, $\overline{\boldsymbol{\theta}}\left(\cdot\right)$
is well-defined. Part (v) of Lemma \ref{lem:: properties of Q star}
plays the most important role. Heuristically, for any given $\boldsymbol{\vartheta}$,
the value of function $\overline{\boldsymbol{\theta}}\left(\boldsymbol{\vartheta}\right)$
is obtained by solving $\frac{\partial Q^{\dagger}\left(\boldsymbol{\theta},\boldsymbol{\vartheta}\right)}{\partial\boldsymbol{\theta}}=\mathbf{0}$.
At the same time, function $\frac{\partial Q^{\dagger}\left(\boldsymbol{\theta},\boldsymbol{\vartheta}\right)}{\partial\boldsymbol{\theta}}$
relates to the first order derivative of $L_{C\mid V,M}\left(\boldsymbol{\theta}\right)$,
the population objective function defined in Section \ref{subsec:Maximum-Likelihood-Estimator}.
By the identification assumption and the convexity of $-L_{C\mid V,M}\left(\boldsymbol{\theta}\right)$,
only the true value $\boldsymbol{\theta}^{\ast}$ satisfies that $\frac{\partial}{\partial\boldsymbol{\theta}}L_{C\mid V,M}\left(\boldsymbol{\theta}\right)=\mathbf{0}$,
which implies that $\frac{\partial Q^{\dagger}\left(\boldsymbol{\theta},\boldsymbol{\vartheta}\right)}{\partial\boldsymbol{\theta}}\mid_{\boldsymbol{\vartheta}=\boldsymbol{\theta}}=\mathbf{0}$
holds only at $\left(\boldsymbol{\theta},\boldsymbol{\vartheta}\right)=\left(\boldsymbol{\theta}^{\ast},\boldsymbol{\theta}^{\ast}\right)$.
If a consistent initial estimator is used, such as $\check{\boldsymbol{\theta}}$,
then Lemma \ref{lem:: properties of Q star} is sufficient for the
consistency of the IDC estimator as stated below.
\begin{theorem}[Consistent initial value] \label{thm:: consistency-fixed S}
Suppose Assumptions \ref{Assm:: consistency 0} and \ref{Assm:: consistency 1}
hold. If $\widehat{\boldsymbol{\theta}}^{\left(0\right)}\stackrel{p}{\rightarrow}\boldsymbol{\theta}^{\ast}$
as $n\rightarrow\infty$, then $\widehat{\boldsymbol{\theta}}^{\left(S\right)}\stackrel{p}{\rightarrow}\boldsymbol{\theta}^{\ast}$
as $n\rightarrow\infty$ for any $S$. \end{theorem}
On the other hand, if the initial estimator is consistent only under
extra assumptions, such as $\widehat{\boldsymbol{\theta}}_{T}$ and
$\widehat{\boldsymbol{\theta}}_{P}$, or even inconsistent, then we
need a contraction mapping assumption on $\overline{\boldsymbol{\theta}}\left(\cdot\right)$.
\begin{assumption}[Contraction Mapping] \label{ass:: contraction mapping assumption}
For any $\boldsymbol{\vartheta}\in\Theta$, there exists a constant
$C<1$ such that
\[
\left\Vert \overline{\boldsymbol{\theta}}\left(\overline{\boldsymbol{\theta}}\left(\boldsymbol{\vartheta}\right)\right)-\overline{\boldsymbol{\theta}}\left(\boldsymbol{\vartheta}\right)\right\Vert \leq C\left\Vert \overline{\boldsymbol{\theta}}\left(\boldsymbol{\vartheta}\right)-\boldsymbol{\vartheta}\right\Vert .
\]
\end{assumption}
Admittedly, Assumption \ref{ass:: contraction mapping assumption}
is a high-level assumption. Based on the evidence from the simulation,
the assumption holds for various values of $\boldsymbol{\theta}^{\ast}$
and distributions of $\boldsymbol{V}$ and $M$. Assumption \ref{ass:: contraction mapping assumption}
relates to the contraction mapping assumption (Assumption 6) in \citet{pastorello2003iterative}
but is weaker, the reason being that the true $\boldsymbol{\theta}^{\ast}$
is the unique fixed point, see Lemma \ref{lem:: properties of Q star}
(v). Specifically, Assumption \ref{ass:: contraction mapping assumption}
only requires that the distance between $\overline{\boldsymbol{\theta}}\left(\boldsymbol{\vartheta}\right)$
and $\boldsymbol{\vartheta}$ get smaller after both being mapped
by $\overline{\boldsymbol{\theta}}\left(\cdot\right)$. Instead, the
contraction mapping assumption (Assumption 6) in \citet{pastorello2003iterative}
requires that the distance between two arbitrary $\boldsymbol{\vartheta}^{1}$
and $\boldsymbol{\vartheta}^{2}$ get smaller after both being mapped
by $\overline{\boldsymbol{\theta}}\left(\cdot\right)$.
\begin{theorem}[Inconsistent initial value] \label{thm:: consistency}
Under Assumptions \ref{Assm:: consistency 0}-\ref{ass:: contraction mapping assumption},
$\widehat{\boldsymbol{\theta}}^{\left(S\right)}\stackrel{p}{\rightarrow}\boldsymbol{\theta}^{\ast}$
as $n\rightarrow\infty$ if $S\rightarrow\infty$.\end{theorem}
\subsection{Asymptotic Distributions and Inference \label{subsec:Asymptotic-Distributions-Inference} }
Under Assumptions \ref{Assm:: consistency 0} and \ref{Assm:: consistency 1},
the MLE $\widetilde{\boldsymbol{\theta}}$ is asymptotically normally
distributed with asymptotic variance given by the Fisher information
matrix
\[
\mathcal{I}\left(\boldsymbol{\theta}^{\ast}\right)\equiv\frac{\partial^{2}}{\partial\boldsymbol{\theta}\partial\boldsymbol{\theta}^{\prime}}Q^{\ast}\left(\boldsymbol{\theta}^{\ast}\right),
\]
where $Q^{\ast}\left(\boldsymbol{\theta}^{\ast}\right)\equiv p\lim_{n\rightarrow\infty}\frac{1}{n}Q_{n}^{\ast}\left(\boldsymbol{\theta}^{\ast}\right)$.
In this section, we show that our IDC estimator has the same asymptotic
distribution as $\widetilde{\boldsymbol{\theta}}$, based on which
we introduce a valid bootstrap inference procedure.
The conditions required for proving the asymptotic distribution result
depend on the initial estimator $\widehat{\boldsymbol{\theta}}^{\left(0\right)}$.
If a consistent initial estimator is used, then the following assumption
is sufficient. For any matrix $A$, denote $\left\Vert A\right\Vert $
as its spectral norm.
\begin{assumption}[Information Dominance] \label{ass::Information Dominance}
It holds that
\[
\left\Vert \left(\frac{\partial^{2}Q^{\dagger}\left(\boldsymbol{\theta},\boldsymbol{\vartheta}\right)}{\partial\boldsymbol{\theta}\partial\boldsymbol{\theta}^{\prime}}\mid_{\boldsymbol{\theta}=\boldsymbol{\theta}^{\ast},\boldsymbol{\vartheta}=\boldsymbol{\theta}^{\ast}}\right)^{-1}\frac{\partial^{2}Q^{\dagger}\left(\boldsymbol{\theta},\boldsymbol{\vartheta}\right)}{\partial\boldsymbol{\theta}\partial\boldsymbol{\vartheta}^{\prime}}\mid_{\boldsymbol{\theta}=\boldsymbol{\theta}^{\ast},\boldsymbol{\vartheta}=\boldsymbol{\theta}^{\ast}}\right\Vert <1.
\]
\end{assumption}
The detailed expressions of $\frac{\partial^{2}Q^{\dagger}\left(\boldsymbol{\theta},\boldsymbol{\vartheta}\right)}{\partial\boldsymbol{\theta}\partial\boldsymbol{\theta}^{\prime}}$
and $\frac{\partial^{2}Q^{\dagger}\left(\boldsymbol{\theta},\boldsymbol{\vartheta}\right)}{\partial\boldsymbol{\theta}\partial\boldsymbol{\vartheta}^{\prime}}$
can be found in Appendix \ref{sec:Notations-and-Equalities}. Assumption
\ref{ass::Information Dominance} is often called the information
dominance condition and is tantamount to the local contraction mapping
condition. It is weaker than Assumption \ref{ass:: contraction mapping assumption}.
Because we have an initial consistent estimator of $\boldsymbol{\theta}^{\ast}$,
Assumption \ref{ass::Information Dominance} can be verified. Additionally,
because the matrix $\frac{\partial^{2}Q^{\dagger}\left(\boldsymbol{\theta},\boldsymbol{\vartheta}\right)}{\partial\boldsymbol{\theta}\partial\boldsymbol{\theta}^{\prime}}$
is block diagonal with each block having dimensions $p\times p$,
computing its inverse is feasible.
The following theorem shows that when $S$ is sufficiently large,
the IDC estimator $\widehat{\boldsymbol{\theta}}^{\left(S\right)}$
is equal to $\widetilde{\boldsymbol{\theta}}$ up to a term of order
smaller than $n^{-1/2}$.
\begin{theorem} \label{Thm:: equivalence IDC and theta tilda} (i)
Suppose Assumptions \ref{Assm:: consistency 0}, \ref{Assm:: consistency 1},
and \ref{ass::Information Dominance} hold. If $\widehat{\boldsymbol{\theta}}^{\left(0\right)}\stackrel{p}{\rightarrow}\boldsymbol{\theta}^{\ast}$
as $n\rightarrow\infty$, then $\widehat{\boldsymbol{\theta}}^{\left(S\right)}-\widetilde{\boldsymbol{\theta}}=o_{p}\left(n^{-1/2}\right)$
if $S\geq\log\left(n\right)$. (ii) Under Assumptions \ref{Assm:: consistency 0}-\ref{ass:: contraction mapping assumption},
$\widehat{\boldsymbol{\theta}}^{\left(S\right)}-\widetilde{\boldsymbol{\theta}}=o_{p}\left(n^{-1/2}\right)$
if $S>n^{\delta}$ for some $\delta>0$. \end{theorem}
Theorem \ref{Thm:: equivalence IDC and theta tilda} shows that we
do not lose efficiency when employing the proposed IDC estimator as
long as $S$ is large enough. A direct implication of the theorem
is that $\widehat{\boldsymbol{\theta}}^{\left(S\right)}$ has the
same asymptotic distribution as $\widetilde{\boldsymbol{\theta}}$
for sufficiently large $S$.
\begin{corollary} \label{Cor:: Asymptotics} (i) Suppose Assumptions
\ref{Assm:: consistency 0}, \ref{Assm:: consistency 1}, and \ref{ass::Information Dominance}
hold. If $\widehat{\boldsymbol{\theta}}^{\left(0\right)}\stackrel{p}{\rightarrow}\boldsymbol{\theta}^{\ast}$
as $n\rightarrow\infty$ and $S\geq\log\left(n\right)$, then $\sqrt{n}\left(\widehat{\boldsymbol{\theta}}^{\left(S\right)}-\boldsymbol{\theta}^{\ast}\right)\stackrel{d}{\rightarrow}\mathcal{N}\left(\mathbf{0},\mathcal{I}^{-1}\left(\boldsymbol{\theta}^{\ast}\right)\right)$
as $n\rightarrow\infty$. (ii) Under Assumptions \ref{Assm:: consistency 0}-\ref{ass:: contraction mapping assumption},
if $S>n^{\delta}$ for some $\delta>0$, then $\sqrt{n}\left(\widehat{\boldsymbol{\theta}}^{\left(S\right)}-\boldsymbol{\theta}^{\ast}\right)\stackrel{d}{\rightarrow}\mathcal{N}\left(\mathbf{0},\mathcal{I}^{-1}\left(\boldsymbol{\theta}^{\ast}\right)\right)$
as $n\rightarrow\infty$. \end{corollary}
To conduct inference on $\boldsymbol{\theta}^{\ast}$ based upon $\widehat{\boldsymbol{\theta}}^{\left(S\right)}$,
we need to consistently estimate the Fisher information matrix and
compute its inverse. Because the dimension of $\mathcal{I}\left(\boldsymbol{\theta}^{\ast}\right)$
is $dp\times dp$, calculating the inverse of its estimator is not
only time-consuming but also unreliable when $d$ is large. As a result,
we proceed by applying the following parametric bootstrap, which is
feasible thanks to the fact that the IDC estimator is fast to compute.
Given $\left\{ \left(\boldsymbol{V}_{i},M_{i}\right)\right\} _{i=1}^{n}$
of the original sample, we draw the bootstrap sample $\boldsymbol{C}_{in}^{\ddagger}$
for $i=1,\ldots,n$ from the multinomial logistic regression model
with the conditional probability mass function given by
\[
\mathrm{MNL}\left(\boldsymbol{C}_{in}^{\ddagger};\widehat{\boldsymbol{\eta}}_{i},M_{i}\right)\textrm{, where }\widehat{\eta}_{ik}=\boldsymbol{V}_{i}^{\prime}\widehat{\boldsymbol{\theta}}_{k}^{\left(S\right)}\textrm{ for }k=1,\ldots d.
\]
The bootstrap version of the iterative estimator $\widehat{\boldsymbol{\theta}}^{\ddagger\left(S\right)}$
is obtained by applying the algorithm introduced in Section \ref{sec_sub:Iterative-Distributed-Computing}
with bootstrap sample $\left\{ \left(\boldsymbol{C}_{in}^{\ddagger},\boldsymbol{V}_{i},M_{i}\right)\right\} _{i=1}^{n}$.
Assume that we start with a consistent initial estimator. Based on
Theorem \ref{Thm:: equivalence IDC and theta tilda} (i), we have
that for $S\geq\log\left(n\right)$,
\begin{equation}
\widehat{\boldsymbol{\theta}}^{\left(S\right)}-\boldsymbol{\theta}^{\ast}=\mathcal{I}^{-1}\left(\boldsymbol{\theta}^{\ast}\right)\frac{1}{n}\frac{d}{d\boldsymbol{\theta}}l_{C\mid V,M}\left(\boldsymbol{\theta}^{\ast}\right)+o_{p}\left(n^{-1/2}\right).\label{eq:Theta S influence function}
\end{equation}
Define the score function for $\boldsymbol{\theta}$ as
\[
\dot{\mathbf{l}}\left(\boldsymbol{\theta}\mid\boldsymbol{c},\boldsymbol{v},m\right)\equiv\frac{d}{d\boldsymbol{\theta}}\log\mathrm{MNL}\left(\boldsymbol{c};\boldsymbol{\eta},m\right)\textrm{, where }\eta_{k}\equiv\boldsymbol{v}^{\prime}\boldsymbol{\theta}_{k}.
\]
The iterative estimator $\widehat{\boldsymbol{\theta}}^{\left(S\right)}$
is asymptotically linear with influence function $\widetilde{\mathbf{l}}\left(\boldsymbol{\theta}^{\ast}\mid\boldsymbol{c},\boldsymbol{v},m\right)$:
\[
\sqrt{n}\left(\widehat{\boldsymbol{\theta}}^{\left(S\right)}-\boldsymbol{\theta}^{\ast}\right)=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\widetilde{\mathbf{l}}\left(\boldsymbol{\theta}^{\ast}\mid\boldsymbol{C}_{i},\boldsymbol{V}_{i},M_{i}\right)+o_{p}\left(1\right),
\]
where $\widetilde{\mathbf{l}}\left(\boldsymbol{\theta}^{\ast}\mid\boldsymbol{c},\boldsymbol{v},m\right)\equiv\mathcal{I}^{-1}\left(\boldsymbol{\theta}^{\ast}\right)\dot{\mathbf{l}}\left(\boldsymbol{\theta}^{\ast}\mid\boldsymbol{c},\boldsymbol{v},m\right)$.
Applying the same derivation, we can show that the bootstrap version
of the estimator is also asymptotically linear with the influence
function evaluated at $\widehat{\boldsymbol{\theta}}^{\left(S\right)}$:
\[
\sqrt{n}\left(\widehat{\boldsymbol{\theta}}^{\ddagger\left(S\right)}-\widehat{\boldsymbol{\theta}}^{\left(S\right)}\right)=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\widetilde{\mathbf{l}}\left(\widehat{\boldsymbol{\theta}}^{\left(S\right)}\mid\boldsymbol{C}_{in}^{\ddagger},\boldsymbol{V}_{i},M_{i}\right)+o_{p}\left(1\right).
\]
The Lindeberg-Feller central limit theorem proves the bootstrap consistency.
The proof for the case of an inconsistent initial estimator is analogous.
Let $\stackrel{d^{\ddagger}}{\rightarrow}$ denote the convergence
in bootstrap distribution.
\begin{theorem} \label{thm:: Bootstrap} (i) Suppose Assumptions
\ref{Assm:: consistency 0}, \ref{Assm:: consistency 1}, and \ref{ass::Information Dominance}
hold. If $\widehat{\boldsymbol{\theta}}^{\left(0\right)}\stackrel{p}{\rightarrow}\boldsymbol{\theta}^{\ast}$
as $n\rightarrow\infty$ and $S\geq\log\left(n\right)$, then $\sqrt{n}\left(\widehat{\boldsymbol{\theta}}^{\ddagger\left(S\right)}-\widehat{\boldsymbol{\theta}}^{\left(S\right)}\right)\stackrel{d^{\ddagger}}{\rightarrow}\mathcal{N}\left(\mathbf{0},\mathcal{I}^{-1}\left(\boldsymbol{\theta}^{\ast}\right)\right)$
as $n\rightarrow\infty$. (ii) Under Assumptions \ref{Assm:: consistency 0}-\ref{ass:: contraction mapping assumption},
if $S>n^{\delta}$ for some $\delta>0$, then $\sqrt{n}\left(\widehat{\boldsymbol{\theta}}^{\ddagger\left(S\right)}-\widehat{\boldsymbol{\theta}}^{\left(S\right)}\right)\stackrel{d^{\ddagger}}{\rightarrow}\mathcal{N}\left(\mathbf{0},\mathcal{I}^{-1}\left(\boldsymbol{\theta}^{\ast}\right)\right)$
as $n\rightarrow\infty$. \end{theorem}
\section{Monte Carlo Simulation\label{sec:Monte-Carlo-Simulation}}
In this section, we evaluate the performance of our IDC estimator
from various perspectives. We present the finite sample performance
of the IDC estimator by looking at the effects of separately increasing
$d$ and $n$ on mean squared error (MSE) and running time. We include
the maximum likelihood estimator to show that the IDC estimator performs
similarly to the MLE in terms of MSE and is always feasible even in
cases where MLE is intractable. Lastly, we study the finite sample
size and power of our bootstrap inference procedure.
\subsection{Estimation}
In what follows, we present results on the finite sample performance
of the IDC estimator in terms of MSE and running time in four tables.
The reported running times in the tables are obtained from a cluster
of 25 AWS EC2 instances with 12 vCPUs and 16GB memory. Such a configuration
can be formed on commonly used cloud computing platforms within minutes.
We employ this basic configuration to illustrate that our IDC estimator
achieves superior performance compared to existing estimators, even
when computational resources are suboptimal for distributed computing.
All of the results presented in this section are repeated five hundred
times and averaged.
We first study the finite sample performance of the IDC estimator
with consistent $\check{\boldsymbol{\theta}}$ initialization, specifically
what happens to MSE ($n$ increasing, $d$ fixed) and running time
($d$ increasing, $n$ fixed). We consider the following data generating
process (DGP):
\textbf{DGP-A {[}MNL{]}}: We set $p=5$. The covariate vector $\boldsymbol{V}$
follows the standard normal distribution; $M$ follows the discrete
uniform distribution on $\left[20,30\right]$; and the values of $\mathbf{\boldsymbol{\theta}}^{*}$
are obtained by random draws from the standard normal distribution.
\begin{table}[H]
\centering
\begin{tabular}{ccccccccccccccccc}
\toprule
& & & & & & \multicolumn{3}{c}{$n=500$} & & \multicolumn{3}{c}{$n=1000$} & & \multicolumn{3}{c}{$n=2000$}\tabularnewline
\midrule
$\widehat{\boldsymbol{\theta}}^{\left(0\right)}$ & & $S$ & & $d$ & & MSE & & Time & & MSE & & Time & & MSE & & Time\tabularnewline
\midrule
$\check{\boldsymbol{\theta}}$ & & $10$ & & $10$ & & $.0030$ & & $45$s & & $.0020$ & & $97$s & & $.0012$ & & $252$s\tabularnewline
& & & & $20$ & & $.0083$ & & $52$s & & $.0052$ & & $113$s & & $.0020$ & & $277$s\tabularnewline
& & & & $50$ & & $.0373$ & & $85$s & & $.0182$ & & $174$s & & $.0094$ & & $453$s\tabularnewline
& & & & $100$ & & $.0793$ & & $157$s & & $.0381$ & & $349$s & & $.0171$ & & $756$s\tabularnewline
& & & & $150$ & & $.1749$ & & $211$s & & $.0585$ & & $455$s & & $.0223$ & & $1007$s\tabularnewline
& & & & & & & & & & & & & & & & \tabularnewline
& & $40$ & & $10$ & & $.0037$ & & $168$s & & $.0021$ & & $320$s & & $.0012$ & & $672$s\tabularnewline
& & & & $20$ & & $.0081$ & & $196$s & & $.0040$ & & $352$s & & $.0019$ & & $739$s\tabularnewline
& & & & $50$ & & $.0365$ & & $320$s & & $.0185$ & & $576$s & & $.0091$ & & $1142$s\tabularnewline
& & & & $100$ & & $.0661$ & & $590$s & & $.0336$ & & $890$s & & $.0158$ & & $1948$s\tabularnewline
& & & & $150$ & & $.1815$ & & $770$s & & $.0577$ & & $1184$s & & $.0211$ & & $2822$s\tabularnewline
\bottomrule
\end{tabular}\caption{Finite sample performance of IDC estimator with $\check{\boldsymbol{\theta}}$
initialization}
\label{tab: IDC-Consistent}
\end{table}
From Table \ref{tab: IDC-Consistent}, we observe that the MSE of
the IDC estimator decreases as the sample size increases. When the
number of iterations $S$ increases, we see an improvement in the
MSE. However, the improvement is marginal, suggesting that the IDC
estimator with $\check{\boldsymbol{\theta}}$ initialization stabilizes
with only a few iterations. We also see from Table \ref{tab: IDC-Consistent}
that the running time is approximately a linear function of $d$.
When the number of cores available does not exceed the number of choices,
the additional computational cost of increasing $d$ is very small.
After the cores are fully occupied by the number of processes, the
running time becomes approximately linear. The nominal value of running
time depends on the hardware specifications.
\begin{table}[H]
\centering
\begin{tabular}{ccccccccccccc}
\toprule
{} & & \multicolumn{3}{c}{$n=500$} & & \multicolumn{3}{c}{$n=1000$} & & \multicolumn{3}{c}{$n=2000$}\tabularnewline
\midrule
$d$ & & MSE & & Time & & MSE & & Time & & MSE & & Time\tabularnewline
\midrule
$10$ & & $.0040$ & & $36$s & & $.0020$ & & $105$s & & $.0015$ & & $350$s\tabularnewline
$20$ & & $.0101$ & & $54$s & & $.0045$ & & $152$s & & $.0026$ & & $652$s\tabularnewline
$50$ & & $.0263$ & & $96$s & & $.0187$ & & $308$s & & $.0098$ & & $1523$s\tabularnewline
$100$ & & $.0852$ & & $250$s & & $.0405$ & & $862$s & & $.0151$ & & $4375$s\tabularnewline
$150$ & & $.1179$ & & $457$s & & $.0536$ & & $1523$s & & $.0291$ & & $9352$s\tabularnewline
\bottomrule
\end{tabular}\caption{MSE and running time of MLE $\widetilde{\boldsymbol{\theta}}$}
\label{tab:CMLE}
\end{table}
To compare the performance of our IDC estimator with the MLE $\widetilde{\boldsymbol{\theta}}$,
we simulate the MSE and running time of $\widetilde{\boldsymbol{\theta}}$
from the same DGP and present the result in Table \ref{tab:CMLE}.\footnote{We also write code to try estimators in \citet{bohning1988monotonicity},
\citet{bohning1992multinomial}, and \citet{simon2013blockwise} for
comparison. However, our simulation result suggests that their performance
depends crucially on the number of iterations.} It can be seen that the MSE of the IDC estimator with $\check{\boldsymbol{\theta}}$
initialization is very close to that of the MLE even when the number
of iterations is only $10$. Note that the main advantage of the parallel
estimator is best observed for high enough $d$ because for low $d$,
the communication between parallel processes is unnecessary and hence
parallel computing increases the running time unnecessarily. The superior
performance of the IDC estimator is apparent when $d$ is large. For
instance, when $d=150$ and $n=2000$, the IDC estimator with $S=10$
achieves a similar MSE as the MLE with only about one-tenth of the
running time. For higher-dimensional cases, such as when $d$ exceeds
$150$, computing the MLE becomes computationally intensive and may
not be practical for many applications. In comparison, the IDC estimator
with $S=10$ demonstrates more efficient computation times, requiring
approximately 5, 10, and 20 minutes for sample sizes $n=500$, $1000$,
and $2000$, respectively. Moreover, the running time for the IDC
estimator can be further decreased if more compute instances are used.
For example, using 96 instances, the running time of the IDC estimator
for $d=150$ can be further reduced to $34$, $51$, and $188$ seconds
for $n=500$, $1000$, $2000$ respectively even for $S=40$. Compared
to the corresponding running time of MLE, the running time of the
IDC estimator using 96 instances is more than 10, 30, and 50 times
shorter.
\begin{table}[H]
\centering
\begin{tabular}{ccccccccccccccccc}
\toprule
{} & & & & & & \multicolumn{3}{c}{$n=500$} & & \multicolumn{3}{c}{$n=1000$} & & \multicolumn{3}{c}{$n=2000$}\tabularnewline
\midrule
$\widehat{\boldsymbol{\theta}}^{\left(0\right)}$ & & $S$ & & $d$ & & MSE & & Time & & MSE & & Time & & MSE & & Time\tabularnewline
\midrule
$\widehat{\boldsymbol{\theta}}_{T}$ & & $10$ & & $10$ & & $.0048$ & & $42$s & & $.0023$ & & $95$s & & $.0017$ & & $247$s\tabularnewline
& & & & $20$ & & $.0098$ & & $49$s & & $.0058$ & & $106$s & & $.0038$ & & $278$s\tabularnewline
& & & & $50$ & & $.0466$ & & $83$s & & $.0197$ & & $175$s & & $.0100$ & & $442$s\tabularnewline
& & & & $100$ & & $.0809$ & & $151$s & & $.0381$ & & $339$s & & $.0184$ & & $754$s\tabularnewline
& & & & $150$ & & $.1725$ & & $206$s & & $.0590$ & & $451$s & & $.0242$ & & $998$s\tabularnewline
& & & & & & & & & & & & & & & & \tabularnewline
& & $40$ & & $10$ & & $.0051$ & & $160$s & & $.0025$ & & $311$s & & $.0017$ & & $667$s\tabularnewline
& & & & $20$ & & $.0102$ & & $187$s & & $.0052$ & & $337$s & & $.0022$ & & $728$s\tabularnewline
& & & & $50$ & & $.0464$ & & $309$s & & $.0198$ & & $564$s & & $.0094$ & & $1128$s\tabularnewline
& & & & $100$ & & $.0867$ & & $581$s & & $.0407$ & & $876$s & & $.0179$ & & $1938$s\tabularnewline
& & & & $150$ & & $.1808$ & & $758$s & & $.0588$ & & $1177$s & & $.0219$ & & $2801$s\tabularnewline
& & & & & & & & & & & & & & & & \tabularnewline
& & & & & & \multicolumn{3}{c}{$n=500$} & & \multicolumn{3}{c}{$n=1000$} & & \multicolumn{3}{c}{$n=2000$}\tabularnewline
\midrule
$\widehat{\boldsymbol{\theta}}^{\left(0\right)}$ & & $S$ & & $d$ & & MSE & & Time & & MSE & & Time & & MSE & & Time\tabularnewline
\midrule
\multirow{1}{*}{$\widehat{\boldsymbol{\theta}}_{P}$} & & $10$ & & $10$ & & $.0068$ & & $52$s & & $.0023$ & & $97$s & & $.0016$ & & $244$s\tabularnewline
& & & & $20$ & & $.0118$ & & $59$s & & $.0058$ & & $111$s & & $.0038$ & & $290$s\tabularnewline
& & & & $50$ & & $.0466$ & & $83$s & & $.0197$ & & $177$s & & $.0100$ & & $438$s\tabularnewline
& & & & $100$ & & $.0808$ & & $155$s & & $.0382$ & & $341$s & & $.0184$ & & $751$s\tabularnewline
& & & & $150$ & & $.1739$ & & $209$s & & $.0588$ & & $457$s & & $.0258$ & & $1008$s\tabularnewline
& & & & & & & & & & & & & & & & \tabularnewline
& & $40$ & & $10$ & & $.0050$ & & $160$s & & $.0025$ & & $322$s & & $.0017$ & & $651$s\tabularnewline
& & & & $20$ & & $.0102$ & & $198$s & & $.0052$ & & $355$s & & $.0022$ & & $728$s\tabularnewline
& & & & $50$ & & $.0464$ & & $318$s & & $.0198$ & & $570$s & & $.0095$ & & $1140$s\tabularnewline
& & & & $100$ & & $.0866$ & & $577$s & & $.0407$ & & $881$s & & $.0179$ & & $1957$s\tabularnewline
& & & & $150$ & & $.1926$ & & $761$s & & $.0587$ & & $1151$s & & $.0245$ & & $2811$s\tabularnewline
\bottomrule
\end{tabular}\caption{Finite sample performance of IDC estimator with $\widehat{\boldsymbol{\theta}}_{T}$
and $\widehat{\boldsymbol{\theta}}_{P}$ initialization}
\label{tab:IDC Taddy}
\end{table}
In Table \ref{tab:IDC Taddy}, we present the MSE and running time
of IDC estimators with $\widehat{\boldsymbol{\theta}}_{T}$ and $\widehat{\boldsymbol{\theta}}_{P}$
as initial estimators, respectively, for two different numbers of
iterations $S$. Comparing MSEs of three IDC estimators with different
initial values: $\check{\boldsymbol{\theta}}$ (Table \ref{tab: IDC-Consistent})
and $\widehat{\boldsymbol{\theta}}_{T}$ (Table \ref{tab:IDC Taddy})
or $\widehat{\boldsymbol{\theta}}_{P}$ (Table \ref{tab:IDC Taddy}),
we observe that the IDC estimator with the consistent initial estimator
reduces the MSEs for the same number of iterations.
\begin{table}[H]
\centering
\begin{tabular}{cccccccccccccc}
\toprule
& & $d$ & & $n$ & & $\widetilde{\boldsymbol{\theta}}$ & & $\widehat{\boldsymbol{\theta}}_{PB}^{I}$ & & $\widehat{\boldsymbol{\theta}}_{P}^{I}$ & & $\widehat{\boldsymbol{\theta}}_{T}^{I}$ & \tabularnewline
\midrule
DGP-A & & $20$ & & $500$ & & $.0102$ & & $.0083$ & & $.0102$ & & $.0102$ & \tabularnewline
& & $20$ & & $1000$ & & $.0049$ & & $.0048$ & & $.0052$ & & $.0052$ & \tabularnewline
& & $50$ & & $1000$ & & $.0155$ & & $.0158$ & & $.0198$ & & $.0198$ & \tabularnewline
& & & & & & & & & & & & & \tabularnewline
DGP-B & & $20$ & & $500$ & & $.0007$ & & $.0008$ & & $.0006$ & & $.0303$ & \tabularnewline
& & $20$ & & $1000$ & & $.0002$ & & $.0002$ & & $.0003$ & & $.0403$ & \tabularnewline
& & $50$ & & $1000$ & & $.0005$ & & $.0005$ & & $.0006$ & & $.0121$ & \tabularnewline
& & & & & & & & & & & & & \tabularnewline
DGP-C & & $20$ & & $500$ & & $.0060$ & & $.0062$ & & $.0072$ & & $.0061$ & \tabularnewline
& & $20$ & & $1000$ & & $.0032$ & & $.0032$ & & $.0032$ & & $.0045$ & \tabularnewline
& & $50$ & & $1000$ & & $.0315$ & & $.0318$ & & $.0305$ & & $.0323$ & \tabularnewline
\bottomrule
\end{tabular}\caption{MSE comparison of competing estimators. Number of iterations $S=20$.}
\label{tab:Comparison}
\end{table}
Table \ref{tab:Comparison} presents MSEs of different estimators
for three $\left(d,n\right)$ pairs each. We set the largest $d$
be $50$ so that MLE can be computed in a reasonable time. $\widehat{\boldsymbol{\theta}}_{PB}^{I}$,
$\widehat{\boldsymbol{\theta}}_{T}^{I}$, and $\widehat{\boldsymbol{\theta}}_{P}^{I}$
denote the IDC estimators with $\check{\boldsymbol{\theta}}$, $\widehat{\boldsymbol{\theta}}_{T}$,
and $\widehat{\boldsymbol{\theta}}_{P}$ as the initial estimators
respectively. Besides DGP-A, we consider two additional DGPs to study
the performance of the IDC estimator under different data settings.
In all DGPs, we let $p=5$.
\textbf{DGP-B {[}Poisson{]}}: The random variable $\boldsymbol{V}$
follows a standard normal distribution; $C_{ik}$ follows a Poisson
distribution with mean $e^{\eta_{ik}^{*}}$; and $M$ is obtained
by summing up realizations of the Poisson draws for different choices.
\textbf{DGP-C {[}Mixture{]}}: We let $\boldsymbol{V}$ follow a mixture
of Gaussian distributions with means $0$ and $4$ with standard deviations
$1$ for both distributions. $M$ is also set to follow a mixture
of Gaussian distributions with means $10$ and $60$ and rounded to
the closest integer. The standard deviations are $1$ and $5$ respectively.
We have made these modifications so that some choices are rarely selected
and ensure the robustness of our estimator in those cases.
Based on the simulation result, our IDC algorithm is successfully
executed for all DGPs and exhibits stability. In contrast, we encounter
errors for DGP-C when computing the estimators in \citet{bohning1988monotonicity},
\citet{bohning1992multinomial}, and \citet{simon2013blockwise}.
We see from Table \ref{tab:Comparison} that $\widehat{\boldsymbol{\theta}}_{PB}^{I}$
performs close to $\widetilde{\boldsymbol{\theta}}$ for all DGPs
and $\left(d,n\right)$ pairs. In DGP-B, the initial estimator $\widehat{\boldsymbol{\theta}}_{P}$
is the maximum likelihood estimator. As a result, $\widehat{\boldsymbol{\theta}}_{P}^{I}$
starts with not only a consistent but asymptotically efficient initial
estimator. Even in this case, $\widehat{\boldsymbol{\theta}}_{PB}^{I}$
has comparable MSEs.
In summary, the IDC estimators with all three initial estimators have
finite sample performance similar to the MLE for the DGPs studied
in this section. They are much faster to compute than the MLE for
large $d$ and are feasible even when the MLE might be intractable.
Moreover, if the IDC estimator starts with the consistent initial
estimator $\check{\boldsymbol{\theta}}$, its finite sample performance
will be further improved and is almost the same as the MLE.
\subsection{Inference \label{subsec:Simulation Inference}}
In this section, we illustrate the bootstrap inference procedure introduced
in Section \ref{subsec:Asymptotic-Distributions-Inference}. We investigate
the finite sample performance of the procedure including the size
and power. All the results are based on one thousand Monte Carlo repetitions,
where the number of bootstrap repetitions is five hundred.
We consider the null hypothesis that some element of $\boldsymbol{\theta}^{\ast}$
equals to a specific value. Data are generated from DGP-A introduced
in the previous section. Let the null and the alternative hypotheses
be that $H_{0}:\theta_{11}^{\ast}=0$ and $H_{1}:\theta_{11}^{\ast}\neq0$.
The test statistic is computed as $\left|\frac{\widehat{\theta}_{11}^{I}}{\widehat{se}_{b}\left(\widehat{\theta}_{11}^{I}\right)/\sqrt{n}}\right|$,
where $\widehat{se}_{b}\left(\widehat{\theta}_{11}^{I}\right)$ is
the bootstrap estimate of the standard error of $\widehat{\theta}_{11}^{I}$.
The number of iterations is 10 when computing the IDC estimator. We
set the nominal size as $5\%$ and use the $97.5\%$ quantile of the
standard normal distribution as the critical value.
\begin{table}[H]
\centering \resizebox{\textwidth}{!}{
\begin{tabular}{llcccccccccccccc}
\toprule
& Dev. & & $-0.2$ & & $-0.1$ & & $-0.05$ & & $0$ & & $0.05$ & & $0.1$ & & $0.2$\tabularnewline
\midrule
$d=20$ & $n=250$ & & $.397$ & & $.137$ & & $.067$ & & $.041$ & & $.060$ & & $.139$ & & $.455$\tabularnewline
& $n=500$ & & $.676$ & & $.216$ & & $.102$ & & $.058$ & & $.115$ & & $.234$ & & $.704$\tabularnewline
& $n=1000$ & & $.952$ & & $.395$ & & $.128$ & & $.055$ & & $.167$ & & $.473$ & & $.963$\tabularnewline
& & & & & & & & & & & & & & & \tabularnewline
& Dev. & & $-0.3$ & & $-0.2$ & & $-0.1$ & & $0$ & & $0.1$ & & $0.2$ & & $0.3$\tabularnewline
\midrule
$d=50$ & $n=250$ & & $.333$ & & $.189$ & & $.099$ & & $.077$ & & $.091$ & & $.186$ & & $.313$ \tabularnewline
& $n=500$ & & $.662$ & & $.424$ & & $.172$ & & $.066$ & & $.198$ & & $.391$ & & $.658$ \tabularnewline
& $n=1000$ & & $.895$ & & $.723$ & & $.247$ & & $.063$ & & $.212$ & & $.697$ & & $.924$ \tabularnewline
\bottomrule
\end{tabular}} \caption{Finite sample rejection probabilities for different values of $\theta_{11}^{\ast}$,
$n$, and $d$}
\label{tab: size and power}
\end{table}
In Table \ref{tab: size and power}, we report the finite sample rejection
probabilities of our test for different values of $\theta_{11}^{\ast}$.
Values in the first row of the table indicate the deviation of $\theta_{11}^{\ast}$
from the null hypothesis. When the deviation is zero, the null hypothesis
is true. It can be seen from the table that the finite sample rejection
rates get closer to the nominal size when the sample size increases.
And when the true value $\theta_{11}^{\ast}$ deviates more from the
null hypothesis, the rejection probabilities increase. The same pattern
appears for both $d=20$ and $d=50$. The finite sample performance
of the test when $d=50$ is not as good as that when $d=20$. This
is predictable because there are many more unknown parameters in the
model when $d=50$ than when $d=20$. We expect the results to
improve as the sample size increases for any fixed $d$.
\section{Conclusion \label{sec:Conclusion} }
In this paper, we propose an iterative distributed computing estimator
for the multinomial logistic model that is fast to compute even when
the number of choices is large. When the number of iterations goes
to infinity, we show that our estimator is both consistent and asymptotically
efficient. Based on the simulation study, the computational time of
our estimator increases linearly with the number of choices. Moreover,
our estimator has comparable finite sample performance to MLE when
the latter is computationally feasible.
Extensions abound. First, our IDC estimator can be combined with several
existing algorithms to accommodate more complex settings. For example,
when minimizing $Q_{kn}\left(\boldsymbol{\theta}_{k},\overline{\boldsymbol{\mu}}_{n}\left(\widehat{\boldsymbol{\theta}}^{\left(s-1\right)}\right)\right)$
for $k=1,\ldots,d$, we can replace the gradient of the objective
function with its stochastic approximation calculated from a randomly
selected subset of the data. Such an algorithm is an online algorithm
and might reduce the running time especially when $n$ is large. We
can also employ a one-step Newton-Raphson to compute $\arg\min_{\boldsymbol{\theta}_{k}\in\varTheta_{k}}Q_{kn}\left(\boldsymbol{\theta}_{k},\overline{\boldsymbol{\mu}}_{n}\left(\widehat{\boldsymbol{\theta}}^{\left(s-1\right)}\right)\right)$
in each iteration or even replace the Hessian matrix with its dominant.
These modifications to the IDC estimator have the potential to further enhance computational efficiency, depending on the application. However, their theoretical properties require further investigation. Second, the asymptotic properties
of the IDC estimator in this paper are established for a large but fixed
number of choices. Asymptotic theory allowing for the number of choices
to diverge with the sample size is yet to be established. Third, in
cases where the number of covariates is also large, $\ell_{1}$ or
$\ell_{2}$ regularization could be adopted. In a companion paper,
we develop an asymptotic theory for a regularized iterative distributed
computing estimator.
\newpage{}