# Speeding up Fourier transform and eigenvector calculation

I have a matrix $$M$$ in Bloch basis (momentum space). Then I have to take Fourier transform of $$M$$ to go to real space. The size of this lattice in real space can be fixed to some value $$N_{unit}$$. Furthermore, I have to determine the eigenvectors of this real space matrix and study dynamics so it is very slow for bigger matrix sizes.

1. Let us say $$M$$ is a square lattice with periodic boundaries can be reduced to a matrix hsquare (kx, ky) by taking Fourier transform, where $$k_x$$ and $$k_y$$ are the conjugate variables to $$x$$ and $$y$$.
1. Then we take the inverse Fourier transform of hsquare (kx, ky) to go back to real space, i.e. matrix hsquareFin(x,y).
2. Furthermore, I have to find the eigenvectors of hsquareFin(x,y). That I plot on top of underlying lattice which is square lattice here.
3. Lastly, I am interested in the dynamics of how a pulse excited at a point in the real space will propagate according to the Schrödinger equation.

My MWE:

Nunit = 400; (* Number of unitcells *)

cond = {J1x -> 1/1, J2x -> 1, J1y -> 1/1};

hsquare[kx_, ky_] = (2 J1x Cos[kx] + 2 J1y Cos[ky] + x1) /. cond

square[kx_, ky_] = 2 Cos[kx] + 2 Cos[ky];

(*The Hamiltonian matrix for the Finite Case*)

DistributeDefinitions[Nunit, hsquare];
hsquareFinx[ky_] =
ParallelTable[#, {i, Nunit}, {j, Nunit}] &@(InverseFourierTransform[
hsquare[kx, ky]/Sqrt[2 \[Pi]], kx,
x] /. {DiracDelta -> KroneckerDelta, x -> j - i, x1 -> i}) //
ArrayFlatten;

DistributeDefinitions[hsquareFinx];
hsquareFin =
ParallelTable[#, {i, Nunit}, {j, Nunit}] &@(InverseFourierTransform[
hsquareFinx[ky]/Sqrt[2 \[Pi]], ky,
y] /. {DiracDelta -> KroneckerDelta, y -> j - i}) //
ArrayFlatten;

graphsq[ky_] =
Table[#, {i, Nunit}, {j, Nunit}] &@(InverseFourierTransform[
square[kx, ky]/Sqrt[2 \[Pi]], kx,
x] /. {DiracDelta -> KroneckerDelta, x -> j - i}) //
ArrayFlatten;

Table[#, {i, Nunit}, {j, Nunit}] &@(InverseFourierTransform[
graphsq[ky]/Sqrt[2 \[Pi]], ky,
y] /. {DiracDelta -> KroneckerDelta, y -> j - i}) //
ArrayFlatten, VertexLabels -> Automatic];


The eigenvectors are calculated

graphcoordsq = GraphEmbedding[sqg]; (*Generating underlying lattice*)

