  We propose a tensor network encoding the set of all eigenstates of a fully
many-body localized system in one dimension. Our construction, conceptually
based on the ansatz introduced in Phys. Rev. B 94, 041116(R) (2016), is built
from two layers of unitary matrices which act on blocks of $\ell$ contiguous
sites.
  We argue this yields an exponential reduction in computational time and
memory requirement as compared to all previous approaches for finding a
representation of the complete eigenspectrum of large many-body localized
systems with a given accuracy. Concretely, we optimize the unitaries by
minimizing the magnitude of the commutator of the approximate integrals of
motion and the Hamiltonian, which can be done in a local fashion. This further
reduces the computational complexity of the tensor networks arising in the
minimization process compared to previous work. We test the accuracy of our
method by comparing the approximate energy spectrum to exact diagonalization
results for the random field Heisenberg model on 16 sites. We find that the
technique is highly accurate deep in the localized regime and maintains a
surprising degree of accuracy in predicting certain local quantities even in
the vicinity of the predicted dynamical phase transition. To demonstrate the
power of our technique, we study a system of 72 sites and we are able to see
clear signatures of the phase transition. Our work opens a new avenue to study
properties of the many-body localization transition in large systems.
