  We develop a multilevel Monte Carlo (MLMC)-FEM algorithm for linear, elliptic
diffusion problems in polytopal domain $\mathcal D\subset \mathbb R^d$, with
Besov-tree random coefficients. This is to say that the logarithms of the
diffusion coefficients are sampled from so-called Besov-tree priors, which have
recently been proposed to model data for fractal phenomena in science and
engineering. Numerical analysis of the fully discrete FEM for the elliptic PDE
includes quadrature approximation and must account for a) nonuniform pathwise
upper and lower coefficient bounds, and for b) low path-regularity of the
Besov-tree coefficients. Admissible non-parametric random coefficients
correspond to random functions exhibiting singularities on random fractals with
tunable fractal dimension, but involve no a-priori specification of the fractal
geometry of singular supports of sample paths. Optimal complexity and
convergence rate estimates for quantities of interest and for their second
moments are proved. A convergence analysis for MLMC-FEM is performed which
yields choices of the algorithmic steering parameters for efficient
implementation. A complexity (``error vs work'') analysis of the MLMC-FEM
approximations is provided.
