Repository: Qiskit/qiskit-terra (python)

Issue:
Add forward and backward options to finite difference gradient
### What should we add?

Currently, `FiniteDiffEstimatorGradient` and `FiniteDiffSamplerGradient` support only a central difference. A forward difference or backward difference require less circuit than the central difference, so they are faster. Also, there are cases that forward difference or backward difference can work better than central difference. 
https://en.wikipedia.org/wiki/Finite_difference

Pull request description:
Add gradient methods to FinDiffGradients
<!--
⚠️ If you do not respect this template, your pull request will be closed.
⚠️ Your pull request title should be short detailed and understandable for all.
⚠️ Also, please add a release note file using reno if the change needs to be
  documented in the release notes.
⚠️ If your pull request fixes an open issue, please link to the issue.

- [ ] I have added the tests to cover my changes.
- [ ] I have updated the documentation accordingly.
- [ ] I have read the CONTRIBUTING document.
-->

### Summary

Resolves #9096

### Details and comments




Patch:
diff --git a/qiskit/algorithms/gradients/finite_diff_estimator_gradient.py b/qiskit/algorithms/gradients/finite_diff_estimator_gradient.py
index 5ab14150b..60ff8d9a7 100644
--- a/qiskit/algorithms/gradients/finite_diff_estimator_gradient.py
+++ b/qiskit/algorithms/gradients/finite_diff_estimator_gradient.py
@@ -14,6 +14,7 @@
 
 from __future__ import annotations
 
+import sys
 from typing import Sequence
 
 import numpy as np
@@ -28,14 +29,27 @@ from qiskit.quantum_info.operators.base_operator import BaseOperator
 from .base_estimator_gradient import BaseEstimatorGradient
 from .estimator_gradient_result import EstimatorGradientResult
 
+if sys.version_info >= (3, 8):
+    # pylint: disable=no-name-in-module, ungrouped-imports
+    from typing import Literal
+else:
+    from typing_extensions import Literal
+
 
 class FiniteDiffEstimatorGradient(BaseEstimatorGradient):
     """
     Compute the gradients of the expectation values by finite difference method.
     """
 
-    def __init__(self, estimator: BaseEstimator, epsilon: float, options: Options | None = None):
-        """
+    def __init__(
+        self,
+        estimator: BaseEstimator,
+        epsilon: float,
+        options: Options | None = None,
+        *,
+        method: Literal["central", "forward", "backward"] = "central",
+    ):
+        r"""
         Args:
             estimator: The estimator used to compute the gradients.
             epsilon: The offset size for the finite difference gradients.
@@ -43,14 +57,27 @@ class FiniteDiffEstimatorGradient(BaseEstimatorGradient):
                 The order of priority is: options in ``run`` method > gradient's
                 default options > primitive's default setting.
                 Higher priority setting overrides lower priority setting
+            method: The computation method of the gradients.
+
+                    - ``central`` computes :math:`\frac{f(x+e)-f(x-e)}{2e}`,
+                    - ``forward`` computes :math:`\frac{f(x+e) - f(x)}{e}`,
+                    - ``backward`` computes :math:`\frac{f(x)-f(x-e)}{e}`
+
+                where :math:`e` is epsilon.
 
         Raises:
             ValueError: If ``epsilon`` is not positive.
