  The Bethe-Salpeter equation plays a crucial role in understanding the physics
of correlated fermions, relating to optical excitations in solids as well as
resonances in high-energy physics. Yet, it is notoriously difficult to control
numerically, typically requiring an effort that scales polynomially with energy
scales and accuracy. This puts many interesting systems out of computational
reach. Using the intermediate representation and sparse modelling for
two-particle objects on the Matsubara axis, we develop an algorithm that solves
the Bethe-Salpeter equation in $O(L^8)$ time with $O(L^4)$ memory, where $L$
grows only logarithmically with inverse temperature, bandwidth, and desired
accuracy, This opens the door for computations in hitherto inaccessible
regimes. We benchmark the method on the Hubbard atom and on the multi-orbital
weak-coupling limit, where we observe the expected exponential convergence to
the analytical results. We then showcase the method for a realistic impurity
problem.
