We present an efficient separation of variables algorithm for the evaluation of imaginary time Feynman diagrams appearing in the bold pseudo-particle strong coupling expansion of the Anderson impurity model. The algorithm uses a fitting method based on AAA rational approximation and numerical optimization to obtain a sum-of-exponentials expansion of the hybridization function, which is then used to decompose the diagrams. A diagrammatic formulation of the algorithm leads to an automated procedure for diagrams of arbitrary order and topology. We also present methods of stabilizing the self-consistent solution of the pseudo-particle Dyson equation. The result is a low-cost and high-order accurate impurity solver for quantum embedding methods using general multi-orbital hybridization functions at low temperatures, appropriate for low-tointermediate expansion orders. In addition to other benchmark examples, we use our solver to perform a dynamical mean-field theory study of a minimal model of the strongly correlated compound Ca2RuO4, describing the anti-ferromagnetic transition and the inand out-of-plane anisotropy induced by spin-orbit coupling.