  Network reconstruction consists in determining the unobserved pairwise
couplings between $N$ nodes given only observational data on the resulting
behavior that is conditioned on those couplings -- typically a time-series or
independent samples from a graphical model. A major obstacle to the scalability
of algorithms proposed for this problem is a seemingly unavoidable quadratic
complexity of $\Omega(N^2)$, corresponding to the requirement of each possible
pairwise coupling being contemplated at least once, despite the fact that most
networks of interest are sparse, with a number of non-zero couplings that is
only $O(N)$. Here we present a general algorithm applicable to a broad range of
reconstruction problems that significantly outperforms this quadratic baseline.
Our algorithm relies on a stochastic second neighbor search (Dong et al., 2011)
that produces the best edge candidates with high probability, thus bypassing an
exhaustive quadratic search. If we rely on the conjecture that the
second-neighbor search finishes in log-linear time (Baron & Darling, 2020;
2022), we demonstrate theoretically that our algorithm finishes in subquadratic
time, with a data-dependent complexity loosely upper bounded by $O(N^{3/2}\log
N)$, but with a more typical log-linear complexity of $O(N\log^2N)$. In
practice, we show that our algorithm achieves a performance that is many orders
of magnitude faster than the quadratic baseline -- in a manner consistent with
our theoretical analysis -- allows for easy parallelization, and thus enables
the reconstruction of networks with hundreds of thousands and even millions of
nodes and edges.
