I'm going to answer my own question since the following method seems to be very effective. I'm making it an answer so people can upvote or downvote it independently of the question if they think it is good or bad.
Answer: use randomized matrix probing applied to the diagonal of the matrix.
Let $A$ be the operator we wish to find the diagonal of, and let $\omega_1,\omega_2,..,\omega_k$ be a small number of gaussian random vectors. Then apply $A$ to the random vectors to get $A \omega_1, A \omega_2, ..., A \omega_k$ and solve the following least-squares minimization problem,
$$\min_{\text{diagonal }D} ||D \omega_1-A\omega_1||^2 + ||D \omega_2-A\omega_2||^2 + ... + ||D \omega_k-A\omega_k||^2.$$
The minimum has the exact formula,
$$d_i=\frac{\omega_1^i A \omega_1^1 + \omega_2^i A \omega_2^i... + \omega_k^i A \omega_k^i}{(\omega_1^i)^2 + (\omega_2^i)^2 + ... + (\omega_k^i)^2}.$$
Matlab code, for example:
omegas = randn(16,3);
dprobe=sum(omegas.*(A*omegas),2)./sum(omegas.^2,2);
In my example matrix, with 3 probing vectors, the exact diagonal and probed diagonal compare as follows:
[dprobe, diag(A)]
ans =
1.0e+04 *
2.3297 2.4985
0.4596 0.4921
0.1322 0.0897
0.2838 0.1764
0.0989 0.0999
0.0106 0.0071
0.0068 0.0068
0.0469 0.0571
0.0070 0.0070
0.0355 0.0372
0.0059 0.0060
0.0071 0.0064
0.0067 0.0067
0.0026 0.0021
0.0012 0.0012
0.0015 0.0013
Update: I've been experimenting applying these ideas to symmetric block matrices, since a matrix i'm working with is almost block diagonal in a wavelet-like basis. It seems to work pretty well for building preconditioners so long as the matrix is "block-diagonally-dominant" (definition is a little tricky), and so long as you symmetrize the least-squares reconstructed blocks.
Recall that a matrix partitioned into blocks $A_{i,j}$ is block-diagonally-dominant if
$$||A_{i,i}^{-1}||^{-1} \ge \sum_j ||A_{i,j}||.$$
Given gaussian random $\omega$'s as above, we seek to find the following least-squares block diagonal reconstruction:
$$\min_{\mathrm{block~diagonals~}B} ||B \omega_1-A\omega_1||^2 + ||B \omega_2-A\omega_2||^2 + ... + ||B \omega_k-A\omega_k||^2.$$
After some tensor product manipulations, you can find the exact formula for the $l$'th block $\tilde{B}^{(l)}$ by solving the local problems:
$$\tilde{B}^{(l)} = [(A\omega_1)^{(l)}\omega_1^{(l)T} + ... + (A\omega_k)^{(l)}\omega_k^{(l)T}][\omega_1^{(l)}\omega_1^{(l)T} + ... + \omega_k^{(l)}\omega_k^{(l)T}]^{-1},$$
where $(A\omega_i)^{(l)}$ and $\omega_i^{(l)}$ are the portions of $A\omega_i$ and $\omega_i$ corresponding to the indices of the $l$'th block.
If I just use these $\tilde{B}$'s, the preconditioning seems to be pretty bad, but if I symmetrize as follows,
$$B^{(l)} = (\tilde{B}^{(l)} + \tilde{B}^{(l)T})/2,$$
in my experiments it becomes almost as good as if I had used the true diagonal blocks (often better!). Here is an example matrix in pictures,