+            TypeError: If ``method`` is invalid.
         """
         if epsilon <= 0:
             raise ValueError(f"epsilon ({epsilon}) should be positive.")
         self._epsilon = epsilon
         self._base_parameter_values_dict = {}
+        if method not in ("central", "forward", "backward"):
+            raise TypeError(
+                f"The argument method should be central, forward, or backward: {method} is given."
+            )
+        self._method = method
         super().__init__(estimator, options)
 
     def _run(
@@ -74,12 +101,25 @@ class FiniteDiffEstimatorGradient(BaseEstimatorGradient):
             metadata_.append({"parameters": [circuit.parameters[idx] for idx in indices]})
 
             offset = np.identity(circuit.num_parameters)[indices, :]
-            plus = parameter_values_ + self._epsilon * offset
-            minus = parameter_values_ - self._epsilon * offset
-            n = 2 * len(indices)
-            job = self._estimator.run(
-                [circuit] * n, [observable] * n, plus.tolist() + minus.tolist(), **options
-            )
+            if self._method == "central":
+                plus = parameter_values_ + self._epsilon * offset
+                minus = parameter_values_ - self._epsilon * offset
+                n = 2 * len(indices)
+                job = self._estimator.run(
+                    [circuit] * n, [observable] * n, plus.tolist() + minus.tolist(), **options
+                )
+            elif self._method == "forward":
+                plus = parameter_values_ + self._epsilon * offset
+                n = len(indices) + 1
+                job = self._estimator.run(
+                    [circuit] * n, [observable] * n, [parameter_values_] + plus.tolist(), **options
+                )
+            elif self._method == "backward":
+                minus = parameter_values_ - self._epsilon * offset
+                n = len(indices) + 1
+                job = self._estimator.run(
+                    [circuit] * n, [observable] * n, [parameter_values_] + minus.tolist(), **options
+                )
             jobs.append(job)
 
         # combine the results
@@ -90,8 +130,13 @@ class FiniteDiffEstimatorGradient(BaseEstimatorGradient):
 
         gradients = []
         for result in results:
-            n = len(result.values) // 2  # is always a multiple of 2
-            gradient_ = (result.values[:n] - result.values[n:]) / (2 * self._epsilon)
+            if self._method == "central":
+                n = len(result.values) // 2  # is always a multiple of 2
+                gradient_ = (result.values[:n] - result.values[n:]) / (2 * self._epsilon)
+            elif self._method == "forward":
+                gradient_ = (result.values[1:] - result.values[0]) / self._epsilon
+            elif self._method == "backward":
+                gradient_ = (result.values[0] - result.values[1:]) / self._epsilon
             gradients.append(gradient_)
         opt = self._get_local_options(options)
         return EstimatorGradientResult(gradients=gradients, metadata=metadata_, options=opt)
diff --git a/qiskit/algorithms/gradients/finite_diff_sampler_gradient.py b/qiskit/algorithms/gradients/finite_diff_sampler_gradient.py
index c5c9fdc00..784ffc427 100644
--- a/qiskit/algorithms/gradients/finite_diff_sampler_gradient.py
+++ b/qiskit/algorithms/gradients/finite_diff_sampler_gradient.py
@@ -14,6 +14,7 @@
 
 from __future__ import annotations
 
+import sys
 from typing import Sequence
 
 import numpy as np
@@ -26,6 +27,12 @@ from qiskit.providers import Options
 from .base_sampler_gradient import BaseSamplerGradient
 from .sampler_gradient_result import SamplerGradientResult
 
+if sys.version_info >= (3, 8):
+    # pylint: disable=no-name-in-module, ungrouped-imports
+    from typing import Literal
+else:
+    from typing_extensions import Literal
+
 
 class FiniteDiffSamplerGradient(BaseSamplerGradient):
     """Compute the gradients of the sampling probability by finite difference method."""
@@ -35,8 +42,10 @@ class FiniteDiffSamplerGradient(BaseSamplerGradient):
         sampler: BaseSampler,
         epsilon: float,
         options: Options | None = None,
+        *,
+        method: Literal["central", "forward", "backward"] = "central",
     ):
-        """
+        r"""
         Args:
             sampler: The sampler used to compute the gradients.
             epsilon: The offset size for the finite difference gradients.
@@ -44,13 +53,26 @@ class FiniteDiffSamplerGradient(BaseSamplerGradient):
                 The order of priority is: options in ``run`` method > gradient's
                 default options > primitive's default setting.
                 Higher priority setting overrides lower priority setting
+            method: The computation method of the gradients.
+
+                    - ``central`` computes :math:`\frac{f(x+e)-f(x-e)}{2e}`,
+                    - ``forward`` computes :math:`\frac{f(x+e) - f(x)}{e}`,
+                    - ``backward`` computes :math:`\frac{f(x)-f(x-e)}{e}`
+
+                where :math:`e` is epsilon.
 
         Raises:
             ValueError: If ``epsilon`` is not positive.
+            TypeError: If ``method`` is invalid.
         """
         if epsilon <= 0:
             raise ValueError(f"epsilon ({epsilon}) should be positive.")
         self._epsilon = epsilon
