  Multilevel estimators aim at reducing the variance of Monte Carlo statistical
estimators, by combining samples generated with simulators of different costs
and accuracies. In particular, the recent work of Schaden and Ullmann (2020) on
the multilevel best linear unbiased estimator (MLBLUE) introduces a framework
unifying several multilevel and multifidelity techniques. The MLBLUE is
reintroduced here using a variance minimization approach rather than the
regression approach of Schaden and Ullmann. We then discuss possible extensions
of the scalar MLBLUE to a multidimensional setting, i.e. from the expectation
of scalar random variables to the expectation of random vectors. Several
estimators of increasing complexity are proposed: a) multilevel estimators with
scalar weights, b) with element-wise weights, c) with spectral weights and d)
with general matrix weights. The computational cost of each method is
discussed. We finally extend the MLBLUE to the estimation of second-order
moments in the multidimensional case, i.e. to the estimation of covariance
matrices. The multilevel estimators proposed are d) a multilevel estimator with
scalar weights and e) with element-wise weights. In large-dimension
applications such as data assimilation for geosciences, the latter estimator is
computationnally unaffordable. As a remedy, we also propose f) a multilevel
covariance matrix estimator with optimal multilevel localization, inspired by
the optimal localization theory of M\'en\'etrier and Aulign\'e (2015). Some
practical details on weighted MLMC estimators of covariance matrices are given
in appendix.
