  The approximative calculation of iterated nested expectations is a recurring
challenging problem in applications. Nested expectations appear, for example,
in the numerical approximation of solutions of backward stochastic differential
equations (BSDEs), in the numerical approximation of solutions of semilinear
parabolic partial differential equations (PDEs), in statistical physics, in
optimal stopping problems such as the approximative pricing of American or
Bermudan options, in risk measure estimation in mathematical finance, or in
decision-making under uncertainty. Nested expectations which arise in the above
named applications often consist of a large number of nestings. However, the
computational effort of standard nested Monte Carlo approximations for iterated
nested expectations grows exponentially in the number of nestings and it
remained an open question whether it is possible to approximately calculate
multiply iterated high-dimensional nested expectations in polynomial time. In
this article we tackle this problem by proposing and studying a new class of
full-history recursive multilevel Picard (MLP) approximation schemes for
iterated nested expectations. In particular, we prove under suitable
assumptions that these MLP approximation schemes can approximately calculate
multiply iterated nested expectations with a computational effort growing at
most polynomially in the number of nestings $ K \in \mathbb{N} = \{1, 2, 3,
\ldots \} $, in the problem dimension $ d \in \mathbb{N} $, and in the
reciprocal $\frac{1}{\varepsilon}$ of the desired approximation accuracy $
\varepsilon \in (0, \infty) $.
