We describe a formulation of multi-reference perturbation theory that obtains a rigorous upper bound to the second order energy by minimizing the Hylleraas functional in the space of matrix product states (MPS). The first order wavefunctions so obtained can also be used to compute the third order energy with little overhead. Our formulation has several advantages including (i) flexibility with respect to the choice of zeroth order Hamiltonian, (ii) recovery of the exact uncontracted multi-reference perturbation theory energies in the limit of large MPS bond dimension, (iii) no requirement to compute high body density matrices, (iv) an embarrassingly parallel algorithm (scaling up to the number of virtual orbitals, squared, processors). Preliminary numerical examples show that the MPS bond dimension required for accurate first order wavefunctions scales sub-linearly with the size of the basis.