  We introduce a unitary matrix-product operator (UMPO) based variational
method that approximately finds all the eigenstates of fully many-body
localized (fMBL) one-dimensional Hamiltonians. The computational cost of the
variational optimization scales linearly with system size for a fixed bond
dimension of the UMPO ansatz. We demonstrate the usefulness of our approach by
considering the Heisenberg chain in a strongly disordered magnetic field for
which we compare the approximation to exact diagonalization results.
