This work presents a novel realization approach to Quantum Boltzmann Machines (QBMs). The preparation of the required Gibbs states, as well as the evaluation of the loss function’s analytic gradient is based on Variational Quantum Imaginary Time Evolution, a technique that is typically used for ground state computation. In contrast to existing methods, this implementation facilitates near-term compatible QBM training with gradients of the actual loss function for arbitrary parameterized Hamiltonians which do not necessarily have to be fully-visible but may also include hidden units. The variational Gibbs state approximation is demonstrated with numerical simulations and experiments run on real quantum hardware provided by IBM Quantum. Furthermore, we illustrate the application of this variational QBM approach to generative and discriminative learning tasks using numerical simulation. The evaluation of the first- as well as higher-order quantum gradients and that are needed to compute the loss function’s analytic gradient are intregrated with Qiskit’s gradient framework.