Eigenvalues and plot are very slow

Question

I'mt trying to find eigenvalues (and plot them) but the evaluation takes way too long. I don't think it supposed to be a complicated computations. What am I doing wrong?

n = 3;
mat1 = Table[E^(I α (x - 1)) KroneckerDelta[x, y], {x, 1, n}, {y, 1, n}];
mat2 = Table[E^(-I α (x - 1)) KroneckerDelta[x, y], {x, 1, n}, {y, 1, n}];
mat3 = Table[KroneckerDelta[x - 1, y], {x, 1, n}, {y, 1, n}];
mat4 = Table[KroneckerDelta[x + 1, y], {x, 1, n}, {y, 1, n}];
sigmap = {{0, 1}, {0, 0}};
sigmam = {{0, 0}, {1, 0}};
H = KroneckerProduct[mat1, sigmap] + KroneckerProduct[mat2, sigmam] 
    + E^(I k α) KroneckerProduct[mat2, sigmap]
    + E^(-I k α) KroneckerProduct[mat1, sigmam] 
    + KroneckerProduct[mat3, sigmap] + KroneckerProduct[mat4, sigmam];
energies = Eigenvalues[H] /. α -> Pi/701;
Plot[energies, {k, 0, 2}]

Answer 1

This answer is edited in order to provide a summary of the problem, the various strategies employed in tackling it and supply some statistics.

The original problem had to do with time, ie it took too long to produce output:

energies = Eigenvalues[H] /. α -> Pi/701;

even for modest values of the underlying parameters (mainly n).

People in the comments suggested that the bottleneck was due to precision issues. Namely, @KraZug suggested to replace 701 with 701. (see code segment above) while @Thies Heidecke opted for moving the ReplaceAll statement inside Eigenvalues[]. These suggestions aimed to 'force' the Eigenvalues bult-in to treat its input numerically, instead of symbolically.

At this point, it is worth the time to take a step back and contemplate on the nature of H. It is a $2n$ x $2n$ parametric matrix; its entries are, for the most part, complex (parametric) numbers and integer quantities (1 and 0). Indeed, forcing the entries of H to be numeric quantities would allow the built-in support for numeric calculations to kick-in. But neither suggestion fully achieved that end.

Consider the following example:

f[x_] := {x}
f[a] /. a -> 1.

The output of the last evaluation is {1} as we'd might probably expect but that is not what we needed. Tracing out the evaluation steps, ReplaceAll[] matches a only after the definition of f has been evaluated. That means, that when f acts upon its parameter (f[a]), a's value is still symbolic not numeric.

Also, evaluating f[a /. a -> 1] again produces {1},only this time, the evaluation of f's argument is performed in integer terms. Perhaps the output of f[ a / 2 /. a -> 1] and f[ a / 2 /. a -> 1.] is more enlightening. In the later case, f acts upon its input,after it has been transformed into a numeric quantity ({0.5}) while in the former it returns {1/2}.

On a separate note, @J.M. suggested in the comments to use SparseArray[] in order to tackle the problem of high dimensionality (n=200). That was the intended use-case of the original question's code. In my orignal answer, and for my original solution, after using SparseAray[], I stumbled upon Mathematica's notice that the calculation of Eigenvalues[] "...is likely to be faster with dense matrix methods, the sparse input matrix will be converted".

