# Why is this iterative process with matrix calculation so slow?

I am trying something similar to the following code with Ntime as large as 100 or so. Now it's very slow as shown by the Ntime=3 case that costs 40s on my laptop. Each time the matrix elements are updated by u and also used to calculate fun and new u.

Is there any way to speed up? I recently learned that Developer`ToPackedArray could be a nice trick. But I failed to apply it here with much difference.

Ntime = 3; Nband = 500; bias = 50;
t1[u_, n_] := 1 + u[[Mod[n + 1 - 1, Nband] + 1]] - u[[n]];
mat[u_] :=
bias IdentityMatrix[Nband] +
DiagonalMatrix[Table[t1[u, n], {n, 1, Nband - 1}], 1] +
DiagonalMatrix[Table[t1[u, n], {n, 1, Nband - 1}], -1] +
DiagonalMatrix[Table[t1[u, n], {n, Nband, Nband}], Nband - 1] +
DiagonalMatrix[Table[t1[u, n], {n, Nband, Nband}], -(Nband - 1)];
fun0[u_] :=
Total[Eigenvalues[mat[u]][[1 ;; Nband/2]]] +
Sum[(u[[n]] - u[[n + 1]])^2, {n, Nband - 1}];
fun[u_] :=
Table[fun0[ReplacePart[u, i -> u[[i]] + 1]], {i, Nband}] - fun0[u];

u = Table[Cos[0.2 i], {i, Nband}];
ulist = {u};
Table[u = u + fun[u]; AppendTo[ulist, u];, {i, 1, Ntime}]; //
AbsoluteTiming // First
• When you call Eigenvalues[mat[u]][[1 ;; Nband/2]]], are you looking for the larger half of the spectrum? Currently, you take that part of the spectrum that consists of the eigenvalues that have largest absolute values. – Henrik Schumacher Apr 22 at 16:23
• @HenrikSchumacher Ah, it probably doesn’t matter so much. I’m just giving the MWE with some calculations similar to what I really want to do. But yes, it’s more like the larger half. – xiaohuamao Apr 22 at 19:47