We derive a hierarchy of stochastic evolution equations for pure states (quantum trajectories) to efficiently solve open quantum system dynamics with non-Markovian structured environments. From this hierarchy of pure states (HOPS) the exact reduced density operator is obtained as an ensemble average. We demonstrate the power of HOPS by applying it to the Spin-Boson model, the calculation of absorption spectra of molecular aggregates and energy transfer in a photosynthetic pigment-protein complex.The treatment of the dynamics of realistic open quantum systems still poses both conceptual and computational challenges. These arise from non-Markovian behavior due to a structured environment or strong systemenvironment interaction [1,2]. Severe assumptions, like weak-coupling or Markov approximation, are often made for practical reasons. However, they fail for many systems of interest. In these situations one relies on computationally demanding numerical methods. Among these are path integral approaches [3,4] or hierarchical equations of motion [5,6] for the system's reduced density matrix.In this Letter we follow a different strategy and derive a hierarchy of stochastic differential equations for pure states in the system Hilbert space (quantum trajectories). From this hierarchy of pure states (HOPS) the exact reduced density operator is obtained as an ensemble average. Our approach is based upon non-Markovian Quantum State diffusion (NMQSD), derived in its general form in Refs. [7][8][9][10]. NMQSD has been applied to various physical problems including the description of energy transfer in photosynthesis [11,12]. On a more fundamental side, NMQSD has been studied in the context of continuous measurement theory [13,14] and spontaneous wavefunction localization [15]. Other stochastic approaches, with various levels of applicability have been suggested [16][17][18].Although the NMQSD approach is formally exact, it seemed numerically difficult to handle, because of the appearance of a functional derivative with respect to a stochastic process. Only a few exactly solvable models are known (see e.g. [19][20][21][22]). In previous works we have replaced that functional derivative by an operator ansatz and dealt with it in the so called ZOFE approximation [11,[23][24][25] that allows for a very efficient numerical solution and agrees remarkably well with established results for a large number of problems. However, in certain cases this method is known to fail (see e.g. [11,25]). In Ref.[26] a hierarchical approach is applied to the operator ansatz of the functional derivative. Our new HOPS presented here is not based on the previously assumed ansatz, it is numerically exact, converges rapidly and offers a systematic way to check for convergence by increasing the number of equations taken into account. In addition, it offers the advantages of stochastic Schrödinger equations, e.g. one deals with pure states (and not large density matrices) and the calculation of independent realizations can be parallelized trivially.In the ...