Optimising eigenvalue matrix calculation

Question

I'm quite new to mathematica and am trying to figure out how to optimise my code as currently it takes a very long time to compile and makes it hard to bug test. The first part that creates the graph of the eigenvalues against k is relatively fast, but the second part takes a large chunk of time so i'm wondering what i should change about this. Any help/advice would be really appreciated.

t4 = 0.3
phi = Pi/2

haldzigzaghamiltonian[n_, k_, delta_] := Table[
   If[Mod[u, 2] != 0 && v == (u + 1), -2*Cos[k/2],
        If[Mod[u, 2] == 0 && v == (u - 1), -2*Cos[k/2],
        If[Mod[u + v - 1, 4] == 0 && Abs[u - v] < 2, -1,
      If[
       Mod[u, 2] != 0 && (v == (u + 2) || v == (u - 2)), -t4*2*
        Cos[phi + k/2],
       If[
        Mod[u, 2] == 0 && (v == (u + 2) || v == (u - 2)), -t4*2*
         Cos[k/2 - phi], 
        If[Mod[u, 2] != 0 && u == v, -t4*2*Cos[phi + k] + delta, 
         If[Mod[u, 2] == 0 && u == v, -t4*2*Cos[k + phi] - delta, 
          0]]]]]]],
   {u, 2*n}, {v, 2*n}];

vx[n_, k_] := Table[If[u == v, -t4*Sin[k + phi],
        If[(u + v > 3) && (Abs[u - v] == 2), -t4*Sin[k/2 - phi],
         If[(Mod[v, 2] == 0) && Abs[u - v] == 1, -Sin[k/2], 0]]],
   {u, 2*n}, {v, 2*n}];

oka[n_] := 
  Table[If[Mod[u, 2] != 0 && v == u, 1, 0], {v, 2*n}, {u, 2*n}];
okb[n_] := 
  Table[If[Mod[u, 2] == 0 && v == u, 1, 0], {v, 2*n}, {u, 2*n}];

nvalue = 3;
k = Range[0, 2*Pi, 2*Pi/99];
delta = 0.3;
eigvals = 
  Table[N[
    Re[Eigenvalues[
      haldzigzaghamiltonian[nvalue, k[[i]], delta]]]], {i, 100}] ;

ListPlot[
 Table[{k[[i]], eigvals[[i]][[j]]}, {i, 0, 100}, {j, 0, 2*nvalue}]]

nvalue = 20;
kvalue = 7*Pi/9;
del = Table[
   eigvecsk = 
    Eigenvectors[{haldzigzaghamiltonian[nvalue, k[[i]], delta]}];
   eigvalsk = 
    Eigenvalues[haldzigzaghamiltonian[nvalue, k[[i]], delta]];
   f[n_] := If[eigvalsk[[n]] <= 0, 1, 0];
   e[n_, nmax_, k_] := eigvalsk[[n]];
   onm[n_, m_, nmax_] := eigvecsk[[n]]\[Conjugate] . eigvecsk[[m]];
   vxnm[n_, m_, nmax_, k_] := 
    eigvecsk[[n]]\[Conjugate] . vx[nmax, k] . eigvecsk[[m]];
   func[n_, m_, nmax_, 
     k_] := (f[n] - f[m])/(e[n, nmax, k] - e[m, nmax, k])^2*
     onm[n, m, nmax]*vxnm[m, n, nmax, k];
   N[Re[Sum[
      If[p != q, func[p, q, nvalue, k[[i]]], 0], {p, 1, 2*nvalue}, {q,
        1, 2*nvalue}]]], {i, 100}];
ListLinePlot[Table[{k[[b]], del[[b]]}, {b, 0, 100}], 
 PlotRange -> Full]

Answer 1

There are some tricks that can vastly shorten the computation time for generating del. However, I want to mention that f in your code is always 0 due to the fact that $$\text{Eigenvectors of a symmetric matrix are mutually orthogonal.}$$

The matrix haldzigzaghamiltonian in your code is symmetric and real (unless you set delta to something else having non-zero imaginary part), making onm always 0 (so are f and the Sum).

--

Anyways, vxnm is taking up the most computation time (~99%) because it's evaluating the same vx[nvalue, k[[i]]] for (2*nvalue)*(2*nvalue) times during the Sum, only to generate one entity in del. To prevent this, define vx0 = vx[nvalue, k[[i]]] before the Sum and change vx[nmax, k] to vx0 in vxnm:

