I am interested in calculating correlation functions of the form \langle \psi |O_i O_j |\psi\rangle. A naive way to do this is to create the state |\psi_{ij}\rangle = O_iO_j|\psi\rangle using the apply() function, and to then take inner(psi',psi). This however takes much longer than running correlation_matrix in the case where O_{i,j} are single-site operators, which makes me wonder why the above procedure for calculating the correlators is inefficient (unfortunately I was not able to completely grok the source code for correlation_matrix). Things can be sped up by moving the orthogonality center to sites i and j in between applications of apply, but the result is still several times slower than directly calling correlation_matrix.
Is there a way I can get similarly fast correlators without using correlation_matrix? I ask because in my applications, O_i are multi-site operators (e.g. O_i = S^+_i S^-_{i+1}), which (to my knowledge) correlation_matrix does not support.
It’s a good question about why correlation_matrix is faster than using MPOs, which is something we find users are often wanting to do (most likely because MPOs are easy to make using our OpSum system and can handle very general operators).
The reason for the speed of correlation_matrix is that is uses the algorithm described here. We used to ask users to code this algorithm themselves, which is why we have that page in the C++ documentation. If you want that speed for more general operators, such as two-site operators, we’d ask you to code that yourself as we have no plans to add that ability to the correlation_matrix function (the reason being that at some point we need to limit the scope of what that function does, otherwise it would become a “measure anything people want” function and would become too complex - it’s already rather complex & I wasn’t surprised you found it hard to read).
In words, the reason the algorithm linked above is faster is that it doesn’t do a full pass over the MPS for every value of i and j. Instead it does a kind of nested double for loop where for each value of i it does the minimal computations needed to get the correlator for all j > i.