I am trying to use the code from pennylane’s QNSPSA tutorial linked here to run a variational quantum eigensolver.
However, if I use a template for my ansatz such as qml.templates.AllSinglesDoubles
I get an error: AdjointUndefinedError.
I have stepped through my code, and I have found that the issue is that BasisState does not have a defined adjoint.
Is there a good way to fix this?
For reference, this is my full code:
from types import FunctionType
from typing import Optional
import torch
from torch import Tensor
import math
import random
import pennylane as qml
from pennylane import numpy as np
from scipy.linalg import sqrtm
import warnings
from functools import partial
class QNSPSA:
"""Quantum natural SPSA optimizer. Refer to https://quantum-journal.org/papers/q-2021-10-20-567/
for a detailed description of the methodology. When disable_metric_tensor
is set to be True, the metric tensor estimation is disabled, and QNSPSA is
reduced to be a SPSA optimizer.
From: https://github.com/PennyLaneAI/qml/blob/master/demonstrations/qnspsa.py
Args:
stepsize (float): The learn rate.
regularization (float): Regularitzation term to the Fubini-Study
metric tensor for numerical stability.
finite_diff_step (float): step size to compute the finite difference
gradient and the Fubini-Study metric tensor.
resamplings (int): The number of samples to average for each parameter
update.
blocking (boolean): When set to be True, the optimizer only accepts
updates that lead to a loss value no larger than the loss value
before update, plus a tolerance. The tolerance is set with the
parameter history_length.
history_length (int): When blocking is True, the tolerance is set to be
the average of the cost values in the last history_length steps.
disable_metric_tensor (boolean): When set to be True, the optimizer is
reduced to be a (1st-order) SPSA optimizer.
seed (int): Seed for the random sampling.
"""
def __init__(
self,
stepsize=1e-3,
regularization=1e-3,
finite_diff_step=1e-2,
resamplings=1,
blocking=True,
history_length=5,
disable_metric_tensor=False,
seed=None,
):
self.stepsize = stepsize
self.reg = regularization
self.finite_diff_step = finite_diff_step
self.metric_tensor = None
self.k = 1
self.resamplings = resamplings
self.blocking = blocking
self.last_n_steps = np.zeros(history_length)
self.history_length = history_length
self.disable_metric_tensor = disable_metric_tensor
random.seed(seed)
return
def step(self, cost, params):
"""Update trainable arguments with one step of the optimizer.
.. warning::
When blocking is set to be True, use step_and_cost instead, as loss
measurements are required for the updates for the case.
Args:
cost (qml.QNode): the QNode wrapper for the objective function for
optimization
params (np.array): Parameter before update.
Returns:
np.array: The new variable values after step-wise update.
"""
if self.blocking:
warnings.warn(
"step_and_cost() instead of step() is called when "
"blocking is turned on, as the step-wise loss value "
"is required by the algorithm.",
stacklevel=2,
)
return self.step_and_cost(cost, params)[0]
if self.disable_metric_tensor:
return self.__step_core_first_order(cost, params)
return self.__step_core(cost, params)
def step_and_cost(self, cost, params):
"""Update trainable parameters with one step of the optimizer and return
the corresponding objective function value after the step.
Args:
cost (qml.QNode): the QNode wrapper for the objective function for
optimization
params (np.array): Parameter before update.
Returns:
tuple[np.array, float]: the updated parameter and the objective
function output before the step.
"""
params_next = (
self.__step_core_first_order(cost, params)
if self.disable_metric_tensor
else self.__step_core(cost, params)
)
if not self.blocking:
loss_curr = cost(params)
return params_next, loss_curr
params_next, loss_curr = self.__apply_blocking(cost, params, params_next)
return params_next, loss_curr
def __step_core(self, cost, params):
# Core function that returns the next parameters before applying blocking.
grad_avg = np.zeros(params.shape)
tensor_avg = np.zeros((params.size, params.size))
for i in range(self.resamplings):
grad_tapes, grad_dir = self.__get_spsa_grad_tapes(cost, params)
metric_tapes, tensor_dirs = self.__get_tensor_tapes(cost, params)
raw_results = qml.execute(grad_tapes + metric_tapes, cost.device, None)
grad = self.__post_process_grad(raw_results[:2], grad_dir)
metric_tensor = self.__post_process_tensor(raw_results[2:], tensor_dirs)
grad_avg = grad_avg * i / (i + 1) + grad / (i + 1)
tensor_avg = tensor_avg * i / (i + 1) + metric_tensor / (i + 1)
self.__update_tensor(tensor_avg)
return self.__get_next_params(params, grad_avg)
def __step_core_first_order(self, cost, params):
# Reduced core function that returns the next parameters with SPSA rule.
# Blocking not applied.
grad_avg = np.zeros(params.shape)
for i in range(self.resamplings):
grad_tapes, grad_dir = self.__get_spsa_grad_tapes(cost, params)
raw_results = qml.execute(grad_tapes, cost.device, None)
grad = self.__post_process_grad(raw_results, grad_dir)
grad_avg = grad_avg * i / (i + 1) + grad / (i + 1)
return params - self.stepsize * grad_avg
def __post_process_grad(self, grad_raw_results, grad_dir):
# With the quantum measurement results from the 2 gradient tapes,
# compute the estimated gradient. Returned gradient is a tensor
# of the same shape with the input parameter tensor.
loss_forward, loss_backward = grad_raw_results
grad = (loss_forward - loss_backward) / (2 * self.finite_diff_step) * grad_dir
return grad
def __post_process_tensor(self, tensor_raw_results, tensor_dirs):
# With the quantum measurement results from the 4 metric tensor tapes,
# compute the estimated raw metric tensor. Returned raw metric tensor
# is a tensor of shape (d x d), d being the dimension of the input parameter
# to the ansatz.
tensor_finite_diff = (
tensor_raw_results[0][0][0]
- tensor_raw_results[1][0][0]
- tensor_raw_results[2][0][0]
+ tensor_raw_results[3][0][0]
)
metric_tensor = (
-(
np.tensordot(tensor_dirs[0], tensor_dirs[1], axes=0)
+ np.tensordot(tensor_dirs[1], tensor_dirs[0], axes=0)
)
* tensor_finite_diff
/ (8 * self.finite_diff_step * self.finite_diff_step)
)
return metric_tensor
def __get_next_params(self, params, gradient):
grad_vec, params_vec = gradient.reshape(-1), params.reshape(-1)
new_params_vec = np.linalg.solve(
self.metric_tensor,
(-self.stepsize * grad_vec + np.matmul(self.metric_tensor, params_vec)),
)
return new_params_vec.reshape(params.shape)
def __get_perturbation_direction(self, params):
param_number = len(params) if isinstance(params, list) else params.size
sample_list = random.choices([-1, 1], k=param_number)
direction = np.array(sample_list).reshape(params.shape)
return direction
def __get_spsa_grad_tapes(self, cost, params):
# Returns the 2 tapes along with the sampled direction that will be
# used to estimate the gradient per optimization step. The sampled
# direction is of the shape of the input parameter.
direction = self.__get_perturbation_direction(params)
cost.construct([params + self.finite_diff_step * direction], {})
tape_forward = cost.tape.copy(copy_operations=True)
cost.construct([params - self.finite_diff_step * direction], {})
tape_backward = cost.tape.copy(copy_operations=True)
return [tape_forward, tape_backward], direction
def __update_tensor(self, tensor_raw):
tensor_avg = self.__get_tensor_moving_avg(tensor_raw)
tensor_regularized = self.__regularize_tensor(tensor_avg)
self.metric_tensor = tensor_regularized
self.k += 1
def __get_tensor_tapes(self, cost, params):
# Returns the 4 tapes along with the 2 sampled directions that will be
# used to estimate the raw metric tensor per optimization step. The sampled
# directions are 1d vectors of the length of the input parameter dimension.
dir1 = self.__get_perturbation_direction(params)
dir2 = self.__get_perturbation_direction(params)
perturb1 = dir1 * self.finite_diff_step
perturb2 = dir2 * self.finite_diff_step
dir_vecs = dir1.reshape(-1), dir2.reshape(-1)
tapes = [
self.__get_overlap_tape(cost, params, params + perturb1 + perturb2),
self.__get_overlap_tape(cost, params, params + perturb1),
self.__get_overlap_tape(cost, params, params - perturb1 + perturb2),
self.__get_overlap_tape(cost, params, params - perturb1),
]
return tapes, dir_vecs
def __get_overlap_tape(self, cost, params1, params2):
op_forward = self.__get_operations(cost, params1)
op_inv = self.__get_operations(cost, params2)
with qml.tape.QuantumTape() as tape:
for op in op_forward:
qml.apply(op)
for op in reversed(op_inv):
print(op)
print(op.adjoint())
op.adjoint()
qml.probs(wires=cost.tape.wires.labels)
return tape
def __get_operations(self, cost, params):
# Given a QNode, returns the list of operations before the measurement.
cost.construct([params], {})
return cost.tape.operations
def __get_tensor_moving_avg(self, metric_tensor):
# For numerical stability: averaging on the Fubini-Study metric tensor.
if self.metric_tensor is None:
self.metric_tensor = np.identity(metric_tensor.shape[0])
return self.k / (self.k + 1) * self.metric_tensor + 1 / (self.k + 1) * metric_tensor
def __regularize_tensor(self, metric_tensor):
# For numerical stability: Fubini-Study metric tensor regularization.
tensor_reg = np.real(sqrtm(np.matmul(metric_tensor, metric_tensor)))
return (tensor_reg + self.reg * np.identity(metric_tensor.shape[0])) / (1 + self.reg)
def __apply_blocking(self, cost, params_curr, params_next):
# For numerical stability: apply the blocking condition on the parameter update.
cost.construct([params_curr], {})
tape_loss_curr = cost.tape.copy(copy_operations=True)
cost.construct([params_next], {})
tape_loss_next = cost.tape.copy(copy_operations=True)
loss_curr, loss_next = qml.execute([tape_loss_curr, tape_loss_next], cost.device, None)
# self.k has been updated earlier.
ind = (self.k - 2) % self.history_length
self.last_n_steps[ind] = loss_curr
tol = (
2 * self.last_n_steps.std()
if self.k > self.history_length
else 2 * self.last_n_steps[: self.k - 1].std()
)
if loss_curr + tol < loss_next:
params_next = params_curr
return params_next, loss_curr
geo_file = "h2.xyz"
symbols, coordinates = qml.qchem.read_structure(geo_file)
hamiltonian, qubits = qml.qchem.molecular_hamiltonian(symbols, coordinates)
print(coordinates)
print(hamiltonian)
print("Number of qubits = ", qubits)
dev = qml.device("default.qubit", wires=qubits, shots = 1000)
def AllSinglesDoubles(init_state, singles = [], doubles =[]):
def circuit(params, wires):
qml.templates.AllSinglesDoubles(weights = params, wires = wires,
hf_state = init_state, singles = singles, doubles = doubles)
return circuit, (len(singles) + len(doubles))
hf_state = np.array([1, 1, 0, 0], requires_grad=True)
singles, doubles = qml.qchem.excitations(2, 4)
ansatz, num_params = AllSinglesDoubles(init_state = hf_state, singles = singles, doubles = doubles)
print(num_params)
parameters = np.random.uniform(low=0, high=2 * np.pi, size=num_params, requires_grad=True)
print(parameters)
@qml.qnode(dev)
def cost(params):
ansatz(params, wires = range(qubits))
return qml.expval(hamiltonian)
opt = QNSPSA()
qngd_cost = []
parameters, energy = opt.step_and_cost(cost, parameters)
qngd_cost.append(energy)
max_iterations = 200
conv_tol = 10e-6
exact_value = -1.136189454088
with qml.Tracker(cost.device) as tracker:
for n in range(max_iterations):
parameters, energy = opt.step_and_cost(cost, parameters)
conv = np.abs(energy - qngd_cost[-1])
qngd_cost.append(energy)
if n % 4 == 0:
print(
"Iteration = {:}, Energy = {:.8f} Ha".format(n, energy)
)
if conv <= conv_tol:
break
print("\nFinal convergence parameter = {:.8f} Ha".format(conv))
print("Number of iterations = ", n)
print("Final value of the ground-state energy = {:.8f} Ha".format(energy))
print("Accuracy with respect to the FCI energy: {:.8f} Ha ({:.8f} kcal/mol)".format(
np.abs(energy - exact_value), np.abs(energy - exact_value) * 627.503))
print()
print("Final circuit parameters = \n", parameters)
print(tracker.totals)'
I can spot that your call for 