  The Bethe-Salpeter equation (BSE) is currently the state of the art in the
description of neutral electron excitations in both solids and large finite
systems. It is capable of accurately treating charge-transfer excitations that
present difficulties for simpler approaches. We present a local basis set
formulation of the BSE for molecules where the optical spectrum is computed
with the iterative Haydock recursion scheme, leading to a low computational
complexity and memory footprint. Using a variant of the algorithm we can go
beyond the Tamm-Dancoff approximation (TDA). We rederive the recursion
relations for general matrix elements of a resolvent, show how they translate
into continued fractions, and study the convergence of the method with the
number of recursion coefficients and the role of different terminators. Due to
the locality of the basis functions the computational cost of each iteration
scales asymptotically as $O(N^3)$ with the number of atoms, while the number of
iterations is typically much lower than the size of the underlying
electron-hole basis. In practice we see that , even for systems with thousands
of orbitals, the runtime will be dominated by the $O(N^2)$ operation of
applying the Coulomb kernel in the atomic orbital representation
