  We present a deterministic algorithm for the efficient evaluation of
imaginary time diagrams based on the recently introduced discrete Lehmann
representation (DLR) of imaginary time Green's functions. In addition to the
efficient discretization of diagrammatic integrals afforded by its
approximation properties, the DLR basis is separable in imaginary time,
allowing us to decompose diagrams into linear combinations of nested sequences
of one-dimensional products and convolutions. Focusing on the strong coupling
bold-line expansion of generalized Anderson impurity models, we show that our
strategy reduces the computational complexity of evaluating an $M$th-order
diagram at inverse temperature $\beta$ and spectral width $\omega_{\max}$ from
$\mathcal{O}((\beta \omega_{\max})^{2M-1})$ for a direct quadrature to
$\mathcal{O}(M (\log (\beta \omega_{\max}))^{M+1})$, with controllable
high-order accuracy. We benchmark our algorithm using third-order expansions
for multi-band impurity problems with off-diagonal hybridization and spin-orbit
coupling, presenting comparisons with exact diagonalization and quantum Monte
Carlo approaches. In particular, we perform a self-consistent dynamical
mean-field theory calculation for a three-band Hubbard model with strong
spin-orbit coupling representing a minimal model of Ca$_2$RuO$_4$,
demonstrating the promise of the method for modeling realistic strongly
correlated multi-band materials. For both strong and weak coupling expansions
of low and intermediate order, in which diagrams can be enumerated, our method
provides an efficient, straightforward, and robust black-box evaluation
procedure. In this sense, it fills a gap between diagrammatic approximations of
the lowest order, which are simple and inexpensive but inaccurate, and those
based on Monte Carlo sampling of high-order diagrams.
