Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
    
I wonder if anybody's already tried SparseArray[Band[{1, 1}] -> 1/Diagonal[A]]... –  J. M. 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
add comment

2 Answers

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 *)
share|improve this answer
    
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
add comment

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 *)
share|improve this answer
    
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
2  
@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
show 3 more comments

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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