1
$\begingroup$

I have following code. The function TraceSystem is to take partial trace of matrix and I think this part is fine. The main code is below in dd and I think the main problem is in here for Eigenvector. I tried to run L=10, and it took forever. People told me that python can run L=18 in 5 min. Is there anyway to improve the speed of this Mathematica code? Thanks!

 SwapParts[expr_, pos1_, pos2_] := 
 ReplacePart[#, #, {pos1, pos2}, {pos2, pos1}] &[expr]
 TraceSystem[D_, s_] := (
  
  Qubits = Reverse[Sort[s]];
  TrkM = D;
  
  z = (Dimensions[Qubits][[1]] + 1);
  
  For[q = 1, q < z, q++,
   n = Log[2, (Dimensions[TrkM][[1]])];
   M = TrkM;
   k = Qubits[[q]];
   If[k == n,
    TrkM = {};
    For[p = 1, p < 2^n + 1, p = p + 2,
     TrkM = 
       Append[TrkM, 
        Take[M[[p, All]], {1, 2^n, 2}] + 
         Take[M[[p + 1, All]], {2, 2^n, 2}]];
      ],
    For[j = 0, j < (n - k), j++,
      b = {0};
      For[i = 1, i < 2^n + 1, i++,
       If[(Mod[(IntegerDigits[i - 1, 2, n][[n]] + 
              IntegerDigits[i - 1, 2, n][[n - j - 1]]), 2]) == 1 && 
         Count[b, i]  == 0, 
        Permut = {i, (FromDigits[
             SwapParts[(IntegerDigits[i - 1, 2, n]), {n}, {n - j - 
                1}], 2] + 1)};
        b = 
         Append[b, (FromDigits[
             SwapParts[(IntegerDigits[i - 1, 2, n]), {n}, {n - j - 
                1}], 2] + 1)];
        c = Range[2^n];
        perm = 
         SwapParts[
          c, {i}, {(FromDigits[
              SwapParts[(IntegerDigits[i - 1, 2, n]), {n}, {n - j - 
                 1}], 2] + 1)}];
        
        M = M[[perm, perm]];
        
         ]    
       ]         ;
      TrkM = {};
      For[p = 1, p < 2^n + 1, p = p + 2,
       TrkM = 
         Append[TrkM, 
          Take[M[[p, All]], {1, 2^n, 2}] + 
           Take[M[[p + 1, All]], {2, 2^n, 2}]];
       ]
         ];
    ]
   
   ]
  
  ; Return[TrkM])

dd = Table[Flatten[{h, L = 10;
      H = -KroneckerProduct[
         KroneckerProduct[PauliMatrix[3], IdentityMatrix[2^(L - 2)]], 
         PauliMatrix[3]];
      For[i = 1, i <= L - 1, i++, 
       H = H - KroneckerProduct[IdentityMatrix[2^(i - 1)], 
          KroneckerProduct[
           KroneckerProduct[PauliMatrix[3], PauliMatrix[3]], 
           IdentityMatrix[2^(L - i - 1)]]]];
      For[j = 1, j <= L, j++, 
       H = H - h*
          KroneckerProduct[IdentityMatrix[2^(j - 1)], 
           KroneckerProduct[PauliMatrix[1], 
            IdentityMatrix[2^(L - j)]]]];
      H = 
       N[SparseArray[
         H - 0.01*
           KroneckerProduct[PauliMatrix[3], 
            IdentityMatrix[2^(L - 1)]]]];
      LES = 
       Flatten[Eigenvectors[-H, 1, 
         Method -> {"Arnoldi", "Criteria" -> "RealPart"}]]; 
      DenM = KroneckerProduct[LES, LES]; 
      RedDenM = N[SparseArray[TraceSystem[DenM, {1, 2}]]];
      aa = Chop[Eigenvalues[RedDenM]];
      Tr[(1/aa[[1]])*aa]}], {h, 0, 3, 0.01}]; // AbsoluteTiming
ListPlot[Thread[{First@#, Rest@#}] & /@ dd, PlotRange -> Full, 
 AxesLabel -> {h, Z}]
$\endgroup$
2
  • 3
    $\begingroup$ You posted a wall of code with no comments or explanation. Unless the problem is glaringly obvious, it's unlikely that people will read through your code, understand what it is supposed to do, then try to come up with a solution for you. Start by measuring timings and finding the bottleneck, then reduce down the code to the problematic section (but still include all definitions needed to run! i.e. make an MWE). Explain what you were trying to do and how. $\endgroup$ – MarcoB Dec 2 '20 at 14:54
  • 3
    $\begingroup$ Also look into switching away from For loops and Append to more idiomatic approaches (see Alternatives to procedural loops and iterating over lists in Mathematica). $\endgroup$ – MarcoB Dec 2 '20 at 14:54