The $\Lambda$CDM model for the Universe is highly successful in explaining cosmological observations to date, and its parameters tightly constrained by Cosmic Microwave Background (CMB) experiments such as Planck. Higher-order statistics, like the three-point correlation function or bispectrum in Fourier space, will be indispensable for furthering our understanding of the Universe. While these methodologies have been developed over the years and applied to CMB analyses, similar work on large-scale structure is still in its infancy. Additionally, information from future galaxy surveys such as LSST and Euclid will soon exceed that available from the CMB, demonstrating a pressing need for such tools.

The theoretical modelling of non-linear gravitational interactions is difficult beyond the perturbative regime, necessitating large, expensive $N$-body dark matter simulations to understand the small-scale dynamics. Additionally, the direct numerical computation of the matter bispectrum is intractable due to the multiplicity of triangular configurations. In this Thesis, we make breakthroughs in both of these problems. First, we present the newly rewritten MODAL-LSS formalism that enables efficient and optimal estimation of the full bispectrum for any matter density field to unprecedented accuracy, as well as demonstrating rapid convergence which makes it ideal for the analysis of large datasets. This has allowed us to benchmark fast dark matter codes (e.g. particle-mesh or L-PICOLA) against GADGET-3 using the bispectrum, showing quantitatively how the mismatch at large $k$ can be improved with a simple boosting technique in the power spectrum. We have also estimated the non-Gaussian contribution to the dark matter bispectrum covariance, which cannot be computed analytically in the non-linear regime. This will be vital for the extraction of cosmological parameters from data in the future.

In preparation for the analysis of future galaxy datasets we have also investigated the non-trivial problem of linking the underlying dark matter density field to the observed galaxy distribution. As an important milestone we have investigated the effects of the halo profile, the Halo Occupation Distribution (HOD) model, and multivariate assembly bias models of the halo occupation and concentration on the power spectrum and full bispectrum of a subhalo catalogue derived from the ROCKSTAR halo finder. These fast, phenomenological methods allow us to pave the way for the efficient generation of mock galaxy catalogues.