On the naive contraction implementation

Hello everyone,
In the current implementation of the naive contraction algorithm of MPS with MPO, a new MPS is computed from contracting each tensor of the MPS with the corresponding tensor of the MPO, then this new state is left orthogonalized and after that truncated sweeping from right to left (actually the second part is true only if truncate argument is true, but I think than when it is false the state should be right ortoghonalized in order to have always the same gauge in output from a contraction, just for clean programming).
Computing the new MPS can be memory expensive if the link indices of the MPO have a significant size. I was thinking to change the implementation in favor to a new one that works as follow:

  • Compute the tensor T_2 = \psi [1]* A[1]
  • for j = 2:N
  • T_1= T_2
  • T_2 = \psi[j]*A[j]
  • Combine the link indices
  • \psi[j-1], \psi[j] = T1,T2
  • orthogonalize!(\psi,j)
  • T_2 = \psi[j] #take the modified tensor back
  • end of the for block

This should be equivalent in terms of stability, since at each iteration I’m going to orthogonalize on only the indices that have already been contracted while leaving untouched the non contracted ones. But since for simply algebraic reasons, \chi_j<\min\{2^{j+1}, 2^{N-j}\} performing the orthogonalization at each step instead that at the end would allow for memory saving (unfortunately the computation could take longer, but the difference in time is much smaller than the difference in memory, especially for large bonds of the MPS and of the MPO).

NOTE: As described, the resulting MPS will be right orthogonal, which would be ideal if we are going to start a truncation after, but if we are not, then the procedure should start from the last tensor of the MPS, obtaining at the end a right orthogonalized MPS