# Speeding up construction of simple tridiagonal matrix

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; // AbsoluteTiming
banded; // AbsoluteTiming This occurs when Band is used to insert the (packed) bdiag list into the diagonal.

• 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, 2012 at 20:31
• @rm-rf a man dare not have a life to keep pace with you ;-p Oct 13, 2012 at 22:05

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.

• 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, 2012 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"}] & • Good addition. +1 Oct 14, 2012 at 1:57