The next generation of experiments measuring the anomalous magnetic moment of the electron are expected to exceed the precision of 1 ppb (1 part per billion); to match it on the theoretical side, the four loop QED contribution to the anomaly must be evaluated with a precision of a few parts per mill, i.e. a factor ten better than the latest numerical value. In this paper we present an approach to the problem which relies on the systematic exploitation of the integration by parts identities among Feynman amplitudes for reducing the very many integrals appearing in the calculation to a much smaller number of master integrals, discussing the algebraic algorithms used for solving the integration by part identities and the implementation of the algorithm in the computer algebra program FORM.