# Cost functions: Multiple wire measurements + backends for Hermitians

#1

Hey!

I would like to construct a cost function to minimise that looks like this

C = < z1 > + < z2 > + < z1z2 >

Where is the expectation over the Hermitian observable
zz_observable = (1/2)*np.array([[1,0,0,0],[0,-1,0,0],[0,0,-1,0],[0,0,0,1]])

1. The Hermitian observable is not supported by the qiskit.basicaer backend. Do you know what might be a good replacement? I was hoping to use the noise functionality of the qiskit simulator.

2. If inside the qnode you perform the < z1 > and < z1z2 > computations the error returned is
QuantumFunctionError: Each wire in the quantum circuit can only be measured once.
Does this mean that we can’t compute both < z1 > and < z1z2 >?

Full code + errors reproduced below.
CODE

# Conversion between problem idxs (with signs) to qubits idxs (without signs)

sign = lambda x: np.sign(x)

def sign_qubit(term): # s denotes the sign and a and b refer to the idx. Minus 1 from idx to return to standard indexing
a = abs(term[0]) - 1
b = abs(term[1]) - 1
sa = sign(term[0])
sb = sign(term[1])
return sa, a, sb, b

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

zz_observable = (1/2)*np.array([[1,0,0,0],[0,-1,0,0],[0,0,-1,0],[0,0,0,1]])

problem  = [[2, 1]]
qubits  = ['zz1', 'zz2']


# QAOA phases

def mixer(beta):
for q in range(N): qml.RX(beta, wires=q)
return None

def separator(gamma, problem):
for term in problem:
sa, a, sb, b = sign_qubit(term)
qml.RZ(sa * gamma, wires=a)
qml.RZ(sb * gamma, wires=b)
## ZZ operation
qml.CNOT(wires=[a, b])
qml.RZ(sb * gamma * 2., wires=b)
qml.CNOT(wires=[a, b])
return None

@qml.qnode(dev)
def qaoa(params):

gamma = params[0]
beta = params[1]
separator(gamma, problem)
mixer(beta)

# Cost
single_qubit_expectations = [qml.expval.PauliZ(i) for i in range(N)]
two_qubit_expectations = [qml.expval.Hermitian(zz_observable,
wires=[abs(a)-1,abs(b)-1]) for a,b in problem]
single_qubit_expectations.extend(two_qubit_expectations)
return single_qubit_expectations


# initialise the optimizer

opt = qml.GradientDescentOptimizer(stepsize=0.4)

def cost(params):
return np.sum(qaoa(params))


# set the number of steps

steps = 100


# set the initial parameter values

params = np.array([0.01, 0.01])

for i in range(steps):
print('here')
# update the circuit parameters
params = opt.step(cost, params)
print('here')

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

print('Optimized rotation angles: {}'.format(params))


## ERROR here

QuantumFunctionError Traceback (most recent call last)
in
64 print(‘here’)
65 # update the circuit parameters
—> 66 params = opt.step(cost, params)
67 print(‘here’)
68

61 “”"
62
64

85 else:
88 return g
89

18 else:
19 x = tuple(args[i] for i in argnum)
—> 20 return unary_operator(unary_f, x, *nary_op_args, **nary_op_kwargs)
21 return nary_f
22 return nary_operator