egvecsq = Eigenvectors[hsquareFin // N];

Manipulate[
Show[ListDensityPlot[
Table[
Append[Chop@graphcoordsq[[i]], Chop@(Abs[egvecsq[[n, i]]]^2)] //
N, {i, Length[egvecsq[[n]]]}], PlotLegends -> Automatic,
PlotLabel -> "Probability density"], sqg], {n, 1, Length[egvecsq],
1}]


Is there a way to speed up the process? Because I also have to study the time dynamics of the eigenstates later.

• It seems as if your matrixes have many 0 entries. I have not waited for the computations to finish. But I think the resulting matrices will be banded with a very simple structure. You should seek to exploit that, e.g. by using SparseArray. Commented Jan 21, 2022 at 18:52
• @Shamina It could be better to explain you problem as mathematical problem as well. Commented Jan 22, 2022 at 5:16
• Yes. The Table in the function hsquareFin should not loop over all {i,j} pairs. Try to loop over those pairs that lead to a nonzero entry. Afterwards, you can use the SparseArray constructor to construct the sparse matrix from the nonzero positions and the nonzero values. Commented Jan 24, 2022 at 17:35
• @CATrevillian Thanks a lot! Let me have a look into that. I don’t think I’ve fully got it but I am going to give it a try. Commented Jan 28, 2022 at 21:38
• @Shamina Do you understand that we have to compute matrix $Nunit^2\times Nunit^2$, therefore for Nunit=400 it gonna be $160000\times 160000$. What is the reason to compute so big matrix? For Nunit<=100 there is algorithm faster then your code. Commented Jan 29, 2022 at 4:37

This code shows how fast SparseArray compare to ParallelTable for moderate $$N_{unit}\le 10^2$$. Note, that during computation code generates $$N_{unit}^2\times N_{unit}^2$$ matrixes hsquareFin1=hsquareFin and solutions egvecsq1=egvecsq. For $$N_{unit}>100$$ the computation time is limited by time of Eigenvectors evaluation, and therefore we can't decrease computation time for Nunit=400 as it required. To compare computation time we can organize Module as follows

Clear["Global*"]

cond = {J1x -> 1, J2x -> 1, J1y -> 1};

hsquare[kx_, ky_] = (2 J1x Cos[kx] + 2 J1y Cos[ky] + x1) /. cond;

square[kx_, ky_] = 2 Cos[kx] + 2 Cos[ky];

fcomp[nn_] :=
Module[{Nunit = nn},
t1 = AbsoluteTiming[DistributeDefinitions[Nunit, hsquare];
hsquareFinx[ky_] =
ParallelTable[#, {i, Nunit}, {j,
Nunit}] &@(InverseFourierTransform[
hsquare[kx, ky]/Sqrt[2 \[Pi]], kx,
x] /. {DiracDelta -> KroneckerDelta, x -> j - i,
x1 -> i}) // ArrayFlatten;

DistributeDefinitions[hsquareFinx];
hsquareFin =
ParallelTable[#, {i, Nunit}, {j,
Nunit}] &@(InverseFourierTransform[
hsquareFinx[ky]/Sqrt[2 \[Pi]], ky,
y] /. {DiracDelta -> KroneckerDelta, y -> j - i}) //
ArrayFlatten;
graphsq[ky_] =
Table[#, {i, Nunit}, {j, Nunit}] &@(InverseFourierTransform[
square[kx, ky]/Sqrt[2 \[Pi]], kx,
x] /. {DiracDelta -> KroneckerDelta, x -> j - i}) //
ArrayFlatten;

sqg =
Table[#, {i, Nunit}, {j, Nunit}] &@(InverseFourierTransform[
graphsq[ky]/Sqrt[2 \[Pi]], ky,
y] /. {DiracDelta -> KroneckerDelta, y -> j - i}) //
ArrayFlatten, VertexLabels -> Automatic];
(*The eigenvectors are calculated*)

