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.
hamil
andinitS
from your code? $\endgroup$hamil
and other variables will be created automatically. By the way, my problem is solved and it was just my mistake....... Thank you anyway~ $\endgroup$tf <= ti
, thenstep
will get the valueNull
, which is returned byPrint
inside theIf
. This is probably not what you want, as theNestList
will fail ungracefully because its third argument isNull
. $\endgroup$MatrixExp
instead of doing the explicit diagonalizations and time steps. $\endgroup$MatrixExp
costs about 45% more time than explicit diagonalization. Thank you. $\endgroup$