We present a new strategy for contracting tensor networks in arbitrary geometries. This method is designed to follow as strictly as possible the renormalization group philosophy, by first contracting tensors in an exact way and, then, performing a controlled truncation of the resulting tensor. We benchmark this approximation procedure in two dimensions against an exact contraction. We then apply the same idea to a three dimensional system. The underlying rational for emphasizing the exact coarse graining renormalization group step prior to truncation is related to monogamy of entanglement.PACS numbers: 03.67. Mn, 05.10.Cc, 05.50.+q Computing analytically the properties of a quantum model is, in general, not possible. It is then necessary to resort to classical simulations that rely on approximation techniques. From a quantum information perspective, the presence of a small amount of entanglement in the system has been identified as the key ingredient that allows for efficient classical simulations. This has been exploited in one dimension (1D) in a series of algorithms based on a description of the system as a tensor network. Following this approach, condensed matter systems can be simulated using Matrix Product States (MPS, [1,2]) as the testbed for the variational procedure of the density-matrix renormalization group (DMRG,[3]) to study ground states of local 1D Hamiltonians. Beyond the 1D setting, the computation of physical magnitudes using tensor networks is limited by the numerical effort necessary to perform the contraction of the tensor network, i.e. to sum all its indices. The efficiency for performing this task is limited by the area law scaling of the entanglement entropy in the system [4][5][6][7]. To overcome this problem, several strategies aim at finding the best possible approximation to the contraction of tensor networks after identifying the relevant degrees of freedom of the system [8][9][10].Let us briefly recall the key elements of the tensor network representation. Given a quantum state of N particles |ψ = c i1,...,iN |i 1 . . . i N its coefficients can be represented as a contraction of local tensors c i1,...,iN = tr(A 1,i1 . . . A N,iN ), where each local tensor A j carries a physical index i j and ancillary indices (which are not written) that get contracted according to a prescribed geometry. The rank of these ancillary indices, that we shall call χ, controls the amount of entanglement which is captured by the tensor representation. If the tensors are simple matrices on a line, the tensor network is called MPS [1,2]. Other possible geometries are regular squared grids in any dimensions that correspond to PEPS [11], and tree-like structures that go under the name of TTN [12] and MERA [13].In this letter we propose a new strategy to contract tensor networks in general geometries, that we shall illustrate in detail for PEPS in 2D and 3D. The method is based on following as strictly as possible the renormalization group (RG) philosophy. First, an exact contraction of a set of local tenso...