Half a year ago, John Wright and I were considering almost the same question, in connection with the quantum tomography problem; we even asked a few people, including Suvrit and fedja. The only very slight difference is we were hoping to show that $M$ is close to its own diagonal, rather than to any diagonal matrix. (Note that since we have the hypothesis that $M$'s diagonal is close to its spectrum, our question is equivalent to asking if $M$ is close to the diagonal matrix made from its spectrum.)
Now if you insist on this, that $M$ be close to its own diagonal, then I think the answer to your question would be "no". I'm not 100% certain, but here's my reasoning. Let $A$ be any diagonal matrix of nonnegative reals (sorted, say) and let $H$ be any Hermitian matrix. Now suppose $M = \exp(-i t H) A \exp(i t H)$, so $M$ is PSD with spectrum given by $A$'s diagonal. Think of $t$ as a positive real tending to $0$. Now let's consider $\Delta := D(M) - A = D(M-A)$. By Taylor expansion,
\begin{align*}
\Delta &= D(A + i t (AH - HA) - \tfrac{t^2}{2}(AH - 2H^2 + HA) + O(t^3) - A) \\
&= \tfrac{t^2}{2} D(AH - 2H^2 + HA) + O(t^3),
\end{align*}
where we used that $AH - HA$ has $0$ diagonal. Thus, thinking of $t \to 0$, we will have $|D(M) - A|_{\mathrm{tr}}$ proportional to $t^2$ (or smaller).
That means you're hoping that $M$ is $\Theta(t^2)$-close to diagonal. Now if, like us, you are further hoping $M$ is $\Theta(t^2)$-close to its own diagonal then you're out of luck: it's equivalent to showing that $|M - A|_{\mathrm{tr}} = \Theta(t^2)$, but as we saw, $M - A = it(AH - HA) + O(t^2)$, and hence $|M - A|_{\mathrm{tr}} = \Theta(t)$ (unless $AH- HA = 0$, but this need not be the case).
If you restrict $M$ to have trace 1, then we felt the desired result was true but with a weaker conclusion of $O(\sqrt{\epsilon})$ rather than $2\epsilon$.
--
In case this is related to our upcoming tomography papers, I'd be happy to continue the discussion over email :)
Best,
Ryan
