Variational tensor network optimization has become a powerful tool for studying classical statistical models in two dimensions. However, its application to three-dimensional systems remains limited, primarily due to the high computational cost associated with evaluating the free energy density and its gradient. This process requires contracting a triple-layer tensor network composed of a projected entangled pair operator and projected entangled pair states. In this paper, we employ a split corner-transfer renormalization group scheme tailored for the contraction of such triple-layer networks, which reduces the computational complexity from O (D^12) in conventional methods to O (D^9). Through numerical benchmarks on the three-dimensional classical Ising model, we demonstrate that the proposed scheme achieves numerical results comparable to the most recent Monte Carlo simulations, while also providing a substantial speedup over previous variational tensor network approaches. This makes this method well-suited for efficient gradient-based optimization in three-dimensional tensor network simulations.
Xu et al. (Tue,) studied this question.