Update SPSA demo to count circuit executions manually.

demonstrations/tutorial_spsa.py

@@ -14,7 +14,7 @@
    tutorial_vqe_qng Accelerating VQEs with quantum natural gradient
    qnspsa Quantum natural SPSA optimizer
 
-*Authors: Antal Szava & David Wierichs — Posted: 19 March 2021. Last updated: 10 February 2023.*
+*Authors: Antal Szava & David Wierichs — Posted: 19 March 2021. Last updated: 23 February 2023.*
 
 In this tutorial, we investigate using a stochastic optimizer called
 the Simultaneous Perturbation Stochastic Approximation (SPSA) algorithm to optimize quantum
@@ -191,52 +191,38 @@ def circuit(param):
 
 ##############################################################################
 # We will execute a few optimizations in this demo, so let's prepare a convenience
-# function that runs an optimizer instance and records the cost values and the
-# number of device executions along the way. The latter will be an interesting
-# quantity to evaluate the optimization cost on hardware!
+# function that runs an optimizer instance and records the cost values
+# along the way. Together with the number of device executions, this will be an
+# interesting quantity to evaluate the optimization cost on hardware!
 
 
-def run_optimizer(optimizer, cost_function, init_param, num_steps, interval):
+def run_optimizer(opt, cost_function, init_param, num_steps, interval):
     # Copy the initial parameters to make sure they are never overwritten
     param = init_param.copy()
 
     # Obtain the device used in the cost function
     dev = cost_function.device
 
-    # Initialize the memory for cost values and device executions during the optimization
+    # Initialize the memory for cost values during the optimization
     cost_history = []
-    exec_history = [0]
     # Monitor the initial cost value
     cost_history.append(cost_function(param))
-    execs_per_cost_eval = dev.num_executions
 
-    print(
-        f"\nRunning the {optimizer.__class__.__name__} optimizer for {num_steps} iterations."
-    )
+    print(f"\nRunning the {opt.__class__.__name__} optimizer for {num_steps} iterations.")
     for step in range(num_steps):
         # Perform an update step
-        param = optimizer.step(cost_function, param)
-
-        # Monitor the device executions, deducting the executions for cost monitoring
-        exec_history.append(dev.num_executions - (step + 1) * execs_per_cost_eval)
+        param = opt.step(cost_function, param)
 
         # Monitor the cost value
         cost_history.append(cost_function(param))
 
         # Print out the status of the optimization
         if step % interval == 0:
-            print(
-                f"Iteration = {step:3d}, "
-                f"Device executions = {exec_history[step]:4d}, "
-                f"Cost = {cost_history[step]}"
-            )
+            print(f"Iteration = {step:3d}, Cost = {cost_history[step]}")
 
-    print(
-        f"Iteration = {num_steps:3d}, Device executions = {exec_history[-1]:4d}, "
-        f"Cost = {cost_history[-1]}"
-    )
+    print(f"Iteration = {num_steps:3d}, Cost = {cost_history[-1]}")
 
-    return cost_history, exec_history
+    return cost_history
 
 
 ##############################################################################
@@ -269,9 +255,10 @@ def run_optimizer(optimizer, cost_function, init_param, num_steps, interval):
 
 num_steps_spsa = 200
 opt = qml.SPSAOptimizer(maxiter=num_steps_spsa, c=0.15, a=0.2)
-cost_history_spsa, exec_history_spsa = run_optimizer(
-    opt, cost_function, init_param, num_steps_spsa, 20
-)
+cost_history_spsa = run_optimizer(opt, cost_function, init_param, num_steps_spsa, 20)
+# We have 201 energy recordings, including the initial value, and spent 2 circuit
+# evaluations in each step
+exec_history_spsa = np.arange(num_steps_spsa + 1) * 2
 
 ##############################################################################
 # Now let's perform the same optimization using gradient descent. We set the
@@ -285,16 +272,16 @@ def run_optimizer(optimizer, cost_function, init_param, num_steps, interval):
 
 num_steps_grad = 15
 opt = qml.GradientDescentOptimizer(stepsize=0.3)
-cost_history_grad, exec_history_grad = run_optimizer(
-    opt, cost_function, init_param, num_steps_grad, 3
-)
+cost_history_grad = run_optimizer(opt, cost_function, init_param, num_steps_grad, 3)
+# We have 16 energy recordings, including the initial value, and spent 2 circuit
+# evaluations *per parameter* in each step
+exec_history_grad = np.arange(num_steps_grad + 1) * (2 * np.prod(param_shape))
 
 ##############################################################################
 # SPSA and gradient descent comparison
 # ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 #
-# At this point, nothing else remains but to check which of these approaches did
-# better!
+# At this point, nothing else remains but to check which of these approaches did better!
 import matplotlib.pyplot as plt
 
 plt.figure(figsize=(10, 6))
@@ -324,7 +311,7 @@ def run_optimizer(optimizer, cost_function, init_param, num_steps, interval):
 
 ##############################################################################
 # This means that SPSA found the minimum up to an absolute accuracy of 0.01 while
-# using ten times fewer device executions than gradient descent! That's a huge
+# using multiple times fewer device executions than gradient descent! That's an important
 # saving, especially when running the algorithm on actual quantum hardware.
 #
 # SPSA and the variational quantum eigensolver
@@ -398,15 +385,17 @@ def circuit(param):
 # This random seed was used in the original VQE demo and is known to allow the
 # gradient descent algorithm to converge to the global minimum.
 np.random.seed(0)
-init_param = np.random.normal(0, np.pi, (2, num_qubits, 3), requires_grad=True)
+param_shape = (2, num_qubits, 3)
+init_param = np.random.normal(0, np.pi, param_shape, requires_grad=True)
 
 # Initialize the optimizer - optimal step size was found through a grid search
 opt = qml.GradientDescentOptimizer(stepsize=2.2)
 
 # Run the optimization
-cost_history_grad, exec_history_grad = run_optimizer(
-    opt, cost_function, init_param, num_steps_grad, 3
-)
+cost_history_grad = run_optimizer(opt, cost_function, init_param, num_steps_grad, 3)
+# We have 16 energy recordings, including the initial value, and spent 2 * 15 circuit
+# evaluations *per parameter* in each step, as there are 15 Hamiltonian terms
+exec_history_grad = np.arange(num_steps_grad + 1) * (2 * np.prod(param_shape) * 15)
 
 final_energy = cost_history_grad[-1]
 print(f"\nFinal estimated value of the ground state energy = {final_energy:.8f} Ha")
@@ -463,9 +452,10 @@ def circuit(param):
 opt = qml.SPSAOptimizer(maxiter=num_steps_spsa, c=0.3, a=1.5)
 
 # Run the SPSA algorithm
-cost_history_spsa, exec_history_spsa = run_optimizer(
-    opt, cost_function, init_param, num_steps_spsa, 20
-)
+cost_history_spsa = run_optimizer(opt, cost_function, init_param, num_steps_spsa, 20)
+# We have 161 energy recordings, including the initial value, and spent 2 * 15 circuit
+# evaluations in each step, as there are 15 terms in the Hamiltonian
+exec_history_spsa = np.arange(num_steps_spsa + 1) * (2 * 15)
 final_energy = cost_history_spsa[-1]
 
 print(f"\nFinal estimated value of the ground state energy = {final_energy:.8f} Ha")