22 arguments as fun, but returns the gradient instead. The function fun
23 should be scalar-valued. The gradient has the same type as the argument."""
—> 24 vjp, ans = _make_vjp(fun, x)
25 if not vspace(ans).size == 1:
26 raise TypeError("Grad only applies to real scalar-output functions. "

8 def make_vjp(fun, x):
9 start_node = VJPNode.new_root(x)
—> 10 end_value, end_node = trace(start_node, fun, x)
11 if end_node is None:
12 def vjp(g): return vspace(x).zeros()

8 with trace_stack.new_trace() as t:
9 start_box = new_box(x, t, start_node)
—> 10 end_box = fun(start_box)
11 if isbox(end_box) and end_box._trace == start_box._trace:
12 return end_box._value, end_box._node

13 else:
14 subargs = subvals(args, zip(argnum, x))
—> 15 return fun(*subargs, **kwargs)
16 if isinstance(argnum, int):
17 x = args[argnum]

in cost(params)
54
55 def cost(params):
—> 56 return np.sum(qaoa(params))
57
58 # set the number of steps

~/anaconda3/envs/vqa/lib/python3.6/site-packages/PennyLane-0.3.1-py3.6.egg/pennylane/decorator.py in wrapper(*args, **kwargs)
151 def wrapper(*args, **kwargs):
152 “”“Wrapper function”""
–> 153 return qnode(*args, **kwargs)
154
155 # bind the jacobian method to the wrapped function

~/anaconda3/envs/vqa/lib/python3.6/site-packages/PennyLane-0.3.1-py3.6.egg/pennylane/qnode.py in call(self, *args, **kwargs)
455 # pylint: disable=no-member
456 args = autograd.builtins.tuple(args) # prevents autograd boxed arguments from going through to evaluate
–> 457 return self.evaluate(args, **kwargs) # args as one tuple
458
459 @ae.primitive

42 parents = tuple(box._node for _ , box in boxed_args)
43 argnums = tuple(argnum for argnum, _ in boxed_args)
—> 44 ans = f_wrapped(*argvals, **kwargs)
45 node = node_constructor(ans, f_wrapped, argvals, kwargs, argnums, parents)
46 return new_box(ans, trace, node)

46 return new_box(ans, trace, node)
47 else:
—> 48 return f_raw(*args, **kwargs)
49 f_wrapped.fun = f_raw

~/anaconda3/envs/vqa/lib/python3.6/site-packages/PennyLane-0.3.1-py3.6.egg/pennylane/qnode.py in evaluate(self, args, **kwargs)
497 m_wires = list(w for ex in self.ev for w in ex.wires)
498 if len(m_wires) != len(set(m_wires)):
–> 499 raise QuantumFunctionError(‘Each wire in the quantum circuit can only be measured once.’)
500
501 def check_op(op):

QuantumFunctionError: Each wire in the quantum circuit can only be measured once.

#2

The Hermitian observable is not supported by the qiskit.basicaer backend. Do you know what might be a good replacement? I was hoping to use the noise functionality of the qiskit simulator.

We are currently in the process of refactoring the Qiskit plugin, due to the recent release of Qiskit v0.9 and 0.10 this week, so this is definitely on the list of features we will be adding!

In the meantime, you could check out the PennyLane-Forest plugin, which also allows qubit simulations with noise, and supports all core PennyLane expectation values.

Note that the Forest plugin currently must be installed by cloning from the GitHub repository directly:

git clone https://github.com/rigetti/pennylane-forest.git
cd pennylane-forest
pip install -e .


If inside the qnode you perform the \langle z_1\rangle and \langle z_1z_2\rangle computations the error returned is

QuantumFunctionError: Each wire in the quantum circuit can only be measured once.


Does this mean that we can’t compute both \langle z_1\rangle and \langle z_1z_2\rangle?

Yes - you cannot compute both \langle z_1\rangle and \langle z_1z_2\rangle from within the same QNode. The reason for this is that hardware is a first-class citizen in PennyLane — and since measurement affects the state of a qubit in hardware, it cannot be measured twice at the end of a computation.

One way around this is to define two QNodes, that both use the same ansatz, but with different measurements:

def ansatz(params):
for i in range(N):

gamma = params[0]
beta = params[1]
separator(gamma, problem)
mixer(beta)

@qml.qnode(dev)
def qaoa_one_qubit(params):
ansatz(params)
return [qml.expval.PauliZ(i) for i in range(N)]

@qml.qnode(dev)
def qaoa_two_qubit(params):
ansatz(params)
wires =[[abs(a)-1, abs(b)-1] for a, b in problem]
return [qml.expval.Hermitian(zz_observable, wires=w) for w in wires]


You can then sum the results from the two QNodes using a classical cost function:

def cost(params):
one_qubit_sum = np.sum(qaoa_one_qubit(params))
two_qubit_sum = np.sum(qaoa_two_qubit(params))
return one_qubit_sum + two_qubit_sum


Hope that helps! Let me know if you have any other questions.

#3

Hey Josh, thanks for the good info & the ansatz tip.

Does the pennylane-forest plugin support >1 qubit Hermitian observables?
If no, is there a plugin that supports >1 qubit Hermitian observables and noisy simulation?

The pennylane-forest plugin is now complaining about a ‘2x2 matrix required’. Code + Error below. I think this is the limitation on zz observables, but I have installed the master branch of pennylane which accepts >2d Hermitians. This is still related to ‘backends for Hermitians’ but I can open a new topic if you prefer.

PENNYLANE VERSION
import pennylane
pennylane.version
‘0.3.1’

CODE

dev = qml.device(‘forest.numpy_wavefunction’, wires=N)

zz_observable = (1/2)*np.array([[1,0,0,0],[0,-1,0,0],[0,0,-1,0],[0,0,0,1]])

# QAOA phases

def mixer(beta):
for q in range(N): qml.RX(beta, wires=q)
return None

def separator(gamma, problem):
for term in problem:
sa, a, sb, b = sign_qubit(term)
qml.RZ(sa * gamma, wires=a)
qml.RZ(sb * gamma, wires=b)
## ZZ operation
qml.CNOT(wires=[a, b])
qml.RZ(sb * gamma * 2., wires=b)
qml.CNOT(wires=[a, b])
return None

def ansatz(params):
for i in range(N):

gamma = params[0]
beta = params[1]
separator(gamma, problem)
mixer(beta)


@qml.qnode(dev)
def qaoa_one_qubit(params):
ansatz(params)
return [qml.expval.PauliZ(i) for i in range(N)]

@qml.qnode(dev)
def qaoa_two_qubit(params):
ansatz(params)
wires =[[abs(a)-1, abs(b)-1] for a, b in problem]
return [qml.expval.Hermitian(zz_observable, wires=w) for w in wires]

# initialise the optimizer

def cost(params):
return np.sum(qaoa_two_qubit(params)) #+ np.sum(qaoa_two_qubit(params))

steps = 10000

# set the initial parameter values

params = np.array([0.01, 0.01])

for i in range(steps):
# update the circuit parameters
params = opt.step(cost, params)

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


print(‘Optimized rotation angles: {}’.format(params))

## ERROR

ValueError Traceback (most recent call last)
in
56 for i in range(steps):
57 # update the circuit parameters
—> 58 params = opt.step(cost, params)
59
60 if (i+1) % 5 == 0:

61 “”"
62
64

85 else:
88 return g
89

18 else:
19 x = tuple(args[i] for i in argnum)
—> 20 return unary_operator(unary_f, x, *nary_op_args, **nary_op_kwargs)
21 return nary_f
22 return nary_operator

22 arguments as fun, but returns the gradient instead. The function fun
23 should be scalar-valued. The gradient has the same type as the argument."""
—> 24 vjp, ans = _make_vjp(fun, x)
25 if not vspace(ans).size == 1:
26 raise TypeError("Grad only applies to real scalar-output functions. "

8 def make_vjp(fun, x):
9 start_node = VJPNode.new_root(x)
—> 10 end_value, end_node = trace(start_node, fun, x)
11 if end_node is None:
12 def vjp(g): return vspace(x).zeros()

8 with trace_stack.new_trace() as t:
9 start_box = new_box(x, t, start_node)
—> 10 end_box = fun(start_box)
11 if isbox(end_box) and end_box._trace == start_box._trace:
12 return end_box._value, end_box._node

13 else:
14 subargs = subvals(args, zip(argnum, x))
—> 15 return fun(*subargs, **kwargs)
16 if isinstance(argnum, int):
17 x = args[argnum]

in cost(params)
47
48 def cost(params):
—> 49 return np.sum(qaoa_two_qubit(params)) #+ np.sum(qaoa_two_qubit(params))
50
51 # set the number of steps

~/anaconda3/envs/vqa/lib/python3.6/site-packages/PennyLane-0.3.1-py3.6.egg/pennylane/decorator.py in wrapper(*args, **kwargs)
151 def wrapper(*args, **kwargs):
152 “”“Wrapper function”""
–> 153 return qnode(*args, **kwargs)
154
155 # bind the jacobian method to the wrapped function

~/anaconda3/envs/vqa/lib/python3.6/site-packages/PennyLane-0.3.1-py3.6.egg/pennylane/qnode.py in call(self, *args, **kwargs)
455 # pylint: disable=no-member
456 args = autograd.builtins.tuple(args) # prevents autograd boxed arguments from going through to evaluate
–> 457 return self.evaluate(args, **kwargs) # args as one tuple
458
459 @ae.primitive

42 parents = tuple(box._node for _ , box in boxed_args)
43 argnums = tuple(argnum for argnum, _ in boxed_args)
—> 44 ans = f_wrapped(*argvals, **kwargs)
45 node = node_constructor(ans, f_wrapped, argvals, kwargs, argnums, parents)
46 return new_box(ans, trace, node)

46 return new_box(ans, trace, node)
47 else:
—> 48 return f_raw(*args, **kwargs)
49 f_wrapped.fun = f_raw

~/anaconda3/envs/vqa/lib/python3.6/site-packages/PennyLane-0.3.1-py3.6.egg/pennylane/qnode.py in evaluate(self, args, **kwargs)
510 check_op(op)
511
–> 512 ret = self.device.execute(self.queue, self.ev)
513 return self.output_type(ret)
514

~/anaconda3/envs/vqa/lib/python3.6/site-packages/PennyLane-0.3.1-py3.6.egg/pennylane/_device.py in execute(self, queue, expectation)
210
211 self.pre_expval()
–> 212 expectations = [self.expval(e.name, e.wires, e.parameters) for e in expectation]
213 self.post_expval()
214

~/anaconda3/envs/vqa/lib/python3.6/site-packages/PennyLane-0.3.1-py3.6.egg/pennylane/_device.py in (.0)
210
211 self.pre_expval()
–> 212 expectations = [self.expval(e.name, e.wires, e.parameters) for e in expectation]
213 self.post_expval()
214

~/projects/vqa_optimisation/qiskit_pennylane/packages/pennylane-forest/pennylane-forest/pennylane_forest/numpy_wavefunction.py in expval(self, expectation, wires, par)
69 if self.shots == 0:
70 # exact expectation value
—> 71 ev = self.ev(A, wires)
72 else:
73 # estimate the ev

~/projects/vqa_optimisation/qiskit_pennylane/packages/pennylane-forest/pennylane-forest/pennylane_forest/numpy_wavefunction.py in ev(self, A, wires)
91 “”"
92 # Expand the Hermitian observable over the entire subsystem
—> 93 A = self.expand_one(A, wires)
94 return np.vdot(self.state, A @ self.state).real
95

~/projects/vqa_optimisation/qiskit_pennylane/packages/pennylane-forest/pennylane-forest/pennylane_forest/numpy_wavefunction.py in expand_one(self, U, wires)
105 “”"
106 if U.shape != (2, 2):
–> 107 raise ValueError(‘2x2 matrix required.’)
108 if len(wires) != 1:
109 raise ValueError(‘One target subsystem required.’)

ValueError: 2x2 matrix required.

#4

Hi @mxn.wls, thanks for pointing this out.

Since the multi-wire qml.expval.Hermitian support was only added to the main PennyLane library last week, we still need to add support to the plugins.

I have just submitted a PR which adds this to the PennyLane-Forest plugin, see here: https://github.com/rigetti/pennylane-forest/pull/13

Until it is merged, you can clone and install this branch directly to get it working:

git clone https://github.com/rigetti/pennylane-forest.git
cd pennylane-forest
git checkout multi_qubit_observables
pip install -e .


#5

Back to this problem though with a slightly different tack:

1. How does the runtime scale with multiple calls to the circuit in order to find multiple wire measurements of the same qubit? Is is linear i.e. every circuit call is x time so if I run two circuits to measure the same qubit twice the runtime is 2x? Example below.

# Code

import pennylane as qml
from pennylane import numpy as np

N = 4
n_clauses = 3
dev = qml.device(‘forest.numpy_wavefunction’, wires=N)
zz_observable = np.array([[1,0,0,0],[0,-1,0,0],[0,0,-1,0],[0,0,0,1]])

@qml.qnode(dev)
def qaoa_two_qubit(params, problem_arg=None):
return [qml.expval.Hermitian(zz_observable, wires=w) for w in problem_arg]

problem = [[1,2],[2,3]]

def cost(params):
two_qubit_exp = qaoa_two_qubit(params, problem_arg=problem)
return np.sum(two_qubit_exp)

steps = 1000
params = np.random.uniform(0, np.pi*2., 2)

for i in range(steps):
params = opt.step(cost, params)

# Error

QuantumFunctionError: Each wire in the quantum circuit can only be measured once.
(As expected, can’t measure qubit 2 in 2 different measurements)

# Solution multiple calls to the circuit

import pennylane as qml
from pennylane import numpy as np
N = 4
n_clauses = 3
dev = qml.device(‘forest.numpy_wavefunction’, wires=N)
zz_observable = np.array([[1,0,0,0],[0,-1,0,0],[0,0,-1,0],[0,0,0,1]])

@qml.qnode(dev)
def qaoa_two_qubit(params, term=None):
qml.RX(params[0], wires=term[0])
return qml.expval.Hermitian(zz_observable, wires=term)

problem = [[1,2],[2,3]]

def cost(params):
two_qubit_exp = [ ]
for term in problem:
two_qubit_exp.append(qaoa_two_qubit(params, term=term))
return np.sum(two_qubit_exp)

steps = 1000
params = np.random.uniform(0, np.pi*2., 2)

for i in range(steps):
params = opt.step(cost, params)
print(‘cost: %.4f’ % cost(params))

# Returns

cost: -1.7844
cost: -1.7925
cost: -1.8003
cost: -1.8078… qed

1. Is there a way to get around this: If I could sample the quantum circuit irl I can make all my gradient + expectation value calculations (I think?) based on N samples. So the runtime doesn’t scale with the number of same wire measurements, if I wanted to calculate and I don’t have to run the circuit twice.

2. Given this, I’d like to understand the limitation of pennylane & where that comes from, it seems to me given a bunch of samples (from one circuit run) it might be possible to calculate all the the expectation values and gradients.

Thanks!

#6

Hi mxn.wls,

1. Yes, the runtime scales linearly. If you want to evaluate the same circuit, but with a different final measurement setting, the entire circuit is reevaluated. This is because PennyLane is hardware-based, and this is a constraint you have on real hardware. Saying this, we do realize that there are efficiency gains that could be had using simulators, e.g., caching the state right before measurement and re-using that.

2. Not sure I follow here, unfortunately, but I’ll give a shot to answering. Say we want to estimate quantities (expectation values, gradients, etc.) with 1000 samples. How many circuit evaluations are needed for different goals?

• Estimating an expectation value: 1000 circuit evaluations.
• Estimating a different expectation value on the same wire whose measurement operator is linearly independent from other ones you have already measured: also 1000 evaluations.
• Estimating a different expectation value on the same wire whose measurement operator is linearly dependent from other ones you have already measured: 0 evaluations
• Estimating gradients with respect to a single parameter: 2000 evaluations
There are tricks and shortcuts you can leverage when working with simulator, but these are not directly built into PennyLane
1. The limitations are physical: Estimating one expectation value does not tell me anything about another expectation value which is linearly independent from it; also, using the “parameter-shift” gradient rule, which is compatible/scalable with hardware, requires evaluation of two expectation values (at shifted parameters)

Even with all the above, we are thinking about other ways we can make things more efficient if the user is not using physical hardware, and may implement these in the future.

#7

Hey,

There is an issue with calling the same qnode multiple times with different wire measurements.

What we need

• Measure expectation values <z1 z2> <z3 z4> and <z2 z3>

Due to constraint of only measuring a qubit once, must use the same circuit twice.

If we try to do this, it throws up errors. I believe this is to do with the circuit initialization. So the question is: How do I measure multiple expectations without hardcoding a second circuit? (Needs to be able to be done iteratively dependent on the number of ‘collisions’ between variables in the cost function.)

Here’s some examples.
call_1() seems to be that the circuit only measures the first expectation value of the second problem

call_2() doesn’'t seem to have a second index to call (because it was initialized with 2 expectation value calls in this case) so crashed…

 import pennylane as qml
from pennylane import numpy as np

N = 4
dev = qml.device('forest.numpy_wavefunction', wires=N)
zz_observable = np.array([[1., 0., 0., 0.], [0., -1., 0., 0.], [0., 0., -1., 0.], [0., 0., 0., 1.]])

@qml.qnode(dev)
def c1(wires_arg=None):
return [qml.expval.Hermitian(zz_observable, wires=w) for w in wires_arg]

def call_1():
problem = {}
problem[0] = [[0,1]]
problem[1] = [[0,1],[2,3]]

for group, measurements in problem.items():
tmp = c1(wires_arg=measurements)
print('Problem: ', measurements, 'Expectations: ', tmp)
return

call_1()

@qml.qnode(dev) # new circuit
def c2(wires_arg=None):
return [qml.expval.Hermitian(zz_observable, wires=w) for w in wires_arg]

def call_2():
problem = {}
problem[0] = [[0, 1], [2, 3]]
problem[1] = [[0,1]]

for group, measurements in problem.items():
tmp = c2(wires_arg=measurements)
print('Problem: ', measurements, 'Expectations: ', tmp)
return

call_2()


Output

Problem: [[0, 1]] Expectations: [1.]
Problem: [[0, 1], [2, 3]] Expectations: [1.]
Problem: [[0, 1], [2, 3]] Expectations: [1. 1.]

Traceback (most recent call last):
File “/home/maxiwelian/projects/vqa_optimisation/comparison/faulty_hermitian_demo.py”, line 38, in
call_2()
File “/home/maxiwelian/projects/vqa_optimisation/comparison/faulty_hermitian_demo.py”, line 34, in call_2
tmp = c2(wires_arg=measurements)
File “/home/maxiwelian/anaconda3/envs/aqua/lib/python3.7/site-packages/PennyLane-0.3.1-py3.7.egg/pennylane/decorator.py”, line 153, in wrapper
try:
File “/home/maxiwelian/anaconda3/envs/aqua/lib/python3.7/site-packages/PennyLane-0.3.1-py3.7.egg/pennylane/qnode.py”, line 461, in call
return f_raw(*args, **kwargs)
File “/home/maxiwelian/anaconda3/envs/aqua/lib/python3.7/site-packages/PennyLane-0.3.1-py3.7.egg/pennylane/qnode.py”, line 501, in evaluate
File “/home/maxiwelian/anaconda3/envs/aqua/lib/python3.7/site-packages/PennyLane-0.3.1-py3.7.egg/pennylane/qnode.py”, line 501, in
File “/home/maxiwelian/anaconda3/envs/aqua/lib/python3.7/site-packages/PennyLane-0.3.1-py3.7.egg/pennylane/operation.py”, line 374, in wires
File “/home/maxiwelian/anaconda3/envs/aqua/lib/python3.7/site-packages/PennyLane-0.3.1-py3.7.egg/pennylane/operation.py”, line 374, in
File “/home/maxiwelian/anaconda3/envs/aqua/lib/python3.7/site-packages/PennyLane-0.3.1-py3.7.egg/pennylane/variable.py”, line 147, in val
IndexError: index 2 is out of bounds for axis 0 with size 2

#8

Hi @mxn.wls, this is something I’ve been working on fixing

You’re assessment is spot on — currently, the QNode caches the quantum circuit on the first evaluation, and uses this cached circuit from then on (to save time/memory).

I’ve created a new PR here on the PennyLane GitHub that alters this behaviour, and allows you to specify that the QNode should re-evaluate the function and re-construct the circuit on every execution.

To install this PR, follow these steps:

git clone https://github.com/XanaduAI/pennylane.git
cd pennylane && git checkout variables_during_construction
pip install -e .


Once it’s installed, all you need to do is call the function qml.eager_mode() after importing PennyLane:

import pennylane as qml
from pennylane import numpy as np

qml.eager_mode()


and then continue with your existing script as is. I just verified that it works, giving the output

Problem:  [[0, 1]] Expectations:  [1.]
Problem:  [[0, 1], [2, 3]] Expectations:  [1. 1.]
Problem:  [[0, 1], [2, 3]] Expectations:  [1. 1.]
Problem:  [[0, 1]] Expectations:  [1.]


Note: This PR is still going through review, and the name of the function is subject to change once it is merged in to the latest release!

Hope that helps, let me know if you have any questions.

#9

Thanks @nathan and @josh. Both helpful / worked.

Max