# Grad only applies to real scalar-output functions 2

Hi,

Hope everyone’s had a great holiday. I’m using the following code:

``````import pennylane as qml
from pennylane import numpy as np
import time

def local_hadamard_test(weights, problem, l=None, lp=None, j=None, part=None):

ancilla_idx = problem.get_n_qubits()

# First Hadamard gate applied to the ancillary qubit.

# For estimating the imaginary part of the coefficient "mu", we must add a "-i"
# phase gate.
if part == "Im" or part == "im":
qml.PhaseShift(-np.pi / 2, wires=ancilla_idx)

# Variational circuit generating a guess for the solution vector |x>
problem.variational_block(weights)

# Controlled application of the unitary component A_l of the problem matrix A.
problem.CA(ancilla_idx, l)

# Adjoint of the unitary U_b associated to the problem vector |b>.
# In this specific example Adjoint(U_b) = U_b.
problem.U_b()

# Controlled Z operator at position j. If j = -1, apply the identity.
if j != -1:
qml.CZ(wires=[ancilla_idx, j])

# Unitary U_b associated to the problem vector |b>.
problem.U_b()

# In this specific example Adjoint(A_lp) = A_lp.
problem.CA(ancilla_idx, lp)

# Second Hadamard gate applied to the ancillary qubit.

# Expectation value of Z for the ancillary qubit.
return qml.expval(qml.PauliZ(wires=ancilla_idx))

def hadamard_overlap_test(weights, problem, l=None, lp=None, part=None):
n_qubits = problem.get_n_qubits()
ancilla_idx = n_qubits * 2

# H on ancilla index

# Variational circuit generating a guess for the solution vector |x> applied to the top half
problem.variational_block(weights, offset=n_qubits)

# unitary U_b associated to the problem vector |b> applied to the bottom half
# In this specific example Adjoint(U_b) = U_b.
problem.U_b()

# Controlled application of the unitary component A_l of the problem matrix A on the top half.
problem.CA(ancilla_idx, l, offset=n_qubits)

# Controlled application of Adjoint(A_lp) applied to the bottom half
# In this specific example Adjoint(A_lp) = A_lp. #TODO: is it really?
problem.CA(ancilla_idx, lp)

if part == "Im":
qml.RZ(phi=-np.pi/2, wires=ancilla_idx)

# bell basis observable
[qml.CNOT(wires=(i+n_qubits, i)) for i in range(n_qubits)]
[qml.Hadamard(wires=i) for i in range(n_qubits, n_qubits*2 + 1)]

# to get P(0) - P(1) we need to perform linear classical post-processing which involves using the probabilities
return qml.probs(wires=range(n_qubits*2 + 1))

# Computes the mu coefficients
def mu(weights, local_hadamard_test, problem, l=None, lp=None, j=None):
"""Generates the coefficients to compute the "local" cost function C_L."""

mu_real = local_hadamard_test(weights, problem, l=l, lp=lp, j=j, part="Re")
mu_imag = local_hadamard_test(weights, problem, l=l, lp=lp, j=j, part="Im")

return mu_real + 1.0j * mu_imag

"""Returns the normalization constant <psi|psi>, where |psi> = A |x>."""
norm = 0.0

for l in range(0, len(c)):
for lp in range(0, len(c)):
norm = norm + c[l] * np.conj(c[lp]) * mu(weights, local_hadamard_test, problem, l, lp, -1)

return abs(norm)

def get_bin(state: int, n_qubits):
"""
@param
state: a measurement outcome as an int
return: (-1 or 1, corresponding to whether the prob on the bitstring should be added or subtracted)
"""
acc = 1

# if aux qubit is 1
if state & 2**(n_qubits*2):
acc *= -1

for i in range(n_qubits):
if state & 2**i and state & 2**(i + n_qubits):
acc *= -1

return acc

def gamma(weights, hadamard_overlap_test, problem, l=None, lp=None):
n_qubits = problem.get_n_qubits()

probs_real = hadamard_overlap_test(weights, problem, l=l, lp=lp, part="Re")
probs_imag = hadamard_overlap_test(weights, problem, l=l, lp=lp, part="Im")

gamma_real = 0
gamma_imag = 0

for state, prob in enumerate(probs_real):
gamma_real += get_bin(state, n_qubits) * prob