graphcoordsq =
GraphEmbedding[sqg];(*Generating underlying lattice*)
egvecsq = Eigenvectors[hsquareFin // N // Normal];][[1]];

t2 = AbsoluteTiming[
s1 =
With[{n = Nunit},
SparseArray[{Band[{1, 1}] ->
Flatten[Table[Table[i, {i, n}], {n}]],
Band[{1, n + 1}] -> 1, Band[{n + 1, 1}] -> 1,
Band[{1, 2}, {n^2, n^2}] -> Table[If[i == n, 0, 1], {i, n}],
Band[{2, 1}, {n^2, n^2}] ->
Table[If[i == n, 0, 1], {i, n}]}, {n^2, n^2}]];
s2 =
With[{n = Nunit},
SparseArray[{Band[{1, n + 1}] -> 1, Band[{n + 1, 1}] -> 1,
Band[{1, 2}, {n^2, n^2}] -> Table[If[i == n, 0, 1], {i, n}],
Band[{2, 1}, {n^2, n^2}] ->
Table[If[i == n, 0, 1], {i, n}]}, {n^2, n^2}]];
sqg1 =
AdjacencyGraph[s2 // ArrayFlatten, VertexLabels -> Automatic];
graphcoordsq1 =
GraphEmbedding[sqg1];(*Generating underlying lattice*)
hsquareFin1 = (s1 // ArrayFlatten);
egvecsq1 = Eigenvectors[hsquareFin1 // N // Normal];][[1]]; {nn,
t1, t2}]


Here t1 is code proposed by Shamina, and t2 is my code with hand made array. First we call

fcomp[100]

Out[]= {100, 115.673, 82.1941}


We see that computation times t1, t2 not so differ due to large Eigenvectors evaluation time. We also can compare matrixes

{hsquareFin1 == hsquareFin, sqg1 == sqg, egvecsq1 == egvecsq}

Out[]= {True, True, True}


Now we can prepare list for range 5<=Nunit<=40 and plot it separately and together with point Nunit=100

lst = Table[fcomp[i], {i, 5, 40}]

Out[]= {{5, 0.164432, 0.0043509}, {6, 0.205826, 0.0050914}, {7,
0.220041, 0.0064627}, {8, 0.253057, 0.0075893}, {9, 0.295401,
0.0123119}, {10, 0.319457, 0.0125374}, {11, 0.364707,
0.0142791}, {12, 0.384723, 0.016998}, {13, 0.429124,
0.0211146}, {14, 0.493359, 0.0247619}, {15, 0.513203,
0.0269879}, {16, 0.593006, 0.0354736}, {17, 0.648195,
0.0353599}, {18, 0.685608, 0.0465485}, {19, 0.743417,
0.0658679}, {20, 0.844621, 0.0619918}, {21, 0.890146,
0.0794076}, {22, 0.962821, 0.0800831}, {23, 1.04226, 0.108473}, {24,
1.15235, 0.111343}, {25, 1.2167, 0.129602}, {26, 1.32637,
0.158706}, {27, 1.41704, 0.162905}, {28, 1.50065, 0.182941}, {29,
1.67359, 0.218453}, {30, 1.73249, 0.219222}, {31, 1.92085,
0.252828}, {32, 2.05668, 0.287723}, {33, 2.26161, 0.313731}, {34,
2.35782, 0.356584}, {35, 2.57131, 0.425878}, {36, 2.74442,
0.454498}, {37, 2.99826, 0.524451}, {38, 3.26045, 0.578484}, {39,
3.5149, 0.66792}, {40, 3.75821, 0.66004}}


To visualize result we can use this code

fm[nn_] :=
Module[{Nunit = nn},
s1 = With[{n = Nunit},
SparseArray[{Band[{1, 1}] ->
Flatten[Table[Table[i, {i, n}], {n}]], Band[{1, n + 1}] -> 1,
Band[{n + 1, 1}] -> 1,
Band[{1, 2}, {n^2, n^2}] -> Table[If[i == n, 0, 1], {i, n}],
Band[{2, 1}, {n^2, n^2}] ->
Table[If[i == n, 0, 1], {i, n}]}, {n^2, n^2}]];
s2 = With[{n = Nunit},
SparseArray[{Band[{1, n + 1}] -> 1, Band[{n + 1, 1}] -> 1,
Band[{1, 2}, {n^2, n^2}] -> Table[If[i == n, 0, 1], {i, n}],
Band[{2, 1}, {n^2, n^2}] ->
Table[If[i == n, 0, 1], {i, n}]}, {n^2, n^2}]];
VertexLabels -> Automatic]; graphcoordsq1 = GraphEmbedding[sqg1];
hsquareFin1 = (s1 // ArrayFlatten);
egvecsq1 = Eigenvectors[hsquareFin1 // N // Normal]; {egvecsq1,
graphcoordsq1}]

With[{nn = 12}, fm[nn];
Manipulate[
Show[ListDensityPlot[
Table[Append[Chop@graphcoordsq1[[i]],
Chop@(Abs[egvecsq1[[j, i]]]^2)] // N, {i,
Length[egvecsq1[[j]]]}], PlotLegends -> Automatic,
Frame -> False, PlotRange -> {{0, nn - 1/2}, {0, nn - .5}, All},
ColorFunction -> Hue], sqg1], {j, 1, Length[egvecsq1], 1},
FrameLabel -> {"", "", "Probability density"}, ImageMargins -> 0,
FrameMargins -> Small, Frame -> False]]
`