Unexpected inefficiency when calculating matrix product by Dot

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[]}];
ind = Range[Length@base];
(*calculate the tag and sort it*)
unsortT = Map[(Dot[p, #]) &, base];
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[] // Flatten;
element = Select[state[], #[] != 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[]]}] // N;
tag = Map[(Dot[p, #[]]) &, element];
row = Extract[index2, Map[(Position[index1, #]) &, tag]] // Flatten;

(*tunneling Hamiltonian matrix in SparseArray form*)
];
(*Matrix element of whole tunneling Hamiltonian*)
tunnelingH[t_List, base_List, index_List] :=
Block[{length = Length@base[]},
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]];
];
(*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]
];

(*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.

• What do I need to execute to create hamil and initS from your code? – Marius Ladegård Meyer Jun 20 '15 at 16:34
• You can just run this code. The hamil and other variables will be created automatically. By the way, my problem is solved and it was just my mistake....... Thank you anyway~ – AchillesJJ Jun 20 '15 at 16:47
• Yeah, I saw Jens' answer just as I observed the huge imaginary part myself ;) Another comment to your code, if tf <= ti, then step will get the value Null, which is returned by Print inside the If. This is probably not what you want, as the NestList will fail ungracefully because its third argument is Null. – Marius Ladegård Meyer Jun 20 '15 at 17:05
• Might be more efficient to use MatrixExp instead of doing the explicit diagonalizations and time steps. – Daniel Lichtblau Jun 20 '15 at 23:17
• I have tried. MatrixExp costs about 45% more time than explicit diagonalization. Thank you. – AchillesJJ Jun 21 '15 at 3:12

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]]


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

MatrixExp[ -I dt hamil]
`