for state, prob in enumerate(probs_imag):
gamma_imag += get_bin(state, n_qubits) * prob

# print(f"gamma: {time.time() - start:.2f}")

return 2 * (gamma_real + 1.0j * gamma_imag) # see appendix C for the 2x coeff

"""Global version of the cost function. Tends to zero when A|x> is proportional to |b>."""
c, _ = problem.get_coeffs()

norm = 0.0
overlap = 0.0

for l in range(0, len(c)):
for lp in range(0, len(c)):
# start = time.time()
norm = norm + c[l] * np.conj(c[lp]) * mu(weights, local_hadamard_test, problem, l, lp, -1)
# print(f"norm accum ({l*len(c) + lp})")

overlap = overlap + c[l] * np.conj(c[lp]) * gamma(weights, hadamard_overlap_test, problem, l, lp)

norm = abs(norm)

print(norm, overlap)

return 1 - overlap / norm # TODO: double check this expression

# convert matrix A encoded as a string (eg. "IZZ") into qml code
def A_to_code (idx, ancilla_idx, terms, offset=0):

if idx < 0:
raise ValueError("Index of linear combination must be >= 0.")

target_pauli = list(terms[idx])

order_idx = offset

for i in range(len(target_pauli)):
if target_pauli[i] == 'I':
order_idx += 1
None
if target_pauli[i] == 'X':
qml.CNOT(wires = (ancilla_idx, order_idx))
order_idx += 1
if target_pauli[i] == 'Y':
qml.CY(wires = (ancilla_idx, order_idx))
order_idx += 1
if target_pauli[i] == 'Z':
qml.CZ(wires = (ancilla_idx, order_idx))
order_idx += 1

import functools as ft

pauli_dict = {"I": qml.Identity.compute_matrix(), "X": qml.PauliX.compute_matrix(), "Y": qml.PauliY.compute_matrix(), "Z": qml.PauliZ.compute_matrix()}

def A_to_num (n_qubits: int, coefs: np.tensor, terms):
"""
Given an array of coeffs c and an array of A_l formatted as a list of strings, return A
@params
coefs (eg. [1, 0.2, 0.2])
terms (eg. ["III", "XZI", "XII"])

returns an np.array
"""
if len(coefs) != len(terms):
raise ValueError("Number of coefficients does not match number of terms.")

if n_qubits <= 0:
raise ValueError("Number of qubits is not a number greater than 0.")

terms_len = len(terms)
for i in range(terms_len):
if len(terms[i]) != n_qubits:
raise ValueError("Number of terms in each Pauli gate combination must be the same as number of qubits.")

dim = 2**n_qubits
mat = np.zeros((dim, dim), dtype=np.complex64)

for (c, pauli) in zip(coefs, terms):
pauli = [pauli_dict[key] for key in pauli]
if pauli == ["I"]:
mat += c * ft.reduce(np.kron, pauli)
else:
mat += c * ft.reduce(np.kron, pauli)

return mat

# these classes encode the linear system we are trying to solve
from abc import ABC, abstractmethod
class Problem(ABC):
def __init__(self, n_qubits, c, A_terms) -> None:
super().__init__()
self.n_qubits = n_qubits
self.A_num = A_to_num(n_qubits, c, A_terms)
self.A_terms = A_terms

# normalize c
self.c = np.array(c) / np.linalg.norm(self.A_num, ord=2)

# Total number of qubits; here we add an ancillary qubit
self.tot_qubits = self.n_qubits + 1
# Index of ancillary qubit (Python lists are 0-indexed)
self.ancilla_idx = self.n_qubits

@abstractmethod
def get_coeffs():
"""gets c, A_l"""
pass

@abstractmethod
def get_n_qubits():
"""gets number of qubits of your problem"""
pass

@abstractmethod
def U_b():
"""Unitary matrix rotating the ground state to the problem vector |b> = U_b |0>."""
pass

@abstractmethod
def CA(idx):
pass

@abstractmethod
def variational_block(weights):
pass

class ToyProblem(Problem):
def __init__(self, n_qubits):
c = [1, 0.25]
A_terms = ["III", "IIZ"]

super().__init__(n_qubits, c, A_terms)

self.param_shape = n_qubits

def get_coeffs(self):
return self.c, self.A_terms

