We give a quasi-polynomial time classical algorithm for estimating the ground state energy and for computing low energy states of quantum impurity models. Such models describe a bath of free fermions coupled to a small interacting subsystem called an impurity. The full system consists of n fermionic modes and has a Hamiltonian H= H0+ Himp, where H0 is quadratic in creation–annihilation operators and Himp is an arbitrary Hamiltonian acting on a subset of O(1) modes. We show that the ground energy of H can be approximated with an additive error 2 -b in time n3exp [ O(b3) ]. Our algorithm also finds a low energy state that achieves this approximation. The low energy state is represented as a superposition of exp [ O(b3) ] fermionic Gaussian states. To arrive at this result we prove several theorems concerning exact ground states of impurity models. In particular, we show that eigenvalues of the ground state covariance matrix decay exponentially with the exponent depending very mildly on the spectral gap of H0. A key ingredient of our proof is Zolotarev’s rational approximation to the x function. We anticipate that our algorithms may be used in hybrid quantum-classical simulations of strongly correlated materials based on dynamical mean field theory. We implemented a simplified practical version of our algorithm and benchmarked it using the single impurity Anderson model.