9
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ I wonder if anybody's already tried SparseArray[Band[{1, 1}] -> 1/Diagonal[A]]... $\endgroup$ – J. M. will be back soon Jun 14 '12 at 13:44
  • $\begingroup$ @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. $\endgroup$ – Mr.Wizard Jun 14 '12 at 13:48
10
$\begingroup$

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 *)
$\endgroup$
  • $\begingroup$ 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). $\endgroup$ – Fazlollah Jun 14 '12 at 14:00
3
$\begingroup$

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 *)
$\endgroup$
  • $\begingroup$ Mr. Wizard, what could be fatser than using Band[] in constructing sparse matrices of the large scale? $\endgroup$ – Fazlollah Jun 14 '12 at 14:01
  • 2
    $\begingroup$ @FazlollahSoleymani In fact, Band can be quite slow. Have a look here for some alternatives. $\endgroup$ – Leonid Shifrin Jun 14 '12 at 14:16
  • $\begingroup$ 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? $\endgroup$ – rm -rf Jun 14 '12 at 14:57
  • $\begingroup$ @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. $\endgroup$ – Mr.Wizard Jun 14 '12 at 15:00
  • $\begingroup$ @R.M This is a known bug. You may also find this and this discussion interesting. $\endgroup$ – Leonid Shifrin Jun 14 '12 at 21:02
0
$\begingroup$

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.

$\endgroup$

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.