  We introduce a new family of numerical algorithms for approximating solutions
of general high-dimensional semilinear parabolic partial differential equations
at single space-time points. The algorithm is obtained through a delicate
combination of the Feynman-Kac and the Bismut-Elworthy-Li formulas, and an
approximate decomposition of the Picard fixed-point iteration with multilevel
accuracy. The algorithm has been tested on a variety of semilinear partial
differential equations that arise in physics and finance, with very
satisfactory results. Analytical tools needed for the analysis of such
algorithms, including a semilinear Feynman-Kac formula, a new class of
semi-norms and their recursive inequalities, are also introduced. They allow us
to prove for semilinear heat equations with gradient-independent nonlinearity
that the computational complexity of the proposed algorithm is bounded by
$O(d\,\varepsilon^{-(4+\delta)})$ for any $\delta \in (0,\infty)$ under
suitable assumptions, where $d\in \mathbb{N}$ is the dimensionality of the
problem and $\varepsilon\in(0,\infty)$ is the prescribed accuracy.
