# 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, If the matrix is banded with bandwidth $$b$$, you can get the diagonal with $$b+1$$ matrix-vector products. Sum the 1st, (b+2)th, (2b+3)th, ... columns of the identity and multiply it against the matrix. Then dot the result with each of those columns of the identity to get the corresponding diagonal entries.
• $$[1, 0, 0, 1, 0, 0, 1, 0]A$$ - gives the 1st, 4th, and 7th diagonal values in the 1st, 4th, and 7th entries
• $$[0, 1, 0, 0, 1, 0, 0, 1]A$$ - gives the 2nd, 5th, and 8th diagonal values in the 2nd, 5th, and 8th entries
• $$[0, 0, 1, 0, 0, 1, 0, 0]A$$ - gives the 3rd and 6th diagonal values in the 3rd and 6th entries
If the matrix isn't exactly banded but the out-of-band elements are small enough, this should be a reasonable approximation with the choice of $$b$$ giving a trade off between cost and accuracy.