  Mixed state entanglement measures can act as a versatile probes of many-body
systems. However, they are generally hard to compute, often relying on tricky
optimizations. One measure that is straightforward to compute is the
logarithmic negativity, yet done naively even this is still limited to small
system sizes. Here, we introduce a method to compute the logarithmic negativity
for arbitrary subsystems of a densely represented state, as well as block
subsystems of matrix product states. The method combines lazily evaluated,
tensor network representations of the partially transposed density matrix with
stochastic Lanczos quadrature, and is easily extendible to other quantities and
classes of many-body states. As examples, we compute the entanglement within
random pure states for density matrices of up to 30 qubits, explore scrambling
in a many-body quench, and match the results of conformal field theory in the
ground-state of the Heisenberg model for density matrices of up to 1000 spins.
An implementation of the algorithm has been made available in the open-source
library \textit{quimb}.
