We describe a numerical algorithm for approximating the equilibrium-reduced
density matrix and the effective (mean force) Hamiltonian for a set of system
spins coupled strongly to a set of bath spins when the total system
(system+bath) is held in canonical thermal equilibrium by weak coupling with a
"super-bath". Our approach is a generalization of now standard typicality
algorithms for computing the quantum expectation value of observables of bare
quantum systems via trace estimators and Krylov subspace methods. In
particular, our algorithm makes use of the fact that the reduced system
density, when the bath is measured in a given random state, tends to
concentrate about the corresponding thermodynamic averaged reduced system
density. Theoretical error analysis and numerical experiments are given to
validate the accuracy of our algorithm. Further numerical experiments
demonstrate the potential of our approach for applications including the study
of quantum phase transitions and entanglement entropy for long-range
interaction systems