  We present a quantum Monte Carlo method capable of sampling the full density
matrix of a many-particle system at finite temperature. This allows arbitrary
reduced density matrix elements and expectation values of complicated non-local
observables to be evaluated easily. The method resembles full configuration
interaction quantum Monte Carlo but works in the space of many-particle
operators instead of the space of many-particle wave functions. One simulation
provides the density matrix at all temperatures simultaneously, from $T=\infty$
to $T=0$, allowing the temperature dependence of expectation values to be
studied. The direct sampling of the density matrix also allows the calculation
of some previously inaccessible entanglement measures. We explain the theory
underlying the method, describe the algorithm, and introduce an
importance-sampling procedure to improve the stochastic efficiency. To
demonstrate the potential of our approach, the energy and staggered
magnetization of the isotropic antiferromagnetic Heisenberg model on small
lattices, the concurrence of one-dimensional spin rings, and the Renyi $S_2$
entanglement entropy of various sublattices of the $6 \times 6$ Heisenberg
model are calculated. The nature of the sign problem in the method is also
investigated.
