  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.
