# Source code for pennylane.optimize.spsa

# Copyright 2018-2022 Xanadu Quantum Technologies Inc.

# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at

# Unless required by applicable law or agreed to in writing, software
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
"""SPSA optimizer"""

import numpy as np
from pennylane.measurements import Shots

[docs]class SPSAOptimizer:
r"""The Simultaneous Perturbation Stochastic Approximation method (SPSA)
is a stochastic approximation algorithm for optimizing cost functions whose evaluation may involve noise.

While other gradient-based optimization methods usually attempt to compute
evaluating the cost function twice in each iteration step. This cost may result in
a significant decrease in the overall cost of function evaluations for the entire optimization.
It is based on an approximation of the unknown gradient :math:\hat{g}(\hat{\theta}_{k})
through a simultaneous perturbation of the input parameters:

.. math::
\hat{g}_k(\hat{\theta}_k) = \frac{y(\hat{\theta}_k+c_k\Delta_k)-
y(\hat{\theta}_k-c_k\Delta_k)}{2c_k} \begin{bmatrix}
\Delta_{k1}^{-1} \\
\Delta_{k2}^{-1} \\
\vdots \\
\Delta_{kp}^{-1}
\end{bmatrix}\text{,}

where

* :math:k is the current iteration step,
* :math:\hat{\theta}_k are the input parameters at iteration step :math:k,
* :math:y is the objective function,
* :math:c_k=\frac{c}{k^\gamma} is the gain sequence corresponding to evaluation step size
and it can be controlled with

* scaling parameter :math:c and
* scaling exponent :math:\gamma

* :math:\Delta_{ki}^{-1} \left(1 \leq i \leq p \right) are the inverted elements of
random pertubation vector :math:\Delta_k.

:math:\hat{\theta}_k is updated to a new set of parameters with

.. math::
\hat{\theta}_{k+1} = \hat{\theta}_{k} - a_k\hat{g}_k(\hat{\theta}_k)\text{,}

where the gain sequences :math:a_k=\frac{a}{(A+k)^\alpha} controls parameter update step size.

The gain sequence :math:a_k can be controlled with

* scaling parameter :math:a,
* scaling exponent :math:\alpha and
* stability constant :math:A

For more details, see Spall (1998a)
<https://www.jhuapl.edu/SPSA/PDF-SPSA/Spall_An_Overview.PDF>_.

.. note::

* One SPSA iteration step of a cost function that involves computing the expectation value of
a Hamiltonian with M terms requires :math:2*M quantum device executions.
* The forward-pass value of the cost function is not computed when stepping the optimizer.
Therefore, in case of using step_and_cost method instead of step, the number
of executions will include the cost function evaluations.

**Examples:**

For VQE/VQE-like problems, the objective function can be the following:

>>> coeffs = [0.2, -0.543, 0.4514]
>>> obs = [qml.X(0) @ qml.Z(1), qml.Z(0) @ qml.Hadamard(2),
...             qml.X(3) @ qml.Z(1)]
>>> H = qml.Hamiltonian(coeffs, obs)
>>> num_qubits = 4
>>> dev = qml.device("default.qubit", wires=num_qubits)
>>> @qml.qnode(dev)
... def cost(params, num_qubits=1):
...     qml.BasisState(np.array([1, 1, 0, 0]), wires=range(num_qubits))
...     for i in range(num_qubits):
...         qml.Rot(*params[i], wires=0)
...         qml.CNOT(wires=[2, 3])
...         qml.CNOT(wires=[2, 0])
...         qml.CNOT(wires=[3, 1])
...     return qml.expval(H)
...
>>> params = np.random.normal(0, np.pi, (num_qubits, 3), requires_grad=True)

Once constructed, the cost function can be passed directly to the
step or step_and_cost function of the optimizer:

>>> max_iterations = 100
>>> opt = qml.SPSAOptimizer(maxiter=max_iterations)
>>> for _ in range(max_iterations):
...     params, energy = opt.step_and_cost(cost, params, num_qubits=num_qubits)
>>> print(energy)
-0.4294539602541956

The algorithm provided by SPSA does not rely on built-in automatic differentiation capabilities of the interface being used
and therefore the optimizer can be used in more complex hybrid classical-quantum workflow with any of the interfaces:

>>> n_qubits = 1
>>> max_iterations = 20
>>> dev = qml.device("default.qubit", wires=n_qubits)
>>> @qml.qnode(dev, interface="tf")
... def layer_fn_spsa(inputs, weights):
...     qml.AngleEmbedding(inputs, wires=range(n_qubits))
...     qml.BasicEntanglerLayers(weights, wires=range(n_qubits))
...     return qml.expval(qml.Z(0))
...
>>> opt = qml.SPSAOptimizer(maxiter=max_iterations)
... def fn(params, tensor_in, tensor_out):
...     with tf.init_scope():
...             for _ in range(max_iterations):
...                     # Some classical steps before the quantum computation
...                     params_a, layer_res = opt.step_and_cost(layer_fn_spsa,
...                                     np.tensor(params))
...                     params = params_a[1]
...                     tensor_out = layer_res
...                     # Some classical steps after the quantum computation
...     return layer_res
...
>>> tensor_in = tf.Variable([0.27507603])
>>> tensor_out = tf.Variable([0])
>>> params = tf.Variable([[3.97507603],
...     [3.12950603],
...     [1.00854038],
...     [1.25907603]])
>>> loss = fn(params, tensor_in, tensor_out)
>>> print(loss)
tf.Tensor(-0.9995854230771829, shape=(), dtype=float64)

Keyword Args:
maxiter (int): the maximum number of iterations expected to be performed.
Used to determine :math:A, if :math:A is not supplied, otherwise ignored.
alpha (float): A hyperparameter to calculate :math:a_k=\frac{a}{(A+k+1)^\alpha}
for each iteration. Its asymptotically optimal value is 1.0.
gamma (float): An hyperparameter to calculate :math:c_k=\frac{c}{(k+1)^\gamma}
for each iteration. Its asymptotically optimal value is 1/6.
c (float): A hyperparameter related to the expected noise. It should be
approximately the standard deviation of the expected noise of the cost function.
A (float): The stability constant; if not provided, set to be 10% of the maximum number
of expected iterations.
a (float): A hyperparameter expected to be small in noisy situations,
its value could be picked using A, :math:\alpha and :math:\hat{g_0} (\hat{\theta_0}).
For more details, see Spall (1998b)
<https://www.jhuapl.edu/spsa/PDF-SPSA/Spall_Implementation_of_the_Simultaneous.PDF>_.
"""

