# Extracting diagonal of an approximately diagonal matrix when we don't know its entries

What is a good way to extract the diagonal from a symmetric matrix that is already almost diagonal, but where you don't have the matrix elements (only the ability to apply it to vectors)?

Further constraints are, (1) applying the n-by-n matrix n-times to explicitly construct the diagonal would be prohibitively costly, and (2) the small elements of the diagonal are important in addition to the large elements.

Here is an example picture of the type of matrix I want to extract the diagonal from (in a small-scale test case): 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, 