How to construct this matrix fast?

Question

I write a function which is a matrix as below:

mat[kx_, n_] := (
   Clear[a, b, h];
   h = Table[0, {i, 1, 2 n}, {j, 1, 2 n}];
   a[m_] := 2 m - 1;
   b[m_] := 2 m;
   t1 = 1;
   aa = 1;
   Do[
   h[[b[i], a[i]]] = 
    t1 (E^(I kx ((Sqrt[3] aa)/2)) + E^(-I kx ((Sqrt[3] aa)/2)));
   h[[a[i], b[i]]] = 
    t1 (E^(I kx ((Sqrt[3] aa)/2)) + E^(-I kx ((Sqrt[3] aa)/2)));
   If[i + 1 <= n, h[[b[i + 1], a[i]]] = t1, Null];
   If[i - 1 > 0, h[[a[i - 1], b[i]]] = t1, Null];
   , {i, 1, n}];
   h)

I test the speed of the matrix construction

In[90]:= Table[mat[kx, 7], {kx, 0., 1., 1/3000}]; // AbsoluteTiming

Out[90]= {1.971113, Null}

the dimension of the mat[kx,7] is 14, so I test the following randomreal matrix

In[91]:= Table[
   RandomReal[{0, 1}, {14, 14}], {kx, 0., 1., 
    1/3000}]; // AbsoluteTiming

Out[91]= {0.027002, Null}

it is 100 times faster.

So I wonder if my mat function could be 100 times faster. I know my programming is poor, but I don't know how to code it elegent and efficient.So can somebody help me with it and point it out which is the most important factor that make my mat function so slow?

Answer 1

Here are a few methods. I am focusing on large values of n. A couple of these return SparseArray objects but that may even be desirable depending on your application.

mat2[kx_, n_] :=
 With[{x = 2 Cos[Sqrt[3] kx/2] // TrigToExp},
  Riffle[
   Array[PadRight[{x, 0, 1}, 2 n, 0, 2 # - 1] &, n],
   Array[PadRight[{1, 0, x}, 2 n, 0, 2 # - n] &, n]
  ]
 ]

mat4[kx_, n_] :=
 With[{x = 2 Cos[Sqrt[3] kx/2] // TrigToExp},
  With[{a = SparseArray@{x, 0, 1}, b = SparseArray@{1, 0, x}},
   Riffle[
    Array[PadRight[a, 2 n, 0, 2 # - 1] &, n],
    Array[PadRight[b, 2 n, 0, 2 # - n] &, n]
   ]
  ]
 ]

mat6[kx_, n_] := With[{x = 2 Cos[Sqrt[3] kx/2] // TrigToExp},
  SparseArray[{
    Band[{1, 2}, Automatic, {2, 2}] -> x,
    Band[{2, 1}, Automatic, {2, 2}] -> x,
    Band[{1, 4}, Automatic, {2, 2}] -> 1,
    Band[{4, 1}, Automatic, {2, 2}] -> 1
    }, {2 n, 2 n}]
  ]

Timings for n = 1500: (search site for timeAvg):

timeAvg @ #[1, 1500] & /@ {mat, mat2, mat4, mat6}

{0.3338, 0.1278, 0.011608, 0.01684}

---

Since I have learned that you intend to apply the function for a fixed n and a changing kx it does make sense to compile the generated array and substitute kx values afterward. I would do it like this:

mem : matC[n_Integer?Positive] :=
  mem = Block[{kx}, Compile[{kx}, mat2[kx, n] // Evaluate]]

Now:

Table[matC[7][kx], {kx, 0`, 1`, 1`/3000}] // timeAvg
matC[7] /@ Range[0`, 1`, 1`/3000]         // timeAvg

0.03496

0.01684

The memoization keeps the form with Table from being very slow, but it is still not as fast as Map. Even though matC[7] will evaluate to a memoized compiled function after the first call it is still called 3,000 times in Table, whereas with Map it evaluates once, and then this compiled function is directly applied to each value in the range.