Optimizing state overlap / multi-qubit qml.expval.Hermitian()

Hi there,

First, I’d like to congratulate you all on making an awesome Python package! I recently saw @jmarrazola give a talk on quantum neural networks and I’ve been interested to try PennyLane for some new research I’ve recently started.

I am working on a proposed photonic architecture to prepare arbitrary multi-qubit (discrete, not CV) states, and I am trying to use PennyLane to optimize the trainable parameters to maximize the circuit fidelity <psi_out|psi_target>. Ideally, what I would like to do is something like this:

dev = qml.device('default.qubit', wires=2)

@qml.qnode(dev)
def circuit(trainable_params):
    my_configurable_circuit(trainable_params)
    target_state = 1/np.sqrt(2) * np.array([1,0,0,1])
    return qml.expval.Overlap(target_state)

circuit()
> 0.92523 (value of <psi|target>)

Since there’s no Overlap() method available, I tried improvising:

@qml.qnode(dev)
def circuit(trainable_params):
    my_configurable_circuit(trainable_params)
    target_state = 1/np.sqrt(2) * np.array([1,0,0,1])
    target_herm_op = np.outer(target_state.conj(), target_state)
    return qml.expval.Hermitian(target_herm_op, [0,1]) # gives <psi|target>^2

This gives the error ValueError: Hermitian: wrong number of wires. 2 wires given, 1 expected. (As a side note, the documentation is slightly ambiguous to a non-careful reader whether Hermitian() could be a (M<=N)-qubit operator which gets “padded” with identity operators. Perhaps change the wires keyword to wire for operators and expectations which only act on one qubit?)

Is there a way to do multi-qubit expectations like this in PennyLane? I understand this doesn’t necessarily correspond to a physically observable value that could be implemented on quantum hardware, but it’s easy to implement in TensorFlow and would be a nice feature to have in the default simulated qubit/qumode backend.

Hi @bencbartlett! Warning, this is a bit of a long post :slight_smile:

Since PennyLane treats hardware as a first class citizen when creating QNodes, we depend on the targeted device or plugin implementing the requested operations, and as such we’re a bit limited in the expectation values we can return.

We can be creative though if we want — for example, the PennyLane-Forest plugin allows for any one-qubit expectation value (even though pyQuil only supports Pauli-Z measurements!) through some nifty change of basis tricks.

In this case, there is no reason why we can’t support multi-mode Hermitian observables — to test it out, I made the following local changes to my PennyLane installation:

diff --git a/pennylane/expval/qubit.py b/pennylane/expval/qubit.py
index a47dc1e..a84029b 100644
--- a/pennylane/expval/qubit.py
+++ b/pennylane/expval/qubit.py
@@ -156,7 +156,7 @@ class Hermitian(Expectation):
         A (array): square hermitian matrix.
         wires (Sequence[int] or int): the wire the operation acts on
     """
-    num_wires = 1
+    num_wires = 0
     num_params = 1
     par_domain = 'A'
     grad_method = 'F'
diff --git a/pennylane/plugins/default_qubit.py b/pennylane/plugins/default_qubit.py
index 993e918..81803e7 100644
--- a/pennylane/plugins/default_qubit.py
+++ b/pennylane/plugins/default_qubit.py
@@ -367,10 +367,13 @@ class DefaultQubit(Device):
         Returns:
           float: expectation value :math:`\expect{A} = \bra{\psi}A\ket{\psi}`
         """
-        if A.shape != (2, 2):
-            raise ValueError('2x2 matrix required.')
+        if A.shape == (2, 2):
+            A = self.expand_one(A, wires)
+        elif A.shape == (4, 4):
+            A = self.expand_two(A, wires)
+        else:
+            raise ValueError('Only 1 or 2 wire expectation values supported.')
 
-        A = self.expand_one(A, wires)
         expectation = np.vdot(self._state, A @ self._state)
 
         if np.abs(expectation.imag) > tolerance:

That is:

  1. Modified pennylane/expval/qubit.py so that the Hermitian observable can act on any number of wires, and
  2. Modified the default.qubit simulator to implement two-mode observables.

And everything works! I tried your example (slightly modified):

import pennylane as qml
from pennylane import numpy as np

dev = qml.device('default.qubit', wires=2)

@qml.qnode(dev)
def circuit(x, target_observable=None):
    qml.RX(x[0], wires=0)
    qml.RY(x[1], wires=0)
    qml.RZ(x[2], wires=0)
    qml.CNOT(wires=[0, 1])
    return qml.expval.Hermitian(target_observable, wires=[0, 1])

target_state = 1/np.sqrt(2) * np.array([1,0,0,1])
target_herm_op = np.outer(target_state.conj(), target_state)
weights = np.array([0.5, 0.1, 0.2])

>>> circuit(weights, target_observable=target_herm_op)
0.5905564040875388

