The approximate contraction of a tensor network of projected entangled pair states (PEPS) is a fundamental ingredient of any PEPS algorithm, required for the optimization of the tensors in ground state search or time evolution, as well as for the evaluation of expectation values. An exact contraction is in general impossible, and the choice of the approximating procedure determines the efficiency and accuracy of the algorithm. We analyze different previous proposals for this approximation, and show that they can be understood via the form of their environment, i.e. the operator that results from contracting part of the network. This provides physical insight into the limitation of various approaches, and allows us to introduce a new strategy, based on the idea of clusters, that unifies previous methods. The resulting contraction algorithm interpolates naturally between the cheapest and most imprecise and the most costly and most precise method. We benchmark the different algorithms with finite PEPS, and show how the cluster strategy can be used for both the tensor optimization and the calculation of expectation values. Additionally, we discuss its applicability to the parallelization of PEPS and to infinite systems.In the past few years, tensor network states (TNS) have been revealed as a very promising choice for the numerical simulation of strongly correlated quantum many-body systems. A sustained effort has led to significant conceptual and technical advancement of these methods, e.g. in .In the case of one-dimensional systems, matrix product states (MPS) are the variational class of TNS underlying the density matrix renormalization group (DMRG) [23]. Insight gained from quantum information theory has allowed the understanding of DMRGʼs enormous success at approximating ground states of spin chains, and the extension of the technique to dynamical problems [1-6] and lattices of more complex geometry [7][8][9].Projected entangled pair states (PEPS) [7] constitute a family of TNS that naturally generalizes MPS to spatial dimensions larger than 1 and arbitrary lattice geometry. Like MPS, PEPS incorporate the area law by construction, which makes them a very promising variational ansatz for strongly correlated systems which might not be tractable by other means, e.g. frustrated or fermionic states where quantum Monte Carlo methods suffer from the sign problem. Although originally defined for spin systems, PEPS have been subsequently formulated for fermions [24][25][26][27], and their potential in the numerical simulation of fermionic phases has been demonstrated [27][28][29][30]. But in contrast to the case for MPS, even local expectation values cannot be computed exactly in the case of PEPS. This is because the exact evaluation of the TN that represents the observables has an exponential cost in the system size. The same difficulty affects the contraction of the TN that surrounds a given tensor, the so-called environment, required for the local update operations in the course of optimization algorithms. It is neve...