+        if method not in ("central", "forward", "backward"):
+            raise TypeError(
+                f"The argument method should be central, forward, or backward: {method} is given."
+            )
+        self._method = method
         super().__init__(sampler, options)
 
     def _run(
@@ -70,10 +92,23 @@ class FiniteDiffSamplerGradient(BaseSamplerGradient):
                 indices = [circuit.parameters.data.index(p) for p in parameters_]
             metadata_.append({"parameters": [circuit.parameters[idx] for idx in indices]})
             offset = np.identity(circuit.num_parameters)[indices, :]
-            plus = parameter_values_ + self._epsilon * offset
-            minus = parameter_values_ - self._epsilon * offset
-            n = 2 * len(indices)
-            job = self._sampler.run([circuit] * n, plus.tolist() + minus.tolist(), **options)
+            if self._method == "central":
+                plus = parameter_values_ + self._epsilon * offset
+                minus = parameter_values_ - self._epsilon * offset
+                n = 2 * len(indices)
+                job = self._sampler.run([circuit] * n, plus.tolist() + minus.tolist(), **options)
+            elif self._method == "forward":
+                plus = parameter_values_ + self._epsilon * offset
+                n = len(indices) + 1
+                job = self._sampler.run(
+                    [circuit] * n, [parameter_values_] + plus.tolist(), **options
+                )
+            elif self._method == "backward":
+                minus = parameter_values_ - self._epsilon * offset
+                n = len(indices) + 1
+                job = self._sampler.run(
+                    [circuit] * n, [parameter_values_] + minus.tolist(), **options
+                )
             jobs.append(job)
 
         # combine the results
@@ -84,14 +119,33 @@ class FiniteDiffSamplerGradient(BaseSamplerGradient):
 
         gradients = []
         for i, result in enumerate(results):
-            n = len(result.quasi_dists) // 2
-            gradient_ = []
-            for dist_plus, dist_minus in zip(result.quasi_dists[:n], result.quasi_dists[n:]):
-                grad_dist = np.zeros(2 ** circuits[i].num_qubits)
-                grad_dist[list(dist_plus.keys())] += list(dist_plus.values())
-                grad_dist[list(dist_minus.keys())] -= list(dist_minus.values())
-                grad_dist /= 2 * self._epsilon
-                gradient_.append(dict(enumerate(grad_dist)))
+            if self._method == "central":
+                n = len(result.quasi_dists) // 2
+                gradient_ = []
+                for dist_plus, dist_minus in zip(result.quasi_dists[:n], result.quasi_dists[n:]):
+                    grad_dist = np.zeros(2 ** circuits[i].num_qubits)
+                    grad_dist[list(dist_plus.keys())] += list(dist_plus.values())
+                    grad_dist[list(dist_minus.keys())] -= list(dist_minus.values())
+                    grad_dist /= 2 * self._epsilon
+                    gradient_.append(dict(enumerate(grad_dist)))
+            elif self._method == "forward":
+                gradient_ = []
+                dist_zero = result.quasi_dists[0]
+                for dist_plus in result.quasi_dists[1:]:
+                    grad_dist = np.zeros(2 ** circuits[i].num_qubits)
+                    grad_dist[list(dist_plus.keys())] += list(dist_plus.values())
+                    grad_dist[list(dist_zero.keys())] -= list(dist_zero.values())
+                    grad_dist /= self._epsilon
+                    gradient_.append(dict(enumerate(grad_dist)))
+            elif self._method == "backward":
+                gradient_ = []
+                dist_zero = result.quasi_dists[0]
+                for dist_minus in result.quasi_dists[1:]:
+                    grad_dist = np.zeros(2 ** circuits[i].num_qubits)
+                    grad_dist[list(dist_zero.keys())] += list(dist_zero.values())
+                    grad_dist[list(dist_minus.keys())] -= list(dist_minus.values())
+                    grad_dist /= self._epsilon
+                    gradient_.append(dict(enumerate(grad_dist)))
             gradients.append(gradient_)
 
         opt = self._get_local_options(options)
diff --git a/releasenotes/notes/gradient-methods-b2ec34916b83c17b.yaml b/releasenotes/notes/gradient-methods-b2ec34916b83c17b.yaml
new file mode 100644
index 000000000..fedc7c41d
--- /dev/null
+++ b/releasenotes/notes/gradient-methods-b2ec34916b83c17b.yaml
@@ -0,0 +1,10 @@
+---
+features:
+  - |
+    :class:`.FiniteDiffEstimatorGradient` and :class:`FiniteDiffSamplerGradient`
+    have new argument method.
+    There are three methods, "central", "forward", and "backward".
+    This option changes the gradient calculation methods.
+    "central" calculates :math:`\frac{f(x+e)-f(x-e)}{2e}`, "forward"
+    :math:`\frac{f(x+e) - f(x)}{e}`, and "backward" :math:`\frac{f(x)-f(x-e)}{e}` where
+    :math:`e` is the offset epsilon.


Test patch:
diff --git a/test/python/algorithms/gradients/test_estimator_gradient.py b/test/python/algorithms/gradients/test_estimator_gradient.py
index e2aae25be..74bd83738 100644
--- a/test/python/algorithms/gradients/test_estimator_gradient.py
+++ b/test/python/algorithms/gradients/test_estimator_gradient.py
@@ -35,14 +35,20 @@ from qiskit.quantum_info import Operator, SparsePauliOp, Pauli
 from qiskit.quantum_info.random import random_pauli_list
 from qiskit.test import QiskitTestCase
 
+gradient_factories = [
+    lambda estimator: FiniteDiffEstimatorGradient(estimator, epsilon=1e-6, method="central"),
+    lambda estimator: FiniteDiffEstimatorGradient(estimator, epsilon=1e-6, method="forward"),
+    lambda estimator: FiniteDiffEstimatorGradient(estimator, epsilon=1e-6, method="backward"),
+    ParamShiftEstimatorGradient,
+    LinCombEstimatorGradient,
+]
+
 
 @ddt
 class TestEstimatorGradient(QiskitTestCase):
     """Test Estimator Gradient"""
 
-    @combine(
-        grad=[FiniteDiffEstimatorGradient, ParamShiftEstimatorGradient, LinCombEstimatorGradient]
-    )
+    @combine(grad=gradient_factories)
     def test_gradient_operators(self, grad):
         """Test the estimator gradient for different operators"""
         estimator = Estimator()
@@ -51,10 +57,7 @@ class TestEstimatorGradient(QiskitTestCase):
         qc.h(0)
         qc.p(a, 0)
         qc.h(0)
-        if grad is FiniteDiffEstimatorGradient:
-            gradient = grad(estimator, epsilon=1e-6)
-        else:
-            gradient = grad(estimator)
+        gradient = grad(estimator)
         op = SparsePauliOp.from_list([("Z", 1)])
         correct_result = -1 / np.sqrt(2)
         param = [np.pi / 4]
@@ -67,9 +70,7 @@ class TestEstimatorGradient(QiskitTestCase):
         value = gradient.run([qc], [op], [param]).result().gradients[0]
         self.assertAlmostEqual(value[0], correct_result, 3)
 
-    @combine(
-        grad=[FiniteDiffEstimatorGradient, ParamShiftEstimatorGradient, LinCombEstimatorGradient]
-    )
+    @combine(grad=gradient_factories)
     def test_gradient_p(self, grad):
         """Test the estimator gradient for p"""
         estimator = Estimator()
@@ -78,10 +79,7 @@ class TestEstimatorGradient(QiskitTestCase):
         qc.h(0)
         qc.p(a, 0)
         qc.h(0)
-        if grad is FiniteDiffEstimatorGradient:
-            gradient = grad(estimator, epsilon=1e-6)
-        else:
-            gradient = grad(estimator)
+        gradient = grad(estimator)
         op = SparsePauliOp.from_list([("Z", 1)])
         param_list = [[np.pi / 4], [0], [np.pi / 2]]
         correct_results = [[-1 / np.sqrt(2)], [0], [-1]]
@@ -90,9 +88,7 @@ class TestEstimatorGradient(QiskitTestCase):
             for j, value in enumerate(gradients):
                 self.assertAlmostEqual(value, correct_results[i][j], 3)
 
-    @combine(
-        grad=[FiniteDiffEstimatorGradient, ParamShiftEstimatorGradient, LinCombEstimatorGradient]
-    )
+    @combine(grad=gradient_factories)
     def test_gradient_u(self, grad):
         """Test the estimator gradient for u"""
         estimator = Estimator()
@@ -103,10 +99,7 @@ class TestEstimatorGradient(QiskitTestCase):
         qc.h(0)
         qc.u(a, b, c, 0)
         qc.h(0)
-        if grad is FiniteDiffEstimatorGradient:
-            gradient = grad(estimator, epsilon=1e-6)
-        else:
-            gradient = grad(estimator)
+        gradient = grad(estimator)
         op = SparsePauliOp.from_list([("Z", 1)])
 
         param_list = [[np.pi / 4, 0, 0], [np.pi / 4, np.pi / 4, np.pi / 4]]
@@ -116,18 +109,13 @@ class TestEstimatorGradient(QiskitTestCase):
             for j, value in enumerate(gradients):
                 self.assertAlmostEqual(value, correct_results[i][j], 3)
 
-    @combine(
-        grad=[FiniteDiffEstimatorGradient, ParamShiftEstimatorGradient, LinCombEstimatorGradient]
-    )
+    @combine(grad=gradient_factories)
     def test_gradient_efficient_su2(self, grad):
         """Test the estimator gradient for EfficientSU2"""
         estimator = Estimator()
         qc = EfficientSU2(2, reps=1)
         op = SparsePauliOp.from_list([("ZI", 1)])
-        if grad is FiniteDiffEstimatorGradient:
-            gradient = grad(estimator, epsilon=1e-6)
-        else:
-            gradient = grad(estimator)
+        gradient = grad(estimator)
         param_list = [
             [np.pi / 4 for param in qc.parameters],
             [np.pi / 2 for param in qc.parameters],
@@ -149,9 +137,7 @@ class TestEstimatorGradient(QiskitTestCase):
             gradients = gradient.run([qc], [op], [param]).result().gradients[0]
             np.testing.assert_allclose(gradients, correct_results[i], atol=1e-3)
 
-    @combine(
-        grad=[FiniteDiffEstimatorGradient, ParamShiftEstimatorGradient, LinCombEstimatorGradient],
-    )
+    @combine(grad=gradient_factories)
     def test_gradient_2qubit_gate(self, grad):
         """Test the estimator gradient for 2 qubit gates"""
         estimator = Estimator()
@@ -165,10 +151,7 @@ class TestEstimatorGradient(QiskitTestCase):
             for i, param in enumerate(param_list):
                 a = Parameter("a")
                 qc = QuantumCircuit(2)
-                if grad is FiniteDiffEstimatorGradient:
-                    gradient = grad(estimator, epsilon=1e-6)
-                else:
-                    gradient = grad(estimator)
+                gradient = grad(estimator)
 
                 if gate is RZZGate:
                     qc.h([0, 1])
@@ -179,9 +162,7 @@ class TestEstimatorGradient(QiskitTestCase):
                 gradients = gradient.run([qc], [op], [param]).result().gradients[0]
                 np.testing.assert_allclose(gradients, correct_results[i], atol=1e-3)
 
-    @combine(
-        grad=[FiniteDiffEstimatorGradient, ParamShiftEstimatorGradient, LinCombEstimatorGradient]
-    )
+    @combine(grad=gradient_factories)
     def test_gradient_parameter_coefficient(self, grad):
         """Test the estimator gradient for parameter variables with coefficients"""
         estimator = Estimator()
@@ -191,10 +172,7 @@ class TestEstimatorGradient(QiskitTestCase):
         qc.u(qc.parameters[0], qc.parameters[1], qc.parameters[3], 1)
         qc.p(2 * qc.parameters[0] + 1, 0)
         qc.rxx(qc.parameters[0] + 2, 0, 1)
-        if grad is FiniteDiffEstimatorGradient:
-            gradient = grad(estimator, epsilon=1e-6)
-        else:
-            gradient = grad(estimator)
+        gradient = grad(estimator)
         param_list = [[np.pi / 4 for _ in qc.parameters], [np.pi / 2 for _ in qc.parameters]]
         correct_results = [
             [-0.7266653, -0.4905135, -0.0068606, -0.9228880],
@@ -205,9 +183,7 @@ class TestEstimatorGradient(QiskitTestCase):
             gradients = gradient.run([qc], [op], [param]).result().gradients[0]
             np.testing.assert_allclose(gradients, correct_results[i], atol=1e-3)
 
-    @combine(
-        grad=[FiniteDiffEstimatorGradient, ParamShiftEstimatorGradient, LinCombEstimatorGradient]
-    )
+    @combine(grad=gradient_factories)
     def test_gradient_parameters(self, grad):
         """Test the estimator gradient for parameters"""
         estimator = Estimator()
@@ -216,10 +192,7 @@ class TestEstimatorGradient(QiskitTestCase):
         qc = QuantumCircuit(1)
         qc.rx(a, 0)
         qc.rx(b, 0)
-        if grad is FiniteDiffEstimatorGradient:
-            gradient = grad(estimator, epsilon=1e-6)
-        else:
-            gradient = grad(estimator)
+        gradient = grad(estimator)
         param_list = [[np.pi / 4, np.pi / 2]]
         correct_results = [
             [-0.70710678],
@@ -229,9 +202,7 @@ class TestEstimatorGradient(QiskitTestCase):
             gradients = gradient.run([qc], [op], [param], parameters=[[a]]).result().gradients[0]
             np.testing.assert_allclose(gradients, correct_results[i], atol=1e-3)
 
-    @combine(
-        grad=[FiniteDiffEstimatorGradient, ParamShiftEstimatorGradient, LinCombEstimatorGradient]
-    )
+    @combine(grad=gradient_factories)
     def test_gradient_multi_arguments(self, grad):
         """Test the estimator gradient for multiple arguments"""
         estimator = Estimator()
@@ -241,10 +212,7 @@ class TestEstimatorGradient(QiskitTestCase):
         qc.rx(a, 0)
         qc2 = QuantumCircuit(1)
         qc2.rx(b, 0)
-        if grad is FiniteDiffEstimatorGradient:
-            gradient = grad(estimator, epsilon=1e-6)
-        else:
-            gradient = grad(estimator)
+        gradient = grad(estimator)
         param_list = [[np.pi / 4], [np.pi / 2]]
         correct_results = [
             [-0.70710678],
diff --git a/test/python/algorithms/gradients/test_sampler_gradient.py b/test/python/algorithms/gradients/test_sampler_gradient.py
index bfd4453c9..6bdc120cb 100644
--- a/test/python/algorithms/gradients/test_sampler_gradient.py
+++ b/test/python/algorithms/gradients/test_sampler_gradient.py
@@ -34,12 +34,20 @@ from qiskit.primitives import Sampler
 from qiskit.result import QuasiDistribution
 from qiskit.test import QiskitTestCase
 
+gradient_factories = [
+    lambda sampler: FiniteDiffSamplerGradient(sampler, epsilon=1e-6, method="central"),
+    lambda sampler: FiniteDiffSamplerGradient(sampler, epsilon=1e-6, method="forward"),
+    lambda sampler: FiniteDiffSamplerGradient(sampler, epsilon=1e-6, method="backward"),
+    ParamShiftSamplerGradient,
+    LinCombSamplerGradient,
+]
+
 
 @ddt
 class TestSamplerGradient(QiskitTestCase):
     """Test Sampler Gradient"""
 
-    @combine(grad=[FiniteDiffSamplerGradient, ParamShiftSamplerGradient, LinCombSamplerGradient])
+    @combine(grad=gradient_factories)
     def test_gradient_p(self, grad):
         """Test the sampler gradient for p"""
         sampler = Sampler()
@@ -49,10 +57,7 @@ class TestSamplerGradient(QiskitTestCase):
         qc.p(a, 0)
         qc.h(0)
         qc.measure_all()
-        if grad is FiniteDiffSamplerGradient:
-            gradient = grad(sampler, epsilon=1e-6)
-        else:
-            gradient = grad(sampler)
+        gradient = grad(sampler)
         param_list = [[np.pi / 4], [0], [np.pi / 2]]
         correct_results = [
             [{0: -0.5 / np.sqrt(2), 1: 0.5 / np.sqrt(2)}],
@@ -65,7 +70,7 @@ class TestSamplerGradient(QiskitTestCase):
                 for k in quasi_dist:
                     self.assertAlmostEqual(quasi_dist[k], correct_results[i][j][k], 3)
 
-    @combine(grad=[FiniteDiffSamplerGradient, ParamShiftSamplerGradient, LinCombSamplerGradient])
+    @combine(grad=gradient_factories)
     def test_gradient_u(self, grad):
         """Test the sampler gradient for u"""
         sampler = Sampler()
@@ -77,10 +82,7 @@ class TestSamplerGradient(QiskitTestCase):
         qc.u(a, b, c, 0)
         qc.h(0)
         qc.measure_all()
-        if grad is FiniteDiffSamplerGradient:
-            gradient = grad(sampler, epsilon=1e-6)
-        else:
-            gradient = grad(sampler)
+        gradient = grad(sampler)
         param_list = [[np.pi / 4, 0, 0], [np.pi / 4, np.pi / 4, np.pi / 4]]
         correct_results = [
             [{0: -0.5 / np.sqrt(2), 1: 0.5 / np.sqrt(2)}, {0: 0, 1: 0}, {0: 0, 1: 0}],
@@ -92,16 +94,13 @@ class TestSamplerGradient(QiskitTestCase):
                 for k in quasi_dist:
                     self.assertAlmostEqual(quasi_dist[k], correct_results[i][j][k], 3)
 
-    @combine(grad=[FiniteDiffSamplerGradient, ParamShiftSamplerGradient, LinCombSamplerGradient])
+    @combine(grad=gradient_factories)
     def test_gradient_efficient_su2(self, grad):
         """Test the sampler gradient for EfficientSU2"""
         sampler = Sampler()
         qc = EfficientSU2(2, reps=1)
         qc.measure_all()
-        if grad is FiniteDiffSamplerGradient:
-            gradient = grad(sampler, epsilon=1e-6)
-        else:
-            gradient = grad(sampler)
+        gradient = grad(sampler)
         param_list = [
             [np.pi / 4 for param in qc.parameters],
             [np.pi / 2 for param in qc.parameters],
@@ -189,7 +188,7 @@ class TestSamplerGradient(QiskitTestCase):
                 for k in quasi_dist:
                     self.assertAlmostEqual(quasi_dist[k], correct_results[i][j][k], 3)
 
-    @combine(grad=[FiniteDiffSamplerGradient, ParamShiftSamplerGradient, LinCombSamplerGradient])
+    @combine(grad=gradient_factories)
     def test_gradient_2qubit_gate(self, grad):
         """Test the sampler gradient for 2 qubit gates"""
         sampler = Sampler()
@@ -211,16 +210,13 @@ class TestSamplerGradient(QiskitTestCase):
             qc = QuantumCircuit(2)
             qc.append(gate(a), [qc.qubits[0], qc.qubits[1]], [])
             qc.measure_all()
-            if grad is FiniteDiffSamplerGradient:
-                gradient = grad(sampler, epsilon=1e-6)
-            else:
-                gradient = grad(sampler)
+            gradient = grad(sampler)
             gradients = gradient.run([qc], [param]).result().gradients[0]
             for j, quasi_dist in enumerate(gradients):
                 for k in quasi_dist:
                     self.assertAlmostEqual(quasi_dist[k], correct_results[i][j][k], 3)
 
-    @combine(grad=[FiniteDiffSamplerGradient, ParamShiftSamplerGradient, LinCombSamplerGradient])
+    @combine(grad=gradient_factories)
     def test_gradient_parameter_coefficient(self, grad):
         """Test the sampler gradient for parameter variables with coefficients"""
         sampler = Sampler()
@@ -231,10 +227,7 @@ class TestSamplerGradient(QiskitTestCase):
         qc.p(2 * qc.parameters[0] + 1, 0)
         qc.rxx(qc.parameters[0] + 2, 0, 1)
         qc.measure_all()
-        if grad is FiniteDiffSamplerGradient:
-            gradient = grad(sampler, epsilon=1e-6)
-        else:
-            gradient = grad(sampler)
+        gradient = grad(sampler)
         param_list = [[np.pi / 4 for _ in qc.parameters], [np.pi / 2 for _ in qc.parameters]]
         correct_results = [
             [
@@ -297,7 +290,7 @@ class TestSamplerGradient(QiskitTestCase):
                 for k in quasi_dist:
                     self.assertAlmostEqual(quasi_dist[k], correct_results[i][j][k], 2)
 
-    @combine(grad=[FiniteDiffSamplerGradient, ParamShiftSamplerGradient, LinCombSamplerGradient])
+    @combine(grad=gradient_factories)
     def test_gradient_parameters(self, grad):
         """Test the sampler gradient for parameters"""
         sampler = Sampler()
@@ -307,10 +300,7 @@ class TestSamplerGradient(QiskitTestCase):
         qc.rx(a, 0)
         qc.rz(b, 0)
         qc.measure_all()
-        if grad is FiniteDiffSamplerGradient:
-            gradient = grad(sampler, epsilon=1e-6)
-        else:
-            gradient = grad(sampler)
+        gradient = grad(sampler)
         param_list = [[np.pi / 4, np.pi / 2]]
         correct_results = [
             [{0: -0.5 / np.sqrt(2), 1: 0.5 / np.sqrt(2)}],
@@ -321,7 +311,7 @@ class TestSamplerGradient(QiskitTestCase):
                 for k in quasi_dist:
                     self.assertAlmostEqual(quasi_dist[k], correct_results[i][j][k], 3)
 
-    @combine(grad=[FiniteDiffSamplerGradient, ParamShiftSamplerGradient, LinCombSamplerGradient])
+    @combine(grad=gradient_factories)
     def test_gradient_multi_arguments(self, grad):
         """Test the sampler gradient for multiple arguments"""
         sampler = Sampler()
@@ -333,10 +323,7 @@ class TestSamplerGradient(QiskitTestCase):
         qc2 = QuantumCircuit(1)
         qc2.rx(b, 0)
         qc2.measure_all()
-        if grad is FiniteDiffSamplerGradient:
-            gradient = grad(sampler, epsilon=1e-6)
-        else:
-            gradient = grad(sampler)
+        gradient = grad(sampler)
         param_list = [[np.pi / 4], [np.pi / 2]]
         correct_results = [
             [{0: -0.5 / np.sqrt(2), 1: 0.5 / np.sqrt(2)}],