  Fast estimation of the single-particle density matrix is key to many
applications in quantum chemistry and condensed matter physics. The best
numerical methods leverage the fact that the density matrix elements
$f(H)_{ij}$ decay rapidly with distance $r_{ij}$ between orbitals. This decay
is usually exponential. However, for the special case of metals at zero
temperature, algebraic decay of the density matrix appears and poses a
significant numerical challenge. We introduce a gradient-based probing method
to estimate all local density matrix elements at a computational cost that
scales linearly with system size. For zero-temperature metals the stochastic
error scales like $S^{-(d+2)/2d}$, where $d$ is the dimension and $S$ is a
prefactor to the computational cost. The convergence becomes exponential if the
system is at finite temperature or is insulating.
