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.

I have the following code to construct a tridiagonal matrix:

ClearAll[th];
th[nwells_ /; EvenQ@nwells] := Module[
  {size = nwells, bdiag},
  bdiag = RandomReal[{0, 99}, size - 1];
  SparseArray[
   {
    Band[{1, 1}] -> bdiag,
    Band[{1, 2}] -> -0.5,
    Band[{2, 1}] -> -0.5},
   {size - 1, size - 1}
   ]
  ]   

This will be executed millions of times (bdiag is actually something that will change each time, so this is unavoidable). I'd like to speed it up as much as possible. Any ideas? I am interested in values of nwells of the order of 100 to 1000.

EDIT: Let us compare the time taken by the Band version, by MrW's version and Rojo's version for varying sizes:

{
  Table[{i, Do[th[i], {100}] // AbsoluteTiming // First}, {i, 100, 
    5000, 200}],
  Table[{i, Do[banded[i], {100}] // AbsoluteTiming // First}, {i, 100,
     5000, 200}],
Table[{i, Do[banded2[i], {100}] // AbsoluteTiming // First}, {i, 100,
         5000, 200}]
  } // ListLogPlot[#, AxesLabel -> {"size", "t"}] &

Mathematica graphics

(the slowest one is mine). Note the logarithmic axis. Evidently, the Band method falls behind more and more with larger system sizes.

Also, using Band unpacks:

On["Packing"]
th[3000]; // AbsoluteTiming
banded[3000]; // AbsoluteTiming

Mathematica graphics

This occurs when Band is used to insert the (packed) bdiag list into the diagonal.

share|improve this question
4  
The fastest solution is: th[nwells_ /; EvenQ@nwells, theta_] := Module[{size = nwells, bdiag, boffdiag}, boffdiag = SparseArray[{}, size - 2, -0.5]; bdiag = RandomReal[99, size - 1]; SparseArray[{{i_, i_} :> bdiag[[i]]}, size - {1, 1}] + Sum[DiagonalMatrix[boffdiag, i], {i, {-1, 1}}] ]... hurr hurr :D –  rm -rf Oct 13 '12 at 20:31
    
@rm-rf a man dare not have a life to keep pace with you ;-p –  Mr.Wizard Oct 13 '12 at 22:05

2 Answers 2

up vote 10 down vote accepted

This question came up in Chat the other day. Here is the solution I proposed.

banded[n_Integer?EvenQ] :=
  With[
   {main = RandomReal[99, n - 1],
    side = SparseArray[{}, n - 2, -0.5]},
   SparseArray[{i_, i_} :> main[[i]], n - 1] +
    Sum[side ~DiagonalMatrix~ i, {i, {-1, 1}}]
  ]

This uses several tricks and observations. Credit for the first goes to Norbert Pozar, who showed me that this is a very fast way to construct a diagonal SparseArray from a list x:

SparseArray[{i_, i_} :> x[[i]], Length @ x]

The second is my own observation that a DiagonalMatrix made from a SparseArray list is also created quickly. This in enhanced by creating the list with SparseArray[{}, n - 2, -0.5], where the "background" of the array is the element to be repeated, rather than ConstantArray or Table. One can see below that only minimal evaluation takes place:

SparseArray[{}, 10^6, -0.5] // InputForm
SparseArray[Automatic, {1000000}, -0.5, {1, {{0, 0}, {}}, {}}]

DiagonalMatrix is particularly fast with this input form.

These "tricks" are combined with the knowledge that adding sparse SparseArrays is fast, and that Band can be rather slow (again thanks to Norbert Pozar) to create a solution that is about fifty times faster than the original.

share|improve this answer
    
Maybe part of the speed difference is related to unpacking? Also, it appears that how faster this is (faster than Band) depends on the system size--see update in my question. I wonder why. –  acl Oct 13 '12 at 22:45

Given the millions of times that code will be run, for a size range of about 1000, all I could think to speed up @Mr.Wizard 's code is to memoize what stays fixed

i : diags[n_] := 
  i = With[{side = SparseArray[{}, n - 2, -0.5]}, 
    Sum[side~DiagonalMatrix~i, {i, {-1, 1}}]];

banded2[n_Integer?EvenQ] := 
 With[{main = RandomReal[99, n - 1]}, 
  SparseArray[{i_, i_} :> main[[i]], n - 1] + diags[n]]

Comparing the two

{Table[{i, Do[banded2[i], {100}] // AbsoluteTiming // First}, {i, 100,
     5000, 200}], 
  Table[{i, Do[banded[i], {100}] // AbsoluteTiming // First}, {i, 100,
     5000, 200}]} // ListPlot[#, AxesLabel -> {"size", "t"}] &

Mathematica graphics

share|improve this answer
    
Good addition. +1 –  Mr.Wizard Oct 14 '12 at 1:57

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.