Multi-Site Expectation Value

[Umberto Borla]

I am developing the method that, given a MPS in canonical form, calculates the expectation value of multi-site operators that are tensor product of one-site operators. This is what I implemented so far: it works but does not handle operator strings (like the Jordan-Wigner string that is needed in the presence of fermionic operators) yet.

Python: Select all

    def expectation_value_multi_site(self, ops, sites, opstr=None, str_on_first=True):
        r"""Expectation value  ``<psi|op_1 op_2 ... op_n |psi>/<psi|psi>`` of
        tensor products of single site operators.

        Given the MPS in canonical form, it calculates the expectation value of
        a tensor product of single-site operators.
        For example the contraction of three one-site operators on sites `i`
        `i+1` `i+2` would look like::

            |          .--S--B[i]---B[i+1]--B[i+2]--.
            |          |     |      |       |       |
            |          |    op1    op2     op3      |
            |          |     |      |       |       |
            |          |     |      |       |       |
            |          .--S--B*[i]--B*[i+1]-B*[i+2]-.
    Parameters
            ----------
            ops : List of { :class:`~tenpy.linalg.np_conserved.Array` | str }
                sList of one-site operators. This method calculates the
                expectation value of the n-sites operator given by their tensor
                product.
            sites : List of int
                List of site indices.
                ``ops[x]`` acts on site ``sites[x]``.
                Is sorted before use, i.e. the order is ignored.
            opstr : None | (list of) { :class:`~tenpy.linalg.np_conserved.Array` | str }
                Ignored by default (``None``).
        """
        ops = npc.to_iterable_arrays(ops)
        sites = np.sort(sites)

        initial_site = sites[0]
        m = len(sites)
        op = self.get_op(ops, 0)
        theta = self.get_theta(initial_site, 1)
        C = npc.tensordot(op, theta, axes=['p*', 'p0'])
        C = npc.tensordot(theta.conj(), C, axes = [['p0*', 'vL*'],['p', 'vL']])
        for i in range(m-1):
            B = self.get_B(initial_site+i+1, form ='B')
            C = npc.tensordot(C, B,  axes=['vR', 'vL'])
            op = self.get_op(ops, i+1)
            C = npc.tensordot(op, C, axes=['p*','p'])
            if i == m-2:
                C_exp = npc.inner(B.conj(), C, axes=[['vL*','p*','vR*'],['vR*','p','vR']])
                break
            else:
                C = npc.tensordot(B.conj(), C, axes=[['vL*','p*'],['vR*','p']])
        return np.real_if_close(C_exp)

I am also looking for suggestions on what to do with the "sites" argument. Right now it is not very useful as it only indicates the consecutive sites where the operators are, but it is still necessary to specify all the operators between the initial and final site, and put identities if we don't want anything to happen on that site. One option would be to have the operators [op_1.. op_n] on the sites [site_1... site_n], with the sites being ordered but not consecutive, and have an identity automatically inserted on the unspecified sites. What do you think?

[Leon]

What exactly is the difference between this method and the native MPS.expectation_value()? The latter can already take in a sequence of operators.

[Johannes]

It takes multiple onsite operators acting on n neighboring sites. To do that with MPS.expecation_value(), you need to generate a full n-site operator before by Hand, which becomes impossible for n > 15. This function contracts in a different order than the MPS.expectation_value(), so it stays "efficient" if even if you have a product of 20 onsite operators to be measured. Regarding the structure of the internal contraction, it's more like the MPS.correlation_function() with a single "site1", but allowing for more than two sites to act on with non-trivial operators.

[mention]Umberto Borla[/mention] You say "``ops[x]`` acts on site ``sites[x]``", where x is the list index in the given arguments. For that reason it's a bad idea to sort the sites, if you don't apply the same reordering to ops, otherwise it's very confusing.
It would be convenient, if one can give arguments ops=['Sx', 'Sz', 'Sy'], sites=[5, 2, 3] (which I think is what you tried to say in the documentation). You can then zip those and sort them together:

Python: Select all

