  Cosmological probes pose an inverse problem where the measurement result is
obtained through observations, and the objective is to infer values of model
parameters which characterize the underlying physical system -- our Universe.
Modern cosmological probes increasingly rely on measurements of the small-scale
structure, and the only way to accurately model physical behavior on those
scales, roughly 65 Mpc/h or smaller, is via expensive numerical simulations. In
this paper, we provide a detailed description of a novel statistical framework
for obtaining accurate parameter constraints by combining observations with a
very limited number of cosmological simulations. The proposed framework
utilizes multi-output Gaussian process emulators that are adaptively
constructed using Bayesian optimization methods. We compare several approaches
for constructing multi-output emulators that enable us to take possible
inter-output correlations into account while maintaining the efficiency needed
for inference. Using Lyman alpha forest flux power spectrum, we demonstrate
that our adaptive approach requires considerably fewer --- by a factor of a few
in Lyman alpha P(k) case considered here --- simulations compared to the
emulation based on Latin hypercube sampling, and that the method is more robust
in reconstructing parameters and their Bayesian credible intervals.
