This is an interesting topic, for sure, and there has been a lot of theoretical progress on it.
The main reference for the system currently used by OpSum is:
“Matrix Product Operators, Matrix Product States, and ab initio Density Matrix Renormalization Group algorithms” Chan, Keselman, Nakatani, Li, White
https://arxiv.org/abs/1605.02611.
For a brief summary of the construction there, I’m attaching the pdf notes below that I wrote for one of the Flatiron postdocs helping us to generalize that system. A key thing to realize about the method in that paper (or something I was initially confused about) is that one places all of the coefficients of the operators in the sum into a ‘bond matrix’ on each bond. But this representation is not simultaneous: these are different representations of the same operator, where one chooses a different bond to ‘expand around’. Then for each such representation, you compute the SVD of the coefficient matrix, then just the “V” matrix from the SVD is sufficient to compute an optimally compressed representation of the MPO, as detailed in the paper and notes.
Note that while the above algorithm gives optimal bond dimensions (up to some minor corner cases), it does not have optimal complexity for systems with arbitrary, long-range interactions, because the SVD of the coefficient matrix scales as the third power of the system size, (or square if you use a technique like randomized SVD), whereas it should be possible to find a general, linear-scaling algorithm.
Last but not least, Parker and Zaletel have developed a general theory of “canonical forms” of MPOs representing sums of operators (i.e. MPOs where there are special ‘settings’ of the bond indices that make strings of identity operators).
https://journals.aps.org/prb/abstract/10.1103/PhysRevB.102.035147
It is worth a read in tandem with the Chan, Keselman, et al. paper above because it deals with the concepts but in a slightly more general (and I think also more readable) way.
