Tensor product in Pennylane

[AlaricCheng]

Hello,

I'm confused by the output of qml.operation.Tensor. It doesn't seem like a usual tensor product that I would expect. Surely, something like Tensor(qml.PauliX(0), qml.PauliX(1)) works perfectly well. But if the wires is not so ordered, then the output is wired. Below are several example I tries out.

• The output of Tensor(qml.PauliX(0), qml.PauliX(0)) is actually $X*X = I$.

• The output of Tensor(qml.PauliX(1), qml.PauliY(0)) is $X_0 \otimes Y_1$ instead of $Y_0 \otimes X_1$.

• The output of Tensor(qml.PauliX(0), qml.PauliY(1), qml.PauliX(0)) is $X_0 \otimes Y_1 \otimes X_2$,.

• But that of Tensor(qml.PauliX(1), qml.PauliY(1), qml.PauliX(0)) is $(X_0 * Y_0) \otimes X_1$.

From these examples, it seems that Tensor will reindex the wires of each observables, and if two consecutive observables have the same wire, then they will be first multiplied. I don't know whether the rules I summurized are complete.

It would be more clear if you could give a detailed description in the documentation, with explicit mathematical expressions.

But this 'wiring strategy' is not so intuitive. I think a better strategy would be to reorder the wire in the numerical order, and multiply the observable with the same wires. For example, (1) the output of Tensor(qml.PauliX(0), qml.PauliY(1)) and Tensor(qml.PauliY(1), qml.PauliX(0)) would be the same $X_0 \otimes Y_1$, as the wires in observables have already indicated the order; (2) the output of Tensor(qml.PauliX(1), qml.PauliY(0), qml.PauliX(2), qml.PauliZ(1)) would be $Y_0 \otimes X_1 Z_1 \otimes X_2$, where the $X_1$ and $Z_1$ are multiplied. I hope I have clarified the idea.

Thank you for your time and your good work :)

[josh146]

HI @AlaricCheng! This is a known issue with the Tensor class; during implementation, we made the assumptions that

• It would always be used with operators on different wires, and

• The operators would always be provided in strictly increasing wire order.

As a result, any deviation from this (multiple operators on the same wire, or decreasing non-strictly increasing wire order) will lead to undefined behaviour.

One option would simply be to include stricter input validation, so that only operator inputs satisfying the above conditions are allowed. However, we are planning to instead extend the behaviour of the tensor class to cover the cases described in your post.

[mariaschuld]

Just a comment: The first point, that operators act on different wires, should be taken care after we merge a refactor of the way PL handles wires soon. From then on, the Tensor class will complain if your observables act on the same wires.

I'm not sure if that also resolves the second point though.

[mariaschuld]

Maybe my reply in Issue #591 will also helps to illuminate the issue, @AlaricCheng? Basically, we cannot in future infer the wire ordering just from the wire labels, because users are soon allowed to define non-ascending or non-numeric wire labels...

But you are right, we need to see if we can "catch out" the most obvious cases where a user goes wrong.

[dwierichs]

Update on this (@josh146, @mariaschuld)
Given the custom wires policy that does not allow PL to infer a reasonable ordering of the wires, I only see the following two improvements to be made:

1. Tensor products with terms that act on the same (set of) wires should be multiplied independently from the ordering in the Tensor call or a multiplication via tensor1 @ tensor2 @ ...
   This can easily be achieved by introducing a sorted call into qml.operation.Tensor.matrix before using itertools.groupby:

    @property
    def matrix(self):
        # group the observables based on what wires they act on
        U_list = []
        key = lambda x: x.wires.labels # new line added
        _obs = sorted(self.obs, key=key) # new line added
        for _, g in itertools.groupby(self.obs, key): # changed line, recycled the key here
            # extract the matrices of each diagonalizing gate
            mats = [i.matrix for i in g]

            if len(mats) > 1:
                # multiply all unitaries together before appending
                mats = [multi_dot(mats)]

            # append diagonalizing unitary for specific wire to U_list
            U_list.append(mats[0])

        # Return the Hermitian matrix representing the observable
        # over the defined wires.
        return functools.reduce(np.kron, U_list)

Note that itertool's groupby "kind of" assumes the input to be sorted:

Generally, the iterable needs to already be sorted on the same key function.

2. If operations have partial overlap of the wires they act on, like Tensor(qml.Hermitian(A, [0, 1]), qml.PauliZ(0)), matrix should either multiply these correctly together as well (feels quite non-trivial to do in all cases with multiple wires overlaps) or it should throw an error, as matrix multiplication is not the main intention for this class. Note that with the current behaviour, the matrix does not have the shape (2**len(self.wires), 2**len(self.wires)) in such scenarios.

The advantage of these two changes is that they only happen in the matrix function, so that no overhead is added before calling matrix, which anyway is more expensive than the symbolic manipulations performed in the other functions/properties.

[mariaschuld]

@dwierichs I wonder if these planned changes will offer a more thorough way to solve the issue once and for all:

• We delete the Tensor class in favour of the global arithmetic operation of multiplication (which would be defined in __matmul__, so @ will be repurposed).
• When two operators are multiplied, a new operation is created that stores the two in decomposition. In other words, the new op has a representation as the product of other ops, just like templates do now.
• The matrix representation is getting a new kwarg wire_order with which one can specify the overall order of the wires. If it is not given, the wire order will be inferred according to a specific, but well-documented rule from the operations (in this case, it could just be using the wires of the first operation, concatenate those of the second that haven't appeared yet).

This logic would hopefully cover the intuition:

• PauliX(0) @ PauliZ(0) is a multiplication in the same subspace
• PauliX(0) @ PauliZ(1) is a tensor product of two different spaces, and we infer the default order [0, 1] for that "space"
• PauliX(1) @ PauliZ(0) is a tensor product of two different spaces, and we infer the default order [1, 0]
• Hermitian(A, [0, 1]) @ PauliZ(0) would understand PauliZ(0) as PauliZ(0) @ Identity(1) and multiply it to the Hermitian

I hope this will be a set of clear and consistent rules?

The deeper issue here is that how we use the tensor symbol in maths always implies that the order of chaining determines the subspaces. However, in our operator abstraction, the subspaces are determined by the wires. The Tensor class is a strange mix of the two, which leads to the issues I think.

[mariaschuld]

I am quite happy to add an error for the overlapping wires case in the meantime though! But I am not sure what the expected behaviour actually is for PauliX(1) @ PauliZ(0) -> the position contradicts the intuitive chronological sorting of wires...

[dwierichs]

Thanks for the pointer, @mariaschuld !
I guess you're right, the operator overhaul will likely cover this more elegantly!
Also, it indeed is not really clear what to do for PauliX(1)@PauliZ(0), I agree.
I would be in favor of adding some sort of warning, because there isn't really a scenario in which the output for partially overlapping wires is how one would expect it to be 😅

[albi3ro]

I believe this issue can also be closed with the replacement of Tensor with Prod.