  We develop a multilevel approach to compute approximate solutions to backward
differential equations (BSDEs). The fully implementable algorithm of our
multilevel scheme constructs sequential martingale control variates along a
sequence of refining time-grids to reduce statistical approximation errors in
an adaptive and generic way. We provide an error analysis with explicit and
non-asymptotic error estimates for the multilevel scheme under general
conditions on the forward process and the BSDE data. It is shown that the
multilevel approach can reduce the computational complexity to achieve
precision $\epsilon$, ensured by error estimates, essentially by one order (in
$\epsilon^{-1}$) in comparison to established methods, which is substantial.
Computational examples support the validity of the theoretical analysis,
demonstrating efficiency improvements in practice.