# pylint: disable-msg=too-many-arguments
def __init__(self, maxiter=None, alpha=0.602, gamma=0.101, c=0.2, A=None, a=None):
self.a = a
self.A = A
if not maxiter and not A:
raise TypeError("One of the parameters maxiter or A must be provided.")
if not A:
self.A = maxiter * 0.1
if not a:
self.a = 0.05 * (self.A + 1) ** alpha
self.c = c
self.alpha = alpha
self.gamma = gamma
self.k = 1
self.ak = self.a / (self.A + 1) ** self.alpha

[docs]    def step_and_cost(self, objective_fn, *args, **kwargs):
r"""Update the parameter array :math:\hat{\theta}_k with one step of the
optimizer and return the step and the corresponding objective function. The number
of steps stored by the k attribute of the optimizer is counted internally when calling step_and_cost and cost.

Args:
objective_fn (function): the objective function for optimization
*args : variable length argument array for objective function
**kwargs : variable length of keyword arguments for the objective function

Returns:
tuple[list [array], float]: the new variable values :math:\hat{\theta}_{k+1} and the
objective function output prior to the step.
"""

self.k += 1

forward = objective_fn(*args, **kwargs)

# unwrap from list if one argument, cleaner return
if len(new_args) == 1:
return new_args[0], forward
return new_args, forward

[docs]    def step(self, objective_fn, *args, **kwargs):
r"""Update trainable arguments with one step of the optimizer. The number
of steps is being counted through calls to step_and_cost and cost.

Args:
objective_fn (function): the objective function for optimization
*args : variable length argument array for objective function
**kwargs : variable length of keyword arguments for the objective function

Returns:
list [array]: the new variable values :math:\hat{\theta}_{k+1}.
"""

self.k += 1

# unwrap from list if one argument, cleaner return
if len(new_args) == 1:
return new_args[0]

return new_args

[docs]    def compute_grad(self, objective_fn, args, kwargs):
r"""Approximate the gradient of the objective function at the
given point.

Args:
objective_fn (function): The objective function for optimization
args (tuple): tuple of NumPy array containing the current parameters
for objective function
kwargs (dict): keyword arguments for the objective function

Returns:
tuple (array): NumPy array containing the gradient
:math:\hat{g}_k(\hat{\theta}_k)
"""
ck = self.c / self.k**self.gamma

delta = []
thetaplus = list(args)
thetaminus = list(args)

for index, arg in enumerate(args):
# Use the symmetric Bernoulli distribution to generate
# the coordinates of delta. Note that other distributions
# may also be used (they need to satisfy certain conditions).
di = np.random.choice([-1, 1], size=arg.shape)
multiplier = ck * di
thetaplus[index] = arg + multiplier
thetaminus[index] = arg - multiplier
delta.append(di)
yplus = objective_fn(*thetaplus, **kwargs)
yminus = objective_fn(*thetaminus, **kwargs)
try:
# pylint: disable=protected-access
dev_shots = objective_fn.device.shots
if isinstance(dev_shots, Shots):
shots = dev_shots if dev_shots.has_partitioned_shots else Shots(None)
elif objective_fn.device.shot_vector is not None:
shots = Shots(objective_fn.device._raw_shot_sequence)  # pragma: no cover
else:
shots = Shots(None)
if np.prod(objective_fn.func(*args).shape(objective_fn.device, shots)) > 1:
raise ValueError(
"The objective function must be a scalar function for the gradient "
"to be computed."
)
except AttributeError:
if yplus.size > 1:
raise ValueError(  # pylint: disable=raise-missing-from
"The objective function must be a scalar function for the gradient "
"to be computed."
)
grad = [(yplus - yminus) / (2 * ck * di) for di in delta]

r"""Update the variables to take a single optimization step.

Args:
function at point :math:\hat{\theta}_{k}
args (tuple): the current value of the variables :math:\hat{\theta}_{k}

Returns:
list [array]: the new values :math:\hat{\theta}_{k+1}"""
self.ak = self.a / (self.A + self.k) ** self.alpha
args_new = list(args)
trained_index = 0
for index, arg in enumerate(args):
args_new[index] = arg - self.ak * grad[trained_index]
trained_index += 1

return args_new


Using PennyLane

Development

API

Internals