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 that 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.