# How to extract and compute on the diagonal entities of a sparse matrix very fast?

As could be seen in the following code:

AbsoluteTiming[
n  = 100000;
A  = SparseArray[{
Band[{1, 120}] -> -2., Band[{950, 1}] -> -1.,
Band[{1, 1}] -> 20., Band[{1, 100}] -> 2.,
Band[{6, 800}] -> 1.1}, {n, n}, 0.];
b  = SparseArray[Table[1., {i, n}]];
DA = Diagonal[A];

(* I think constructing B is time-consuming. My ParallelTable[] does
not work or show any improvement herein! *)

B = SparseArray[Table[(1/DA[[i]]), {i, 1, n}]];
V = DiagonalMatrix[SparseArray[B]];
]


I'm trying to extract the diagonal entries of a very large sparse matrix and to compute $1/a_{ii}$ to make my new large sparse diagonal matrix $V$. This process takes around 18 seconds, and I would like to accelerate this process.

• I wonder if anybody's already tried SparseArray[Band[{1, 1}] -> 1/Diagonal[A]]... Commented Jun 14, 2012 at 13:44
• @J.M. yes, but it's much slower than DiagonalMatrix[1/Diagonal[A]] -- I cannot remember who first put me onto it but Band often is not fast(est). -- I found a record: it was Norbert Pozar who first showed me that Band can be much slower than alternatives. Commented Jun 14, 2012 at 13:48

You can create your new diagonal matrix V in a single step as:

V = DiagonalMatrix@SparseArray[1/Normal[Diagonal[A]]];


On my machine, this takes 0.05 seconds, compared to 9 seconds for your code above (excluding time taken to construct A).

You can verify that they're both the same:

DiagonalMatrix[SparseArray[B]] == DiagonalMatrix@SparseArray[1/Normal[Diagonal[A]]]
(* True *)

• Thanks a lot R.M., your tip works very well. In my computer, now it takes 3.5 seconds, compared to 19 seconds at the beginning (including the time for constructing A).
– Faz
Commented Jun 14, 2012 at 14:00

I'm probably missing something important here, but it seems to me that one does not have convert back and forth to Normal form, meaning that DiagonalMatrix[ 1/Diagonal[A] ] works:

DiagonalMatrix[ 1/Diagonal[A] ] == DiagonalMatrix[ SparseArray[B] ]

(* True *)

• Mr. Wizard, what could be fatser than using Band[] in constructing sparse matrices of the large scale?
– Faz
Commented Jun 14, 2012 at 14:01
• @FazlollahSoleymani In fact, Band can be quite slow. Have a look here for some alternatives. Commented Jun 14, 2012 at 14:16
• This was the first thing I tried, but I get Power::infy: "Infinite expression 1/0. encountered" even though the diagonals are all 20... know why?
– rm -rf
Commented Jun 14, 2012 at 14:57
• @R.M. 0 is the background of the SparseArray. If you look at the object afterward you will see that ComplexInfinity is the new background. You could manipulate the background directly if needed. Commented Jun 14, 2012 at 15:00
• @R.M This is a known bug. You may also find this and this discussion interesting. Commented Jun 14, 2012 at 21:02

Let me join this old thread. One can set the default element to 1. and inverse the array as you want

HoldPattern@setDef[SparseArray@s___, x_] := SparseArray[#, #2, x, #4] &@s;
V == DiagonalMatrix[1/setDef[Diagonal[A], 1.]]
(* True *)


Timings:

Do[DiagonalMatrix@SparseArray[1/Normal[Diagonal[A]]], {1000}] // AbsoluteTiming
(* {12.273275, Null} *)

Do[DiagonalMatrix[1/setDef[Diagonal[A], 1.]], {1000}]; // AbsoluteTiming
(* {9.803133, Null} *)


It is nice to have a bit faster solution. However, I don't know the simpler way to define the default element than Leonid's technique.