expval ignores wire labels and uses position implicitly for tensor product ops

[Jaybsoni]

Context:
Consider the following circuit:

dev = qml.device("default.qubit", wires=3)

@qml.qnode(dev)
def circ():
    qml.Hadamard(wires=0)
    qml.CNOT(wires=[0, 2])
    return qml.expval(qml.Identity(wires=0) @ qml.PauliZ(wires=1))      # -->  0.9999999999999998
    # return qml.expval(qml.Identity(wires=1) @ qml.PauliZ(wires=0))   # -->  0.9999999999999998
    # return qml.expval(qml.PauliZ(wires=0) @ qml.Identity(wires=1))   # -->  0.0
    # return qml.expval(qml.PauliZ(wires=1) @ qml.Identity(wires=0))   # -->  0.0

But the expected behaviour should be:

... 
    return qml.expval(qml.Identity(wires=0) @ qml.PauliZ(wires=1))      # -->  0.9999999999999998
    # return qml.expval(qml.Identity(wires=1) @ qml.PauliZ(wires=0))   # -->  0.0
    # return qml.expval(qml.PauliZ(wires=0) @ qml.Identity(wires=1))   # -->  0.0
    # return qml.expval(qml.PauliZ(wires=1) @ qml.Identity(wires=0))   # -->  0.9999999999999998

It seems as though the wire labels are being ignored and instead the position of the operations in the tensor product are being used to determine how to compute the expectation value.

Additionally:
This PR also addresses the following bug which was discovered during the testing of this feature:

    dev = qml.device("default.qubit", wires=3)
    theta = 0.432
    phi = 0.123
    varphi = -0.543
 
    @qml.qnode(dev)
    def circuit(o):
        qfunc()
        return qml.expval(qml.apply(o))

    obs0 = qml.PauliZ(wires=2)

    obs1 = qml.Identity(wires=0) @ qml.Identity(wires=1) @ qml.PauliZ(wires=2)
    obs2 = qml.Identity(wires=1) @ qml.Identity(wires=0) @ qml.PauliZ(wires=2)

    obs3 = qml.Identity(wires=0) @ qml.PauliZ(wires=2) @ qml.Identity(wires=1)
    obs4 = qml.Identity(wires=1) @ qml.PauliZ(wires=2) @ qml.Identity(wires=0)

    obs5 = qml.PauliZ(wires=2) @ qml.Identity(wires=0) @ qml.Identity(wires=1)
    obs6 = qml.PauliZ(wires=2) @ qml.Identity(wires=1) @ qml.Identity(wires=0)

    obs_lst = [obs0, obs1, obs2, obs3, obs4, obs5, obs6]

    for i, o in enumerate(obs_lst):
        print(f"------ {i} ------- \n{circuit(o)}")

Which produces the following result:

------ 0 -------
0.8561624160163042
------ 1 -------
0.8561624160163042
------ 2 -------
0.8561624160163042
------ 3 -------
0.8561624160163042
------ 4 -------
0.9081301906949608
------ 5 -------
0.9924450321351934
------ 6 -------
0.8561624160163042

In the above result we notice that the expectation values computed for observables 4, 5 do not match the others.

Description of the Change:
The expval(self, observable, shot_range=None, bin_size=None) function defined in _qubit_device.py uses the eigenvalues of the provided tensor product observable and the probability vector of the state to compute the expectation value when shots=None.

The first bug arises because of lexicographic sorting of the eigvals based on the wires the observables act on in the function qml.Tensor.get_eigvals(). This results in the eigvals of the two observables PauliZ @ Identity and Identity @ PauliZ being the same.

This is addressed by removing the sorting in the get_eigvals() function.

The second bug arises because of the way we permute the probability vector before computing the expectation value. The wire order of the observables is passed into the device.probability() function which then computes the marginal probability distribution over the observable wires and permutes the distribution according to the wire order of the observable wires. This permutation is incorrect and instead we should be using the inverse of this permutation. The reason such a bug was not caught sooner is because in most cases the tensor product observables usually contain 2 terms at most and/or there is only one qubit flipped out of position. In these common cases, the wire order permutation is its own inverse, and thus produces the correct expected value. But in situations with more then one qubit flip, it breaks.

This is addressed by computing the permutation before hand, then permuting the wires before they are passed into the probability function.

Related GitHub Issues:
This PR closes this issue

[github-actions]

Hello. You may have forgotten to update the changelog!
Please edit doc/releases/changelog-dev.md with:

• A one-to-two sentence description of the change. You may include a small working example for new features.
• A link back to this PR.
• Your name (or GitHub username) in the contributors section.

[Jaybsoni]

[sc-15486]

[Jaybsoni]

The reason some of the tests are failing is because of the following edge case:

A = np.array( ... some hermitian 8x8 matrix ... )
H = qml.Hermitian(A, wires=[2,1,0])

def circ(): 
... some ops ...
    return qml.expval(H)   # now that we sort the wires for the probability vector but not for the
                           # Eigen_values, we get the "wrong" answer! 

The answer would be "correct" if we interpreted the matrix as already being transformed under the provided wire order.

[antalszava]

Looks good 🥇 Amazing solution, goes a long way with understanding what's happening here.

Codefactor seems to cause some issues still and left come comments, but once those are addressed, this is good to be merged from my side.