  We derive a stochastic hierarchy of matrix product states (HOMPS) for
non-Markovian dynamics in open quantum system at finite temperature, which is
numerically exact and efficient. HOMPS is obtained from the recently developed
stochastic hierarchy of pure states (HOPS) by expressing HOPS in terms of
formal creation and annihilation operators. The resulting stochastic first
order differential equation is then formulated in terms of matrix product
states and matrix product operators. In this way the exponential complexity of
HOPS can be reduced to scale polynomial with the number of particles. The
validity and efficiency of HOMPS is demonstrated for the spin-boson model and
long chains where each site is coupled to a structured, strongly non-Markovian
environment.