del = Table[
 vx0 = vx[nvalue, k[[i]]];
 eigvecsk = Eigenvectors[{haldzigzaghamiltonian[nvalue, k[[i]], delta]}];
 eigvalsk = Eigenvalues[haldzigzaghamiltonian[nvalue, k[[i]], delta]];
 f[n_] := If[eigvalsk[[n]] <= 0, 1, 0];
 e[n_] := eigvalsk[[n]];
 onm[n_, m_] := eigvecsk[[n]]\[Conjugate].eigvecsk[[m]];
 vxnm[n_, m_, nmax_, k_] := eigvecsk[[n]]\[Conjugate].vx0.eigvecsk[[m]];
 func[n_, m_, nmax_, k_] := (f[n] - f[m])/(e[n] - e[m])^2*onm[n, m]*vxnm[m, n, nmax, k];
 Re@Sum[If[p != q, func[p, q, nvalue, k[[i]]], 0], {p, 1, 2*nvalue}, {q, 1, 2*nvalue}]
, {i, 100}] // AbsoluteTiming;
First[%]

(* 18.8918 *)

The following graph is a benchmark, showing that this makes it about 44 times faster, asymptotically:

Answer 2

Being an interpreted language, Mathematica is very slow with all loop constructs such as Table. So one should try to use matrix constructors whose backend is alread compiled. Moreover, the return haldzigzaghamiltonian is a sparse matrix with only few bands. Hence we can employ a SparseArray which will also safe quite some memory if n gets larger.

haldzigzaghamiltonian2[n_, k_, delta_] := With[{
    diag0 = -2. t4 Cos[phi + k] - (-1.)^Range[2 n] delta,
    diag1 = 
     Riffle[ConstantArray[-2. Cos[0.5 k], n], 
      ConstantArray[-1., n - 1]],
    diag2 = 2. t4 Cos[0.5 k + (-1.)^Range[2 n - 2] phi]
    },
   SparseArray[
    {
     Band[{1, 1}] -> diag0,
     Band[{1, 2}] -> diag1,
     Band[{2, 1}] -> diag1,
     Band[{1, 3}] -> diag2,
     Band[{3, 1}] -> diag2
     },
    {2 n, 2 n}
    ]
   ];
oka2[n_] := DiagonalMatrix[SparseArray[Mod[Range[2 n], 2]]];
okb2[n_] := DiagonalMatrix[SparseArray[1 - Mod[Range[2 n], 2]]];

nvalue = 3;
k = Subdivide[0., 2. Pi, 99];
delta = 0.3;
eigvals = Table[Sort@ Eigenvalues[haldzigzaghamiltonian2[nvalue, k[[i]], delta]], {i, 100}];
ListLinePlot[Thread[{k, #}] & /@ Transpose[eigvals]]

Now to the rest of your code: Having had a closer look on it, I am quite confident that a more precise and significantly faster way to get its result would be

del = ConstantArray[0., 100];

;)

The reason is that onm[n_, m_, nmax_] is essentially equivalent to KroneckerDelta because Eigenvectors orthonormalizes the eigenvectors for machine precision inputs. But since you evaluate func only if m and n are distint, func always returns 0. (or something very close to it due to floating point errors).

Just for completeness, here the way I would rewrite your code:

vx2[n_, k_] := SparseArray[{
    Band[{1, 1}] -> -t4*Sin[k + phi],
    Band[{1, 2}] -> -Sin[k/2] Mod[Range[1, 2 n - 1], 2],
    Band[{2, 1}] -> -Sin[k/2] (1 - Mod[Range[1, 2 n - 1], 2]),
    Band[{1, 3}] -> -t4*Sin[k/2 - phi],
    Band[{3, 1}] -> -t4*Sin[k/2 - phi]
    },
   {2 n, 2 n}
   ];
id = IdentityMatrix[2 nvalue, SparseArray, WorkingPrecision -> MachinePrecision];
nvalue = 20;
del = Table[
   {eigvalsk, eigvecsk} = Eigensystem[haldzigzaghamiltonian2[nvalue, k[[i]], delta]];
   fvec = 1. - UnitStep[eigvalsk];
   fdiff = (# - Transpose[#] &[ConstantArray[fvec, 2 nvalue]]);
   ediff = # - Transpose[#] &[ConstantArray[eigvalsk, 2 nvalue]];
   V = ConjugateTranspose[eigvecsk].vx[nvalue, k[[i]]].eigvecsk;
   W = ConjugateTranspose[eigvecsk].eigvecsk;
   B = fdiff (ediff + id)^-2 V W;
   Total[B, 2] - Total[Diagonal[B]]
   , {i, 100}];

It gets the job done in half a second.