  We introduce an optimal strategy to sample quantum outcomes of local
measurement strings for isometric tensor network states. Our method generates
samples based on an exact cumulative bounding function, without prior
knowledge, in the minimal amount of tensor network contractions. The algorithm
avoids sample repetition and, thus, is efficient at sampling distribution with
exponentially decaying tails. We illustrate the computational advantage
provided by our optimal sampling method through various numerical examples,
involving condensed matter, optimization problems, and quantum circuit
scenarios. Theory predicts up to an exponential speedup reducing the scaling
for sampling the space up to an accumulated unknown probability $\epsilon$ from
$\mathcal{O}(\epsilon^{-1})$ to $\mathcal{O}(\log(\epsilon^{-1}))$ for a
decaying probability distribution. We confirm this in practice with over one
order of magnitude speedup or multiple orders improvement in the error
depending on the application. Our sampling strategy extends beyond local
observables, e.g., to quantum magic.
