  We propose a control variate multilevel Monte Carlo method for the kinetic
BGK model of the Boltzmann equation subject to random inputs. The method
combines a multilevel Monte Carlo technique with the computation of the optimal
control variate multipliers derived from local or global variance minimization
problems. Consistency and convergence analysis for the method equipped with a
second-order positivity-preserving and asymptotic-preserving scheme in space
and time is also performed. Various numerical examples confirm that the
optimized multilevel Monte Carlo method outperforms the classical multilevel
Monte Carlo method especially for problems with discontinuities.
