Tensor operations

[quantshah]

In case of multi-partite quantum systems, we would ideally like to define operators acting on the larger Hilbert space for more than one qubit. I am starting this thread to discuss this and possibly implement it in Pennylane if it is interesting.

With the addition of new features following #175 , we can now calculate expectation values of multi-wire operators by defining the Hermitian matrix corresponding to the operator in the larger Hilbert space as pennylane.expval.Hermitian(operator_matrix, wires=[0, 1]).

• Such operators could be constructed as matrices from single qubit operations e.g.,

  operator_matrix = tensor(PauliX, PauliZ)
  pennylane.expval.Hermitian(operator_matrix, wires=[0, 1])

• Or, we could have a something like pennnylane.expval.Tensor(PauliX(0), PauliY(2))

I am guessing that a hardware device should also be able to measure/apply such multiwire operators? Then, do we need something like a TensoredQNode to deal with such operators?

Josh explained a trick where we can do such measurements by converting the basis to a |+> just before measurement. This is discussed in this StackExchange as trading off entangling measurement with entangling gates:

https://quantumcomputing.stackexchange.com/questions/1588/do-multi-qubit-measurements-make-a-difference-in-quantum-circuits

In the default_qubit backend, we then need to implement this trick perhaps which detects that the measurement is multi-wire and uses an entangling gate operation before local measurement to overcome the lack of a joint measurement, hence a TensoredQnode.

Let me know your comments or suggestions. @josh146 @co9olguy

[josh146]

Or, we could have a something like pennnylane.expval.Tensor(PauliX(0), PauliY(2))

I really like this syntax! We could even overload the __add__/__mul__ method, so that PauliX(0) + PauliY(2) or PauliX(0) * PauliY(2) is equivalent to qml.expval.Tensor(PauliX(0), PauliY(2)) :laughing:

I am guessing that a hardware device should also be able to measure/apply such multiwire operators? Then, do we need something like a TensoredQNode to deal with such operators?

This is correct. The way to think about the process is as follows:

1. Make a PR to PennyLane that adds the functionality to operations.py and expectatations.py, and a way for the Tensor operation/expectation to 'store' the names of the operations/wires that are tensored as a string/int.

   This information will be passed by the QNode to the device - no matrices should be passed. If the user wants to use explicit matrices, they can use the existing qml.expval.Hermitian.

   When this is merged, it will not work with any existing plugins, but this is okay - they will just raise a Not supported operation exception.

2. Add support to the plugins. Since all the frameworks differ somehow, the way this should be implemented will differ as well. Some frameworks have explicit tensor product support that we would simply use, others are may not, and we will have to map it to an explicit matrix within the plugin.

   Step 2 should probably be done in two separate steps:

   • Add support to default.qubit (not sure how this would look for default.gaussian/CV for now). This can probably just be done by mapping to explicit matrices inside the plugin (as we are already doing), and just using np.kron if the operation is a Tensor.

     This will also help us iron out any hidden bugs.

   • Finally, add support for the other external plugins.

[smite]

If you want to use operator overloading to implement the tensoring operation, I'd suggest __mul__, because it seems more appropriate than __add__ (direct sum?) or __matmul__.

[smite]

A few more thoughts:

The syntax Tensor(PauliX, PauliY, wires=[7,4]), which otherwise looks nicer to me, will run into trouble if you want to tensor parametrized Expectations. Maybe you could use functools.partial to handle the parameters but it's a bit complicated.

What about tensoring multiqubit expectations? Tensor(PauliX, TwoQubitExpval, wires=[7, [1,4]])?

On the other hand for Tensor(PauliX(wires=[7]), PauliY(wires=[4])) to work, the two Pauli expectation instances must not be queued into the QNode.ev queue on their own when they are constructed. You'd also need to do additional checking to make sure the tensor product is valid, i.e. the operations act on non-overlapping sets of wires.

[josh146]

The syntax Tensor(PauliX, PauliY, wires=[7,4]), which otherwise looks nicer to me, will run into trouble if you want to tensor parametrized Expectations. Maybe you could use functools.partial to handle the parameters but it's a bit complicated.

One solution would be to restrict this to PauliX, PauliY, PauliZ, and Identity for now.

On the other hand for Tensor(PauliX(wires=[7]), PauliY(wires=[4])) to work, the two Pauli expectation instances must not be queued into the QNode.ev queue on their own when they are constructed. You'd also need to do additional checking to make sure the tensor product is valid, i.e. the operations act on non-overlapping sets of wires.

I assume this logic can be taken into account in the TensorProduct(Operations) subclass?

[quantshah]

