Dyadic data are ubiquitous in our society. International trade, travels, population flows, military alliances, partnerships between firms, research collaboration, and many others can be represented as dyadic data, where each dyad represents a pair of countries, firms, or individuals, depending on the context. Dyadic data analysis is particularly prevalent in the literature of international trade, where regression-based analysis, the so-called gravity model, serves as a primary analytical approach in these fields since the pioneering work by Tinbergen (1962) (see also, e.g., Anderson, 1979, 2011; Head and Mayer, 2014 and references therein). For reviews of recent econometric literature on dyadic data analysis in general, see, for example, Graham (2020a,b). Despite the popularity of dyadic data, there are only a few causal inference methods tailored specifically for dyadic data analysis, with some exceptions such as Baier and Bergstrand (2009), Arpino et al. (2017), and Nagengast and Yotov (2023). This may be due to the non-standard and complex endogeneity structure often encountered in typical applications of dyadic data. For example, suppose we are interested in the impacts of free trade agreements (FTA) on trade flows between countries. The treatment variable, FTA, should be considered endogenous because both the decision to enter into FTA and the trade outcome should be influenced by each country's economic factors and the economic and political relationship between the countries involved. Thus, if one tries to resolve the endogeneity issue by using the instrumental variables (IV) method, for instance, then he/she needs to prepare at least three different types of IVs: those accounting for confounding factors at the "origin" country, those at the "destination", and pair-specific factors.

This paper considers the task of linear regression with shuffled labels, i.e., $\mathbf Y = \mathbf \Pi \mathbf X \mathbf B + \mathbf W$, where $\mathbf Y \in \mathbb R^{n\times m}, \mathbf Pi \in \mathbb R^{n\times n}, \mathbf X\in \mathbb R^{n\times p}, \mathbf B \in \mathbb R^{p\times m}$, and $\mathbf W\in \mathbb R^{n\times m}$, respectively, represent the sensing results, (unknown or missing) corresponding information, sensing matrix, signal of interest, and additive sensing noise. Given the observation $\mathbf Y$ and sensing matrix $\mathbf X$, we propose a one-step estimator to reconstruct $(\mathbf \Pi, \mathbf B)$. From the computational perspective, our estimator's complexity is $O(n^3 + np^2m)$, which is no greater than the maximum complexity of a linear assignment algorithm (e.g., $O(n^3)$) and a least square algorithm (e.g., $O(np^2 m)$). From the statistical perspective, we divide the minimum $snr$ requirement into four regimes, e.g., unknown, hard, medium, and easy regimes; and present sufficient conditions for the correct permutation recovery under each regime: $(i)$ $snr \geq \Omega(1)$ in the easy regime; $(ii)$ $snr \geq \Omega(\log n)$ in the medium regime; and $(iii)$ $snr \geq \Omega((\log n)^{c_0}\cdot n^{{c_1}/{srank(\mathbf B)}})$ in the hard regime ($c_0, c_1$ are some positive constants and $srank(\mathbf B)$ denotes the stable rank of $\mathbf B$). In the end, we also provide numerical experiments to confirm the above claims.