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):

enter image description here


2 Answers 2


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); 

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, enter image description here

  • 3
    $\begingroup$ Don't hesitate to answer your own question. This is exactly the kind of thing we want on SciComp (I've done it). You may want to wait on accepting your own answer, in case a better answer appears. It's much better to answer your own question (when you have an answer) than to leave it unanswered; we encourage all users to follow your example, if possible. $\endgroup$ Mar 17, 2013 at 20:00
  • $\begingroup$ Ok, glad to hear that! I will wait on accepting for a few days in case anyone else has a better answer. $\endgroup$
    – Nick Alger
    Mar 17, 2013 at 22:41

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.

For example, for an 8x8 matrix, A, with bandwidth 2,

  • $[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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.