# How to construct this matrix fast?

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?

• Try a simple change mat2[kx_, n_] := ( h = Table[0, {i, 1, 2 n}, {j, 1, 2 n}]; a[m_] = 2 m - 1; b[m_] = 2 m; t1 = 1; aa = 1; e = (E^(I kx ((Sqrt[3] aa)/2)) + E^(-I kx ((Sqrt[3] aa)/2))) // FullSimplify; Do[h[[b[i], a[i]]] = t1 e; h[[a[i], b[i]]] = t1 e; 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). – b.gates.you.know.what Aug 25 '13 at 8:18
• @b.gatessucks Thank you! Now it is 0.9second compared to the original 1.9 second. Can it be even faster? I hope it can reach 0.027s. – matheorem Aug 25 '13 at 8:21
• You can probably remove at least one of the Ifs and you can always try Compile. – b.gates.you.know.what Aug 25 '13 at 9:01
• @b.gatessucks Compile which part? how to ? I have tried before, not working well. maybe I tried wrong way – matheorem Aug 25 '13 at 9:08
• @b.gatessucks Thank you for reminds me to use Compile. But I got a littel problem using Compile. Can you take a look at the message I drop to Mr.Wizard? – matheorem Aug 25 '13 at 12:23

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.

• Thank you Mr.Wizard! you always have many clever methods. BTW, I was inspired by b.gatessucks to use Compile. If I compile specific dimension(e.g. matcompile = Compile[{kx}, Evaluate[mat[kx, 7]]],then Table[mat[kx, 7], {kx, 0., 1., 1/3000}]only takes 0.05s, and it is 10times faster than mat4 if n=7 )but I tried to make it general depending to dimension. So I write matcompile2[m_] := Compile[{kx}, Evaluate[mat[kx, m]]]. But of course this won't work! So how to do it? Do you have any idea? – matheorem Aug 25 '13 at 12:04
• @matheorem I'll try to take a look at that later today or tomorrow. – Mr.Wizard Aug 25 '13 at 18:08
• Hello, Mr.Wizard. Have you have any idea now? If you're busy, it doesn't matter, I can wait. – matheorem Aug 26 '13 at 8:51
• @matheorem I'm sorry, I forgot to come back to this and I'm just going to bed. Remind me again later. – Mr.Wizard Aug 26 '13 at 9:10
• Hi, Mr.Wizard! just a remind, maybe you have got an idea – matheorem Aug 28 '13 at 7:37