op_i = list(zip(ops, sites)
op_i.sort(key=lambda tup: tup[1])

I think it would be even better to already give the list in a "zipped" form, i.e., as ops=[('Sx', 5), ('Sz', 2), ('Sy',3)] (which is exactly what the list(zip(...)) gives).
Giving the argument ops=[('Sx', 5), ('Sz', 2), ('Sy',3)] should be equivalent to ops=[('Sx', 5), ('Sy',3), ('Sz', 2)].
EDIT: Exception: for fermionic "Cd", "C" operators, changing the order should possibly give an overall minus sign.

Thinking more generally, what we really need (also for a project I'm working on) is a further generalization of this function for measuring many terms and summing them up, based on your function looking like

Python: Select all

def expectation_value_terms_sum(self, operators, prefactors):
    """Calculate expectation values for a bunch of terms and sum them up.

    Parameters
    ----------
    operators : list of terms
        Each `term` should have the form ``[(Op1, site1), (Op2, site2), ...]``.
    prefactors : list of (complex) floats
        Prefactors for the ``total_sum`` to be evaluated.

    Returns
    -------
    terms : list of (complex) float
        ``terms[i]`` contains the value for ``self.expectation_value_multi_site(operators[i])``.
    terms_sum : list of (complex) float
        Equivalent to ``sum([t*f for t, f in zip(terms, prefactors)])``.

    See also
    --------
    :meth:`expectation_value_multi_site`: evaluates a single `term`.
    """
    # TODO: this is a naive implementation. For efficiency, a caching is necessary:
    # We need to evaluate the expectation values of terms starting with the same operator(s)
    # from the left at the same time to avoid doing the same contractions many times.
    terms = [self.expectation_value_multi_site(ops) for ops in operators]
    terms_sum = sum([t*f for t, f in zip(terms, prefactors)])
    return terms, terms_sum

As a side note, I think it would be great if we automatically insert the necessary Jordan Wigner strings.
We should also do this in MPS.correlation_function() ...

[Umberto Borla]

I implemented this function that takes a term containing operators acting on non-neighbouring sites and returns a list of operators, obtained by filling the intermediate sites with "Id" or "JW" where appropriate. The output can be used by the method expectation_value_multi_sites

Python: Select all

    def _term_to_jw(self, term, op_string = "JW"):
        """Given a term that is a tensor product of  several operators acting on
        possibly non-neighboring sites [i_min, ... i_n, ... , i_max],
        this puts Jordan-Wigner strings where appropriate on the remaining
        intermediate sites.

        Parameters
        ----------
        term : list of (str, int)
        List of tuples ``op, i`` where `i` is the site `op` acts on, and `op`
        is given as a string naming an onsite operator.

        Returns
        -------
        modified_term : list of str
            List that contains operators (given as strings) acting on all
            sites from i_min to i_max.
        --------
        """
        i_min = min([term[i][1] for i in range(len(term))])
        i_max = max([term[i][1] for i in range(len(term))])
        new_term = [None for i in range(i_max-i_min+1)]
        for op, i in term[0:]:
            if new_term[i-i_min] is not None:
                new_term[i-i_min] = new_term[i-i_min] + " " + op
            else:
                new_term[i-i_min] = op
            if op in self.sites[i].need_JW_string:
                for j in range(i-1):
                    if new_term[j] is not None:
                            new_term[j] = new_term[j] + " " + op_string
                    else:
                        new_term[j] = op_string
        new_term = ["Id" if op is None else op for op in new_term]
        new_term = [op.replace(op_string + " " + op_string, "").rstrip().lstrip() for op in new_term]
        return new_term

The limitations right now are that it can only be given operators in the form of strings (their names). It can handle string operators different from JW, but always assumes that they square to one. Please let me know if you think I have missed something

[Johannes]

Thanks!
This function does more than just inserting Jordan wigner strings: it maps a representation of a term in the form [('N', 1), ('C', 5), ('Cd', 3)] into a form having operators on each site for some range, ['N', 'Id', 'JW Cd', 'JW', 'C'].

Some minor issues:

• I think it should be for j in range(i-1-i_min):
• use the function Site.op_needs_JW(op) instead of checking op in Site.op_need_JW, then it will also handle operators which are already put together (like "Cd C")
• It's probably safer to replace " JW JW" with "" even before inserting the second "JW". Right now your function might end up with empty strings at some places (e.g. in the given example above)

Do you want to make a pull request for the git repository to clearly mark the changes as your development? Or should I just upload the changes?

[Umberto Borla]

Hi! This is the function modified according to your useful suggestions:

Python: Select all

    def _term_to_jw(self, term, op_string = "JW"):
        """Given a term that is a tensor product of  several operators acting on
        possibly non-neighboring sites [i_min, ... i_n, ... , i_max],
        this puts Jordan-Wigner strings where appropriate on the remaining
        intermediate sites, collects the operators acting on the same site
        in the correct order and returns a list of operators (as strings) that can be used
        by the expectation_value_multi_site method.

        Parameters
        ----------
        term : list of (str, int)
        List of tuples ``op, i`` where `i` is the site `op` acts on, and `op`
        is given as a string naming an onsite operator.

        Returns
        -------
        modified_term : list of str
            List that contains operators (given as strings) acting on all
            sites from i_min to i_max.
        --------
        """
        i_min = min([term[i][1] for i in range(len(term))])
        i_max = max([term[i][1] for i in range(len(term))])
        new_term = [None for i in range(i_max-i_min+1)]
        for op, i in term[0:]:
            if new_term[i-i_min] is not None:
                new_term[i-i_min] = new_term[i-i_min] + " " + op
            else:
                new_term[i-i_min] = op
            if self.sites[i-i_min].op_needs_JW(op):
                for j in range(i-i_min):
                    if new_term[j] is not None:
                            new_term[j] = new_term[j] + " " + op_string
                    else:
                        new_term[j] = op_string
        new_term = ["Id" if op is None else op for op in new_term]
        new_term = [op.replace(op_string + " " + op_string, "").rstrip().lstrip() for op in new_term]
        new_term = ["Id" if not op else op for op in new_term]
        return new_term

• I updated the documentation, but feel free to rephrase it if you find a better way to explain what the function does.
• It was indeed for j in range(i-i_min)
• I did not find an easier way to replace "JW_JW" with empty strings. To circumvent the problem you pointed out, I added a line in the end that checks for empty strings and replaces them with "Id".

I tested it for several cases and everything seems to work as expected. I don't mind if you upload the changes yourself.

[Johannes]

thank you very much, I've now included these new MPS methods into the master branch.