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.
qml.Hadamard(wires=ancilla_idx)
# 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()
# Controlled application of Adjoint(A_lp).
# In this specific example Adjoint(A_lp) = A_lp.
problem.CA(ancilla_idx, lp)
# Second Hadamard gate applied to the ancillary qubit.
qml.Hadamard(wires=ancilla_idx)
# 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
qml.Hadamard(ancilla_idx)
# 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
def psi_norm(weights, c, local_hadamard_test, problem):
"""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
def cost_global(problem, weights, local_hadamard_test, hadamard_overlap_test):
"""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>."""
[qml.Hadamard(wires=i) for i in [0,1]]
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)
local_hadamard_test = qml.QNode(local_hadamard_test, dev_mu, interface="autograd")
hadamard_overlap_test = qml.QNode(hadamard_overlap_test, dev_gamma, interface="autograd")
problem = ToyProblem(3)
opt = qml.GradientDescentOptimizer(0.1)
w = np.random.randn(problem.param_shape, requires_grad=True)
start = time.time()
for it in range(10):
w, cost = opt.step_and_cost(lambda w: cost_global(problem, w, local_hadamard_test, hadamard_overlap_test), w)
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",
"message": "Grad only applies to real scalar-output functions. Try jacobian, elementwise_grad or holomorphic_grad.",
"stack": "---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
/Users/bigsad/Downloads/Algorithm-Research/Student-Hub/Indy-Ng/cost_global_error.ipynb セル 1 line 2
<a href='vscode-notebook-cell:/Users/bigsad/Downloads/Algorithm-Research/Student-Hub/Indy-Ng/cost_global_error.ipynb#W0sZmlsZQ%3D%3D?line=292'>293</a> start = time.time()
<a href='vscode-notebook-cell:/Users/bigsad/Downloads/Algorithm-Research/Student-Hub/Indy-Ng/cost_global_error.ipynb#W0sZmlsZQ%3D%3D?line=294'>295</a> for it in range(10):
--> <a href='vscode-notebook-cell:/Users/bigsad/Downloads/Algorithm-Research/Student-Hub/Indy-Ng/cost_global_error.ipynb#W0sZmlsZQ%3D%3D?line=295'>296</a> w, cost = opt.step_and_cost(lambda w: cost_global(problem, w, local_hadamard_test, hadamard_overlap_test), w)
<a href='vscode-notebook-cell:/Users/bigsad/Downloads/Algorithm-Research/Student-Hub/Indy-Ng/cost_global_error.ipynb#W0sZmlsZQ%3D%3D?line=297'>298</a> print(\"Step {:3d} Cost_L = {:9.7f}\".format(it, cost), flush=True)
<a href='vscode-notebook-cell:/Users/bigsad/Downloads/Algorithm-Research/Student-Hub/Indy-Ng/cost_global_error.ipynb#W0sZmlsZQ%3D%3D?line=299'>300</a> it += 1
File ~/Downloads/Algorithm-Research/Student-Hub/Indy-Ng/.venv/lib/python3.11/site-packages/pennylane/optimize/gradient_descent.py:59, in GradientDescentOptimizer.step_and_cost(self, objective_fn, grad_fn, *args, **kwargs)
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 \"\"\"
---> 59 g, forward = self.compute_grad(objective_fn, args, kwargs, grad_fn=grad_fn)
60 new_args = self.apply_grad(g, args)
62 if forward is None:
File ~/Downloads/Algorithm-Research/Student-Hub/Indy-Ng/.venv/lib/python3.11/site-packages/pennylane/optimize/gradient_descent.py:117, in GradientDescentOptimizer.compute_grad(objective_fn, args, kwargs, grad_fn)
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 \"\"\"
116 g = get_gradient(objective_fn) if grad_fn is None else grad_fn
--> 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)
File ~/Downloads/Algorithm-Research/Student-Hub/Indy-Ng/.venv/lib/python3.11/site-packages/pennylane/_grad.py:118, in grad.__call__(self, *args, **kwargs)
115 self._forward = self._fun(*args, **kwargs)
116 return ()
--> 118 grad_value, ans = grad_fn(*args, **kwargs) # pylint: disable=not-callable
119 self._forward = ans
121 return grad_value
File ~/Downloads/Algorithm-Research/Student-Hub/Indy-Ng/.venv/lib/python3.11/site-packages/autograd/wrap_util.py:20, in unary_to_nary.<locals>.nary_operator.<locals>.nary_f(*args, **kwargs)
18 else:
19 x = tuple(args[i] for i in argnum)
---> 20 return unary_operator(unary_f, x, *nary_op_args, **nary_op_kwargs)
File ~/Downloads/Algorithm-Research/Student-Hub/Indy-Ng/.venv/lib/python3.11/site-packages/pennylane/_grad.py:139, in grad._grad_with_forward(fun, x)
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. \"
141 \"Try jacobian, elementwise_grad or holomorphic_grad.\"
142 )
144 grad_value = vjp(vspace(ans).ones())
145 return grad_value, ans
TypeError: Grad only applies to real scalar-output functions. Try jacobian, elementwise_grad or holomorphic_grad."
}
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:
License: Apache License 2.0
Location: /Users/bigsad/Downloads/Algorithm-Research/Student-Hub/Indy-Ng/.venv/lib/python3.11/site-packages
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.autograd (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)
Yep that can happen. Glad you figured it out!