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.
73,915 characters · 17 sections · 29 citation commands
Inference From Random Restarts
In many areas of economics, equilibria are used as predictions of observable outcomes and as inputs to counterfactual and policy analysis. When equilibria are not unique, however, the predictive content of a model is weakened, as often models do not provide guidance on which equilibrium will arise in practice. Counterfactual analysis and policy evaluation then become problematic, since different equilibria may respond differently to changes in the environment. For these reasons, applied researchers often seek reassurance that the equilibrium they compute is unique, or that alternative equilibria are “empirically irrelevant.”
Theoretical sufficient conditions for equilibrium uniqueness are often difficult to verify in applied research. As a result, even for well-studied models such as the differentiated-products demand framework of Berry, Levinsohn, and Pakes BerryLevinsohnPakes1995, formal guarantees of equilibrium uniqueness are generally unavailable. In the absence of such results, applied researchers frequently rely on numerical evidence to assess whether multiplicity is a concern in practice.
Iterative numerical methods used to compute equilibria in economics typically depend on the choice of initial conditions. Typically, convergence from a particular initialization does not imply other initial conditions will not end up at a different solution. Nevertheless, a common empirical practice is to restart equilibrium solvers from multiple initial values and to assess whether they converge to the same outcome. Repeated convergence under different starting values is then interpreted as evidence that either the solution is unique or that the solution they found has a dominant basin of attraction.
This reasoning is articulated explicitly in the seminal work of Pakes and McGuire PakesMcGuire1994. While noting that convergence of their algorithm does not theoretically guarantee uniqueness, they report computing equilibria from multiple initial conditions and always converging to the same fixed point, concluding that nonuniqueness “does not seem to be a problem” for the specifications they study. This restart-based diagnostic has since become standard in applied dynamic industrial organization and related areas.
Subsequent work shows that such numerical evidence can be misleading. Besanko et al. besanko2010learning, for example, study a model closely related to the Ericson--Pakes ericson1995markov framework and find up to nine distinct equilibria for some parameterizations, despite using similar computational methods. Taken together, these findings highlight an existing gap in the literature: a lack of formal understanding of the relationship between convergence to the same solution from multiple starting points and the belief over uniqueness of the solution.
Related restart-based reasoning also appears outside economics, in statistics, optimization, and operations research, where nonconvex problems are routinely solved using local or iterative algorithms. In these settings, it is standard to run solvers from multiple starting points or random initializations and to assess robustness based on whether the resulting solutions coincide or achieve similar objective values. For example, in operations research, multistart methods are a canonical heuristic for global optimization, in which a local solver is repeatedly initialized from different starting points and repeated convergence to the same solution is treated as empirical evidence that a global optimum has been identified ugray2007scatter. Similar practices appear in applied statistical estimation and scientific computing, where agreement across starting values is used as a diagnostic for solution stability or dominance.
Across these fields, restart-based arguments play an important practical role, yet they are typically heuristic in nature and lack a formal inferential interpretation. Observing repeated convergence raises natural questions: What does such evidence imply about the distribution of basins of attraction? Under what conditions does repeated convergence provide meaningful support for uniqueness, and how many runs are required before such evidence becomes persuasive? This paper provides a probabilistic framework for answering these questions and for interpreting numerical evidence from repeated solver runs in a principled way.
This paper sits at the intersection of Bayesian inference for discrete latent structure, geometric analysis of basins of attraction, and the interpretation of algorithmic randomness as data. While each of these literatures contains tools relevant to our setting, none directly addresses the inferential problem induced by repeated random initializations of a general solver. Our contribution is to clarify this gap and provide a framework that is both probabilistically coherent and compatible with algorithmic geometry.
\paragraph{Bayesian Inference For Discrete Latent Distributions} At a formal level, repeated runs of a randomly initialized solver produce categorical observations drawn from an unknown discrete distribution over outcomes. A natural starting point is therefore Bayesian inference for discrete distributions, including Dirichlet and Dirichlet process priors (ferguson1973bayesian; ghosal2017fundamentals). These priors are widely used to model uncertainty over unknown probability masses and to allow for an unbounded number of latent categories.
However, directly applying nonparametric Bayesian machinery in our setting is not straightforward. Dirichlet process priors place mass on infinitely many potential atoms, most of which are never observed. In the present context, solver outcomes are not arbitrary labels: they correspond to basins of attraction induced by the geometry of an objective function and the dynamics of an algorithm. As a result, the parameters of a nonparametric prior are not free, but constrained by geometric and dynamical considerations. Our framework makes explicit which aspects of the outcome distribution are identifiable from solver outputs and what assumptions need to be imposed on the solver for obtaining them.
In the particular case of minimization of scalar functions, there exists a literature exploring optimal stopping rules when it comes to random restarts of local solvers boender1987bayesian, zielinski1981statistical, Boender1995. This literature assumes that the function being optimized has an unknown number of finitely many local optima each with positive basin size. In this setting, the optimal stopping rule is to decide on when to stop depending on the observed values for the local optima. Our work differs from this literature as it can be applied to more general problems. Specifically, we develop a rigorous framework that establishes when the solution outcomes can be interpreted as random variables and when there are only finitely many solutions.
\paragraph{Basins of Attraction and Algorithmic Geometry} A separate literature studies basins of attraction in dynamical systems and optimization. Classical results characterize the stability, regularity, and measure-theoretic properties of basins under gradient flows and related dynamics (hirsch2013differential). In nonconvex optimization, recent work has analyzed how properties of the objective function and algorithm—such as smoothness, step size, and curvature—determine convergence behavior and basin structure (e.g., lee2016gradient, ge2015escaping, bhojanapalli2016global).
This literature provides essential intuition for why basins exist and why their sizes may vary dramatically across solutions. However, it is almost entirely deterministic in nature. Basin geometry is treated as an object to be characterized given full knowledge of the function and dynamics, rather than inferred from finite observations. In contrast, our setting assumes that the researcher observes only the terminal outputs of a solver under random initialization, not the underlying geometry itself. Our contribution is not to refine geometric characterizations of basins, but to use the minimal geometric fact—that basins induce a probability distribution over outcomes under random initialization—as the basis for statistical inference.
\paragraph{The Gap Between Geometry and Inference} Bridging the two perspectives above is nontrivial. While geometry determines basin sizes, solver outputs alone do not reveal the geometric features that generated them. Conversely, standard Bayesian nonparametric priors require specifying or learning infinitely many parameters corresponding to both observed and unobserved outcomes, without reference to the constraints imposed by algorithmic dynamics. This mismatch creates a gap: many prior specifications do not agree with the geometry of the problem.
Our approach can be viewed as addressing precisely this gap. Rather than attempting to infer full geometric structure or to endow the outcome distribution with unrestricted nonparametric flexibility, we prove when finiteness of solutions holds and use that to build small yet expressive parametric models. We focus on inferential questions that are well posed given the data—such as posterior beliefs about dominance or uniqueness—and characterize how evidence accumulates as the number of solver runs grows. In this sense, our framework can be seen as deliberately intermediate: richer than purely heuristic arguments, but more limited and thus more disciplined than unconstrained nonparametric modeling.
Assume we are given a problem, e.g., minimizing a function $f$ or finding roots of a vector field $F$. Assume we have a solver $g$ that takes in an initial guess and potentially produces a solution.\footnote{Later we will use $g_Q$ where $g$ is the solver and $Q$ its dynamics.}
Assume we have run the solver with random independent initial conditions $x_1,\dots,x_n \in X$ be i.i.d.\ draws from some initial distribution and gotten the same solution $g(x_1) = g(x_2) = \dots = g(x_n)=x^*$. What should be our belief about $x^*$ being the only solution? Or even simpler, how big is the set of initial conditions that result in $x^*$? A reasonable answer to these questions hinges on two essential requirements:
We acknowledge there is a difference between such properties of a solver $g$ and intrinsic properties of the problem itself. For any problem there may be many solvers and not all solvers are made equal. For example, some solvers may have a hard time converging to certain solutions while others may not.\footnote{For example, besanko2010learning prove that the Pakes-McGuire algorithm may systematically fail to compute certain classes of equilibria.} The question of choosing a solver, even though important, is beyond the scope of this paper. Still, we will revisit this problem in Section (ref) of the paper and provide some remarks.
This section introduces Dynamic Solvers. Intuitively, dynamic solvers are solvers that start with an initial guess and update it continuously according to a fixed updating rule. As we will show, many common solvers including gradient descent and fixed point iterations fall into this category. We take advantage of their fixed dynamic structure to prove measurability and finiteness of solutions under mild conditions.
Given a dynamic solver, we are interested in its possible long-run outcomes, which we call reachable terminal states (RTS), and the sizes of their basins of attraction.
To stay formal while eschewing unnecessary numerical issues, we model "ideal" solvers who do not suffer from numerical imprecision. We then assume that the practitioner only sees the results through an observation map that captures numerical imprecision.
Dynamic solvers are modeled as projected autonomous dynamical systems. That is, dynamical systems that start from an initial solution guess and update it continuously over the space of possible solutions according to a fixed projected updating rule. This continuous-time idealization is a tool: the inference in Section (ref) depends only on the induced terminal-outcome map $g_Q$ and the resulting basin partition, not on whether the solver is implemented in discrete or continuous time. Subsection (ref) provides a formal discussion of this.
Formally, let $X\subseteq \mathbb{R}^d$ be the space of possible solutions and let $\mu$ be a probability measure on $(X,\mathcal{B}_X)$, where $\mathcal{B}_X$ is the Borel $\sigma$-algebra on $X$. Initial conditions are drawn i.i.d.\ from $\mu$. Let $\lambda_X$ denote $d$-dimensional Lebesgue measure restricted to $X$. We assume $\lambda_X \ll \mu$. Equivalently, for any measurable $A\subseteq X$, if $\lambda_X(A)>0$ then $\mu(A)>0$. In particular, any set of initial conditions with positive Lebesgue measure has positive probability under $\mu$, and therefore will be visited almost surely under infinitely many random restarts.
As algorithms may not converge to any valid solution (they may oscillate, etc.), we introduce an absorbing graveyard state $\dagger$ to represent non-convergence/non-existence of a valid terminating solution. Thus the space of all possible outputs of a solver is $\tilde X = X\sqcup\{\dagger\}$ (equipped with the disjoint union $\sigma$-algebra). Before finally defining a solver, we need the following definitions:
One can think of $\operatorname{TC}_X(x)$ as the set of directions that point "inside" $X$ at $x$ and $N_X(x)$ as the normal vectors of all the possible supporting hyperplanes at $x$.
Using this, we can now define projected dynamics:
Our end goal is to work with maps that are limits of such solutions as $t \to \infty$. We need the following well-behavedness assumption to hold for the solution.
The following is a sufficient condition for this to hold:
Now we can define the dynamic solver induced by $Q$ as the limit of the solution to the projected autonomous system.
The strength of this formulation lies in how easy and explicit it makes showing measurability of a solver in spite of its generality. Here, proving measurability of $g$ is reduced to verifying regularity of $Q$, which is often much easier to do.
Now that we have addressed measurability, we turn to finiteness of reachable terminal states. As we only know a solver $g_Q$ through its dynamics $Q$, an ideal condition for finiteness of reachable terminal states should be stated in terms of properties of $Q$.
Let $\mathcal{R}(f)$ be the roots of $f$ on its domain. The following lemma connects reachable terminal states $g_Q(X)$ to the roots of the dynamic $Q^{\Pi_{\operatorname{TC}_X}}$.
Using this, we have the following proposition, a necessary and sufficient condition for finiteness of reachable terminal states:
For example, in case of gradient flows on $f:\mathbb{R}^d \to \mathbb{R}$, the reachable terminal states of $g_Q$ are contained in the set of stationary points of $f$. So $f$ having finitely many stationary points is a sufficient condition for $g_{-\nabla f}$ to have finitely many reachable terminal states.
Often times, we cannot directly verify finiteness of roots of $Q^{\Pi_{\operatorname{TC}_X}}$ for different reasons: $Q$ may be a complicated function or we may only have access to an oracle that lets us evaluate it. To rescue, we provide conditions under which finiteness of roots holds for almost all dynamic solvers belonging to a class of dynamic solvers in a measure-theoretic sense.
For technical reasons that will become clear later, we need the following extra assumption on the domain of $Q$:
When $\mathcal{Q}$ is finite-dimensional, the standard notion for almost all is full Lebesgue measure. For infinite-dimensional $\mathcal{Q}$, Lebesgue measure is not informative, as it will be either zero or infinite on all open sets. Prevalence (hunt1992prevalence) and Haar-null sets christensen1972sets are two generalizations of "almost all" that agree with full Lebesgue measure in finite-dimensional spaces and still apply to infinite-dimensional ones. As it is more widely used, we use prevalence here.
This gives us a way to argue that almost all solvers in a class have finitely many reachable terminal states, so assuming finiteness of reachable terminal states for an unknown solver of that class is a sensible assumption.
We want sufficient conditions on $\mathcal{Q}$ such that the set of $Q$ such that $g_Q$ has finitely many reachable terminal states is prevalent in $\mathcal{Q}$. Our strategy for proving finiteness of reachable terminal states is to invoke the following lemma.
In the previous section, we showed that $g_Q(X)$ is finite iff $Q^{\Pi_{\operatorname{TC}_X}}$ has finitely many roots. The nice thing about roots of a function is that they form a closed set if the function is continuous. Though $Q^{\Pi_{\operatorname{TC}_X}}$ is not continuous, we show that its roots belong to the union of roots of two separate continuous functions.
What remains is having a condition that implies the discreteness of these sets in $X$. Once we have that, we can invoke Lemma (ref) on $\mathcal{R}(Q)$ and $\mathcal{R}(Q^{\Pi_{\operatorname{TC}_X}}_{|\partial X})$. When a function $f$ is $C^1$, a sufficient condition for roots of $f$ to be isolated is transversality of $f$ to 0: $f \pitchfork 0$. As transversality works well with prevalence, we would like to phrase discreteness in terms of transversality to 0. This works well for $Q$ as it is $C^1$ by our definition of classes of dynamic solvers. $Q^{\Pi_{\operatorname{TC}_X}}_{|\partial X}$, however, is not $C^1$ as it involves $\langle Q(x), n(x) \rangle_+$.
To overcome this, we first show that roots of $Q^{\Pi_{\operatorname{TC}_X}}_{|\partial X}$ can be characterized as roots of a related $C^1$ map. Define: \[ \tilde{Q}:\partial X \times [0,\infty) \to \mathbb{R}^{d}, \quad \tilde{Q}(x, \alpha) = Q(x) - \alpha n(x). \] The following lemma shows why $\tilde{Q} \pitchfork 0$ implies discreteness of $\mathcal{R}(Q^{\Pi_{\operatorname{TC}_X}}_{|\partial X})$.
So to prove discreteness of $\mathcal{R}(Q)$ and $\mathcal{R}(Q^{\Pi_{\operatorname{TC}_X}}_{|\partial X})$, it suffices to prove $Q \pitchfork 0$ and $\tilde{Q} \pitchfork 0$. The following lemma connects prevalence and transversality.
A good candidate condition for prevalence is the perturbed evaluation map being a submersion.
We can use this to prove the following sufficient condition for prevalence of transversality.
As the evaluation maps being submersions (requirements of Theorem (ref)) can feel intangible, we provide the following sufficient condition that is easier to grasp. This condition will be used in section (ref) to verify prevalence of finitely many RTS.
As mentioned before, solvers in practice may not produce the exact output of $g_Q$ due to numerical imprecision. Moreover, in practice the notion of “same solution” is operational: two outputs are identified if they are indistinguishable at a prescribed numerical tolerance. To model this, we introduce a measurable observation/identification map $\mathcal O_\varepsilon:\tilde X \to \tilde X$ where $\varepsilon>0$ represents the effective numerical resolution (e.g., floating-point precision, stopping criteria, etc). The observed output from an initial condition $x\in X$ is $\mathcal O_\varepsilon(g_Q(x))$. In particular, decreasing $\varepsilon$ may refine the observed basin partition by splitting sets of initial conditions whose terminal points agreed previously.
With measurable $g_Q$ that has finitely many reachable final states, $\mathcal O_\varepsilon\circ g_Q$ is also measurable and has finitely many reachable final states. Hence the analysis of this paper applies verbatim to $\mathcal O_\varepsilon\circ g_Q$.
The statistical analysis of the next section depends only on the partition of $X$ into basins induced by $\mathcal O_\varepsilon\circ g_Q$. Thus $\mathcal O_\varepsilon$ matters only insofar as it changes this partition. Varying $\varepsilon$ corresponds to coarsenings/refinements of the observed partition, which is exactly the setting of Assumption (ref). In practice, small numerical errors do not significantly change the induced partition on $X$. Hence we work with $\mathcal O_\varepsilon=\mathrm{Id}$. When numerical imprecision changes the partition, the results of next section should be read as inference about the observed partition induced by $\mathcal O_\varepsilon\circ g_Q$ rather than the idealized one induced by $g_Q$.
It turns out that dealing with and making inference about the size of a basin of attraction is much easier than dealing with the number of reachable terminal states (RTS). This is not surprising as it is not easy to rule out the existence of RTS with infinitesimally small basins of attraction.
Let $x_1, x_2, \dots, x_n$ be i.i.d. samples from $\mu$ and define $Y_i := g_{Q}(x_i)$ for some dynamic solver $g_Q$. Define the event of getting the same RTS $x^* \in \tilde{X}$ for all $n$ runs as: \[ H_n(x^*) = \{Y_1 = Y_2 = \dots = Y_n=x^*\}. \]
We show how conditioning on $H_n$ changes our beliefs about the size of $x^*$'s basin of attraction $s_{x^*}$ and the number of RTS. We start with a simpler model that allows us to make inference about the basin size. We then tackle the harder problem of making inference about the number of RTS using heavier machinery.
Given a prior $\Pi$ on $s_{x^*}$, we are interested in computing the posterior: \[ \Pi(s_{x^*} \mid H_n(x^*)) = \frac{\mathbb{P}(H_n(x^*) \mid s_{x^*}) \Pi(s_{x^*})}{\mathbb{P}(H_n(x^*))}. \] Note that the likelihood $\mathbb{P}(H_n(x^*) \mid s_{x^*}) = s_{x^*}^n$ as each independent run has probability $s_{x^*}$ of converging to $x^*$.
We are interested in how this posterior concentrates around 1 as $n$ increases. Dually, that is how $\Pi(A_\varepsilon\mid H_n(x^*))$ concentrates around $0$ where $A_\varepsilon = \{s_{x^*} < 1 - \varepsilon\}$ for some $\varepsilon > 0$.
A simple choice is to use a Beta prior $\operatorname{Beta}(\alpha, \beta)$ for some $\alpha, \beta > 0$. This prior has full support on $[0,1]$ and is conjugate to the Bernoulli likelihood that will arise from our observations. Under this prior, we have:
This result formalizes the intuition that observing repeated convergence to the same RTS increases our belief that its basin of attraction is large. The exponential concentration rate shows that this increase in belief is quite rapid, hence hundreds of repeated runs can provide strong evidence for a dominant basin of attraction.
Though this result is derived under a Beta prior, we generalize this to a wide class of priors. Essentially, we need to impose some regularity conditions on the prior density near 1, ensuring that the prior does not vanish too quickly as we approach 1. We start by the following bound:
Combining these two lemmas, we obtain the following generalization of the exponential concentration result:
So under the mild condition that the prior assigns positive density to a neighborhood of 1, we again obtain exponential concentration bounds of the posterior around 1. As this bound depends on $\eta$ which is a primitive of our model, it may not be easy to interpret. So far, we assumed positive density near 1. If we can further strengthen this assumption, we can free the inequality from $\eta$.
A necessary condition for uniqueness of an RTS is that its basin of attraction has size one. Yet under any prior $s_x$ that admits a density, this event has zero prior probability. The reason is that under the prior, the event $s_x=1$ has probability zero, and Bayesian updating will not give a positive probability to an event that has zero prior probability. So to talk about uniqueness we need more. We introduce two approaches: a simpler spike and slab prior and a heavier Mixture of Finite Models prior.
So far, we have focused on the sizes of basins of attraction. Bridging the gap between basin sizes and number of RTS is easy if all RTS are $\mu$-observable. Then having a unique RTS is equivalent to it having basin size of one.
In practice, however, this assumption may be too strong. For example, we may want to know whether the local minima of a Morse function $f$ found by GD are global. The Morse condition rules out flatness, hence all local minima have basins of positive size. The problem is that saddle points are also RTS but they have basins of size zero.
To overcome this, we reformulate our question about uniqueness. Instead of asking whether there is a unique RTS, we ask whether there is a unique $\mu$-observable RTS. In the GD example, this corresponds to asking whether there is a unique local minimum rather than whether there is a unique RTS.
Now we can go back to our previous analysis about basin sizes. Remember that any event that has probability zero under the prior will have probability zero under the posterior. So any prior over a basin size that admits density such as the ones discussed so far will never give a positive probability for uniqueness of a $\mu$-observable RTS. Hence we need a prior that assigns positive probability to the event $s_{x^*} = 1$.
Spike and slab priors are a popular choice for such situations where we want to assign positive probability to a specific value (the spike) while allowing for uncertainty around other values (the slab). They are a mixture of a point mass at a specific value and a continuous distribution over the remaining values.
Let $\delta_{x}$ be the Dirac delta measure centered at $x$, and let $W$ be a measure over $[0,1]$ that admits density. Define a random variable $Z$ that determines whether the spike is true or the slab. The following is our spike and slab prior: \[ s_{x^*} | Z=0 \sim W \quad , \quad s_{x^*} | Z=1 = 1 \] If $\Pi(Z=1) = p$, then the marginal prior on $s_{x^*}$ is: \[ \Pi = p \delta_{1} + (1-p) W \]
Now we are interested in $\Pi(s_{x^*}=1\mid H_n(x^*))$. Notice that $\Pi(s_{x^*} = 1\mid H_n(x^*)) = \Pi(Z=1\mid H_n(x^*))$. Using Bayes' rule, we have:
To talk about the rate of convergence, we need to understand the rates of convergence for $\mathbb{E}_W[s_{x^*}^{n}]$.
Having this bound, we can speak to the rate of convergence of $\Pi(s_{x^*}=1\mid H_n(x^*))$.
The spike and slab prior is simple and intuitive. However, it only lets us test whether there is a unique RTS. One way to extend this is to model $(Y_1,\dots,Y_n)$ as i.i.d. draws from a categorical distribution whose number of components is unknown. This is a Mixture of Finite Models (MFM) (miller2018mixture) specialized to categorical data (component labels observed), with no component-specific parameters.
MFMs are a class of Bayesian models that allow for uncertainty in the number of components in a mixture model. They are particularly useful when the true number of components is unknown and we want to infer it from the data.
In our context, an MFM constitutes a prior over the number of components $\mathbb{P}(K)$ and a distribution of basin sizes given $K$:
As is, this model is underspecified. It depends heavily on the choice of $P_K$ and $\alpha_i^j$ without clear guidelines on how to choose them. We will argue that given certain symmetries, there are natural choices for these parameters.
As we have no reason to believe that any RTS has a larger basin of attraction compared to others, the Dirichlet prior should be symmetric. That is ensured by setting $\alpha_k^i = \alpha_k$ for all $1 \leq i \leq k$.
To further tie $\alpha_k$s together, we rely on the following standard principle:
This is reminiscent of the observation map $\mathcal{O}$ discussed at the end of Section (ref). In our setting, we impose invariance under refinement--coarsening: if each RTS is split into $m$ indistinguishable sub-RTS and we then coarsen by summing the $m$ sub-basins, the induced law on the coarsened basin sizes should agree with the $K$-RTS prior.
Thus far, we have reduced infinitely many parameters $\alpha_k^i$ to a single parameter $\alpha$. $\alpha$ captures how concentrated the distribution of basins of attraction is: smaller $\alpha$ gives spikier distributions while larger $\alpha$ gives more uniform distributions. Hence, if $\alpha$ is small, $H_n(x^*)$ is not as informative about $K$ compared to when $\alpha$ is large. So our inference is contingent on the choice of $\alpha$.
Since the prior depends on $P_K$ and $\alpha$, we can either adopt conservative defaults (e.g.\ a heavy-tailed $P_K$ and small $\alpha$), or calibrate these hyperparameters on benchmark instances from the same problem class. To do so, parametrize $P_K(\cdot;\theta)$ by a finite-dimensional parameter $\theta$. Sample independent instances $Q_i\in\mathcal Q$ from an application-specific instance distribution, and for each instance run the solver from $m$ i.i.d.\ initial conditions $x^i_t\sim\mu$ to obtain terminal states $Y^i_t := g_{Q_i}(x^i_t)$, $t=1,\dots,m$. We then choose $(\hat\theta,\hat\alpha)$ by empirical Bayes, i.e.\ by maximizing the marginal likelihood of the observed outcomes $\{(Y^i_1,\dots,Y^i_m)\}_{i=1}^R$ under the model (integrating out basin sizes and mixing over $K$).
Using this lemma, we can now derive rate-tight concentration bounds for $\mathbb{P}(K=1\mid \alpha,H_n(x^*))$.
Empirical industrial organization provides natural applications for our methods. An incomplete list includes static models of supply and demand BerryLevinsohnPakes1995, vertical relations crawford2012welfare,grennan2013price,crawford2018welfare, and dynamic oligopoly collard2013demand,fowlie2016market. In these models, equilibrium uniqueness confers significant theoretic and econometric advantages BerryLevinsohnPakes1995,bajari2007estimating,pakes2007simple. These models have no proofs of equilibrium uniqueness, however, so the random-restart heuristic is frequently employed cussen2026nash,PakesMcGuire1994,crawford2012welfare. In this section, we show that our framework for inference can be applied to industrial organization.
Due to space constraints we focus on models with mixed logit demand and multi-product firms, the workhorse of static supply and demand estimation.\footnote{Although multiple Bertrand-Nash equilibria do not hinder estimation of BLP-style models BerryLevinsohnPakes1995, they make counterfactuals hard or impossible to evaluate. This is the problem created by multiple equilibria that is commonly cited. However, multiple equilibria also constrain the construction of optimal instruments BerryLevinsohnPakes1995. As Conlon and Gortmaker conlon2020best show, the use of optimal instruments greatly enhances the performance of BLP.} Here, we focus on the fixed-point iteration algorithm proposed by Morrow and Skerlos morrow2011fixed.
Let there be $d$ inside goods and an outside option(with price normalized to $0$) in the market, and assume firms have constant returns to scale. $c \in (0,\infty)^d$ denotes the vector of marginal costs and $x \in (0,\infty)^{d}$ the vector of prices. The vector field constructed from firms' first-order conditions can be rewritten as \[ F(x) = x - c - \zeta(x)\label{eqn:ms_zeta_vf} \] where $\zeta$ is the vector of markups whose expression is provided in morrow2011fixed in terms market shares without closed form. They propose to solve (ref) by doing fixed-point iterations on $x \mapsto c + \zeta(x)$. Although the roots of (ref) solve all firms' first-order conditions, they need not be equilibria. Therefore, the authors point out, if multiple solutions may exist, it is important to search for alternative solutions to (ref) by restarting computations from different points in the strategy space. We take this idea further by applying our inference methods.
We first demonstrate how this fits within our framework from Section (ref). Let $X\subseteq\mathbb{R}^d$ denote the strategy set, i.e. prices that firms can choose. skerlos2010fixed prove that, without loss, $X$ can be a box whose lower-bounds are marginal costs and which has finite upper bounds that never constrain firms' choices. This satisfies Assumption (ref) but not Assumption (ref), which is used in Section (ref) to invoke finiteness results for reachable terminal states (RTSs). However, since we know that equilibrium prices always lie inside the box, it is without loss to set $X$ to be any super-set of the box. In particular, we set $X$ to be a superset of the box, in the positive orthant, that satisfies Assumption (ref) by having a $C^2$ boundary.
skerlos2010fixed also prove that, for standard mixed logit models, $\zeta$ is $C^1$. By compactness of $X$ it is Lipschitz, so their algorithm corresponds to a constrained Picard flow with $Q(x) = c + \zeta(x) - x$ and a corresponding dynamic solver $g_Q$ (Remark (ref)). By Proposition (ref), $g_Q$ is measurable.
For finiteness of reachable terminal states, we can not directly check the finiteness of $\mathcal{R}(Q^{\Pi_{\operatorname{TC}_X}})$ as $\zeta$ does not have a closed form expression. We provide two alternative, prevalence-based arguments. The weaker argument leverages the fact that almost all Picard flows on $C^1$ functions have finitely many reachable terminal states (Remark (ref)). The weakness of this argument stems from the fact that all possible $c+\zeta(x)$ make a proper subset of $C^1$ functions. Ideally, we would want to know almost all Picard flows on these particular functions have finitely many reachable terminal states.
The second alternative, not subject to that weakness, is proving that the class dynamic solvers generated by varying marginal costs and markups is a complete vector space. Then, we can invoke Propositions (ref) and (ref) to show that having finitely many reachable terminal states is prevalent.
Canonically, marginal costs are parameterized by $\theta \in \mathbb{R}^{d_C}$ and markups by $\omega \in \mathbb{R}^{d_\zeta}$ conlon2020best. Consider the class of dynamic solvers $\mathcal{Q} = \{ Q_{(\theta,\omega)} \mid (\theta,\omega) \in \mathbb{R}^{d_C}\times\mathbb{R}^{d_\zeta} \}$ where $Q_{(\theta,\omega)}(x) = c_\theta + \zeta_\omega(x) - x$.
$\mathcal{Q}$ may not be closed under standard addition and scalar multiplication, so it may not be a vector space under these operations. This can be problematic as prevalence is defined on vector spaces. The solution is easy. We define addition $\oplus$ and scalar multiplication $\odot$ on $\mathcal{Q}$ that make it a vector space. Intuitively, we put the linear structure induced by the parameters $(\theta,\omega)$ on $\mathcal{Q}$. \[ Q_{(\theta_1, \omega_1)} \oplus Q_{(\theta_2, \omega_2)} := Q_{(\theta_1 + \theta_2, \omega_1 + \omega_2)} \quad , \quad \lambda \odot Q_{(\theta, \omega)} := Q_{(\lambda \theta, \lambda \omega)} \] Consider the perturbations $P=\{Q_{(\theta,0)}|\theta \in \mathbb{R}^{d_C}\}$, that is, perturbations that keep markups as is and move marginal costs. If the structure of $c_\theta$ is rich enough to be a submersion, i.e., we can move the marginal cost in any direction by changing $\theta$, then we can invoke Proposition (ref) to show evaluation maps are submersions and then Proposition (ref). A straightforward choice is to set $d_C=d$ and $c_\theta = e^\theta$ where $e^\theta$ is the element-wise exponential of $\theta$. (This happens to be the standard function for marginal costs in mixed logit models conlon2020best.) In this case, the Jacobian of $c_\theta$ is a diagonal matrix with non-zero entries, so it is a submersion. Hence, having finitely many reachable terminal states is prevalent in $\mathcal{Q}$.
Having put things in context, we can apply results of section (ref). We will show how seeing $n$ random restarts converging to the same RTS (event $H_n$) should change our belief about the size of its basin of attraction and its uniqueness.
Remember that $A_\varepsilon = \{s_{x^*} < 1 - \varepsilon\}$ is the event of the basin of attraction of the observed RTS being smaller than $1-\varepsilon$. Put the uniform prior $\operatorname{Beta}(1,1)$ on $s_{x^*}$. The posterior after observing $H_n$ is $\operatorname{Beta}(n+1,1)$ with CDF $F(x) = x^{n+1}$. Hence, $\Pi(A_\varepsilon \mid H_n) = F(1-\varepsilon) = (1-\varepsilon)^{n+1}$. Table (ref) provides numerical values for its complement: $\Pi(A_\varepsilon^c \mid H_n) = 1 - \Pi(A_\varepsilon \mid H_n)$, probability of $s_{x^*} \geq 1 - \varepsilon$. This agrees with the exponential concentration promises of proposition (ref) as the Beta distribution satisfies assumption (ref).
We first consider the spike and slab model. For that, we need to define two things, $p=\mathbb{P}(Z=1)$ (the spike) and $W$ (the slab). To rule out the problem with having RTS with basins of size zero, instead of asking whether there is a unique RTS, we ask whether there is a unique $\mu$-observable RTS. This corresponds to equilibria that are not fragile to arbitrarily small perturbations of the initial condition.
We set $W=\operatorname{Beta}(1,1)$. We vary $p=\mathbb{P}(Z=1)$ and $n$ and report $\Pi(s_{x^*}=1\mid H_n)$ from (ref). Table (ref) shows the resulting values.
We now turn to the Mixture of Finite Models (MFM) prior of section (ref). We compare two priors on $K$: a geometric prior $\mathbb{P}(K=k)=2^{-k}$ (which has a heavier tail than Poisson) and a (zero-truncated) Poisson(1) prior $\mathbb{P}(K=k)=\frac{e^{-1}}{(1-e^{-1})k!}$ supported on $k\ge 1$. Table (ref) reports the lower bounds implied by Theorem (ref) on $\mathbb{P}(K=1\mid H_n)$ for different values of $\alpha$ and $n$.
We formalize the random-restart heuristic by treating solver outputs under i.i.d.\ initializations as samples from the basin-induced distribution over reachable terminal states (RTS). For a dynamic solver $g_Q$ and restart distribution $\mu$, draw $x_1,\dots,x_n\sim\mu$ and record the terminal outcomes $g_Q(x_i)$. Repeated convergence to a particular RTS $x^*\in\tilde X$ corresponds to the event $H_n(x^*)=\{g_Q(x_1)=\dots=g_Q(x_n)=x^*\}$.
We show that broad classes of dynamic solvers induce measurable terminal-outcome maps and, under mild conditions, admit only finitely many RTS. This makes probabilistic and statistical reasoning about restart outcomes well posed. Conditioning on $H_n(x^*)$, we derive explicit posterior concentration bounds for the basin mass $s_{x^*}=\mu(g_Q^{-1}(x^*))$. We also provide concentration bounds for the posterior on uniqueness of RTS under both a spike-and-slab and an MFM prior. While the posterior on the basin mass concentrates exponentially fast in $n$, the posterior on uniqueness concentrates at a rate polynomial in $n$.
We illustrate how these ideas apply to an industrial organization application. We demonstrate how the problem of finding the equilibrium of mixed logit demand with multi-product firms can be stated as a dynamic solver within our framework. We then report numerical values for posteriors over basin size and uniqueness of solution under different specifications.
The inference of this paper is conditional on event $H_n(x^*)$. In practice, one may be interested in inference under different events, perhaps events involving observation of multiple RTS.
Finally, we address here the choice of solver discussed at the beginning of Section (ref). From a statistical perspective, we prefer all the rts to be $\mu$-observable. However, this is a property of the solver and not the problem itself. As the choice of the best solver depends on the problem at hand, we can only encourage paying attention to this aspect of the problem.