One fundamental limitation of quantum chemical methods is the accuracy of the approximate many-body theoretical framework that is utilized. Accurate many-body formalisms for quantum chemical methods do exist, but these methods are computationally very expensive . Methods also exist that are much less computationally expensive such as Hatree-Fock, Density Functional and the Hybrid Functional theories [2-4], but at a reduced representation of the exact many-body ground state. This severely limits either the system size that can be addressed accurately, or the accuracy of the representation of the many-body ground state. What is essential is a method which represents the many-body ground state accurately, but with a low computational cost. Recently, a method for determining the response, to any order of the perturbation, within the density matrix formalism has been discovered [5-7]. This method is very simple and computationally efficient, and it immediately opens up the possibility of computing the variational many-body ground state to unprecedented accuracy within a simplified computational approach. Within this article, we report on the theoretical development of this methodology, which we refer to as Many Body Density Matrix Perturbation Theory. This theory has many significant advantages over existing methods. One, its computational cost is equivalent to Hartree-Fock or Density Functional theory. Two it is a variational upper bound to the exact many-body ground state energy. Three, like Hartree-Fock, it has no self-interaction. And four, because this theory obeys the Hellman-Feyman theorem, force are easily computable.