  Many properties of current \emph{ab initio} approaches to the quantum
many-body problem, both perturbational or otherwise, are related to the
singularity structure of Rayleigh--Schr\"odinger perturbation theory. A
numerical procedure is presented that in principle computes the complete set of
singularities, including the dominant singularity which limits the radius of
convergence. The method approximates the singularities as eigenvalues of a
certain generalized eigenvalue equation which is solved using iterative
techniques. It relies on computation of the action of the perturbed Hamiltonian
on a vector, and does not rely on the terms in the perturbation series. Some
illustrative model problems are studied, including a Helium-like model with
$\delta$-function interactions for which M{\o}ller--Plesset perturbation theory
is considered and the radius of convergence found.
