8
$\begingroup$

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

$\endgroup$
9
$\begingroup$

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

$\endgroup$
  • 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$ – Geoff Oxberry Mar 17 '13 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 '13 at 22:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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