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.
92,453 characters · 18 sections · 75 citation commands
Batched Adaptive Network Formation
Networks play a central role in many economic and management applications and are often modeled as graphs, with agents as nodes, their connections as edges, and pairwise production or interaction outcomes as edge weights. For example, a manager forms teams of workers to complete tasks, with team output reflected in edge weights. Similarly, a teacher seats students in pairs, or a social platform recommends connections, both aiming to foster friendships or social ties, where the success or failure of these efforts is recorded as binary edge weights. In all these cases, the decision maker---be it a manager, platform, or teacher---faces the same core challenge: how to (re)allocate effective connections that maximize total network production, measured by the sum of edge weights.
The effectiveness of these connections hinges on complementarities between agents, that is, which workers collaborate productively or which students or users are likely to form successful friendships or social ties.\footnote{We use the term complementarity to refer to the interaction effects between heterogeneous agents that drive pairwise outcomes. These effects may reflect productivity complementarities in teams or affinity in social tie formation. We adopt a unified terminology across applications for clarity.} In practice, the decision maker must first learn from data: by forming connections, observing outcomes, and inferring agents' unobserved heterogeneity and complementarities---the exploration phase. Once sufficient information is gathered, a network formation policy can be determined to reallocate connections for improved overall network production---the exploitation phase. Offline policies that carry out the two phases in a single round often perform poorly: updating the policy too early with insufficient data risks inaccurate estimates, while delaying too long may miss valuable opportunities for improvement. Both scenarios reflect an imbalance in the exploration–exploitation trade-off, resulting in suboptimal outcomes.
A more effective approach is to adaptively update the policy as new data become available, dynamically balancing exploration and exploitation. In this paper, we propose such an adaptive, or online, policy designed to form increasingly productive networks. We frame this as a Batched Multi-Armed Bandit (Batched MAB) problem, where each possible pair of agent latent types corresponds to an arm. The algorithm proceeds in batches: after each batch of observed outcomes, it updates beliefs about complementarities and uses these to solve a maximum-weight matching problem to determine the next batch’s network. This updated policy is applied to pair agents for improved outcomes, with the resulting data fed into the next round of inference and optimization. This iterative structure enables the algorithm to progressively refine its decisions, yielding continuous improvement over time.
One practical challenge decision makers often face is that only edge-level outcomes (e.g., the success of teamwork, the development of friendships or social ties) are observable, while individual contributions typically remain unobserved. Moreover, these outcomes depend on the agents' unobserved characteristics (e.g., workers' soft skills, students' or platform users' personalities) often in complex, nonlinear patterns not amenable to parametric specification.\footnote{The importance of unobserved agent heterogeneity on complementarities has been documented across various contexts, including teamwork (weidmann2021team, weigel2024supermodular, and xu2024heterogeneous), medical matching (agarwal2015empirical), and online dating (lee2016effect).} These challenges motivate our adoption of a weighted stochastic block model (WSBM), which captures agents' unobserved heterogeneity through $K$ discrete latent types (see, e.g., bonhomme2021teams and jochmans2024nonparametric). The expected edge-level outcome---reflecting complementarities between agents---depends on the latent types of the agents in each agent pair and is parameterized by $\boldsymbol{\theta}\in\mathbb{R}^{K\times K}$. A key advantage of this WSBM-based framework is that it imposes no functional form restrictions on $\boldsymbol{\theta}$, which allows for more flexible modeling of heterogeneous complementarities compared to production models that that impose structure, such as strongest/weakest-link specifications (see, e.g., JohKamKriLi2018Exploration), linear models (jochmans2019fixed), or nonlinear models with specific functional forms (e.g., the constant elasticity of substitution form adopted by ahmadpoor2019decoding).\footnote{Both jochmans2019fixed and ahmadpoor2019decoding model agents’ unobserved heterogeneity as continuous-valued fixed effects, whereas the WSBM framework captures it through $K$ discrete latent types. Due to this difference, the comparison between our framework and these two works is not entirely direct.}
To rigorously integrate WSBM-based network models into a tractable online decision-making framework, we build on H\'ajek-Le Cam's theory of convergence of statistical experiments, drawing in particular on the recent extension by hirano2023asymptotic for batched adaptive experiments. We establish a batched local asymptotic normality (batched LAN) result for the WSBM-based network production model, assuming that the network formation policies are determined prior to each batch. This result shows that, within each batch, the complex network production experiment can be asymptotically represented by a much simpler limit experiment involving a single draw from a normal distribution. Subsequently, the entire $T$-batch network experiment can be approximated by a $T$-stage Gaussian bandit environment, each stage generates an independent Gaussian signal, with mean determined by $\boldsymbol{\theta}$ and variance depending on the network formation policy employed in the corresponding batch.\footnote{In particular, type pairs that are sampled more frequently within a batch yield more precise parameter estimates (i.e., signals with lower variance); see Section (ref) for details.} Under Gaussianity, these per-batch signals serve as maximum likelihood estimates (MLEs) in the limit.
Nevertheless, computing the finite-sample MLE is generally infeasible for WSBM models due to the intractability of their likelihoods (see graham2020econometric and bonhomme2021teams). To address this, we adopt a variational estimation approach, applied independently to each batch, and show that it is asymptotically equivalent to the MLE. In addition to estimating $\boldsymbol{\theta}$, the variational method also produces local variational approximations (LVAs) as a valuable by-product for inferring agents' latent types. These LVAs, given as categorical distributions over the $K$ possible types, approximate the posterior distribution of the latent type variable $\boldsymbol{Z}$, and are shown via simulation to rapidly concentrate around true latent types. All these results are developed under batched sparse network asymptotics, reflecting practical constraints---for example, that an agent may participate in only a limited number of tasks per batch (e.g., one month), regardless of how many potential collaborators they have.\footnote{Specifically, the batched sparse asymptotics allows the expected degree (i.e., the number of edges per node) in each batch to grow as slowly as logarithmically in the number of nodes.}
These theoretical results convey a simple yet powerful message: batched variational estimates can be treated as Gaussian signals for the production parameter $\boldsymbol{\theta}$, while the batched LVAs can be viewed as categorical signals for agents' latent types $\boldsymbol{Z}$. Building on this insight, we aggregate information on $\boldsymbol{\theta}$ and $\boldsymbol{Z}$ across batches using standard Bayesian methods. From the resulting posteriors, we either compute point estimates such as posterior means in a greedy manner, or sample parameter values following the Thompson sampling principle. These estimates or sampled values are then incorporated into a constrained maximum-weight matching problem, the solution of which determines the network formation policy for the next batch. In particular, we propose a Hybrid Greedy-Thompson (HGT) algorithm, which uses posterior mean as the point estimate for the production parameter (“greedy-in-$\boldsymbol{\theta}$”) and samples agent-type realizations from the LVAs (“Thompson-in-$\boldsymbol{Z}$”). Simulation results show that our HGT algorithm substantially improves overall outcomes within a few batches, while simultaneously refining estimates for both $\boldsymbol{\theta}$ and $\boldsymbol{Z}$ over time. We further evaluate performance in non-stationary environments with agent entry and exit. In these settings, HGT algorithm remains robust, though requiring a few more batches to achieve comparable performance.
Our paper relates to the extensive literature on network models; for a comprehensive review, see graham2020network. In particular, we contribute to the growing body of work on network production models based on the weighted stochastic block model (WSBM) (see, e.g., bonhomme2021teams, jochmans2024nonparametric), which flexibly capture “nonparametric” complementarities among heterogeneous agents. These models have been applied in recent empirical studies such as weigel2024supermodular (tax collector teams) and xu2024heterogeneous (sales-force teams) to estimate heterogeneous complementarities and evaluate offline network formation policies. We extend this literature by developing an online framework in which network formation policies are updated adaptively, enabling continuous improvement in production outcomes over time.
While motivated by the online setting, our theoretical results---namely, the batched local asymptotic normality (LAN) for WSBM-based network production models and the asymptotic equivalence of variational estimators to maximum likelihood estimators under sparse network asymptotics---also readily apply to offline estimation by treating the full sample as a single batch. This work extends the asymptotic theory of feasible variational inference (consistency and asymptotic normality), established for classical (unweighted) stochastic block models (SBMs; see, e.g., bickel2009nonparametric, bickel2013asymptotic), to weighted, exogenously formed network SBMs.
We also contribute to the literature on multi-armed bandit (MAB) algorithms (e.g., thompson1933likelihood, lai1985asymptotically, auer2002using, kasy2021adaptive; and cesa2023adaptive for applications in economics), particularly in the context of team formation and assignment problems (johari2018exploration, johari2021matching, eichhorn2022online). Incorporating WSBM into this setting enables flexible modeling of complementarities between heterogeneous agents whose characteristics are unobserved. However, it introduces two key departures from standard MAB formulations. First, the arms---here, type pairs---are not directly observed, but instead their inferred probabilities through local variational approximations. Second, the network formation policy are solutions to constrained optimization problems, rather than explicit functions of simple descriptive statistics. The first departure is specific to our WSBM-based setting, while the second connects our work to the combinatorial semi-bandit literature, in which feasible policies often involve selecting matchings subject to combinatorial constraints (e.g., audibert2014regret, chen2013combinatorial, and kasy2023matching). To address these challenges, we adopt the batched bandit framework (e.g., perchet2016batched, zhang2020inference) and leverage Le Cam’s asymptotic theory for batched adaptive experiments by hirano2023asymptotic, enabling the integration of complex network econometric modeling into online decision-making.
The rest of the paper is organized as follows. Section (ref) introduces the WSBM-based network model and frames the online network formation problem as a batched MAB. Section (ref) presents our variational approximation estimation strategy, followed by asymptotic analysis in Section (ref). Section (ref) details our algorithm, including Bayesian updating and an agent-level constrained maximum-weight matching problem, making it readily applicable to practical settings. Section (ref) presents a comprehensive Monte Carlo study demonstrating the effectiveness of our algorithm. Section (ref) concludes.
Consider a pool of agents indexed by $1, \dots, n$, represented as nodes in a graph. Each agent belongs to one of $K$ latent types, which capture heterogeneity in agents' unobserved characteristics. The latent type of agent $i$ is recorded in the random variable $Z_i$, which takes values in $\{1,\dots,K\}$. We assume $Z_i$ are independently and identically distributed (i.i.d.) with probabilities
for $a = 1,\dots,K$, where $\pi(1) + \cdots + \pi(K) = 1$.
Agents are assigned into $m$ pairs, represented as edges in the graph, to complete $m$ tasks.\footnote{We use the term “task” broadly to refer to the intended outcome of a pairwise assignment, including both collaborative activities and social interactions---e.g., teaming agents for joint work, seating students to foster friendships, or recommending platform users to follow each other.} Each pair $j \in \{1, \dots, m\}$ is formed by two distinct agents $i_1(j), i_2(j) \in \{1, \dots, n\}$ with $i_1(j) \neq i_2(j)$. Throughout, we focus exclusively on two-agent collaborations and therefore use “network formation” and “agent pairing” interchangeably.\footnote{Extending our framework to collaborations involving more than two agents would require hypergraph theory and generalization of the asymptotic results in Section (ref), which we leave for future work.}
The outcome of the Task $j$, completed by Pair $j$, is modeled as the weight of the corresponding edge in the graph and is captured by a random variable $Y_j$. The weight $Y_j$ can be binary, where $Y_j = 1$ indicates success and $Y_j = 0$ indicates failure, or continuous, where $Y_j \in \mathbb{R}$. We assume that the outcomes $Y_j$ are i.i.d.\ with distribution
which depends on the type pair, $\left(Z_{i_1(j)},Z_{i_2(j)}\right)$, where $\theta = \theta_{(Z_{i_1(j)},Z_{i_2(j)})} \in\mathbb{R}$ is an unknown parameter. For simplicity, in this pairwise production model $p_{\theta}(\cdot|a,b)$, we restrict our analysis to the univariate case, where $\theta_{a b}$ is a scalar for each type pair $(a,b)$, with $a,b = 1,\dots,K$ (abbreviating $\theta_{a b}$ as $\theta$ when used as a subscript). However, our results extend naturally to the multivariate, albeit with more involved notation. We impose no functional form restrictions on $p_{\theta}(\cdot|a,b)$ beyond mild regularity conditions (required for Assumption (ref) to hold).
To motivate the policy design problem, we begin with an oracle setting in which the agents’ latent types and the pairwise production model are known to the decision maker. The decision maker's objective is to form agent pairs that maximize total network output. Specifically, the decision maker's objective is to form agent pairs that maximize total network output. Specifically, she/he aims to solve
where $\Upsilon(a,b) \equiv \mathbb{E}_{\theta}\left(Y_j \,|\, Z_{i_1(j)} = a, Z_{i_2(j)} = b\right)$ is the expected production of a type $(a, b)$ pair, and the pairing policy function $\psi(a,b)$ determines the relative frequency (or probability mass) with which each type pair $(a, b)$ is assigned across $m$ pairs.
The unrestricted policy space $\Psi$, to which $\psi$ belongs, comprises all valid probability distributions over $K\times K$ possible type pairs:
That is, each $\psi(a,b)$ as type-$(a,b)$ pair probability lies in $[0, 1]$, and all probabilities sum to one. In real-world scenarios, the optimization problem often incorporates additional constraints. For instance, in the teamwork example, a key restriction is the workload constraint (see weigel2024supermodular and xu2024heterogeneous) which ensures a fair distribution of tasks across workers. Another commonly imposed constraint in bandit literature is the “clipping constraint,” designed to avoid sampling probabilities that are too small and that would lead to inconsistent parameter estimates. We formally state these two constraints as follows.
Consequently, the optimization problem becomes
where $\Psi^{c} = \left\{\psi \in \Psi \,\big|\, \psi \textrm{~satisfies~} (\ref{eqn:constraint_workload}) \textrm{~and~} (\ref{eqn:constraint_clipping}) \right\}$ denotes the constrained policy space. In this oracle scenario where the agent types $\boldsymbol{Z} = (Z_1,\dots,Z_{n})^\top$ are observed and the pairwise production models $p_{\theta}(\cdot | a,b)$ are known, this optimization problem can be solved readily via linear programming.
In practice, however, neither the agents' types nor the pairwise production models are known to the decision maker. This necessitates a two-phase procedure: first, estimate these parameters using data from sampled agent pairings and their outcomes---the exploration phase; then, use the obtained estimates to solve the network output maximization problem ((ref))---the exploitation phase. When this process is carried out in a single round, it is known as an offline policy, where the decision maker commits to a fixed policy without collecting new data to refine the policy further. A more effective strategy is to run adaptive experiments, in which the agent-pairing policy is updated sequentially using newly collected data. This iterative process, referred to as an online policy, balances exploration (learning parameters) and exploitation (maximizing network output based on current beliefs). In the following subsection, we cast this iterative learn-and-earn procedure as a non-standard Batched Multi-Armed Bandit problem and develop an algorithm specifically tailored for WSBM-based adaptive network formation.
In the classical Multi-Armed Bandit (MAB) framework, a decision maker sequentially selects from a set of treatments---akin to a gambler choosing among slot machine arms---each associated with an unknown reward distribution. The objective is to maximize cumulative rewards over time by balancing exploration (sampling different arms to learn their reward distributions) and exploitation (favoring the arm currently believed to yield the highest expected return). In our setting, each agent-type pair $(a, b)$ represents an arm, and edge-level outcomes serve as stochastic rewards governed by the conditional distribution $p_{\theta}(\cdot|a, b)$. However, the WSBM-based network formation problem is inherently more complex, departing from the standard MAB paradigm in two fundamental ways.
First, in addition to the production model $p_{\theta}$, the agent types $\boldsymbol{Z}$ are themselves unobserved. This introduces a key departure from the standard setting: when assigning a pair, the underlying type combination---and thus the specific arm being pulled---is unknown. Instead, decisions must rely on the posterior distribution of $\boldsymbol{Z}$, obtained via the variational approximation (see Section (ref)), which induces corresponding probabilities over type combinations. As data accumulate, this posterior distribution becomes increasingly concentrated around the true latent types, thereby improving the decision-maker’s understanding of each agent and enabling better-informed pairing decisions.
Second, the development of the sampling policy that governs arm-pulling probabilities differs significantly from that of a standard MAB problem, where policies are often explicit functions of descriptive statistics such as accumulated rewards and arm-pulling frequencies. These functions are typically guided by inferences such as the posterior distribution, point estimates, or upper confidence bounds of the unknown parameters, corresponding to Thompson sampling (thompson1933likelihood), greedy heuristics, or UCB algorithms (lai1985asymptotically, auer2002using), respectively. In contrast, our network formation policy is determined by solving the optimization problem in ((ref)), where the unknown parameters are replaced with their inferred values. This procedure is further complicated by practical restrictions such as the workload constraint.
To tackle these challenges, we formulate our algorithm in a batched bandit setting, where the network formation policy $\psi$ is updated at discrete time intervals. Specifically, we carry out the network formation and production process over $T$ batches for a potentially evolving group of agents or their types. In each batch $t \in \{1,\dots,T\}$, we assign $n_t$ agents into $m_t$ pairs according to a policy $\psi_t(a,b)$; we record the pair assignments in mappings $(i_{t,1}(j), i_{t,2}(j))$ and the corresponding production outcomes in variable $Y_{t,j}$, for $j = 1,\dots,m_t$. Both $T$ and $m_t$ are either predetermined or exogenous, and thus treated as deterministic.\footnote{In real-world applications, both the batched decision making and the number of tasks in each batch, $m_t$, can naturally arise from institutional or market constraints. For instance, the real estate team setting described in xu2024heterogeneous is well suited to a batched implementation, where the manager could periodically update the assignment policy (e.g., monthly or quarterly), and the number of properties $m_t$ to be managed in each period is exogenously determined by market availability. Similar batching structures apply in other domains such as education rohrer2021proximity and online platforms rajkumar2022causal.} Importantly, we require that each policy $\psi_t$ be specified prior to batch $t$, based solely on information available from the previous batches. We formalize this “pre-determined policy” condition in the following assumption.
The pre-determined policy condition in Assumption (ref) is naturally satisfied by any feasible network formation algorithm that determines the pair assignment policy for each batch based solely on data from previous batches, and keeps that policy fixed throughout the batch. This condition allows us to treat the network formation $\boldsymbol{I}_t$ for batch $t$ as fixed, or nonrandom. In settings where agents form pairs themselves, this condition implies the Network Exogeneity assumption (see bonhomme2021teams), which requires that network formation be independent of agent-pair-specific shocks, thereby ruling out cases where agents have prior knowledge of such shocks and use it to influence pair formation.\footnote{Relaxing this condition in such cases requires explicitly modeling the network formation process; see also bonhomme2021teams.} In our setting, however, this concern does not arise, as all pairings are assigned by a decision maker according to pre-specified policies.
We conclude this section with a brief overview of our online network formation policy. The complete procedure is detailed in Algorithm (ref) in Section (ref). Its core components---including variational estimation, asymptotic results, Bayesian updating, and an agent-level optimization problem---are presented in Sections (ref)--(ref).
Online Network Formation Policy: Starting with an initial policy $\psi_1$, the algorithm proceeds iteratively over batches $t = 1, \dots, T$ according to the following procedure:
In this section, we introduce our estimation strategy based on a variational approximation, applied independently to each batch. For notational simplicity, we suppress the batch index $t$ throughout the section (e.g., for the number of agents $n_t$, number of tasks $m_t$, policy function $\psi_{t}$, and outcome variables $Y_{t,j}$ and $\boldsymbol{Y}_t$), unless needed for clarity.
Under the i.i.d.\ assumption for $Y_j$, the conditional likelihood of $\boldsymbol{Y} = (Y_1,\dots,Y_m)^\top$ given agents' latent types $\boldsymbol{Z}$ is
where the vector-valued parameter $\boldsymbol{\theta}\in\mathbb{R}^{K\times K}$ collects all pair-specific production parameters $\theta = \theta_{a b}$ (of pairwise model $p_{\theta}(\cdot|a,b)$), $a,b=1,\dots,K$. The types of the agents assigned to task $j$, $(Z_{i_1(j)},Z_{i_2(j)})$, is from the network formation data $\boldsymbol{I}$. While $\boldsymbol{I}$ determines the structure of the likelihood, it is assumed to be pre-determined and treated as fixed (i.e., non-random), and is thus omitted from the notation. The likelihood of $(\boldsymbol{Y},\boldsymbol{Z})$ is then given by
where $p_{\boldsymbol{\pi}}(\boldsymbol{Z}) = \prod_{i=1}^{n}\pi(Z_i)$ denotes the marginal likelihood of $\boldsymbol{Z}$. Following bickel2013asymptotic, we refer to the model in which $\boldsymbol{Z}$ is observed as the complete graph model.
In the complete graph model, the type probabilities $\boldsymbol{\pi}$ and production parameters $\boldsymbol{\theta}$ can be directly estimated by maximizing the joint likelihood $f_{\boldsymbol{\pi},\boldsymbol{\theta}}(\boldsymbol{Y},\boldsymbol{Z})$. However, when agent types $\boldsymbol{Z}$ are unobserved, maximum likelihood estimation becomes computationally intractable, as it requires working with the marginal likelihood $g_{\boldsymbol{\pi},\boldsymbol{\theta}}(\boldsymbol{Y})$, computed by summing over all possible type assignment vectors $\boldsymbol{e} = (e_1,\dots,e_n)^\top$, where each $e_i\in\{1,\dots,K\}$:
This summation takes $K^{n}$ operations and quickly becomes computationally infeasible as $n$ grows. To circumvent this challenge, we employ a mean-field variational approximation to the marginal likelihood below, a technique widely used in latent variable models to enable tractable inference.
Write the joint likelihood as $f(\boldsymbol{Y},\boldsymbol{Z}) = h(\boldsymbol{Z}|\boldsymbol{Y})g(\boldsymbol{Y})$, we approximate the conditional likelihood $h(\boldsymbol{Z}|\boldsymbol{Y})$ with a factorized term, $q(\boldsymbol{Z}) \equiv \prod_{i}q_i(Z_i)$. Each component $q_i$, referred to as a local variational approximation (LVA), serves as an individual approximation to the posterior distribution for the latent variable $Z_i$ over the support $\{1,\dots,K\}$. The discrepancy between $q(\boldsymbol{Z})$ and $h(\boldsymbol{Z}|\boldsymbol{Y})$ is measured by the Kullback–Leibler (KL) divergence, defined as
Expanding this expression yields
where
and the boldface $\boldsymbol{q} = (q_1,\dots,q_n)$ denotes the collection of local variational approximations. The space where $\boldsymbol{q}$ takes value is given by $$\mathcal{Q} \equiv \left\{ (q_1,\dots,q_n) \in \big([0,1]^{K}\big)^{n} \,\big|\, q_i(1)+\cdots+q_i(K) = 1, ~ \forall i = 1,\dots,n \right\}.$$
The term $J(\boldsymbol{Y})$ is known as the Evidence Lower Bound (ELBO) which provides a lower bound on the log-likelihood $\logg(\boldsymbol{Y})$, with the observed data $\boldsymbol{Y}$ being the “evidence”. This inequality holds because the KL-divergence is always non-negative, implying $J(\boldsymbol{Y}) \leq \logg(\boldsymbol{Y})$. In practice, the ELBO offers a computationally tractable proxy for the generally intractable log-likelihood $\logg(\boldsymbol{Y})$. Specifically, it simplifies to
where computing the summations now requires only $n\times K + m\times K^2$ operations in total.
We aim to make $J_{\boldsymbol{\pi},\boldsymbol{\theta};\boldsymbol{q}}(\boldsymbol{Y})$ as close as possible to the log-likelihood $\logg_{\boldsymbol{\pi},\boldsymbol{\theta}}(\boldsymbol{Y})$, so that it serves as a good approximation. This is achieved by jointly maximizing $J_{\boldsymbol{\pi},\boldsymbol{\theta};\boldsymbol{q}}(\boldsymbol{Y})$ with respect to $\boldsymbol{q}$, which is treated as a newly introduced parameter. In the meantime, we optimize over the model parameters $\boldsymbol{\pi}$ and $\boldsymbol{\theta}$ to obtain their variational estimates. This joint optimization is naturally embedded in the framework of the Expectation-Maximization (EM) algorithm (see, e.g., daudin2008mixture), which proceeds as follows:
This section establishes the asymptotic results of the WSBM and its batched variational estimators above. We begin in Section (ref) with the complete graph model, where the agent types $\boldsymbol{Z}_t$ are observed. In this setting, we show that the complete graph model converges to a Gaussian shift experiment (in H\'ajek-Le Cam's sense), provided that each pairwise production model $p_{\theta}\left(\cdot \,|\, a,b\right)$ satisfies a similar convergence. Section (ref) then demonstrates that the variational estimates (based on $J_{\boldsymbol{\pi},\boldsymbol{\theta};\boldsymbol{q}}(\boldsymbol{Y}_t)$) are asymptotically equivalent to the complete graph MLEs (based on $f_{\boldsymbol{\pi},\boldsymbol{\theta}}(\boldsymbol{Y}_t,\boldsymbol{Z}_t)$). For notational clarity, we vectorize the production parameter $\boldsymbol{\theta}$ as an element of $\mathbb{R}^{K^2}$ throughout this section.
We first establish a batched local asymptotic normality (batched LAN) result for the complete graph model, where the latent types $\boldsymbol{Z}_t$ are observed. This is done under two high-level conditions on the pairwise production models $p_{\theta}(\cdot | a,b)$ for $a,b = 1,\dots,K$, stated in the following assumption.
Assumption (ref)(a) assumes that the production model for each type pair $(a,b)$ is locally asymptotically normal (LAN); see, e.g., van2000asymptotic. That is, the log-likelihood ratio admits a quadratic expansion in which the first-order term converges in distribution to a normal random variable with variance equal to twice the second-order term. This LAN condition is mild and holds for a broad class of smooth parametric models. For continuous $Y_j$, the LAN property follows from the differentiable in quadratic mean (DQM) condition, which requires the existence of a measurable function $\dot{\ell}_{\theta}(Y_j|a,b)$ such that, as $\epsilon \to 0$,
for $a,b=1,\dots,K$. The DQM condition is implied by absolute continuity of the density function $p_{\theta}(\cdot|a,b)$ with a square-integrable first-order derivative $\dot{p}_{\theta}(\cdot|a,b)$. In that case, the score function is given by $\dot{\ell}_{\theta}(\cdot|a,b) = -\dot{p}_{\theta}(\cdot|a,b)/p_{\theta}(\cdot|a,b)$ and the Fisher information is $\mathcal{I}_{\theta}(a,b) = \mathbb{E}_\theta\big[\dot{\ell}_{\theta}^2(Y_j|a,b)\big]$. For binary $Y_j$, the LAN condition arises from a standard second-order Taylor expansion, where the score function is $\dot{\ell}_{\theta}(Y|a,b) = Y - \theta_{a b}$ and the Fisher information is $\mathcal{I}_{\theta}(a,b) = \theta_{a b}(1 - \theta_{a b})$; see bickel2013asymptotic.
Assumption (ref)(b) can be heuristically motivated---though not formally implied---by part (a). Specifically, maximizing the LAN-form log-likelihood in ((ref)) with respect to $h$ yields its maximum likelihood estimate, $\hath_{a b}^\mathrm{MLE} = \mathcal{I}_{\theta}^{-1}(a,b)\frac{1}{\sqrt{m}}\sum_{j=1}^{m}\dot{\ell}_{\theta}(Y_j|a,b)$ which weakly converges to $\mathcal{N}(0,\mathcal{I}_{\theta}^{-1}(a,b))$; in addition, at $\theta_{a b}$, $\sqrt{m}(\hat\theta_{a b}^\mathrm{MLE} - \theta_{a b}) = \hath_{a b}^\mathrm{MLE}$. We nonetheless impose part (b) as a formal assumption, which remains quite mild: For continuous $Y_j$, the asymptotic normality of the corresponding MLE is generally satisfied in smooth parametric models under mild regularity conditions; see, e.g., van2000asymptotic. For binary $Y_j$, the asymptotic normality of the MLE follows directly from the central limit theorem for binomial experiments.
Define $\omega_{a} = \log\frac{\pi(a)}{1 - \sum_{l=1}^{K-1}\pi(l)}$ for $a = 1,\dots,K-1$, and let $\boldsymbol{\omega} = (\omega_{1},\dots,\omega_{K-1})^\top$. This one-to-one transformation allows us to use the notations $f_{\boldsymbol{\pi},\boldsymbol{\theta}}$, $g_{\boldsymbol{\pi},\boldsymbol{\theta}}$, and $J_{\boldsymbol{\pi},\boldsymbol{\theta};\boldsymbol{q}}$ interchangeably with $f_{\boldsymbol{\omega},\boldsymbol{\theta}}$, $g_{\boldsymbol{\omega},\boldsymbol{\theta}}$, and $J_{\boldsymbol{\omega},\boldsymbol{\theta};\boldsymbol{q}}$, respectively. The following proposition establishes the batched LAN result for the complete graph model in which $\boldsymbol{Z}$ is observed.
In establishing this result, the pre-specified policy condition in Assumption (ref) allows the pairing assignments $(i_{t,1}(j),i_{t,2}(j))^\top$ (collected in $\boldsymbol{I}_t$) to be treated as fixed, so the likelihood depends solely on the production outcomes $Y_{t,j}$. The i.i.d.-ness of $Y_{t,j}$ then implies that the joint likelihood of $\boldsymbol{Y}_t$ factorizes into the product of individual likelihoods, from which Proposition (ref) follows readily under the high-level conditions in Assumption (ref). A detailed proof is provided in Appendix (ref) for completeness.
The complete graph LAN result in Proposition (ref), developed at a given point $(\boldsymbol{\pi},\boldsymbol{\theta})$ in the parameter space, characterizes also the asymptotic behavior of statistics in its local neighborhood. Specifically, by Le Cam's third lemma, under the local alternative $f_{\boldsymbol{\omega}+\frac{\boldsymbol{\nu}}{\sqrt{n}},\boldsymbol{\theta}+\frac{\boldsymbol{h}}{\sqrt{m}}}$ and as $m_t\to\infty$ and $n_t\to\infty$,
for $\boldsymbol{\omega}\in\mathbb{R}^{K}$ and $\boldsymbol{h}\in\mathbb{R}^{K^2}$.
More importantly, the LAN result within each batch invokes the Asymptotic Representation Theorem (see le1972limits, hajek1970characterization, van2000asymptotic). This theorem states that a sequence of complex experiments---here, the WSBM-based network production model---can be approximated by a simpler limit experiment, whose likelihood ratio is the limit of the original sequence of likelihood ratios. In our context, the limit experiment corresponds to the so-called Gaussian shift experiment in ((ref)), where we observe single observations, or signals, of random variables $\boldsymbol{\mathcal{Q}}_{\boldsymbol{\omega},t}^{-1}\boldsymbol{\Delta}_{\boldsymbol{\omega},t} \sim \mathcal{N}(\boldsymbol{\nu}, \boldsymbol{\mathcal{Q}}_{\boldsymbol{\omega},t}^{-1})$ and $\boldsymbol{\mathcal{Q}}_{\boldsymbol{\theta},t}^{-1}\boldsymbol{\Delta}_{\boldsymbol{\theta},t} \sim \mathcal{N}(\boldsymbol{h}, \boldsymbol{\mathcal{Q}}_{\boldsymbol{\theta},t}^{-1})$. A direct calculation confirms that this limit experiment has a likelihood ratio equal to $\Lambda_t(\boldsymbol{\nu},\boldsymbol{h})$.
Across batches, the pre-specified policy condition (Assumption (ref)) and the batchwise application of variational estimation ensure that the sequence of batched likelihood ratios converges jointly: $$\big(\Lambda_{f,1}^{(n,m)}(\boldsymbol{\nu},\boldsymbol{h}),\dots,\Lambda_{f,T}^{(n,m)}(\boldsymbol{\nu},\boldsymbol{h})\big) \Rightarrow \big(\Lambda_1(\boldsymbol{\nu},\boldsymbol{h}),\dots,\Lambda_T(\boldsymbol{\nu},\boldsymbol{h})\big).$$ This invokes the multi-stage Asymptotic Representation Theorem (Theorem 2 of hirano2023asymptotic), which states that any statistical decision rule based on a convergent sequence of statistics in the original batched network formation problem admits a corresponding representative rule in the limiting $T$-observation Gaussian shift experiments. We will leverage this result to develop our adaptive network formation policy in Section (ref).
It is worth noting that the Gaussian signals across batches are not identically distributed. In particular, the normal distribution $\mathcal{N}(\boldsymbol{h}, \boldsymbol{\mathcal{Q}}_{\boldsymbol{\theta},t}^{-1})$ has batch-specific variance determined by the network formation policy $\psi_t$ and the production model (via Fisher information $\mathcal{I}_{\theta}(a,b)$). As shown in ((ref)), the variance is given by $\mathcal{Q}_{\boldsymbol{\theta},t}^{-1}(a,b) = \frac{1}{\psi_t(a,b)\mathcal{I}_{\theta}(a,b)}$. This reflects the intuition that, for fixed $\mathcal{I}_{\theta}(a,b)$, assigning more $(a,b)$-type pairs (i.e., a larger $\psi_t(a,b)$ value) leads to more precise estimation of the production parameter $\theta_{a b}$, hence lower variance.
We now return to the real-world setting where the decision maker cannot directly observe the agent types $\boldsymbol{Z}_t$. We show below that the marginal model $g_{\boldsymbol{\pi},\boldsymbol{\theta}}(\boldsymbol{Y}_t)$, and subsequently its variational approximation $J_{\boldsymbol{\omega},\boldsymbol{\theta};\boldsymbol{q}}(\boldsymbol{Y}_t)$, is asymptotically equivalent to the complete graph model $f_{\boldsymbol{\pi},\boldsymbol{\theta}}(\boldsymbol{Y}_t,\boldsymbol{Z}_t)$.
This result is developed under a batched sparse network asymptotic regime, formalized in Assumption (ref). This mirrors the sparse network asymptotics for the classical stochastic block model (see, e.g., bickel2013asymptotic, amini2013pseudo), where the expected number of edges per node, or degree, grows at the rate $O(\logn_t)$. Such a growth rate reflects empirical patterns in large-scale networks, where each node typically forms only a few links. Compared to dense network asymptotics, the sparse regime yields more realistic and accurate approximations in these dateset. In our batched WSBM setting, we similarly allow the average the number of (pre-specified) edges per node, $m_t/n_t$, in each batch $t$ to grow as slowly as $O(\logn_t)$ as the number of nodes $n_t$ increases. This is also motivated by real-world scenarios, such as the sales-force context in xu2024heterogeneous, where each agent can only participate in a limited number of tasks within a given batch period (e.g., one month), regardless of the total number of potential collaborators. We formalize this condition as follows.
To establish the desired equivalence result, we also impose additional regularity conditions on the production model, given in Assumptions (ref) and (ref). These conditions are readily verified once the distribution is specified. For instance, in the case of continuous $Y_j$, they are satisfied when $f$ is normal. In the case of binary outcomes $Y_j$, these conditions, adapted to our weighted SBM setting, can likewise be verified following the approach of bickel2013asymptotic for proving their Theorem 1. These conditions, together with the proof, are organized in Appendix (ref).
Theorem (ref)(i) shows that, under batched sparse network asymptotics, the likelihoods with and without observing the latent agent types $\boldsymbol{Z}_t$ are asymptotically equivalent. This, in turn, implies that the MLEs based on the complete-graph likelihood $f(\boldsymbol{Y}_t,\boldsymbol{Z}_t)$ and those based on the marginal likelihood $g(\boldsymbol{Y}_t)$ are asymptotically equivalent, as formalized in part (ii). While this result primarily serves as a theoretical intermediate step, since maximizing $g(\boldsymbol{Y}_t)$ is generally infeasible due to the computational challenges discussed in Section (ref), it nonetheless provides a critical foundation for analyzing the variational approximation $J(\boldsymbol{Y}_t)$. In Theorem (ref) below, we establish that $J(\boldsymbol{Y}_t)$ is asymptotically equivalent to $\logf(\boldsymbol{Y}_t,\boldsymbol{Z}_t)$ under the same sparse asymptotics. As a result, feasible variational estimates based on $J(\boldsymbol{Y}_t)$ are asymptotically equivalent to the oracle MLEs based on the complete-data likelihood. The proof is provided in Appendix (ref).
Theorem (ref)(i) further implies that the $J_{\boldsymbol{\omega},\boldsymbol{\theta};\boldsymbol{q}}(\boldsymbol{Y})$-based log-likelihood ratio proxy, $\Lambda_{J,t}^{(n,m)}(\boldsymbol{\nu},\boldsymbol{h},\boldsymbol{q}) = \frac{J_{\boldsymbol{\omega}+\boldsymbol{\nu}/\sqrt{n},\boldsymbol{\theta} + \boldsymbol{h}/\sqrt{m};\boldsymbol{q}}(\boldsymbol{Y})}{J_{\boldsymbol{\omega},\boldsymbol{\theta};\boldsymbol{q}}(\boldsymbol{Y})}$, is asymptotically equivalent to the complete-graph log-likelihood ratio, $\Lambda_{f,t}^{(n,m)}(\boldsymbol{\nu},\boldsymbol{h})$; both converge weakly to the same Gaussian shift experiment likelihood ratio $\Lambda_t(\boldsymbol{\nu},\boldsymbol{h})$ (defined in Proposition (ref)). Using arguments similar to those in Proposition (ref), we obtain a LAN-type expansion for $\Lambda_{J,t}^{(n,m)}(\boldsymbol{\nu},\boldsymbol{h},\boldsymbol{q})$ as
where
It is straightforward to verify that $\Lambda_{J,t}^{(n,m)}(\boldsymbol{\nu},\boldsymbol{h},\boldsymbol{Z}) = \Lambda_{f,t}^{(n,m)}(\boldsymbol{\nu},\boldsymbol{h})$, meaning the ELBO function recovers the complete-graph log-likelihood ratio when the variational approximations $q_i(a)$ in the former are replaced by the true class indicators $\mathbbm{1}\{Z_{t,i} = a\}$. Therefore, their asymptotic equivalence implies that each $q_i(a)$ converges to $\mathbbm{1}\{Z_{t,i} = a\}$ within each batch $t$ under the sparse network asymptotics. This aligns with the established consistency results for label estimation of $\boldsymbol{Z}$ in the classical stochastic block models using methods such as the EM algorithm; See, e.g., bickel2009nonparametric, amini2013pseudo, wang2017likelihood for more detailed results. That said, our online network formation algorithm does not aim to explicitly recover the latent labels $\boldsymbol{Z}$ within each batch. Instead, it leverages $\boldsymbol{q}$ in its capacity as batched categorical approximations to the posterior distribution of $\boldsymbol{Z}$, as elaborated in the next section.
This section details the implementation of our online network formation algorithm. Section (ref) describes the Bayesian updating procedure, which refines posterior beliefs about the pairwise production parameters and latent agent types using batched variational estimates as signals. In Section (ref), we incorporate these updated posteriors and translate the policy-level constrained optimization problem from Section (ref) to the agent level, making it readily applicable in practice. Finally, Section (ref) presents the complete algorithmic procedure.
The theoretical results in Section (ref) convey a simple yet powerful insight: the complex batched network formation and production problem can be approximated by a multi-stage Gaussian shift experiment. Specifically, for pairwise production parameter $\theta_{a b}$, where $a, b = 1, \dots, K$, we receive a single Gaussian signal at each stage $t = 1,\dots,T$, with mean $\boldsymbol{\theta}$ and variance determined by the network formation policy $\psi_t(a,b)$ and the production model $p_{\theta}(\cdot|a,b)$. Under Gaussianity, these signals are thus stage-specific maximum likelihood estimates (MLEs) for $\boldsymbol{\theta}$. Since the variational estimates are shown to converge weakly to these MLEs, they serve as feasible finite-sample analogues.
In light of this, we apply Gaussian Bayesian updating pairwise to each $\theta_{a b}$ using their batched variational estimates. For each $(a,b)$ pair with $a,b = 1,\dots,K$, let $\hat\theta_{t,a b}$ denote the variational estimate of $\theta_{a b}$ from batch $t$, and let $\widehatse_{t,a b}$ denote its standard error. Assuming a Gaussian prior $\mathcal{N}(\mu_{0,a b}, \sigma_{0,a b}^2)$, the posterior distribution after round $t$ remains Gaussian, $\mathcal{N}(\mu_{t,a b}, \sigma_{t,a b}^2)$, and is updated recursively via the standard Gaussian Bayesian formula
for $t = 1,\dots, T$ and $a,b = 1,\dots,K$.
Turning to the agents' latent types, recall that for each agent $i$, the variational estimation via EM algorithm in each batch $t$ produces a local variational approximation $\hatq_{t,i}(a)$, a categorical distribution over types $a = 1,\dots,K$, approximating the posterior of the agent's type $Z_i$. These categorical signals can be aggregated recursively for each agent via Bayesian classification. Assuming a categorical prior $\omega_{0,i}(a)$ on $Z_i$, the posterior after batch $t$, denoted by $\omega_{t,i}(a)$, remains categorical and updates as
for $t = 1,\dots,T$, $i = 1,\dots,n$, and $a = 1,\dots,K$.
Beyond a natural choice for sequential learning, Bayesian updating offers an additional practical advantage: it systematically incorporates prior knowledge about unknown parameters, particularly valuable in real-world applications. For instance, managers or teachers may have prior insights into workers’ styles or students’ personalities when forming teams or assigning seatmates. Similarly, users on social platforms often disclose latent traits during registration. When no prior information is available, the framework defaults to uninformative priors: infinite-variance normal distributions for $\theta_{a b}$ and discrete uniform distributions for $Z_i$.
In practical settings, the policy-level constrained optimization problem in ((ref)) must often be reformulated at the agent level for (at least) two key reasons. First, certain practical constraints need to be implemented at the individual level. For instance, the workload constraint reflects the fact that each agent can handle only a limited number of tasks within a given time frame and must therefore be applied individually. Second, reliable inference of latent types requires that each agent be matched frequently enough to ensure the convergence of their estimated type distribution $\hatq_i(a)$). We now reformulate the problem accordingly as a agent-level maximum weight matching problem.
Let $\tilde\theta_{a b}$ and $\postZ_{i}$ denote some realizations (e.g., point estimates, random draws) of the pairwise production parameters and agent latent types, respectively, obtained from their latest posteriors. Based on them, we aim to solve
subject to
where $\tilde\Upsilon(a,b) = \int_{y\in\mathcal{Y}} p_{\tilde\theta}(y|a,b)\mathrm{d} y$ is the expected weight for a type $a$-$b$ pair, and $x_{i i'} = 1$ indicates that agents $i$ and $i'$ are paired (otherwise, $x_{i i'} = 0$). In words, by selecting which edges to form among all possible agent pairs, we seek to maximize the overall expected output ((ref)), subject to the constraint ((ref)) that exactly $m$ edges are formed.
We next reformulate the workload and clipping constraints in ((ref)) and ((ref)) at the agent level. The workload constraint becomes
for some $d_l \leq d_h \in \mathbb{N}^{+}$, ensuring that each agent is assigned a number of tasks within the allowed range. The clipping constraint is similarly recast as
for a given clipping rate $\lambda\in(0, 0.5)$ and all $a,b = 1,\dots,K$.
The maximum weight matching problem in ((ref)), subject to ((ref)) and, when relevant, ((ref)) and ((ref)), can be formulated as an integer linear program, the solution to which yields a valid matching that maximizes total weight.
In the classical MAB problem, the central challenge is balancing exploration (pulling less-explored, potentially suboptimal arms to gain information) and exploitation (pulling arms currently believed to yield the highest reward). This trade-off is typically managed by policies that determine arm-pulling probabilities, which can often be interpreted as functions of the posteriors of the unknown rewards. In particular, the greedy algorithm pulls the arm with the highest average reward, corresponding to the point estimate provided by the Gaussian posterior means. Thompson sampling (thompson1933likelihood) assigns probabilities based on the chance each arm is optimal: it draws a sample from each arm’s posterior and selects the arm with the highest draw. The Upper Confidence Bound (UCB) algorithm (see, e.g., lai1985asymptotically, auer2002using) chooses the arm with the highest upper bound of a confidence (or credible, under Gaussianity) interval, thereby allocating more probability to less-explored arms, which tend to have wider bounds.
Leveraging these insights, we construct our online network formation algorithm using the posteriors $\mathcal{N}(\mu_{t,a b}, \sigma_{t,a b}^2)$ for $\theta_{a b}$ and $\omega_i$ for $Z_i$. Specifically, we sample agent-type realizations $\tildeZ_i$ from categorical posteriors $\omega_i$, following the principle of Thompson sampling. Given the exploration inherent in sampling $\tildeZ_i$, we we adopt a greedy heuristic for $\theta_{a b}$ by simply using the posterior mean as a point estimate, setting $\tilde\theta_{a b} = \mu_{t,a b}$. These sampled types and point estimates, denoted with tildes, are then used in the maximum weight matching optimization described in Section (ref), whose solution defines the sampling policy for the next round. We refer to this “greedy-in-$\boldsymbol{\theta}$, Thompson-in-$\boldsymbol{q}$” strategy as the Hybrid Greedy-Thompson (HGT) algorithm, detailed in Algorithm (ref).
We conduct a series of Monte Carlo experiments to evaluate the finite-sample performance of the Hybrid-Greedy-Thompson (HGT) algorithm introduced in Section (ref). In each experiment, HGT assigns agents into pairwise networks to complete binary-outcome tasks ($Y_j \in \{0, 1\}$). We consider two cases: one with two latent agent types ($K = 2$) and the other with three ($K = 3$), which allows us to assess the algorithm’s robustness to different levels of latent structure. In both cases, key parameters are chosen to broadly reflect real-world settings commonly studied in the empirical literature.
For the $K = 2$ case, we draw on the classroom seating and friendship formation context in rohrer2021proximity to guide our choices of network size and the productivity regime. That study finds that students with similar observable characteristics tend to form friendships more successfully. While our model focuses on unobserved heterogeneity (latent types), we retain the underlying intuition by letting agents with similar latent types have a higher probability of productive interaction. In the $K = 3$ case, we calibrate the agent-to-task ratio and average productivity level to be aligned with xu2024heterogeneous, which studies salesforce team formation and success in real estate deals. A key finding there is that agents with intermediate solo performance tend to have better team productivity when paired with any other type of agent. We retain this feature in our data-generating process by assigning one latent type to consistently deliver strong team performance across all pairings.
For the implementation of the EM algorithm, we initialize from multiple random starting points to reduce sensitivity to local optima, a widely recognized issue for SBM-based models (see, e.g., nowicki2001estimation, bickel2013asymptotic, bonhomme2021teams). All simulation results are based on $100$ repetitions.
We consider a fixed pool of $n = 32$ agents (e.g., students), slightly smaller than the maximum class size studied in rohrer2021proximity. Each student belongs to one of $K = 2$ latent types, representing unobserved characteristics such as personality traits that are unlikely to change over short time horizons. In each replication, student types are independently drawn from a multinomial distribution with equal probabilities $(0.5, 0.5)$. Once generated, these types remain fixed throughout the simulation and are treated as unknown parameters. The probability of $Y_j = 1$ for a given pair $j$---indicating friendship formation---is specified by:
with rows and columns ordered by type $(1, 2)$.
We run our HGT algorithm over $T = 6$ batches. In each batch, the $32$ students are assigned to $384$ seatmate pairs, resulting in an average of $24$ ($= 384\times2/32$) pairings per student. This corresponds to a teacher who reshuffles student seating on a fixed schedule (e.g., weekly), collects data, and updates the pairing policy after each batch. Depending on the frequency of reshuffling, each batch can be interpreted as spanning several months or one or more academic semesters. Since each student participates in exactly the same number of pairings per batch, we impose the agent-level “workload” constraint in ((ref)) with $d_l = d_h = 24$ in this case.
Figure (ref) shows the empirical posterior distributions of each $\theta_{ab}$ for $a, b = 1, 2$, after batches 1, 3, and 6. The histograms become progressively more concentrated and shift toward the true parameter values, marked by grey dotted vertical lines. This indicates that the algorithm effectively refines the parameter estimates via Bayesian updating, though the estimates from each batch (appear to) exhibit small finite-sample biases (which may be attributable to the presence of multiple local optima discussed earlier).\footnote{Unreported simulation results suggest that these biases diminish as the sample size per batch increases.}
Table (ref) reports the false labeling rates, defined as the proportion of incorrectly estimated $Z_i$ values across all $n$ agents: $\sum_{i=1}^{n} \mathbbm{1}\{\hat{Z}_i \neq Z_i\}/n$, where $\hat{Z}_i$ denotes the estimated latent type for agent $i$. Each estimate is assigned using the maximum a posteriori rule: $\hat{Z}_i = \arg\max_{k} \tilde{q}_i(k)$, where $\tilde{q}_i$ is the posterior distribution after batch $t$. From Table (ref), we observe that the false labeling rate declines steadily over successive batches, falling from nearly $22\%$ after the first batch to about $8\%$ the sixth. The steady decline in misclassification rates suggests that the HGT algorithm effectively accumulates information on each agent’s type over time through Bayesian updating based on batched LVAs, echoing the asymptotic insight in Section (ref) that latent type recovery improves as more data accumulate.
We now turn to network production performance---the primary economic outcome of interest. The left panel of Figure (ref) displays the total network productivity (i.e., expected success rate of establishing a friendship) achieved by our HGT algorithm across batches (blue), with error bars indicating $\pm 1$ standard deviation (black) based on 100 replications. For reference, the oracle benchmark (orange) represents the fully informed case where the production parameters $\boldsymbol{\theta}$ and the agents’ latent types $\boldsymbol{Z}$ are known. The results show steady improvement in network productivity across batches, with the average increasing from just above $90$ to near the oracle level by the end of the sixth batch. At the same time, variability across replications declines sharply, as the error bars shrink from nearly $15$ to around $3$, suggesting convergence toward oracle-level performance in most of the 100 simulation runs. The right panel shows regret, defined as the difference in expected network output between the HGT algorithm and the oracle benchmark. This plot essentially mirrors the productivity plot in reverse. We present both plots for the stationary environment and use regret as the primary performance metric in the non-stationary setting discussed next, where total network productivity is affected by changes in agent composition, making regret a more consistent metric for comparison.
We consider a pool of $n = 48$ agents (i.e., workers), each belonging to one of three latent types (e.g., weak-, strong-, or moderate-team players), drawn independently from a multinomial distribution with probabilities $(0.4, 0.3, 0.3)$. As in the previous setting, these latent types are fixed once generated in each simulation repetition and treated as unknown parameters. In each batch, these workers are paired into two-person teams to complete $960$ tasks, yielding an average of $40$ tasks per worker. We impose the agent-level “workload” constraint in ((ref)) with bounds $[d_l, d_r] = [35, 45]$.
The team-level success probability matrix $\boldsymbol{\theta}$ is calibrated to reflect productivity similar to those in xu2024heterogeneous, where agents with medium solo performance tend to be strong team players across all types. We preserve this feature while also aligning the matrix with the target of their average productivity level by setting the success probability for a two-worker team to be
with rows and columns ordered by type $(\text{weak}, \text{strong}, \text{moderate})$. The symmetry of $\boldsymbol{\theta}$ implies that the roles of the two workers in a team are interchangeable.
In this case, we also evaluate our algorithm in a non-stationary environment, where workers may leave and new ones may join over time---a common feature in high-turnover industries. In each batch, we randomly select $\xi$ of the $50$ workers, simulating them being replaced by new hires whose latent types are unknown and drawn from the same population distribution. For each newcomer, the associated $q_i$ values are reset to an uninformative categorical prior with equal probabilities, reflecting the absence of prior knowledge.
We begin by examining network production performance. Figure (ref) reports regret curves under four levels of worker turnover ($\xi = 0, 2, 4, 6$), with panels arranged from top left to bottom right. In each batch, the oracle benchmark is updated in each batch to reflect the current agent composition, and the plotted regret reflects the percentage shortfall in expected successes relative to the oracle. When there is no turnover ($\xi = 0$), the average regret decreases from about $25$ fewer successful tasks to approximately $2.5$ by batch $6$, relative to the oracle benchmark (which achieves around $290$ expected successes out of $960$ tasks). With $\xi = 2$ worker replacements per batch, regret still declines, albeit at a slower rate, and converges to a slightly higher level. This pattern becomes more pronounced at higher turnover levels ($\xi = 4$ and $6$), as shown in the lower panels. This is because the algorithm has no prior knowledge of new hires and must gradually relearn their latent types. Nevertheless, the HGT algorithm continues to yield substantial and consistent gains across all scenarios.
Table (ref) reports the false labeling rates, which decrease steadily across batches for each turnover level, indicating improved accuracy in type estimation over time. As turnover increases, the rate of improvement slows, since the algorithm must relearn the latent types of newly arrived workers.
For completeness, we also plot the empirical posterior distributions of each $\boldsymbol{\theta}$ element for the case of $\xi = 0$ in Figure (ref). The histograms confirm that our algorithm progressively sharpens parameter estimates through Bayesian updating. Results under other turnover levels are similar, as agent replacement does not affect the estimation of $\boldsymbol{\theta}$. We therefore omit those plots.
In this paper, we develop a novel adaptive, or online, network formation algorithm that pairs agents, represented as nodes, into edges to maximize total network output, measured by accumulated edge weights. We adopt a weighted stochastic block model, where each agent belongs to one of $K$ latent types that capture unobserved heterogeneity. Only the pairwise outcomes (edge weights), which depend on the types of the paired agents, are observed, while individual node-level contributions remain unobservable.
We formulate the problem of online network formation for improved production as a Batched Multi-Armed Bandit problem, where each type pair corresponds to an arm. In each batch, we assign agents into pairs (i.e., form the network) according to a pre-specified policy. Based on the observed outcomes, we employ a variational approximation to estimate the pairwise expected production parameter and the latent types of agents. These estimates, as justified by our asymptotic analysis, can be treated as Gaussian or categorical signals, and incorporated into updates of the underlying parameters. We then incorporate the updated posteriors into a maximum-weight matching problem to determine the pairing policy for the next batch. Iterating this procedure across batches leads to stochastically improved network-level production.
Our algorithm offers several key advantages. First, the WSBM framework captures heterogeneous complementarities in a nonparametric manner (i.e., without imposing a functional form on the $K \times K$ matrix of expected outcomes), providing greater flexibility than models like the strongest/weakest link, where outcomes depend solely on the better or worse worker. Second, the batched setting enables within-batch asymptotics, allowing us to invoke the multi-stage asymptotic representation theorem of hirano2023asymptotic, which offers a much simpler approximation for algorithm design. Third, the batch-wise, adaptively updated policy naturally accommodates nonstationary environments, where production parameters or agent types may evolve over time. Finally, by solving a maximum-weight matching problem, our algorithm can incorporate practical constraints (e.g., to balance workload) that commonly arise in organizational settings. Simulations demonstrate that our approach boosts productivity over a small number of batches while simultaneously improving estimates of both production parameters and latent types.