Note that QNodes are a restricted subset of Python functions, in that they must only consist of quantum operations, one per line, finishing with an expectation value. This mimics the behaviour of quantum hardware, and requires that all classical processing be done outside the QNode. In order to pass training data to the QNode, we can instead use keyword arguments, which are never treated as differentiable by PennyLane.


One final note: the default.qubit simulator is more of a reference plugin, and provided simply as a guide for plugin developers. As a result, it is not optimized for high performance, and only supports limited features (i.e., maximum of two-qubit operations/expectations, no mixed state simulations, etc.).

So, I would suggest using either PennyLane-Forest or PennyLane-ProjectQ for production code.

2 Likes

And just for completeness, here is the full optimization:

import pennylane as qml
from pennylane import numpy as np

dev = qml.device('default.qubit', wires=2)

@qml.qnode(dev)
def circuit(x, target_observable=None):
    qml.RX(x[0], wires=0)
    qml.RY(x[1], wires=0)
    qml.RZ(x[2], wires=0)
    qml.CNOT(wires=[0, 1])
    return qml.expval.Hermitian(target_observable, wires=[0, 1])

target_state = 1/np.sqrt(2) * np.array([1,0,0,1])
target_herm_op = np.outer(target_state.conj(), target_state)

weights = np.array([0.5, 0.1, 0.2])

def cost(weights):
    return np.abs(circuit(weights, target_observable=target_herm_op)-1)**2

opt = qml.AdamOptimizer(stepsize=0.4)

for i in range(2000):
    weights = opt.step(cost, weights)

    if (i+1) % 20 == 0:
        print('Cost after step {:5d}: {: .7f}'.format(i+1, cost(weights)))

fidelity = circuit(weights, target_observable=target_herm_op)

print('Optimized rotation angles: {}'.format(weights))
print('Final fidelity: {}'.format(fidelity))

# this is a hidden/undocumented attribute of the default.qubit
# plugin, that is not part of the API.
print('Output state: {}'.format(dev._state))

which produces the output

Optimized rotation angles: [ 1.44393011  1.57000406  1.44393011]
Final fidelity: 0.9999999974878844
Output state: [ 0.70714217-0.00027786j,  0.00000000+0.j  0.00000000+0.j  0.70707128-0.00027786j]

deleted by author…

Hi @josh,

Thanks for your detailed response. I’ve played around with PennyLane a lot in the last few days, and I’ve implemented the necessary subroutines for performing Shor’s algorithm for factoring small numbers. (I took your advice and am using the ProjectQ simulator as the backend – it’s incredibly performant!)

What I’m stuck on currently is a similar issue to my original post. I have N qubits in the first register, and I would like to return the expectation values of the 2^N possible logical values of the register. Ideally, I expect to get a distribution of logical expectation values which looks like this:

However, it seems there is no way to do this in PennyLane, since I am limited to returning single-qubit expectations, which trace out all other degrees of freedom and would return roughly <σZ> ~ 0 for all qubits. If I were doing this directly within ProjectQ, there is a measureGate method which seems like it should allow me to numerically compute the expectations I need, but I haven’t been able to find an equivalent measurement method or multi-qubit expectation in PennyLane. How would you suggest proceeding? Thanks in advance!

Hi @bencbartlett,

It seems like this is an identical issue to your first post, in that you would like to measure the expectation value of the observable

\hat{O}=|x\rangle \langle x|

for each computational basis state |x\rangle of the register, is this correct? The resulting expectation values would then give |\langle x|\psi\rangle|^2 for each x.


This could be done in a similar way as before:

  1. Determine the set of Hermitian matrices corresponding to |x\rangle \langle x|

  2. Modify the PennyLane expectation value operation Hermitian in pennylane/expval/qubit.py to allow it to act on multiple wires (already discussed)

  3. Modifying the ProjectQ plugin so that it knows how to apply the multi-wire expectation value.

The drawback of this approach is that the QNode will have to be redefined and re-run for every Hermitian matrix |x\rangle \langle x|.


The other approach would be a bit more intensive, but would actually involve modifying PennyLane (to remove the restriction that wires can only be measured once), as well as defining the new expectation value |x\rangle \langle x| as its own operator (perhaps qml.expval.BasisState?).

(Note: when we first conceptualized PennyLane, we made the design decision that every QNode should be restricted so that it can always run on quantum hardware, hence the restriction to only measure each wire once.)

I’m curious to hear what your thoughts are — supporting multi-wire expectation values is on our todo list for the next release of PennyLane, if you are interested in helping out with a PR!

Hi, I opened a PR to make the changes for two-qubit expectations. If this is merged then we could go for the CV case and update the changes in the plugin docs to make sure that hardware implementations can tackle this. The basis change trick seems like a neat way to do this for the two-qubit case and perhaps should be documented as a suggestion for hardware plugins. Feel free to comment or suggest changes @bencbartlett

1 Like