# 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]]... –  Ｊ. Ｍ. Jun 14 '12 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. –  Mr.Wizard Jun 14 '12 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). –  Fazlollah Soleymani Jun 14 '12 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? –  Fazlollah Soleymani Jun 14 '12 at 14:01
@FazlollahSoleymani In fact, Band can be quite slow. Have a look here for some alternatives. –  Leonid Shifrin Jun 14 '12 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 Jun 14 '12 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. –  Mr.Wizard Jun 14 '12 at 15:00
@R.M This is a known bug. You may also find this and this discussion interesting. –  Leonid Shifrin Jun 14 '12 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.

-