Unexpected inefficiency when calculating matrix product by Dot

Question

There is an unexpected inefficiency when I calculate matrix product by Dot. The only relevant part of my code is presented as follows:

ClearAll[timeEvol];
timeEvol[hamil_, initS_List, ti_?NumericQ, tf_?NumericQ, dt_?NumericQ] :=
Block[{eigv, matD, step, evolOpt},
{eigv, matD} = Eigensystem[hamil];
matD = Transpose@matD;
step = If[tf > ti, Floor[(tf - ti)/dt], 
 Print["Final time shall be larger than initial time !"]];
evolOpt = matD.DiagonalMatrix[-I*dt*(Exp /@ eigv)].Inverse[matD];
NestList[(Normalize[evolOpt.#]) &, initS, step]
];

The function timeEvol has five variables,

hamil is a sparse matrix generated by SparseArray[] and its dimension is 496$\times$496;

initS is a 1-d list with dimension 496;

ti,tf,dt are positive real number, which are not important for my problem;

Now is the time to present my problem:

Using AbsoluteTiming for evaluation of expression, it shows that Mathematica costs about 130s to calculate the timeEvol function. Most of the time is spent on calculating the expression evolOpt, which is just a simple matrix product problem.(in fact, this part alone costs about over 129s..., calculating the Inverse[matD] and DiagonalMatrix[] just cost less than 0.1s)

It seems Mathematica is not good at matrix product... But the following test shows otherwise. The problem given above(more precisely, the evolOpt part)is equivalent with the following problem

Clear[testMat1, testMat2];
testMat1 = Table[RandomReal[], {500}, {500}];
testMat2 = DiagonalMatrix[Table[RandomReal[] + I*RandomReal[], {500}]];
AbsoluteTiming[testMat1.testMat2.Inverse@testMat1]

It just cost 0.07s...

So, what's the problem here ?

For further usage, I present the whole code below if you want to test on it. You just need to run this code and the output gives the time costed by each part of the expression in function timeEvol.

Thank you.

Needs["Parallel`"];

(*Generate a basis in lexicographic order*)
ClearAll[baseGenerator];
baseGenerator[m_Integer, n_Integer] :=
Reverse@Sort[
Join @@ Permutations /@ IntegerPartitions[n, {m}, Range[n, 0, -1]]]


(*Give index to the basis*)
ClearAll[indexT];
indexT[base_List] :=
Module[{p, ind, unsortT},
p = Table[Sqrt[100.*i + 3.], {i, 1, Length@base[[1]]}];
ind = Range[Length@base];
(*calculate the tag and sort it*)
unsortT = Map[(Dot[p, #]) &, base];
unsortT = Thread[List[unsortT, ind]];
SortBy[unsortT, First]
];


(*Generate Hamiltonian using Reap-Sow technique*)
(*Matrix element of tunneling term Subscript[a, \
i]^\[Dagger]Subscript[a, j]*)
ClearAll[tunnelingHij, tunnelingH];
tunnelingHij[i_Integer, j_Integer, tij_, base_List, index_List] :=
Block[{state, column, element, p, tag, row, index1, index2},

(*calculate Subscript[a, i]^\[Dagger]Subscript[a, j]|Subscript[n, 
1],Subscript[n, 2],...,Subscript[n, 
m]\[RightAngleBracket] and reap the index of non-zero terms*)
state = MapIndexed[{Sqrt[(#1[[i]] + 1)*#1[[j]]],
    If[#1[[j]] > 0,
     Sow[First@#2]; 
     ReplacePart[#1, {i -> #1[[i]] + 1, j -> #1[[j]] - 1}],
     #1 /. #1 -> 0]} &, base] // Reap;

(*column index and matrix element of all non-zero terms*)
column = state[[2]] // Flatten;
element = Select[state[[1]], #[[1]] != 0 &] // N;
index1 = index[[All, 1]];(*tags*)
index2 = index[[All, 2]];(*ind*)

(*calculate tags and find its corresponding row index*)
p = Table[Sqrt[100*k + 3], {k, 1, Length[base[[1]]]}] // N;
tag = Map[(Dot[p, #[[2]]]) &, element];
row = Extract[index2, Map[(Position[index1, #]) &, tag]] // Flatten;

(*tunneling Hamiltonian matrix in SparseArray form*)
Thread[Transpose[List[row, column]] -> tij*element[[All, 1]]]
];
(*Matrix element of whole tunneling Hamiltonian*)
tunnelingH[t_List, base_List, index_List] :=
Block[{length = Length@base[[1]]},
Flatten[Join[
 Table[
  tunnelingHij[i, j, t[[i, j]], base, index], {i, 1, 
   length - 1}, {j, i + 1, length}],
 Table[
  tunnelingHij[j, i, t[[j, i]], base, index], {i, 1, 
   length - 1}, {j, i + 1, length}]], 2]
];


(*Matrix element of density-density interaction*)
ClearAll[ddIntHii, ddInt];
(*Matrix element of d-d interaction Subscript[u, ii]Subscript[a, i]^\
\[Dagger]Subscript[a, i]^\[Dagger]Subscript[a, i]Subscript[a, i]*)
ddIntHii[i_Integer, uii_, base_List] :=
Block[{state, column, element},
state = Flatten[Last[
  MapIndexed[(If[#1[[i]] >= 2, 
       Sow[{First[#2], #1[[i]]*(#1[[i]] - 1)}]]) &, base] // 
   Reap], 1];
column = state[[All, 1]];
element = uii*state[[All, 2]];
Thread[Transpose[List[column, column]] -> element]
];
(*Matrix element of whole d-d interaction*)
ddInt[u_List, base_List] :=
Block[{},
Flatten[Table[ddIntHii[i, u[[i]], base], {i, 1, Length@base[[1]]}],
 1]
];

(*Matrix element of whole Bose-Hubbard Hamiltonian*)
(*m sites, n particle*)
ClearAll[bhHamil];
bhHamil[t_List, u_List, m_Integer, n_Integer] :=
Block[{base, index, dim},
base = baseGenerator[m, n];
index = indexT[base];
dim = (m + n - 1)!/((m - 1)!*n!);
SparseArray[
Join[tunnelingH[t, base, index], ddInt[u, base]], {dim, dim}]
]


(*time evolution*)
ClearAll[timeEvol];
timeEvol[hamil_, initS_List, ti_?NumericQ, tf_?NumericQ, 
dt_?NumericQ] :=
Block[{eigv, matD, step, evolOpt},
Print["time of evaluating EigenSystem is ",AbsoluteTiming[{eigv, matD} =      Eigensystem[hamil];]];
Print["time of evaluating Transpose is ",AbsoluteTiming[matD = Transpose@matD;]];
step = 
If[tf > ti, Floor[(tf - ti)/dt], 
 Print["Final time shall be larger than initial time !"]];
Print["time of evaluting evolOpt is ",AbsoluteTiming[
evolOpt = 
   matD.(DiagonalMatrix[-I*dt*(Exp /@ eigv)].Inverse[matD]);]];
Print["time of evaluating NestList is ",AbsoluteTiming[
 NestList[(Normalize[evolOpt.#]) &, initS, step];]]
];


(*Physical quantities calculation*)
(*Population of the i-th site*)
numSi[i_Integer, state_, base_] := 
Map[(Norm /@ #)^2 &, state].base[[All, i]];

clearAll[m, n];
m = 3;
n = 30;
base = baseGenerator[m, n];
index = indexT[base];
{t, u, initS} = Block[{dim = (m + n - 1)!/((m - 1)!*n!)},
{ConstantArray[-1., {m, m}],
ConstantArray[1., {m}],
Normalize@ConstantArray[1., {dim}]}];
hamil = bhHamil[t, u, m, n];
state = timeEvol[hamil, initS, 0, 1, 1.];

By the way, this code is used for Bose-Hubbard model with exact diagonalization algorithm.

Answer 1

There is a bug in your code which causes one of the factors in the Dot product to have excessively large imaginary parts (of the form 0. +1.32133*10^246 I). This means the multiplication can't be done in machine precision arithmetic, and as a result the calculation slows down considerably.

The mistake is that the exponentiation of the eigenvalues must read:

DiagonalMatrix[Exp[-I*dt*eigv]]

instead of

DiagonalMatrix[-I*dt*(Exp/@eigv)]

The latter contains real exponentials which cause the problem.

As a side note, the whole process of diagonalizing and exponentiation can be compressed into the single command

MatrixExp[ -I dt hamil]