def get_n_qubits(self):
return self.n_qubits

def U_b(self):
"""Unitary matrix rotating the ground state to the problem vector |b> = U_b |0>."""

def CA(self, ancilla_idx, idx, offset=0):
A_to_code(idx, ancilla_idx=ancilla_idx, terms=self.A_terms, offset=offset)

def variational_block(self, weights, offset=0):
[qml.RY(phi=weights[i], wires=i+offset) for i in range(self.n_qubits)]

#############
n_qubits = 3
dev_mu = qml.device("default.qubit", wires=n_qubits+1)
dev_gamma = qml.device("default.qubit", wires=n_qubits*2 + 1)

problem = ToyProblem(3)

start = time.time()

for it in range(10):

print("Step {:3d}       Cost_L = {:9.7f}".format(it, cost), flush=True)

it += 1

print(f"Training time: {time.time() - start}s")
``````

But I’m running into the following error:

``````{
"name": "TypeError",
"stack": "---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)

39 def step_and_cost(self, objective_fn, *args, grad_fn=None, **kwargs):
40     \"\"\"Update trainable arguments with one step of the optimizer and return the corresponding
41     objective function value prior to the step.
42
(...)
56         If single arg is provided, list [array] is replaced by array.
57     \"\"\"
62     if forward is None:

99 r\"\"\"Compute gradient of the objective function at the given point and return it along with
100 the objective function forward pass (if available).
101
(...)
114     will not be evaluted and instead ``None`` will be returned.
115 \"\"\"
--> 117 grad = g(*args, **kwargs)
118 forward = getattr(g, \"forward\", None)
120 num_trainable_args = sum(getattr(arg, \"requires_grad\", False) for arg in args)

115     self._forward = self._fun(*args, **kwargs)
116     return ()
119 self._forward = ans

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

136 vjp, ans = _make_vjp(fun, x)
138 if not vspace(ans).size == 1:
--> 139     raise TypeError(
140         \"Grad only applies to real scalar-output functions. \"
142     )

}
``````

I recognize that my code does output a complex value, but it has worked for a piece of similar code that also theoretically outputs complex values. I’m hesitant to follow the suggestion in the error message because it would complicate and slow down my code by a lot.

`qml.about()`:

``````Name: PennyLane
Version: 0.33.1
Summary: PennyLane is a Python quantum machine learning library by Xanadu Inc.
Home-page: https://github.com/PennyLaneAI/pennylane
Author:
Author-email:
Requires: appdirs, autograd, autoray, cachetools, networkx, numpy, pennylane-lightning, requests, rustworkx, scipy, semantic-version, toml, typing-extensions
Required-by: PennyLane-Lightning, PennyLane-qiskit

Platform info:           macOS-12.6-x86_64-i386-64bit
Python version:          3.11.6
Numpy version:           1.23.5
Scipy version:           1.10.1
Installed devices:
- default.gaussian (PennyLane-0.33.1)
- default.mixed (PennyLane-0.33.1)
- default.qubit (PennyLane-0.33.1)
- default.qubit.jax (PennyLane-0.33.1)
- default.qubit.legacy (PennyLane-0.33.1)
- default.qubit.tf (PennyLane-0.33.1)
- default.qubit.torch (PennyLane-0.33.1)
- default.qutrit (PennyLane-0.33.1)
- null.qubit (PennyLane-0.33.1)
- lightning.qubit (PennyLane-Lightning-0.33.1)
- qiskit.aer (PennyLane-qiskit-0.33.1)
- qiskit.basicaer (PennyLane-qiskit-0.33.1)
- qiskit.ibmq (PennyLane-qiskit-0.33.1)
- qiskit.ibmq.circuit_runner (PennyLane-qiskit-0.33.1)
- qiskit.ibmq.sampler (PennyLane-qiskit-0.33.1)
- qiskit.remote (PennyLane-qiskit-0.33.1)

``````
1 Like

Hey @jkwan314, hope you had a nice holiday as well!

If `cost_global` is giving you anything but a real-valued scalar output, then there will be problems. Is it that the output has a non-zero imaginary component, or is it that it’s a >0D array?

Sorry, I realized I forgot to norm somewhere which is causing an imaginary output. Thanks!

1 Like

Ah! Yep that can happen. Glad you figured it out!