  Quantum dynamical simulations of statistical ensembles pose a significant
computational challenge due to the fact that mixed states need to be
represented. If the underlying dynamics is fully unitary, for example in
ultrafast coherent control at finite temperatures, one approach to approximate
time-dependent observables is to sample the density operator by solving the
Schr\"{o}dinger equation for a set of wave functions with randomized phases. We
show that, on average, random-phase wave functions perform well for ensembles
with high mixedness, whereas at higher purities a deterministic sampling of the
energetically lowest-lying eigenstates becomes superior. We prove that
minimization of the worst-case error for computing arbitrary observables is
uniquely attained by eigenstate-based sampling. We show that this error can be
used to form a qualitative estimate of the set of ensemble purities for which
the sampling performance of the eigenstate-based approach is superior to
random-phase wave functions. Furthermore we present refinements to both schemes
which remove redundant information from the sampling procedure to accelerate
their convergence. Finally, we point out how the structure of low-rank
observables can be exploited to further improve eigenstate-based sampling
schemes.
