Speeding up construction of simple tridiagonal matrix

Question

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"}] &

(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

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

Answer 1

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.

Answer 2

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"}] &