We are currently solving the following issue within our repository. Here is the issue text:
--- BEGIN ISSUE ---
QNSPSA produces irreproducible results
<!-- ⚠️ If you do not respect this template, your issue will be closed -->
<!-- ⚠️ Make sure to browse the opened and closed issues -->

### Information

- **Qiskit Terra version**: qiskit-terra                  0.17.0 
- **Python version**:                                             Python 3.7.3
- **Operating system**: MacOS Big Sur             11.6

### What is the current behavior?                     
Executing a `QNSPSA` optimization with the same random seed gives different results.
We checked that the same optimization gives reproducible results with other optimizers such as `ADAM`.

### Steps to reproduce the problem
Run an arbitrary optimization, e.g., VQE with QNSPSA, and compare the results.


### What is the expected behavior?
The same results for the same seed.


### Suggested solutions




--- END ISSUE ---

Below are some code segments, each from a relevant file. One or more of these files may contain bugs.
--- BEGIN FILES ---
Path: `qiskit/algorithms/optimizers/qnspsa.py`
Content:
```
1 # This code is part of Qiskit.
2 #
3 # (C) Copyright IBM 2021.
4 #
5 # This code is licensed under the Apache License, Version 2.0. You may
6 # obtain a copy of this license in the LICENSE.txt file in the root directory
7 # of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
8 #
9 # Any modifications or derivative works of this code must retain this
10 # copyright notice, and modified files need to carry a notice indicating
11 # that they have been altered from the originals.
12 
13 """The QN-SPSA optimizer."""
14 
15 from typing import Any, Iterator, Optional, Union, Callable, Dict
16 
17 import numpy as np
18 from qiskit.providers import Backend
19 from qiskit.circuit import ParameterVector, QuantumCircuit
20 from qiskit.opflow import StateFn, CircuitSampler, ExpectationBase
21 from qiskit.utils import QuantumInstance
22 
23 from .spsa import SPSA, CALLBACK, TERMINATIONCHECKER, _batch_evaluate
24 
25 # the function to compute the fidelity
26 FIDELITY = Callable[[np.ndarray, np.ndarray], float]
27 
28 
29 class QNSPSA(SPSA):
30     r"""The Quantum Natural SPSA (QN-SPSA) optimizer.
31 
32     The QN-SPSA optimizer [1] is a stochastic optimizer that belongs to the family of gradient
33     descent methods. This optimizer is based on SPSA but attempts to improve the convergence by
34     sampling the **natural gradient** instead of the vanilla, first-order gradient. It achieves
35     this by approximating Hessian of the ``fidelity`` of the ansatz circuit.
36 
37     Compared to natural gradients, which require :math:`\mathcal{O}(d^2)` expectation value
38     evaluations for a circuit with :math:`d` parameters, QN-SPSA only requires
39     :math:`\mathcal{O}(1)` and can therefore significantly speed up the natural gradient calculation
40     by sacrificing some accuracy. Compared to SPSA, QN-SPSA requires 4 additional function
41     evaluations of the fidelity.
42 
43     The stochastic approximation of the natural gradient can be systematically improved by
44     increasing the number of ``resamplings``. This leads to a Monte Carlo-style convergence to
45     the exact, analytic value.
46 
47     Examples:
48 
49         This short example runs QN-SPSA for the ground state calculation of the ``Z ^ Z``
50         observable where the ansatz is a ``PauliTwoDesign`` circuit.
51 
52         .. code-block:: python
53 
54             import numpy as np
55             from qiskit.algorithms.optimizers import QNSPSA
56             from qiskit.circuit.library import PauliTwoDesign
57             from qiskit.opflow import Z, StateFn
58 
59             ansatz = PauliTwoDesign(2, reps=1, seed=2)
60             observable = Z ^ Z
61             initial_point = np.random.random(ansatz.num_parameters)
62 
63             def loss(x):
64                 bound = ansatz.bind_parameters(x)
65                 return np.real((StateFn(observable, is_measurement=True) @ StateFn(bound)).eval())
66 
67             fidelity = QNSPSA.get_fidelity(ansatz)
68             qnspsa = QNSPSA(fidelity, maxiter=300)
69             result = qnspsa.optimize(ansatz.num_parameters, loss, initial_point=initial_point)
70 
71 
72     References:
73 
74         [1] J. Gacon et al, "Simultaneous Perturbation Stochastic Approximation of the Quantum
75         Fisher Information", `arXiv:2103.09232 <https://arxiv.org/abs/2103.09232>`_
76 
77     """
78 
79     def __init__(
80         self,
81         fidelity: FIDELITY,
82         maxiter: int = 100,
83         blocking: bool = True,
84         allowed_increase: Optional[float] = None,
85         learning_rate: Optional[Union[float, Callable[[], Iterator]]] = None,
86         perturbation: Optional[Union[float, Callable[[], Iterator]]] = None,
87         last_avg: int = 1,
88         resamplings: Union[int, Dict[int, int]] = 1,
89         perturbation_dims: Optional[int] = None,
90         regularization: Optional[float] = None,
91         hessian_delay: int = 0,
92         lse_solver: Optional[Callable[[np.ndarray, np.ndarray], np.ndarray]] = None,
93         initial_hessian: Optional[np.ndarray] = None,
94         callback: Optional[CALLBACK] = None,
95         termination_checker: Optional[TERMINATIONCHECKER] = None,
96     ) -> None:
97         r"""
98         Args:
99             fidelity: A function to compute the fidelity of the ansatz state with itself for
100                 two different sets of parameters.
101             maxiter: The maximum number of iterations. Note that this is not the maximal number
102                 of function evaluations.
103             blocking: If True, only accepts updates that improve the loss (up to some allowed
104                 increase, see next argument).
105             allowed_increase: If ``blocking`` is ``True``, this argument determines by how much
106                 the loss can increase with the proposed parameters and still be accepted.
107                 If ``None``, the allowed increases is calibrated automatically to be twice the
108                 approximated standard deviation of the loss function.
109             learning_rate: The update step is the learning rate is multiplied with the gradient.
110                 If the learning rate is a float, it remains constant over the course of the
111                 optimization. It can also be a callable returning an iterator which yields the
112                 learning rates for each optimization step.
113                 If ``learning_rate`` is set ``perturbation`` must also be provided.
114             perturbation: Specifies the magnitude of the perturbation for the finite difference
115                 approximation of the gradients. Can be either a float or a generator yielding
116                 the perturbation magnitudes per step.
117                 If ``perturbation`` is set ``learning_rate`` must also be provided.
118             last_avg: Return the average of the ``last_avg`` parameters instead of just the
119                 last parameter values.
120             resamplings: The number of times the gradient (and Hessian) is sampled using a random
121                 direction to construct a gradient estimate. Per default the gradient is estimated
122                 using only one random direction. If an integer, all iterations use the same number
123                 of resamplings. If a dictionary, this is interpreted as
124                 ``{iteration: number of resamplings per iteration}``.
125             perturbation_dims: The number of perturbed dimensions. Per default, all dimensions
126                 are perturbed, but a smaller, fixed number can be perturbed. If set, the perturbed
127                 dimensions are chosen uniformly at random.
128             regularization: To ensure the preconditioner is symmetric and positive definite, the
129                 identity times a small coefficient is added to it. This generator yields that
130                 coefficient.
131             hessian_delay: Start multiplying the gradient with the inverse Hessian only after a
132                 certain number of iterations. The Hessian is still evaluated and therefore this
133                 argument can be useful to first get a stable average over the last iterations before
134                 using it as preconditioner.
135             lse_solver: The method to solve for the inverse of the Hessian. Per default an
136                 exact LSE solver is used, but can e.g. be overwritten by a minimization routine.
137             initial_hessian: The initial guess for the Hessian. By default the identity matrix
138                 is used.
139             callback: A callback function passed information in each iteration step. The
140                 information is, in this order: the parameters, the function value, the number
141                 of function evaluations, the stepsize, whether the step was accepted.
142             termination_checker: A callback function executed at the end of each iteration step. The
143                 arguments are, in this order: the parameters, the function value, the number
144                 of function evaluations, the stepsize, whether the step was accepted. If the callback
145                 returns True, the optimization is terminated.
146                 To prevent additional evaluations of the objective method, if the objective has not yet
147                 been evaluated, the objective is estimated by taking the mean of the objective
148                 evaluations used in the estimate of the gradient.
149 
150 
151         """
152         super().__init__(
153             maxiter,
154             blocking,
155             allowed_increase,
156             # trust region *must* be false for natural gradients to work
157             trust_region=False,
158             learning_rate=learning_rate,
159             perturbation=perturbation,
160             resamplings=resamplings,
161             callback=callback,
162             second_order=True,
163             hessian_delay=hessian_delay,
164             lse_solver=lse_solver,
165             regularization=regularization,
166             perturbation_dims=perturbation_dims,
167             initial_hessian=initial_hessian,
168             termination_checker=termination_checker,
169         )
170 
171         self.fidelity = fidelity
172 
173     def _point_sample(self, loss, x, eps, delta1, delta2):
174         loss_points = [x + eps * delta1, x - eps * delta1]
175         fidelity_points = [
176             (x, x + eps * delta1),
177             (x, x - eps * delta1),
178             (x, x + eps * (delta1 + delta2)),
179             (x, x + eps * (-delta1 + delta2)),
180         ]
181         self._nfev += 6
182 
183         loss_values = _batch_evaluate(loss, loss_points, self._max_evals_grouped)
184         fidelity_values = _batch_evaluate(self.fidelity, fidelity_points, self._max_evals_grouped)
185 
186         # compute the gradient approximation and additionally return the loss function evaluations
187         gradient_estimate = (loss_values[0] - loss_values[1]) / (2 * eps) * delta1
188 
189         # compute the preconditioner point estimate
190         diff = fidelity_values[2] - fidelity_values[0]
191         diff -= fidelity_values[3] - fidelity_values[1]
192         diff /= 2 * eps ** 2
193 
194         rank_one = np.outer(delta1, delta2)
195         # -0.5 factor comes from the fact that we need -0.5 * fidelity
196         hessian_estimate = -0.5 * diff * (rank_one + rank_one.T) / 2
197 
198         return np.mean(loss_values), gradient_estimate, hessian_estimate
199 
200     @property
201     def settings(self) -> Dict[str, Any]:
202         """The optimizer settings in a dictionary format."""
203         # re-use serialization from SPSA
204         settings = super().settings
205         settings.update({"fidelity": self.fidelity})
206 
207         # remove SPSA-specific arguments not in QNSPSA
208         settings.pop("trust_region")
209         settings.pop("second_order")
210 
211         return settings
212 
213     @staticmethod
214     def get_fidelity(
215         circuit: QuantumCircuit,
216         backend: Optional[Union[Backend, QuantumInstance]] = None,
217         expectation: Optional[ExpectationBase] = None,
218     ) -> Callable[[np.ndarray, np.ndarray], float]:
219         r"""Get a function to compute the fidelity of ``circuit`` with itself.
220 
221         Let ``circuit`` be a parameterized quantum circuit performing the operation
222         :math:`U(\theta)` given a set of parameters :math:`\theta`. Then this method returns
223         a function to evaluate
224 
225         .. math::
226 
227             F(\theta, \phi) = \big|\langle 0 | U^\dagger(\theta) U(\phi) |0\rangle  \big|^2.
228 
229         The output of this function can be used as input for the ``fidelity`` to the
230         :class:~`qiskit.algorithms.optimizers.QNSPSA` optimizer.
231 
232         Args:
233             circuit: The circuit preparing the parameterized ansatz.
234             backend: A backend of quantum instance to evaluate the circuits. If None, plain
235                 matrix multiplication will be used.
236             expectation: An expectation converter to specify how the expected value is computed.
237                 If a shot-based readout is used this should be set to ``PauliExpectation``.
238 
239         Returns:
240             A handle to the function :math:`F`.
241 
242         """
243         params_x = ParameterVector("x", circuit.num_parameters)
244         params_y = ParameterVector("y", circuit.num_parameters)
245 
246         expression = ~StateFn(circuit.assign_parameters(params_x)) @ StateFn(
247             circuit.assign_parameters(params_y)
248         )
249 
250         if expectation is not None:
251             expression = expectation.convert(expression)
252 
253         if backend is None:
254 
255             def fidelity(values_x, values_y):
256                 value_dict = dict(
257                     zip(params_x[:] + params_y[:], values_x.tolist() + values_y.tolist())
258                 )
259                 return np.abs(expression.bind_parameters(value_dict).eval()) ** 2
260 
261         else:
262             sampler = CircuitSampler(backend)
263 
264             def fidelity(values_x, values_y=None):
265                 if values_y is not None:  # no batches
266                     value_dict = dict(
267                         zip(params_x[:] + params_y[:], values_x.tolist() + values_y.tolist())
268                     )
269                 else:
270                     value_dict = {p: [] for p in params_x[:] + params_y[:]}
271                     for values_xy in values_x:
272                         for value_x, param_x in zip(values_xy[0, :], params_x):
273                             value_dict[param_x].append(value_x)
274 
275                         for value_y, param_y in zip(values_xy[1, :], params_y):
276                             value_dict[param_y].append(value_y)
277 
278                 return np.abs(sampler.convert(expression, params=value_dict).eval()) ** 2
279 
280         return fidelity
281
```
--- END FILES ---

Please first localize the bug based on the issue statement, and then generate a patch according to the `git diff` format fenced by three backticks.

Here is an example:

```diff
diff --git a/examples/server_async.py b/examples/server_async.py
--- a/examples/server_async.py
+++ b/examples/server_async.py
@@ -313,4 +313,4 @@
 
 
 if __name__ == "__main__":
-    asyncio.run(run_async_server("."), debug=True)
+    asyncio.run(run_async_server(), debug=True)
diff --git a/examples/server_sync.py b/examples/server_sync.py
--- a/examples/server_sync.py
+++ b/examples/server_sync.py
@@ -313,5 +313,5 @@
 
 
 if __name__ == "__main__":
-    server = run_sync_server(".")
+    server = run_sync_server()
     server.shutdown()

```

diff --git a/qiskit/algorithms/optimizers/qnspsa.py b/qiskit/algorithms/optimizers/qnspsa.py
--- a/qiskit/algorithms/optimizers/qnspsa.py
+++ b/qiskit/algorithms/optimizers/qnspsa.py
@@ -44,6 +44,12 @@
     increasing the number of ``resamplings``. This leads to a Monte Carlo-style convergence to
     the exact, analytic value.
 
+    .. note::
+
+        This component has some function that is normally random. If you want to reproduce behavior
+        then you should set the random number generator seed in the algorithm_globals
+        (``qiskit.utils.algorithm_globals.random_seed = seed``).
+
     Examples:
 
         This short example runs QN-SPSA for the ground state calculation of the ``Z ^ Z``

{"golden_diff": "diff --git a/qiskit/algorithms/optimizers/qnspsa.py b/qiskit/algorithms/optimizers/qnspsa.py\n--- a/qiskit/algorithms/optimizers/qnspsa.py\n+++ b/qiskit/algorithms/optimizers/qnspsa.py\n@@ -44,6 +44,12 @@\n     increasing the number of ``resamplings``. This leads to a Monte Carlo-style convergence to\n     the exact, analytic value.\n \n+    .. note::\n+\n+        This component has some function that is normally random. If you want to reproduce behavior\n+        then you should set the random number generator seed in the algorithm_globals\n+        (``qiskit.utils.algorithm_globals.random_seed = seed``).\n+\n     Examples:\n \n         This short example runs QN-SPSA for the ground state calculation of the ``Z ^ Z``\n", "issue": "QNSPSA produces irreproducible results\n<!-- \u26a0\ufe0f If you do not respect this template, your issue will be closed -->\r\n<!-- \u26a0\ufe0f Make sure to browse the opened and closed issues -->\r\n\r\n### Information\r\n\r\n- **Qiskit Terra version**: qiskit-terra                  0.17.0 \r\n- **Python version**:                                             Python 3.7.3\r\n- **Operating system**: MacOS Big Sur             11.6\r\n\r\n### What is the current behavior?                     \r\nExecuting a `QNSPSA` optimization with the same random seed gives different results.\r\nWe checked that the same optimization gives reproducible results with other optimizers such as `ADAM`.\r\n\r\n### Steps to reproduce the problem\r\nRun an arbitrary optimization, e.g., VQE with QNSPSA, and compare the results.\r\n\r\n\r\n### What is the expected behavior?\r\nThe same results for the same seed.\r\n\r\n\r\n### Suggested solutions\r\n\r\n\r\n\n", "before_files": [{"content": "# This code is part of Qiskit.\n#\n# (C) Copyright IBM 2021.\n#\n# This code is licensed under the Apache License, Version 2.0. You may\n# obtain a copy of this license in the LICENSE.txt file in the root directory\n# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.\n#\n# Any modifications or derivative works of this code must retain this\n# copyright notice, and modified files need to carry a notice indicating\n# that they have been altered from the originals.\n\n\"\"\"The QN-SPSA optimizer.\"\"\"\n\nfrom typing import Any, Iterator, Optional, Union, Callable, Dict\n\nimport numpy as np\nfrom qiskit.providers import Backend\nfrom qiskit.circuit import ParameterVector, QuantumCircuit\nfrom qiskit.opflow import StateFn, CircuitSampler, ExpectationBase\nfrom qiskit.utils import QuantumInstance\n\nfrom .spsa import SPSA, CALLBACK, TERMINATIONCHECKER, _batch_evaluate\n\n# the function to compute the fidelity\nFIDELITY = Callable[[np.ndarray, np.ndarray], float]\n\n\nclass QNSPSA(SPSA):\n    r\"\"\"The Quantum Natural SPSA (QN-SPSA) optimizer.\n\n    The QN-SPSA optimizer [1] is a stochastic optimizer that belongs to the family of gradient\n    descent methods. This optimizer is based on SPSA but attempts to improve the convergence by\n    sampling the **natural gradient** instead of the vanilla, first-order gradient. It achieves\n    this by approximating Hessian of the ``fidelity`` of the ansatz circuit.\n\n    Compared to natural gradients, which require :math:`\\mathcal{O}(d^2)` expectation value\n    evaluations for a circuit with :math:`d` parameters, QN-SPSA only requires\n    :math:`\\mathcal{O}(1)` and can therefore significantly speed up the natural gradient calculation\n    by sacrificing some accuracy. Compared to SPSA, QN-SPSA requires 4 additional function\n    evaluations of the fidelity.\n\n    The stochastic approximation of the natural gradient can be systematically improved by\n    increasing the number of ``resamplings``. This leads to a Monte Carlo-style convergence to\n    the exact, analytic value.\n\n    Examples:\n\n        This short example runs QN-SPSA for the ground state calculation of the ``Z ^ Z``\n        observable where the ansatz is a ``PauliTwoDesign`` circuit.\n\n        .. code-block:: python\n\n            import numpy as np\n            from qiskit.algorithms.optimizers import QNSPSA\n            from qiskit.circuit.library import PauliTwoDesign\n            from qiskit.opflow import Z, StateFn\n\n            ansatz = PauliTwoDesign(2, reps=1, seed=2)\n            observable = Z ^ Z\n            initial_point = np.random.random(ansatz.num_parameters)\n\n            def loss(x):\n                bound = ansatz.bind_parameters(x)\n                return np.real((StateFn(observable, is_measurement=True) @ StateFn(bound)).eval())\n\n            fidelity = QNSPSA.get_fidelity(ansatz)\n            qnspsa = QNSPSA(fidelity, maxiter=300)\n            result = qnspsa.optimize(ansatz.num_parameters, loss, initial_point=initial_point)\n\n\n    References:\n\n        [1] J. Gacon et al, \"Simultaneous Perturbation Stochastic Approximation of the Quantum\n        Fisher Information\", `arXiv:2103.09232 <https://arxiv.org/abs/2103.09232>`_\n\n    \"\"\"\n\n    def __init__(\n        self,\n        fidelity: FIDELITY,\n        maxiter: int = 100,\n        blocking: bool = True,\n        allowed_increase: Optional[float] = None,\n        learning_rate: Optional[Union[float, Callable[[], Iterator]]] = None,\n        perturbation: Optional[Union[float, Callable[[], Iterator]]] = None,\n        last_avg: int = 1,\n        resamplings: Union[int, Dict[int, int]] = 1,\n        perturbation_dims: Optional[int] = None,\n        regularization: Optional[float] = None,\n        hessian_delay: int = 0,\n        lse_solver: Optional[Callable[[np.ndarray, np.ndarray], np.ndarray]] = None,\n        initial_hessian: Optional[np.ndarray] = None,\n        callback: Optional[CALLBACK] = None,\n        termination_checker: Optional[TERMINATIONCHECKER] = None,\n    ) -> None:\n        r\"\"\"\n        Args:\n            fidelity: A function to compute the fidelity of the ansatz state with itself for\n                two different sets of parameters.\n            maxiter: The maximum number of iterations. Note that this is not the maximal number\n                of function evaluations.\n            blocking: If True, only accepts updates that improve the loss (up to some allowed\n                increase, see next argument).\n            allowed_increase: If ``blocking`` is ``True``, this argument determines by how much\n                the loss can increase with the proposed parameters and still be accepted.\n                If ``None``, the allowed increases is calibrated automatically to be twice the\n                approximated standard deviation of the loss function.\n            learning_rate: The update step is the learning rate is multiplied with the gradient.\n                If the learning rate is a float, it remains constant over the course of the\n                optimization. It can also be a callable returning an iterator which yields the\n                learning rates for each optimization step.\n                If ``learning_rate`` is set ``perturbation`` must also be provided.\n            perturbation: Specifies the magnitude of the perturbation for the finite difference\n                approximation of the gradients. Can be either a float or a generator yielding\n                the perturbation magnitudes per step.\n                If ``perturbation`` is set ``learning_rate`` must also be provided.\n            last_avg: Return the average of the ``last_avg`` parameters instead of just the\n                last parameter values.\n            resamplings: The number of times the gradient (and Hessian) is sampled using a random\n                direction to construct a gradient estimate. Per default the gradient is estimated\n                using only one random direction. If an integer, all iterations use the same number\n                of resamplings. If a dictionary, this is interpreted as\n                ``{iteration: number of resamplings per iteration}``.\n            perturbation_dims: The number of perturbed dimensions. Per default, all dimensions\n                are perturbed, but a smaller, fixed number can be perturbed. If set, the perturbed\n                dimensions are chosen uniformly at random.\n            regularization: To ensure the preconditioner is symmetric and positive definite, the\n                identity times a small coefficient is added to it. This generator yields that\n                coefficient.\n            hessian_delay: Start multiplying the gradient with the inverse Hessian only after a\n                certain number of iterations. The Hessian is still evaluated and therefore this\n                argument can be useful to first get a stable average over the last iterations before\n                using it as preconditioner.\n            lse_solver: The method to solve for the inverse of the Hessian. Per default an\n                exact LSE solver is used, but can e.g. be overwritten by a minimization routine.\n            initial_hessian: The initial guess for the Hessian. By default the identity matrix\n                is used.\n            callback: A callback function passed information in each iteration step. The\n                information is, in this order: the parameters, the function value, the number\n                of function evaluations, the stepsize, whether the step was accepted.\n            termination_checker: A callback function executed at the end of each iteration step. The\n                arguments are, in this order: the parameters, the function value, the number\n                of function evaluations, the stepsize, whether the step was accepted. If the callback\n                returns True, the optimization is terminated.\n                To prevent additional evaluations of the objective method, if the objective has not yet\n                been evaluated, the objective is estimated by taking the mean of the objective\n                evaluations used in the estimate of the gradient.\n\n\n        \"\"\"\n        super().__init__(\n            maxiter,\n            blocking,\n            allowed_increase,\n            # trust region *must* be false for natural gradients to work\n            trust_region=False,\n            learning_rate=learning_rate,\n            perturbation=perturbation,\n            resamplings=resamplings,\n            callback=callback,\n            second_order=True,\n            hessian_delay=hessian_delay,\n            lse_solver=lse_solver,\n            regularization=regularization,\n            perturbation_dims=perturbation_dims,\n            initial_hessian=initial_hessian,\n            termination_checker=termination_checker,\n        )\n\n        self.fidelity = fidelity\n\n    def _point_sample(self, loss, x, eps, delta1, delta2):\n        loss_points = [x + eps * delta1, x - eps * delta1]\n        fidelity_points = [\n            (x, x + eps * delta1),\n            (x, x - eps * delta1),\n            (x, x + eps * (delta1 + delta2)),\n            (x, x + eps * (-delta1 + delta2)),\n        ]\n        self._nfev += 6\n\n        loss_values = _batch_evaluate(loss, loss_points, self._max_evals_grouped)\n        fidelity_values = _batch_evaluate(self.fidelity, fidelity_points, self._max_evals_grouped)\n\n        # compute the gradient approximation and additionally return the loss function evaluations\n        gradient_estimate = (loss_values[0] - loss_values[1]) / (2 * eps) * delta1\n\n        # compute the preconditioner point estimate\n        diff = fidelity_values[2] - fidelity_values[0]\n        diff -= fidelity_values[3] - fidelity_values[1]\n        diff /= 2 * eps ** 2\n\n        rank_one = np.outer(delta1, delta2)\n        # -0.5 factor comes from the fact that we need -0.5 * fidelity\n        hessian_estimate = -0.5 * diff * (rank_one + rank_one.T) / 2\n\n        return np.mean(loss_values), gradient_estimate, hessian_estimate\n\n    @property\n    def settings(self) -> Dict[str, Any]:\n        \"\"\"The optimizer settings in a dictionary format.\"\"\"\n        # re-use serialization from SPSA\n        settings = super().settings\n        settings.update({\"fidelity\": self.fidelity})\n\n        # remove SPSA-specific arguments not in QNSPSA\n        settings.pop(\"trust_region\")\n        settings.pop(\"second_order\")\n\n        return settings\n\n    @staticmethod\n    def get_fidelity(\n        circuit: QuantumCircuit,\n        backend: Optional[Union[Backend, QuantumInstance]] = None,\n        expectation: Optional[ExpectationBase] = None,\n    ) -> Callable[[np.ndarray, np.ndarray], float]:\n        r\"\"\"Get a function to compute the fidelity of ``circuit`` with itself.\n\n        Let ``circuit`` be a parameterized quantum circuit performing the operation\n        :math:`U(\\theta)` given a set of parameters :math:`\\theta`. Then this method returns\n        a function to evaluate\n\n        .. math::\n\n            F(\\theta, \\phi) = \\big|\\langle 0 | U^\\dagger(\\theta) U(\\phi) |0\\rangle  \\big|^2.\n\n        The output of this function can be used as input for the ``fidelity`` to the\n        :class:~`qiskit.algorithms.optimizers.QNSPSA` optimizer.\n\n        Args:\n            circuit: The circuit preparing the parameterized ansatz.\n            backend: A backend of quantum instance to evaluate the circuits. If None, plain\n                matrix multiplication will be used.\n            expectation: An expectation converter to specify how the expected value is computed.\n                If a shot-based readout is used this should be set to ``PauliExpectation``.\n\n        Returns:\n            A handle to the function :math:`F`.\n\n        \"\"\"\n        params_x = ParameterVector(\"x\", circuit.num_parameters)\n        params_y = ParameterVector(\"y\", circuit.num_parameters)\n\n        expression = ~StateFn(circuit.assign_parameters(params_x)) @ StateFn(\n            circuit.assign_parameters(params_y)\n        )\n\n        if expectation is not None:\n            expression = expectation.convert(expression)\n\n        if backend is None:\n\n            def fidelity(values_x, values_y):\n                value_dict = dict(\n                    zip(params_x[:] + params_y[:], values_x.tolist() + values_y.tolist())\n                )\n                return np.abs(expression.bind_parameters(value_dict).eval()) ** 2\n\n        else:\n            sampler = CircuitSampler(backend)\n\n            def fidelity(values_x, values_y=None):\n                if values_y is not None:  # no batches\n                    value_dict = dict(\n                        zip(params_x[:] + params_y[:], values_x.tolist() + values_y.tolist())\n                    )\n                else:\n                    value_dict = {p: [] for p in params_x[:] + params_y[:]}\n                    for values_xy in values_x:\n                        for value_x, param_x in zip(values_xy[0, :], params_x):\n                            value_dict[param_x].append(value_x)\n\n                        for value_y, param_y in zip(values_xy[1, :], params_y):\n                            value_dict[param_y].append(value_y)\n\n                return np.abs(sampler.convert(expression, params=value_dict).eval()) ** 2\n\n        return fidelity\n", "path": "qiskit/algorithms/optimizers/qnspsa.py"}], "after_files": [{"content": "# This code is part of Qiskit.\n#\n# (C) Copyright IBM 2021.\n#\n# This code is licensed under the Apache License, Version 2.0. You may\n# obtain a copy of this license in the LICENSE.txt file in the root directory\n# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.\n#\n# Any modifications or derivative works of this code must retain this\n# copyright notice, and modified files need to carry a notice indicating\n# that they have been altered from the originals.\n\n\"\"\"The QN-SPSA optimizer.\"\"\"\n\nfrom typing import Any, Iterator, Optional, Union, Callable, Dict\n\nimport numpy as np\nfrom qiskit.providers import Backend\nfrom qiskit.circuit import ParameterVector, QuantumCircuit\nfrom qiskit.opflow import StateFn, CircuitSampler, ExpectationBase\nfrom qiskit.utils import QuantumInstance\n\nfrom .spsa import SPSA, CALLBACK, TERMINATIONCHECKER, _batch_evaluate\n\n# the function to compute the fidelity\nFIDELITY = Callable[[np.ndarray, np.ndarray], float]\n\n\nclass QNSPSA(SPSA):\n    r\"\"\"The Quantum Natural SPSA (QN-SPSA) optimizer.\n\n    The QN-SPSA optimizer [1] is a stochastic optimizer that belongs to the family of gradient\n    descent methods. This optimizer is based on SPSA but attempts to improve the convergence by\n    sampling the **natural gradient** instead of the vanilla, first-order gradient. It achieves\n    this by approximating Hessian of the ``fidelity`` of the ansatz circuit.\n\n    Compared to natural gradients, which require :math:`\\mathcal{O}(d^2)` expectation value\n    evaluations for a circuit with :math:`d` parameters, QN-SPSA only requires\n    :math:`\\mathcal{O}(1)` and can therefore significantly speed up the natural gradient calculation\n    by sacrificing some accuracy. Compared to SPSA, QN-SPSA requires 4 additional function\n    evaluations of the fidelity.\n\n    The stochastic approximation of the natural gradient can be systematically improved by\n    increasing the number of ``resamplings``. This leads to a Monte Carlo-style convergence to\n    the exact, analytic value.\n\n    .. note::\n\n        This component has some function that is normally random. If you want to reproduce behavior\n        then you should set the random number generator seed in the algorithm_globals\n        (``qiskit.utils.algorithm_globals.random_seed = seed``).\n\n    Examples:\n\n        This short example runs QN-SPSA for the ground state calculation of the ``Z ^ Z``\n        observable where the ansatz is a ``PauliTwoDesign`` circuit.\n\n        .. code-block:: python\n\n            import numpy as np\n            from qiskit.algorithms.optimizers import QNSPSA\n            from qiskit.circuit.library import PauliTwoDesign\n            from qiskit.opflow import Z, StateFn\n\n            ansatz = PauliTwoDesign(2, reps=1, seed=2)\n            observable = Z ^ Z\n            initial_point = np.random.random(ansatz.num_parameters)\n\n            def loss(x):\n                bound = ansatz.bind_parameters(x)\n                return np.real((StateFn(observable, is_measurement=True) @ StateFn(bound)).eval())\n\n            fidelity = QNSPSA.get_fidelity(ansatz)\n            qnspsa = QNSPSA(fidelity, maxiter=300)\n            result = qnspsa.optimize(ansatz.num_parameters, loss, initial_point=initial_point)\n\n\n    References:\n\n        [1] J. Gacon et al, \"Simultaneous Perturbation Stochastic Approximation of the Quantum\n        Fisher Information\", `arXiv:2103.09232 <https://arxiv.org/abs/2103.09232>`_\n\n    \"\"\"\n\n    def __init__(\n        self,\n        fidelity: FIDELITY,\n        maxiter: int = 100,\n        blocking: bool = True,\n        allowed_increase: Optional[float] = None,\n        learning_rate: Optional[Union[float, Callable[[], Iterator]]] = None,\n        perturbation: Optional[Union[float, Callable[[], Iterator]]] = None,\n        last_avg: int = 1,\n        resamplings: Union[int, Dict[int, int]] = 1,\n        perturbation_dims: Optional[int] = None,\n        regularization: Optional[float] = None,\n        hessian_delay: int = 0,\n        lse_solver: Optional[Callable[[np.ndarray, np.ndarray], np.ndarray]] = None,\n        initial_hessian: Optional[np.ndarray] = None,\n        callback: Optional[CALLBACK] = None,\n        termination_checker: Optional[TERMINATIONCHECKER] = None,\n    ) -> None:\n        r\"\"\"\n        Args:\n            fidelity: A function to compute the fidelity of the ansatz state with itself for\n                two different sets of parameters.\n            maxiter: The maximum number of iterations. Note that this is not the maximal number\n                of function evaluations.\n            blocking: If True, only accepts updates that improve the loss (up to some allowed\n                increase, see next argument).\n            allowed_increase: If ``blocking`` is ``True``, this argument determines by how much\n                the loss can increase with the proposed parameters and still be accepted.\n                If ``None``, the allowed increases is calibrated automatically to be twice the\n                approximated standard deviation of the loss function.\n            learning_rate: The update step is the learning rate is multiplied with the gradient.\n                If the learning rate is a float, it remains constant over the course of the\n                optimization. It can also be a callable returning an iterator which yields the\n                learning rates for each optimization step.\n                If ``learning_rate`` is set ``perturbation`` must also be provided.\n            perturbation: Specifies the magnitude of the perturbation for the finite difference\n                approximation of the gradients. Can be either a float or a generator yielding\n                the perturbation magnitudes per step.\n                If ``perturbation`` is set ``learning_rate`` must also be provided.\n            last_avg: Return the average of the ``last_avg`` parameters instead of just the\n                last parameter values.\n            resamplings: The number of times the gradient (and Hessian) is sampled using a random\n                direction to construct a gradient estimate. Per default the gradient is estimated\n                using only one random direction. If an integer, all iterations use the same number\n                of resamplings. If a dictionary, this is interpreted as\n                ``{iteration: number of resamplings per iteration}``.\n            perturbation_dims: The number of perturbed dimensions. Per default, all dimensions\n                are perturbed, but a smaller, fixed number can be perturbed. If set, the perturbed\n                dimensions are chosen uniformly at random.\n            regularization: To ensure the preconditioner is symmetric and positive definite, the\n                identity times a small coefficient is added to it. This generator yields that\n                coefficient.\n            hessian_delay: Start multiplying the gradient with the inverse Hessian only after a\n                certain number of iterations. The Hessian is still evaluated and therefore this\n                argument can be useful to first get a stable average over the last iterations before\n                using it as preconditioner.\n            lse_solver: The method to solve for the inverse of the Hessian. Per default an\n                exact LSE solver is used, but can e.g. be overwritten by a minimization routine.\n            initial_hessian: The initial guess for the Hessian. By default the identity matrix\n                is used.\n            callback: A callback function passed information in each iteration step. The\n                information is, in this order: the parameters, the function value, the number\n                of function evaluations, the stepsize, whether the step was accepted.\n            termination_checker: A callback function executed at the end of each iteration step. The\n                arguments are, in this order: the parameters, the function value, the number\n                of function evaluations, the stepsize, whether the step was accepted. If the callback\n                returns True, the optimization is terminated.\n                To prevent additional evaluations of the objective method, if the objective has not yet\n                been evaluated, the objective is estimated by taking the mean of the objective\n                evaluations used in the estimate of the gradient.\n\n\n        \"\"\"\n        super().__init__(\n            maxiter,\n            blocking,\n            allowed_increase,\n            # trust region *must* be false for natural gradients to work\n            trust_region=False,\n            learning_rate=learning_rate,\n            perturbation=perturbation,\n            resamplings=resamplings,\n            callback=callback,\n            second_order=True,\n            hessian_delay=hessian_delay,\n            lse_solver=lse_solver,\n            regularization=regularization,\n            perturbation_dims=perturbation_dims,\n            initial_hessian=initial_hessian,\n            termination_checker=termination_checker,\n        )\n\n        self.fidelity = fidelity\n\n    def _point_sample(self, loss, x, eps, delta1, delta2):\n        loss_points = [x + eps * delta1, x - eps * delta1]\n        fidelity_points = [\n            (x, x + eps * delta1),\n            (x, x - eps * delta1),\n            (x, x + eps * (delta1 + delta2)),\n            (x, x + eps * (-delta1 + delta2)),\n        ]\n        self._nfev += 6\n\n        loss_values = _batch_evaluate(loss, loss_points, self._max_evals_grouped)\n        fidelity_values = _batch_evaluate(self.fidelity, fidelity_points, self._max_evals_grouped)\n\n        # compute the gradient approximation and additionally return the loss function evaluations\n        gradient_estimate = (loss_values[0] - loss_values[1]) / (2 * eps) * delta1\n\n        # compute the preconditioner point estimate\n        diff = fidelity_values[2] - fidelity_values[0]\n        diff -= fidelity_values[3] - fidelity_values[1]\n        diff /= 2 * eps ** 2\n\n        rank_one = np.outer(delta1, delta2)\n        # -0.5 factor comes from the fact that we need -0.5 * fidelity\n        hessian_estimate = -0.5 * diff * (rank_one + rank_one.T) / 2\n\n        return np.mean(loss_values), gradient_estimate, hessian_estimate\n\n    @property\n    def settings(self) -> Dict[str, Any]:\n        \"\"\"The optimizer settings in a dictionary format.\"\"\"\n        # re-use serialization from SPSA\n        settings = super().settings\n        settings.update({\"fidelity\": self.fidelity})\n\n        # remove SPSA-specific arguments not in QNSPSA\n        settings.pop(\"trust_region\")\n        settings.pop(\"second_order\")\n\n        return settings\n\n    @staticmethod\n    def get_fidelity(\n        circuit: QuantumCircuit,\n        backend: Optional[Union[Backend, QuantumInstance]] = None,\n        expectation: Optional[ExpectationBase] = None,\n    ) -> Callable[[np.ndarray, np.ndarray], float]:\n        r\"\"\"Get a function to compute the fidelity of ``circuit`` with itself.\n\n        Let ``circuit`` be a parameterized quantum circuit performing the operation\n        :math:`U(\\theta)` given a set of parameters :math:`\\theta`. Then this method returns\n        a function to evaluate\n\n        .. math::\n\n            F(\\theta, \\phi) = \\big|\\langle 0 | U^\\dagger(\\theta) U(\\phi) |0\\rangle  \\big|^2.\n\n        The output of this function can be used as input for the ``fidelity`` to the\n        :class:~`qiskit.algorithms.optimizers.QNSPSA` optimizer.\n\n        Args:\n            circuit: The circuit preparing the parameterized ansatz.\n            backend: A backend of quantum instance to evaluate the circuits. If None, plain\n                matrix multiplication will be used.\n            expectation: An expectation converter to specify how the expected value is computed.\n                If a shot-based readout is used this should be set to ``PauliExpectation``.\n\n        Returns:\n            A handle to the function :math:`F`.\n\n        \"\"\"\n        params_x = ParameterVector(\"x\", circuit.num_parameters)\n        params_y = ParameterVector(\"y\", circuit.num_parameters)\n\n        expression = ~StateFn(circuit.assign_parameters(params_x)) @ StateFn(\n            circuit.assign_parameters(params_y)\n        )\n\n        if expectation is not None:\n            expression = expectation.convert(expression)\n\n        if backend is None:\n\n            def fidelity(values_x, values_y):\n                value_dict = dict(\n                    zip(params_x[:] + params_y[:], values_x.tolist() + values_y.tolist())\n                )\n                return np.abs(expression.bind_parameters(value_dict).eval()) ** 2\n\n        else:\n            sampler = CircuitSampler(backend)\n\n            def fidelity(values_x, values_y=None):\n                if values_y is not None:  # no batches\n                    value_dict = dict(\n                        zip(params_x[:] + params_y[:], values_x.tolist() + values_y.tolist())\n                    )\n                else:\n                    value_dict = {p: [] for p in params_x[:] + params_y[:]}\n                    for values_xy in values_x:\n                        for value_x, param_x in zip(values_xy[0, :], params_x):\n                            value_dict[param_x].append(value_x)\n\n                        for value_y, param_y in zip(values_xy[1, :], params_y):\n                            value_dict[param_y].append(value_y)\n\n                return np.abs(sampler.convert(expression, params=value_dict).eval()) ** 2\n\n        return fidelity\n", "path": "qiskit/algorithms/optimizers/qnspsa.py"}]}