Exponential densities on unitary conjugation orbits of Hermitian matrices, known as Harish-Chandra-Itzykson-Zuber (HCIZ) densities, are important in various settings in physics and random matrix theory. However, the basic question of efficient sampling from these densities has remained open. We present two efficient algorithms to sample matrices from distributions that are close to the HCIZ distribution. Both algorithms exploit a natural self-reducible structure that arises from continuous symmetries of the underlying unitary orbit. We will also discuss applications to quantum inference and differentially private rank-k approximation.