I really like this syntax! We could even overload the add/mul method, so that PauliX(0) + PauliY(2) or PauliX(0) * PauliY(2) is equivalent to qml.expval.Tensor(PauliX(0), PauliY(2))

• Could this be confused with a multiplication between the matrices denoting the operators? We can start with the list syntax and then overload the operators.

On the other hand for Tensor(PauliX(wires=[7]), PauliY(wires=[4])) to work, the two Pauli expectation instances must not be queued into the QNode.ev queue on their own when they are constructed. You'd also need to do additional checking to make sure the tensor product is valid, i.e. the operations act on non-overlapping sets of wires.

• I agree with @josh146 on this that the TensorProduct(operations) subclass can do these checks but what if someone does this?

def qfunc1(x):
    qml.RY(x, wires=1)
    op1 = qml.PauliX(wires=1)
    op2 = qml.PauliY(wires=2)
    return qml.expval.Tensor(op1, op2)

• Will the op1 and op2 be already queued? In this case, we do not even reach the qml.expval.Tensor class so if we wish to prevent these ops from getting queued then we need to somehow check that they are not called directly but assigned to a variable op1 and then prevent them from being queued. I am trying to understand when QNode queues the operations. In pennylane.qnode, we have this function which seems to do the queuing,

    def _append_op(self, op):
        """Appends a quantum operation into the circuit queue.

        Args:
            op (:class:`~.operation.Operation`): quantum operation to be added to the circuit
        """
        # EVs go to their own, temporary queue
        if isinstance(op, pennylane.operation.Expectation):
            self.ev.append(op)

So, this prevents anything which is not a Operationfrom getting queued. E.g.,

In [46]: def qfunc(x):
    ...:     qml.RY(x, wires=1)
    ...:     op1 = qml.PauliX(wires=1)
    ...:     params = 2.3
    ...:     return qml.expval.PauliZ(0)

In [47]: node = QNode(qfunc, dev1)
In [48]: node.ops
Out[49]: []
In [50]: node(1)
In [51]: node.ops
Out[51]:
[<pennylane.ops.qubit.RY at 0x127ab0ac8>,
 <pennylane.ops.qubit.PauliX at 0x127ab02e8>,
 <pennylane.expval.qubit.PauliZ at 0x127ab06d8>]

So this is a bit tricky again and needs more discussing.

[quantshah]

The syntax Tensor(PauliX, PauliY, wires=[7,4]), which otherwise looks nicer to me, will run into trouble if you want to tensor parametrized Expectations. Maybe you could use functools.partial to handle the parameters but it's a bit complicated.

I agree. I like Tensor(Gate1(params, wires), Gate2(params, wires)) better for the long run.

[josh146]

Could this be confused with a multiplication between the matrices denoting the operators? We can start with the list syntax and then overload the operators.

I did get a bit overexcited - let's leave operator overloading out for now.

I am trying to understand when QNode queues the operations.

The queuing actually takes place inside the Operation class initialization. See line 297 of operation.py; Operation.init() calls self.queue:

def queue(self):
    """Append the operation to a QNode queue."""
    if QNode._current_context is None:
        raise QuantumFunctionError("Quantum operations can only be used inside a qfunc.")

    QNode._current_context._append_op(self)
    # return self so pre-constructed Expectation instances can be queued
    # and returned in a single statement
    return self  

Because queueing takes place on operation initialization, this doesn't allow a pattern like

op = qml.PauliX(wires=0)

because the queuing will take place on initialization/assignment.

Instead, you could do this (although maybe this is a bit pointless):

op = qml.PauliX
# later on
op(wires=0)

as the initialization is delayed.

---

Taking the above into account, the following would work trvially assuming non-parametrized expectation values as arguments:

return qml.Tensor(qml.expval.PauliX, qml.expval.PauliY, wires=[0, 1])

If we instead want to go with the following:

return qml.Tensor(qml.expval.PauliX(wires=0), qml.expval.PauliY(wires=1))

this would need us to refactor/rethink how operation/expectation queueing works.

---

A final thought; is the following too weird, using @smite's idea of nested wire lists, and where params is an optional argument? I'm thinking probably, but it would be doable with the way queuing works now.

H = qml.expval.Hermitian
Y = qml.expval.PauliY
return qml.Tensor(H, Y, wires=[[0, 1], 2], params=(A, None))

Ideally, the best solution is to get the following to work:

return qml.Tensor(H(A, wires=[0, 1]), Y(0))

[quantshah]

I like this, and will try to formulate what we need to get this working if everyone agrees with this syntax.