Additionaly, @Thies Heidecke pointed out that my using Compile[] was unlikely to provide significant advantages outside from, perhaps, speeding up the construction of the original matrices (inputs). An additional insight was that in the original code, matrix H was dependent on the symbolic parameter k. The workaround was sought in wrapping the calculations of Eigenvalues[H] in a function that allowed only numeric input (eg MyEigenvalues[kk_?NumericQ] - see @Thies Heidecke's answer).

Even though, such an approach allows for numeric input only, such a wrapper does not provide any additional advantage. Consider the folowing definition:

g[y_?NumericQ] := h[x /. x -> y]

Indeed, now x will be 'forced' to take numeric values (the definition of g does not allow symbolic inputs). Also, h will be sure to operate on numeric input. The same result could have been achieved with somethng analogous to f[ a / 2 /. a -> 1.] from the example above.

My original answer was not lightning-fast and the other answer by @Thies Heidecke reported similar issues (especially for n=200) and also insightfully noted in the comments that the time constraints on the problem are imposed by the run-time of Eigenvalues[] function.

The present edit aims to present a solution with more acceptable time contraints. To that extent I opted for compiling and memoizing the input expressions. Also, I compiled the process of calculating the eigenvalues of H (also compiled) in order to take advantage of RuntimeAttributes->{Listable}. The code can be found HERE.

The following graphs display the values of k I tested my code on, along with the time it took to produce the output plot.

(I am using the PlotStyle parameters used in the answer by @Thies Heidecke)

It fails to produce the artifacts that are present in the other solution. This has to do perhaps with the fact that I am sampling parameter k in discrete, non-necessarily integer, intervals. Note, how the new code makes use of Array[] to sample k energies[pts_,lo_,hi_]:=Transpose[eigenvals[Array[#&,pts,{lo,hi}]//N]]

The next graph plots pairs of {k,time} (in minutes, on my machine).

It appears that the energies scale linearly with the number of points we choose to sample k.

original answer :

This is what I tried : n = 200;

mat1[a_, n_] := Array[E^(I a (#1 - 1)) KroneckerDelta[#1, #2] &, {n, n}] mat2[a_, n_] := Array[E^(-I a (#1 - 1)) KroneckerDelta[#1, #2] &, {n, n}]

mat3[n_] := Array[KroneckerDelta[#1 - 1, #2] &, {n, n}] mat4[n_] := Array[KroneckerDelta[#1 + 1, #2] &, {n, n}]

sigmap = {{0, 1}, {0, 0}}; sigmam = {{0, 0}, {1, 0}};

Hfunc = Function[{k, a, n}, KroneckerProduct[mat1[a, n], sigmap] + KroneckerProduct[mat2[a, n], sigmam] + E^(I k a) KroneckerProduct[mat2[a, n], sigmap] + E^(-I k a) KroneckerProduct[mat1[a, n], sigmam] + KroneckerProduct[mat3[n], sigmap] + KroneckerProduct[mat4[n], sigmam]];

H = Function[{k}, Hfunc[k, Pi/701., n]];

energies = Compile[{{k, Real}}, Eigenvalues[H[k]], {{, _Real, 1}}]

ListPlot[Array[energies[#] &, 10, {0, 2}] // Transpose, Joined -> True]

I tried using SparseArray butMathematicacomplained :

Eigenvalues::arh : Because finding 400 out of the 400 eigenvalues and/or eigenvectors is likely to be faster with dense matrix methods, the sparse input matrix will be converted.If fewer eigenvalues and/or eigenvectors would be sufficient, consider restricting this number using the second argument to Eigenvalues. >>

Answer 2

n = 200;
mat1 = Table[E^(I \[Alpha] (x - 1)) KroneckerDelta[x, y], {x, 1, n}, {y, 1, n}];
mat2 = Table[E^(-I \[Alpha] (x - 1)) KroneckerDelta[x, y], {x, 1, n}, {y, 1, n}];
mat3 = Table[KroneckerDelta[x - 1, y], {x, 1, n}, {y, 1, n}];
mat4 = Table[KroneckerDelta[x + 1, y], {x, 1, n}, {y, 1, n}];
sigmap = {{0, 1}, {0, 0}};
sigmam = {{0, 0}, {1, 0}};

H = KroneckerProduct[mat1, sigmap] + KroneckerProduct[mat2, sigmam] + E^(I k \[Alpha]) KroneckerProduct[mat2, sigmap] + E^(-I k \[Alpha]) KroneckerProduct[mat1, sigmam] + KroneckerProduct[mat3, sigmap] + KroneckerProduct[mat4, sigmam];
MyEigenvalues[kk_?NumericQ] := Sort@Eigenvalues[H /. {\[Alpha] -> Pi/701., k -> kk}]

(g = With[{maxk=2},ListLinePlot[Transpose@Table[MyEigenvalues[k], {k, 0, maxk}], DataRange -> {0, maxk}, PlotStyle -> Directive[Opacity[0.15], Blue]]]) //Timing

gives a computation time of 1.19s for the g part on my machine (2010 MBP).

Now increasing the maximum k from 2 to 1000

(data = With[{maxk = 1000}, Transpose@Table[MyEigenvalues[k], {k, 0, maxk}]];) // Timing
(g = With[{maxk = 1000}, ListLinePlot[data, DataRange -> {0, maxk}, PlotStyle -> Directive[Opacity[0.15], Blue]]]) // Timing

gives me 241s and 4s for the computation of data and the plot g respectively, and produces this plot:

I didn't manage to get rid of the artifacts at the multiples of 100, but i think the plot still comes out nice. If anyone has an idea how to get rid of the artifacts, i'm happy to hear and improve the answer.