  The contraction of tensor networks is a central task in the application of
tensor network methods to the study of quantum and classical many body systems.
In this paper, we investigate the impact of gauge degrees of freedom in the
virtual indices of the tensor network on the contraction process, specifically
focusing on boundary matrix product state methods for contracting
two-dimensional tensor networks. We show that the gauge transformation can
affect the entanglement structures of the eigenstates of the transfer matrix
and change how the physical information is encoded in the eigenstates, which
can influence the accuracy of the numerical simulation. We demonstrate this
effect by looking at two different examples. First, we focus on the local gauge
transformation, and analyze its effect by viewing it as an imaginary-time
evolution governed by a diagonal Hamiltonian. As a specific example, we perform
a numerical analysis in the classical Ising model on the square lattice.
Second, we go beyond the scope of local gauge transformations and study the
antiferromagnetic Ising model on the triangular lattice. The partition function
of this model has two tensor network representations connected by a non-local
gauge transformation, resulting in distinct numerical performances in the
boundary MPS calculation.