return qml.expval.Tensor(H(A, wires=[0, 1]), Y(0))

We could assume that

op1 = qml.expval.PauliY(0)

will always queue the operation. It is also not clear why this should not queue the operation anyways.

[smite]

• In principle you could explicitly turn all queuing off using e.g. an Operation class variable, construct the Operations, call Tensor(), and then turn queuing back on, but this would be ugly and error prone so I would advise against it.

• Another not-so-awesome option is to set the do_queue argument of Operation.__init__() to False, which prevents the queuing of that single Operation, but this too is pretty clumsy and easy to forget, which leads into hard-to-catch bugs.

• Yet another way to temporarily turn off the queuing would be to use a special context manager:

  with qml.noqueue: qml.Tensor(RY(0.3, wires=[3]), RY(-6.1, wires=[2]))

  Now the user cannot forget to turn queuing off or back on, but if they construct the gates and assign them into temporary variables outside the context it won't help.

• Finally

  H = qml.expval.Hermitian Y = qml.expval.PauliY return qml.Tensor(H, Y, wires=[[0, 1], 2], params=(A, None))

  is pretty robust against user errors (it leaves the construction of the Operations to the Tensor class or function which can "safely" manipulate the queuing behavior), but the syntax itself may be a bit confusing. Anyway, with the way the QNodes currently work, I think this one makes most sense.

[co9olguy]

I haven't thought about this enough atm to hold any strong opinions, but syntactically it seems less than ideal to have Tensor(expval1, expval1), since the tensor appers outside the expectation values (at least it appears this way upon reading).

The whole point of naming it expval is to reinforce to the users that they will be taking expectation values, not single-shot measurements. Something cleaner would have the form qml.expval(Tensor(op1, op2))

Since we might implement (non-differentiable) measurements sometime soon, maybe it makes sense to think about this as part of a broader context for measurement.

Currently, based on the above, I'm thinking about something like:

def circuit(...):
   Gate(...)
   return qml.expval(op1), qml.expval(Tensor(op2, op3))

The main re-abstraction here is that expval becomes a separate function, which takes in measurement operators (or tensor products). We could include single-shot measurements by having a companion function qml.measure

return qml.measure(op1), qml.Tensor(op2, op3)

[josh146]

The main re-abstraction here is that expval becomes a separate function, which takes in measurement operators (or tensor products). We could include single-shot measurements by having a companion function qml.measure

return qml.measure(op1), qml.Tensor(op2, op3)

Amazingly, this is exactly the same proposal @quantshah and I discussed this morning.

This has a lot of advantages:

1. It makes a lot more sense. We would have a separate submodules qml.ops (for unitaries) and qml.obs (for Hermitian observables). Then, you either specify qml.expval(obs1), qml.var(obs1), or qml.measure(obs1)

2. It makes the line

   return qml.expval(qml.tensor(obs1(wires=1), obs2(A, wires=2)))

   much more intuitive (can read it directly as 'the expectation value of the tensor product of obs1 and obs2').

3. It gets rid of a lot of code duplication

4. It will make the documentation a lot cleaner.

[josh146]

Note: this issue will be addressed once #226 has been addressed.

[co9olguy]

With the refactor of how observables work (#232), we should be able to implement tensor product observables with a simple return qml.expval(op1(wire1), op3(wire3)), qml.expval(op2(wire2)) There isn't even a need for a direct tensor function now, as it reads pretty directly from the above

If we want to get even crazier, we could overload basic arithmetic, and support something like return qml.expval(a*op1(wire1) + b*op2(wire2)) This would be pretty useful for algorithms like VQE.

For devices that don't support arbitrary hermitian observables, we would need some extra logic though.

[quantshah]

Hi, I would like to continue the discussion on the PR. At the moment, let us restrict to only allowing tensors of single-qubit observables. Since we have the new pl.expval() API I think it is natural to have the syntax that @co9olguy suggests, pl.expval(PauliX(0), PauliZ(1)). A hacky way to do this quick now would be to convert a list of observables into a multi-qubit Hermitian observable and then replacing it in the queue. I am sort of playing around with this logic in the PR:

tensored_matrix = get_matrix_rep(list_of_ops)
QNode._current_context._append_op(pl.Hermitian(tensored_matrix, wires=wires))

But would be better to just implement a new queue? tensorop.

[co9olguy]

Hi @quantshah, what you suggest sounds reasonable, I wouldn't even consider it hacky. Basically the expval function should count the number of arguments it receives, and, if more than one, combine them into a new Hermitian observable that is placed in the queue. For two qubits, we might even want to pre-generated all pairs of Pauli matrices