  In this paper we introduce a multilevel Picard approximation algorithm for
semilinear parabolic partial integro-differential equations (PIDEs). We prove
that the numerical approximation scheme converges to the unique viscosity
solution of the PIDE under consideration. To that end, we derive a Feynman-Kac
representation for the unique viscosity solution of the semilinear PIDE,
extending the classical Feynman-Kac representation for linear PIDEs.
Furthermore, we show that the algorithm does not suffer from the curse of
dimensionality, i.e. the computational complexity of the algorithm is bounded
polynomially in the dimension $d$ and the reciprocal of the prescribed accuracy
$\varepsilon$. We also provide a numerical example in up to 10'000 dimensions
to demonstrate its applicability.
