We introduce a tensor renormalization group scheme for coarse graining a two-dimensional tensor network that can be successfully applied to both classical and quantum systems on and off criticality. The key innovation in our scheme is to deform a 2D tensor network into small loops and then optimize the tensors on each loop. In this way, we remove short-range entanglement at each iteration step and significantly improve the accuracy and stability of the renormalization flow. We demonstrate our algorithm in the classical Ising model and a frustrated 2D quantum model. DOI: 10.1103/PhysRevLett.118.110504 Introduction.-In recent years, the tensor network (TN) approach [1,2] has become a powerful theoretical and computational [8, tool for studying condensed matter systems. Many physical quantities, including the partition function of a classical system, the Euclidean path integral of a quantum system, and the expectation value of physical observables, can be expressed in terms of tensor networks. Evaluating these quantities is reduced to the contraction of a multidimensional tensor network. In the two dimensional case, many algorithms [8,32,[37][38][39][40][41]43,[45][46][47][48][49][50][53][54][55][56][57] have been developed to implement the approximate tensor contractions. Among these, the tensor renormalization group approach introduced by Levin and Nave [38] and its generalizations [8,22,39,[43][44][45][46][47]55,56,61] have unique features: the tensor contraction is based on a fully isotropic coarse-graining procedure. Moreover, when applying the method to a system on a finite torus, the computational cost is lower than those based on matrix product states (MPS) [32,37,41,[48][49][50]53,54].However, the Levin-Nave tensor network renormalization (TRG, also referred as LN-TNR here) [38] is based on the singular value decomposition (SVD) of local tensors, which only minimizes the truncation errors of tree tensor networks. Several improvements [45][46][47] have taken into account the effect of the environments, but they are still essentially based on tree tensor networks. These approaches cannot completely remove short-range entanglements during the coarse graining process. For example, in the 2D TN calculation of a partition function (or a path integral) TNR based on simple SVD cannot simplify the corner-double-line (CDL) tensor [38], despite the CDL tensor describing a product state that should be simplified to a one-dimensional tensor. In Ref.[8], this issue was seriously discussed. The authors pointed out that to further remove short-range entanglement, it is crucial to optimize the tensor configurations that contain a loop. However, due to the computational cost, only a crude iterative method is