From 2ee4ffa0a285cba31644767ece671cc4dbfd715e Mon Sep 17 00:00:00 2001
From: a-matsuo <47442626+a-matsuo@users.noreply.github.com>
Date: Wed, 29 Mar 2023 04:50:41 +0900
Subject: [PATCH] Rearrange the gradient result based on the order of the input
 parameters (Qiskit/qiskit-terra#9503)

* finite diff and spsa estimator gradient

* for param shift, lin_comb, and reverse gradients

* reverse estimator gradient

* implemented the new feature

* fix lint

* fix test

* fix and add reno

* fix

* fix style
---
 .../gradients/base_estimator_gradient.py      |  98 +++---
 qiskit_algorithms/gradients/base_qgt.py       |  94 +++---
 .../gradients/base_sampler_gradient.py        |  79 ++---
 .../finite_diff_estimator_gradient.py         |  12 +-
 .../gradients/finite_diff_sampler_gradient.py |  12 +-
 .../gradients/lin_comb_estimator_gradient.py  |  24 +-
 qiskit_algorithms/gradients/lin_comb_qgt.py   |  32 +-
 .../gradients/lin_comb_sampler_gradient.py    |  25 +-
 .../param_shift_estimator_gradient.py         |  22 +-
 .../gradients/param_shift_sampler_gradient.py |  20 +-
 qiskit_algorithms/gradients/qfi.py            |  14 +-
 .../reverse_gradient/reverse_gradient.py      |  22 +-
 .../gradients/spsa_estimator_gradient.py      |  12 +-
 .../gradients/spsa_sampler_gradient.py        |  18 +-
 qiskit_algorithms/gradients/utils.py          |  26 +-
 test/gradients/test_estimator_gradient.py     |  87 +++--
 test/gradients/test_qgt.py                    |  36 +-
 test/gradients/test_sampler_gradient.py       | 319 +++++++++++-------
 .../variational/test_var_qite.py              |  18 +-
 19 files changed, 543 insertions(+), 427 deletions(-)

diff --git a/qiskit_algorithms/gradients/base_estimator_gradient.py b/qiskit_algorithms/gradients/base_estimator_gradient.py
index daf837bd..a1e84b0c 100644
--- a/qiskit_algorithms/gradients/base_estimator_gradient.py
+++ b/qiskit_algorithms/gradients/base_estimator_gradient.py
@@ -35,7 +35,7 @@
     DerivativeType,
     GradientCircuit,
     _assign_unique_parameters,
-    _make_gradient_parameter_set,
+    _make_gradient_parameters,
     _make_gradient_parameter_values,
 )
 
@@ -103,7 +103,8 @@ def run(
             parameters: The sequence of parameters to calculate only the gradients of
                 the specified parameters. Each sequence of parameters corresponds to a circuit in
                 ``circuits``. Defaults to None, which means that the gradients of all parameters in
-                each circuit are calculated.
+                each circuit are calculated. None in the sequence means that the gradients of all
+                parameters in the corresponding circuit are calculated.
             options: Primitive backend runtime options used for circuit execution.
                 The order of priority is: options in ``run`` method > gradient's
                 default options > primitive's default setting.
@@ -127,24 +128,24 @@ def run(
 
         if parameters is None:
             # If parameters is None, we calculate the gradients of all parameters in each circuit.
-            parameter_sets = [set(circuit.parameters) for circuit in circuits]
+            parameters = [circuit.parameters for circuit in circuits]
         else:
             # If parameters is not None, we calculate the gradients of the specified parameters.
             # None in parameters means that the gradients of all parameters in the corresponding
             # circuit are calculated.
-            parameter_sets = [
-                set(parameters_) if parameters_ is not None else set(circuits[i].parameters)
-                for i, parameters_ in enumerate(parameters)
+            parameters = [
+                params if params is not None else circuits[i].parameters
+                for i, params in enumerate(parameters)
             ]
         # Validate the arguments.
-        self._validate_arguments(circuits, observables, parameter_values, parameter_sets)
+        self._validate_arguments(circuits, observables, parameter_values, parameters)
         # The priority of run option is as follows:
         # options in ``run`` method > gradient's default options > primitive's default setting.
         opts = copy(self._default_options)
         opts.update_options(**options)
         # Run the job.
         job = AlgorithmJob(
-            self._run, circuits, observables, parameter_values, parameter_sets, **opts.__dict__
+            self._run, circuits, observables, parameter_values, parameters, **opts.__dict__
         )
         job.submit()
         return job
@@ -155,7 +156,7 @@ def _run(
         circuits: Sequence[QuantumCircuit],
         observables: Sequence[BaseOperator | PauliSumOp],
         parameter_values: Sequence[Sequence[float]],
-        parameter_sets: Sequence[set[Parameter]],
+        parameters: Sequence[Sequence[Parameter]],
         **options,
     ) -> EstimatorGradientResult:
         """Compute the estimator gradients on the given circuits."""
@@ -165,9 +166,9 @@ def _preprocess(
         self,
         circuits: Sequence[QuantumCircuit],
         parameter_values: Sequence[Sequence[float]],
-        parameter_sets: Sequence[set[Parameter]],
+        parameters: Sequence[Sequence[Parameter]],
         supported_gates: Sequence[str],
-    ) -> tuple[Sequence[QuantumCircuit], Sequence[Sequence[float]], Sequence[set[Parameter]]]:
+    ) -> tuple[Sequence[QuantumCircuit], Sequence[Sequence[float]], Sequence[Sequence[Parameter]]]:
         """Preprocess the gradient. This makes a gradient circuit for each circuit. The gradient
         circuit is a transpiled circuit by using the supported gates, and has unique parameters.
         ``parameter_values`` and ``parameters`` are also updated to match the gradient circuit.
@@ -175,7 +176,7 @@ def _preprocess(
         Args:
             circuits: The list of quantum circuits to compute the gradients.
             parameter_values: The list of parameter values to be bound to the circuit.
-            parameter_sets: The sequence of parameters to calculate only the gradients of the specified
+            parameters: The sequence of parameters to calculate only the gradients of the specified
                 parameters.
             supported_gates: The supported gates used to transpile the circuit.
 
@@ -184,10 +185,8 @@ def _preprocess(
             parameter_values and parameters are updated to match the gradient circuit.
         """
         translator = TranslateParameterizedGates(supported_gates)
-        g_circuits, g_parameter_values, g_parameter_sets = [], [], []
-        for circuit, parameter_value_, parameter_set in zip(
-            circuits, parameter_values, parameter_sets
-        ):
+        g_circuits, g_parameter_values, g_parameters = [], [], []
+        for circuit, parameter_value_, parameters_ in zip(circuits, parameter_values, parameters):
             circuit_key = _circuit_key(circuit)
             if circuit_key not in self._gradient_circuit_cache:
                 unrolled = translator(circuit)
@@ -197,15 +196,15 @@ def _preprocess(
             g_parameter_values.append(
                 _make_gradient_parameter_values(circuit, gradient_circuit, parameter_value_)
             )
-            g_parameter_sets.append(_make_gradient_parameter_set(gradient_circuit, parameter_set))
-        return g_circuits, g_parameter_values, g_parameter_sets
+            g_parameters.append(_make_gradient_parameters(gradient_circuit, parameters_))
+        return g_circuits, g_parameter_values, g_parameters
 
     def _postprocess(
         self,
         results: EstimatorGradientResult,
         circuits: Sequence[QuantumCircuit],
         parameter_values: Sequence[Sequence[float]],
-        parameter_sets: Sequence[set[Parameter]],
+        parameters: Sequence[Sequence[Parameter]],
     ) -> EstimatorGradientResult:
         """Postprocess the gradients. This method computes the gradient of the original circuits
         by applying the chain rule to the gradient of the circuits with unique parameters.
@@ -214,17 +213,17 @@ def _postprocess(
             results: The computed gradients for the circuits with unique parameters.
             circuits: The list of original circuits submitted for gradient computation.
             parameter_values: The list of parameter values to be bound to the circuits.
-            parameter_sets: An optional subset of parameters with respect to which the gradients should
-                be calculated.
+            parameters: The sequence of parameters to calculate only the gradients of the specified
+                parameters.
 
         Returns:
             The gradients of the original circuits.
         """
         gradients, metadata = [], []
-        for idx, (circuit, parameter_values_, parameter_set) in enumerate(
-            zip(circuits, parameter_values, parameter_sets)
+        for idx, (circuit, parameter_values_, parameters_) in enumerate(
+            zip(circuits, parameter_values, parameters)
         ):
-            gradient = np.zeros(len(parameter_set))
+            gradient = np.zeros(len(parameters_))
             if (
                 "derivative_type" in results.metadata[idx]
                 and results.metadata[idx]["derivative_type"] == DerivativeType.COMPLEX
@@ -232,18 +231,12 @@ def _postprocess(
                 # If the derivative type is complex, cast the gradient to complex.
                 gradient = gradient.astype("complex")
             gradient_circuit = self._gradient_circuit_cache[_circuit_key(circuit)]
-            g_parameter_set = _make_gradient_parameter_set(gradient_circuit, parameter_set)
+            g_parameters = _make_gradient_parameters(gradient_circuit, parameters_)
             # Make a map from the gradient parameter to the respective index in the gradient.
-            parameter_indices = [param for param in circuit.parameters if param in parameter_set]
-            g_parameter_indices = [
-                param
-                for param in gradient_circuit.gradient_circuit.parameters
-                if param in g_parameter_set
-            ]
-            g_parameter_indices = {param: i for i, param in enumerate(g_parameter_indices)}
+            g_parameter_indices = {param: i for i, param in enumerate(g_parameters)}
             # Compute the original gradient from the gradient of the gradient circuit
             # by using the chain rule.
-            for i, parameter in enumerate(parameter_indices):
+            for i, parameter in enumerate(parameters_):
                 for g_parameter, coeff in gradient_circuit.parameter_map[parameter]:
                     # Compute the coefficient
                     if isinstance(coeff, ParameterExpression):
@@ -261,7 +254,7 @@ def _postprocess(
                         * results.gradients[idx][g_parameter_indices[g_parameter]]
                     )
             gradients.append(gradient)
-            metadata.append([{"parameters": parameter_indices}])
+            metadata.append({"parameters": parameters_})
         return EstimatorGradientResult(
             gradients=gradients, metadata=metadata, options=results.options
         )
@@ -271,7 +264,7 @@ def _validate_arguments(
         circuits: Sequence[QuantumCircuit],
         observables: Sequence[BaseOperator | PauliSumOp],
         parameter_values: Sequence[Sequence[float]],
-        parameter_sets: Sequence[set[Parameter]],
+        parameters: Sequence[Sequence[Parameter]],
     ) -> None:
         """Validate the arguments of the ``run`` method.
 
@@ -279,9 +272,8 @@ def _validate_arguments(
             circuits: The list of quantum circuits to compute the gradients.
             observables: The list of observables.
             parameter_values: The list of parameter values to be bound to the circuit.
-            parameter_sets: The Sequence of parameter sets to calculate only the gradients of
-                the specified parameters. Each set of parameters corresponds to a circuit in
-                ``circuits``.
+            parameters: The sequence of parameters to calculate only the gradients of the specified
+                parameters.
 
         Raises:
             ValueError: Invalid arguments are given.
@@ -292,18 +284,6 @@ def _validate_arguments(
                 f"the number of parameter value sets ({len(parameter_values)})."
             )
 
-        if len(circuits) != len(observables):
-            raise ValueError(
-                f"The number of circuits ({len(circuits)}) does not match "
-                f"the number of observables ({len(observables)})."
-            )
-
-        if len(circuits) != len(parameter_sets):
-            raise ValueError(
-                f"The number of circuits ({len(circuits)}) does not match "
-                f"the number of the specified parameter sets ({len(parameter_sets)})."
-            )
-
         for i, (circuit, parameter_value) in enumerate(zip(circuits, parameter_values)):
             if not circuit.num_parameters:
                 raise ValueError(f"The {i}-th circuit is not parameterised.")
@@ -313,6 +293,12 @@ def _validate_arguments(
                     f"the number of parameters ({circuit.num_parameters}) for the {i}-th circuit."
                 )
 
+        if len(circuits) != len(observables):
+            raise ValueError(
+                f"The number of circuits ({len(circuits)}) does not match "
+                f"the number of observables ({len(observables)})."
+            )
+
         for i, (circuit, observable) in enumerate(zip(circuits, observables)):
             if circuit.num_qubits != observable.num_qubits:
                 raise ValueError(
@@ -321,10 +307,16 @@ def _validate_arguments(
                     f"({observable.num_qubits})."
                 )
 
-        for i, (circuit, parameter_set) in enumerate(zip(circuits, parameter_sets)):
-            if not set(parameter_set).issubset(circuit.parameters):
+        if len(circuits) != len(parameters):
+            raise ValueError(
+                f"The number of circuits ({len(circuits)}) does not match "
+                f"the number of the list of specified parameters ({len(parameters)})."
+            )
+
+        for i, (circuit, parameters_) in enumerate(zip(circuits, parameters)):
+            if not set(parameters_).issubset(circuit.parameters):
                 raise ValueError(
-                    f"The {i}-th parameter set contains parameters not present in the "
+                    f"The {i}-th parameters contains parameters not present in the "
                     f"{i}-th circuit."
                 )
 
diff --git a/qiskit_algorithms/gradients/base_qgt.py b/qiskit_algorithms/gradients/base_qgt.py
index 14bcff1b..2a196cf4 100644
--- a/qiskit_algorithms/gradients/base_qgt.py
+++ b/qiskit_algorithms/gradients/base_qgt.py
@@ -33,7 +33,7 @@
     DerivativeType,
     GradientCircuit,
     _assign_unique_parameters,
-    _make_gradient_parameter_set,
+    _make_gradient_parameters,
     _make_gradient_parameter_values,
 )
 
@@ -145,22 +145,22 @@ def run(
 
         if parameters is None:
             # If parameters is None, we calculate the gradients of all parameters in each circuit.
-            parameter_sets = [set(circuit.parameters) for circuit in circuits]
+            parameters = [circuit.parameters for circuit in circuits]
         else:
             # If parameters is not None, we calculate the gradients of the specified parameters.
             # None in parameters means that the gradients of all parameters in the corresponding
             # circuit are calculated.
-            parameter_sets = [
-                set(parameters_) if parameters_ is not None else set(circuits[i].parameters)
-                for i, parameters_ in enumerate(parameters)
+            parameters = [
+                params if params is not None else circuits[i].parameters
+                for i, params in enumerate(parameters)
             ]
         # Validate the arguments.
-        self._validate_arguments(circuits, parameter_values, parameter_sets)
+        self._validate_arguments(circuits, parameter_values, parameters)
         # The priority of run option is as follows:
         # options in ``run`` method > QGT's default options > primitive's default setting.
         opts = copy(self._default_options)
         opts.update_options(**options)
-        job = AlgorithmJob(self._run, circuits, parameter_values, parameter_sets, **opts.__dict__)
+        job = AlgorithmJob(self._run, circuits, parameter_values, parameters, **opts.__dict__)
         job.submit()
         return job
 
@@ -169,7 +169,7 @@ def _run(
         self,
         circuits: Sequence[QuantumCircuit],
         parameter_values: Sequence[Sequence[float]],
-        parameter_sets: Sequence[Sequence[Parameter]],
+        parameters: Sequence[Sequence[Parameter]],
         **options,
     ) -> QGTResult:
         """Compute the QGTs on the given circuits."""
@@ -179,9 +179,9 @@ def _preprocess(
         self,
         circuits: Sequence[QuantumCircuit],
         parameter_values: Sequence[Sequence[float]],
-        parameter_sets: Sequence[set[Parameter]],
+        parameters: Sequence[Sequence[Parameter]],
         supported_gates: Sequence[str],
-    ) -> tuple[Sequence[QuantumCircuit], Sequence[Sequence[float]], Sequence[set[Parameter]]]:
+    ) -> tuple[Sequence[QuantumCircuit], Sequence[Sequence[float]], Sequence[Sequence[Parameter]]]:
         """Preprocess the gradient. This makes a gradient circuit for each circuit. The gradient
         circuit is a transpiled circuit by using the supported gates, and has unique parameters.
         ``parameter_values`` and ``parameters`` are also updated to match the gradient circuit.
@@ -189,7 +189,7 @@ def _preprocess(
         Args:
             circuits: The list of quantum circuits to compute the gradients.
             parameter_values: The list of parameter values to be bound to the circuit.
-            parameter_sets: The sequence of parameters to calculate only the gradients of the specified
+            parameters: The sequence of parameters to calculate only the gradients of the specified
                 parameters.
             supported_gates: The supported gates used to transpile the circuit.
 
@@ -198,10 +198,8 @@ def _preprocess(
             parameter_values and parameters are updated to match the gradient circuit.
         """
         translator = TranslateParameterizedGates(supported_gates)
-        g_circuits, g_parameter_values, g_parameter_sets = [], [], []
-        for circuit, parameter_value_, parameter_set in zip(
-            circuits, parameter_values, parameter_sets
-        ):
+        g_circuits, g_parameter_values, g_parameters = [], [], []
+        for circuit, parameter_value_, parameters_ in zip(circuits, parameter_values, parameters):
             circuit_key = _circuit_key(circuit)
             if circuit_key not in self._gradient_circuit_cache:
                 unrolled = translator(circuit)
@@ -211,15 +209,20 @@ def _preprocess(
             g_parameter_values.append(
                 _make_gradient_parameter_values(circuit, gradient_circuit, parameter_value_)
             )
-            g_parameter_sets.append(_make_gradient_parameter_set(gradient_circuit, parameter_set))
-        return g_circuits, g_parameter_values, g_parameter_sets
+            g_parameters_ = [
+                g_param
+                for g_param in gradient_circuit.gradient_circuit.parameters
+                if g_param in _make_gradient_parameters(gradient_circuit, parameters_)
+            ]
+            g_parameters.append(g_parameters_)
+        return g_circuits, g_parameter_values, g_parameters
 
     def _postprocess(
         self,
         results: QGTResult,
         circuits: Sequence[QuantumCircuit],
         parameter_values: Sequence[Sequence[float]],
-        parameter_sets: Sequence[set[Parameter]],
+        parameters: Sequence[Sequence[Parameter]],
     ) -> QGTResult:
         """Postprocess the QGTs. This method computes the QGTs of the original circuits
         by applying the chain rule to the QGTs of the circuits with unique parameters.
@@ -228,36 +231,33 @@ def _postprocess(
             results: The computed QGT for the circuits with unique parameters.
             circuits: The list of original circuits submitted for gradient computation.
             parameter_values: The list of parameter values to be bound to the circuits.
-            parameter_sets: An optional subset of parameters with respect to which the QGTs should
-                be calculated.
+            parameters: The sequence of parameters to calculate only the gradients of the specified
+                parameters.
 
         Returns:
             The QGTs of the original circuits.
         """
         qgts, metadata = [], []
-        for idx, (circuit, parameter_values_, parameter_set) in enumerate(
-            zip(circuits, parameter_values, parameter_sets)
+        for idx, (circuit, parameter_values_, parameters_) in enumerate(
+            zip(circuits, parameter_values, parameters)
         ):
             dtype = complex if self.derivative_type == DerivativeType.COMPLEX else float
-            qgt = np.zeros((len(parameter_set), len(parameter_set)), dtype=dtype)
+            qgt = np.zeros((len(parameters_), len(parameters_)), dtype=dtype)
 
             gradient_circuit = self._gradient_circuit_cache[_circuit_key(circuit)]
-            g_parameter_set = _make_gradient_parameter_set(gradient_circuit, parameter_set)
+            g_parameters = _make_gradient_parameters(gradient_circuit, parameters_)
             # Make a map from the gradient parameter to the respective index in the gradient.
-            parameter_indices = [param for param in circuit.parameters if param in parameter_set]
+            # parameters_ = [param for param in circuit.parameters if param in parameters_]
             g_parameter_indices = [
                 param
                 for param in gradient_circuit.gradient_circuit.parameters
-                if param in g_parameter_set
+                if param in g_parameters
             ]
             g_parameter_indices = {param: i for i, param in enumerate(g_parameter_indices)}
-
-            rows, cols = np.triu_indices(len(parameter_indices))
+            rows, cols = np.triu_indices(len(parameters_))
             for row, col in zip(rows, cols):
-                for g_parameter1, coeff1 in gradient_circuit.parameter_map[parameter_indices[row]]:
-                    for g_parameter2, coeff2 in gradient_circuit.parameter_map[
-                        parameter_indices[col]
-                    ]:
+                for g_parameter1, coeff1 in gradient_circuit.parameter_map[parameters_[row]]:
+                    for g_parameter2, coeff2 in gradient_circuit.parameter_map[parameters_[col]]:
                         if isinstance(coeff1, ParameterExpression):
                             local_map = {
                                 p: parameter_values_[circuit.parameters.data.index(p)]
@@ -281,12 +281,13 @@ def _postprocess(
                                 g_parameter_indices[g_parameter1], g_parameter_indices[g_parameter2]
                             ]
                         )
+
             if self.derivative_type == DerivativeType.IMAG:
                 qgt += -1 * np.triu(qgt, k=1).T
             else:
                 qgt += np.triu(qgt, k=1).conjugate().T
             qgts.append(qgt)
-            metadata.append([{"parameters": parameter_indices}])
+            metadata.append([{"parameters": parameters_}])
         return QGTResult(
             qgts=qgts,
             derivative_type=self.derivative_type,
@@ -294,19 +295,19 @@ def _postprocess(
             options=results.options,
         )
 
+    @staticmethod
     def _validate_arguments(
-        self,
         circuits: Sequence[QuantumCircuit],
         parameter_values: Sequence[Sequence[float]],
-        parameter_sets: Sequence[set[Parameter]],
+        parameters: Sequence[Sequence[Parameter]],
     ) -> None:
         """Validate the arguments of the ``run`` method.
 
         Args:
             circuits: The list of quantum circuits to compute the QGTs.
             parameter_values: The list of parameter values to be bound to the circuits.
-            parameter_sets: The sequence of parameter sets with respect to which the QGTs should be
-                computed. Each set of parameters corresponds to a circuit in ``circuits``.
+            parameters: The sequence of parameters with respect to which the QGTs should be
+                computed.
 
         Raises:
             ValueError: Invalid arguments are given.
@@ -314,13 +315,13 @@ def _validate_arguments(
         if len(circuits) != len(parameter_values):
             raise ValueError(
                 f"The number of circuits ({len(circuits)}) does not match "
-                f"the number of parameter value sets ({len(parameter_values)})."
+                f"the number of parameter values ({len(parameter_values)})."
             )
 
-        if len(circuits) != len(parameter_sets):
+        if len(circuits) != len(parameters):
             raise ValueError(
                 f"The number of circuits ({len(circuits)}) does not match "
-                f"the number of the specified parameter sets ({len(parameter_sets)})."
+                f"the number of the specified parameter sets ({len(parameters)})."
             )
 
         for i, (circuit, parameter_value) in enumerate(zip(circuits, parameter_values)):
@@ -332,6 +333,19 @@ def _validate_arguments(
                     f"the number of parameters ({circuit.num_parameters}) for the {i}-th circuit."
                 )
 
+        if len(circuits) != len(parameters):
+            raise ValueError(
+                f"The number of circuits ({len(circuits)}) does not match "
+                f"the number of the list of specified parameters ({len(parameters)})."
+            )
+
+        for i, (circuit, parameters_) in enumerate(zip(circuits, parameters)):
+            if not set(parameters_).issubset(circuit.parameters):
+                raise ValueError(
+                    f"The {i}-th parameters contains parameters not present in the "
+                    f"{i}-th circuit."
+                )
+
     @property
     def options(self) -> Options:
         """Return the union of estimator options setting and QGT default options,
diff --git a/qiskit_algorithms/gradients/base_sampler_gradient.py b/qiskit_algorithms/gradients/base_sampler_gradient.py
index edfd1cd3..b4983fd3 100644
--- a/qiskit_algorithms/gradients/base_sampler_gradient.py
+++ b/qiskit_algorithms/gradients/base_sampler_gradient.py
@@ -31,7 +31,7 @@
 from .utils import (
     GradientCircuit,
     _assign_unique_parameters,
-    _make_gradient_parameter_set,
+    _make_gradient_parameters,
     _make_gradient_parameter_values,
 )
 
@@ -71,7 +71,8 @@ def run(
             parameters: The sequence of parameters to calculate only the gradients of
                 the specified parameters. Each sequence of parameters corresponds to a circuit in
                 ``circuits``. Defaults to None, which means that the gradients of all parameters in
-                each circuit are calculated.
+                each circuit are calculated. None in the sequence means that the gradients of all
+                parameters in the corresponding circuit are calculated.
             options: Primitive backend runtime options used for circuit execution.
                 The order of priority is: options in ``run`` method > gradient's
                 default options > primitive's default setting.
@@ -90,22 +91,22 @@ def run(
             circuits = (circuits,)
         if parameters is None:
             # If parameters is None, we calculate the gradients of all parameters in each circuit.
-            parameter_sets = [set(circuit.parameters) for circuit in circuits]
+            parameters = [circuit.parameters for circuit in circuits]
         else:
             # If parameters is not None, we calculate the gradients of the specified parameters.
             # None in parameters means that the gradients of all parameters in the corresponding
             # circuit are calculated.
-            parameter_sets = [
-                set(parameters_) if parameters_ is not None else set(circuits[i].parameters)
-                for i, parameters_ in enumerate(parameters)
+            parameters = [
+                params if params is not None else circuits[i].parameters
+                for i, params in enumerate(parameters)
             ]
         # Validate the arguments.
-        self._validate_arguments(circuits, parameter_values, parameter_sets)
+        self._validate_arguments(circuits, parameter_values, parameters)
         # The priority of run option is as follows:
         # options in `run` method > gradient's default options > primitive's default options.
         opts = copy(self._default_options)
         opts.update_options(**options)
-        job = AlgorithmJob(self._run, circuits, parameter_values, parameter_sets, **opts.__dict__)
+        job = AlgorithmJob(self._run, circuits, parameter_values, parameters, **opts.__dict__)
         job.submit()
         return job
 
@@ -114,7 +115,7 @@ def _run(
         self,
         circuits: Sequence[QuantumCircuit],
         parameter_values: Sequence[Sequence[float]],
-        parameter_sets: Sequence[set[Parameter]],
+        parameters: Sequence[Sequence[Parameter]],
         **options,
     ) -> SamplerGradientResult:
         """Compute the sampler gradients on the given circuits."""
@@ -124,7 +125,7 @@ def _preprocess(
         self,
         circuits: Sequence[QuantumCircuit],
         parameter_values: Sequence[Sequence[float]],
-        parameter_sets: Sequence[set[Parameter]],
+        parameters: Sequence[Sequence[Parameter]],
         supported_gates: Sequence[str],
     ) -> tuple[Sequence[QuantumCircuit], Sequence[Sequence[float]], Sequence[set[Parameter]]]:
         """Preprocess the gradient. This makes a gradient circuit for each circuit. The gradient
@@ -134,7 +135,7 @@ def _preprocess(
         Args:
             circuits: The list of quantum circuits to compute the gradients.
             parameter_values: The list of parameter values to be bound to the circuit.
-            parameter_sets: The sequence of parameters to calculate only the gradients of the specified
+            parameters: The sequence of parameters to calculate only the gradients of the specified
                 parameters.
             supported_gates: The supported gates used to transpile the circuit.
 
@@ -143,10 +144,8 @@ def _preprocess(
             parameter_values and parameters are updated to match the gradient circuit.
         """
         translator = TranslateParameterizedGates(supported_gates)
-        g_circuits, g_parameter_values, g_parameter_sets = [], [], []
-        for circuit, parameter_value_, parameter_set in zip(
-            circuits, parameter_values, parameter_sets
-        ):
+        g_circuits, g_parameter_values, g_parameters = [], [], []
+        for circuit, parameter_value_, parameters_ in zip(circuits, parameter_values, parameters):
             circuit_key = _circuit_key(circuit)
             if circuit_key not in self._gradient_circuit_cache:
                 unrolled = translator(circuit)
@@ -156,15 +155,15 @@ def _preprocess(
             g_parameter_values.append(
                 _make_gradient_parameter_values(circuit, gradient_circuit, parameter_value_)
             )
-            g_parameter_sets.append(_make_gradient_parameter_set(gradient_circuit, parameter_set))
-        return g_circuits, g_parameter_values, g_parameter_sets
+            g_parameters.append(_make_gradient_parameters(gradient_circuit, parameters_))
+        return g_circuits, g_parameter_values, g_parameters
 
     def _postprocess(
         self,
         results: SamplerGradientResult,
         circuits: Sequence[QuantumCircuit],
         parameter_values: Sequence[Sequence[float]],
-        parameter_sets: Sequence[set[Parameter] | None],
+        parameters: Sequence[Sequence[Parameter] | None],
     ) -> SamplerGradientResult:
         """Postprocess the gradient. This computes the gradient of the original circuit from the
         gradient of the gradient circuit by using the chain rule.
@@ -173,30 +172,24 @@ def _postprocess(
             results: The results of the gradient of the gradient circuits.
             circuits: The list of quantum circuits to compute the gradients.
             parameter_values: The list of parameter values to be bound to the circuit.
-            parameter_sets: The sequence of parameters to calculate only the gradients of the specified
+            parameters: The sequence of parameters to calculate only the gradients of the specified
                 parameters.
 
         Returns:
             The results of the gradient of the original circuits.
         """
         gradients, metadata = [], []
-        for idx, (circuit, parameter_values_, parameter_set) in enumerate(
-            zip(circuits, parameter_values, parameter_sets)
+        for idx, (circuit, parameter_values_, parameters_) in enumerate(
+            zip(circuits, parameter_values, parameters)
         ):
             gradient_circuit = self._gradient_circuit_cache[_circuit_key(circuit)]
-            g_parameter_set = _make_gradient_parameter_set(gradient_circuit, parameter_set)
+            g_parameters = _make_gradient_parameters(gradient_circuit, parameters_)
             # Make a map from the gradient parameter to the respective index in the gradient.
-            parameter_indices = [param for param in circuit.parameters if param in parameter_set]
-            g_parameter_indices = [
-                param
-                for param in gradient_circuit.gradient_circuit.parameters
-                if param in g_parameter_set
-            ]
-            g_parameter_indices = {param: i for i, param in enumerate(g_parameter_indices)}
+            g_parameter_indices = {param: i for i, param in enumerate(g_parameters)}
             # Compute the original gradient from the gradient of the gradient circuit
             # by using the chain rule.
             gradient = []
-            for parameter in parameter_indices:
+            for parameter in parameters_:
                 grad_dist = defaultdict(float)
                 for g_parameter, coeff in gradient_circuit.parameter_map[parameter]:
                     # Compute the coefficient
@@ -215,7 +208,7 @@ def _postprocess(
                         grad_dist[key] += float(bound_coeff) * value
                 gradient.append(dict(grad_dist))
             gradients.append(gradient)
-            metadata.append([{"parameters": parameter_indices}])
+            metadata.append([{"parameters": parameters_}])
         return SamplerGradientResult(
             gradients=gradients, metadata=metadata, options=results.options
         )
@@ -224,17 +217,15 @@ def _postprocess(
     def _validate_arguments(
         circuits: Sequence[QuantumCircuit],
         parameter_values: Sequence[Sequence[float]],
-        parameter_sets: Sequence[set[Parameter]],
+        parameters: Sequence[Sequence[Parameter]],
     ) -> None:
         """Validate the arguments of the ``run`` method.
 
         Args:
             circuits: The list of quantum circuits to compute the gradients.
             parameter_values: The list of parameter values to be bound to the circuit.
-            parameter_sets: The Sequence of Sequence of Parameters to calculate only the gradients of
-                the specified parameters. Each Sequence of Parameters corresponds to a circuit in
-                ``circuits``. Defaults to None, which means that the gradients of all parameters in
-                each circuit are calculated.
+            parameters: The sequence of parameters to calculate only the gradients of the specified
+                parameters.
 
         Raises:
             ValueError: Invalid arguments are given.
@@ -245,12 +236,6 @@ def _validate_arguments(
                 f"the number of parameter value sets ({len(parameter_values)})."
             )
 
-        if len(circuits) != len(parameter_sets):
-            raise ValueError(
-                f"The number of circuits ({len(circuits)}) does not match "
-                f"the number of the specified parameter sets ({len(parameter_sets)})."
-            )
-
         for i, (circuit, parameter_value) in enumerate(zip(circuits, parameter_values)):
             if not circuit.num_parameters:
                 raise ValueError(f"The {i}-th circuit is not parameterised.")
@@ -261,8 +246,14 @@ def _validate_arguments(
                     f"the number of parameters ({circuit.num_parameters}) for the {i}-th circuit."
                 )
 
-        for i, (circuit, parameter_set) in enumerate(zip(circuits, parameter_sets)):
-            if not set(parameter_set).issubset(circuit.parameters):
+        if len(circuits) != len(parameters):
+            raise ValueError(
+                f"The number of circuits ({len(circuits)}) does not match "
+                f"the number of the specified parameter sets ({len(parameters)})."
+            )
+
+        for i, (circuit, parameters_) in enumerate(zip(circuits, parameters)):
+            if not set(parameters_).issubset(circuit.parameters):
                 raise ValueError(
                     f"The {i}-th parameter set contains parameters not present in the "
                     f"{i}-th circuit."
diff --git a/qiskit_algorithms/gradients/finite_diff_estimator_gradient.py b/qiskit_algorithms/gradients/finite_diff_estimator_gradient.py
index 40ba2f65..9d1d751d 100644
--- a/qiskit_algorithms/gradients/finite_diff_estimator_gradient.py
+++ b/qiskit_algorithms/gradients/finite_diff_estimator_gradient.py
@@ -88,21 +88,19 @@ def _run(
         circuits: Sequence[QuantumCircuit],
         observables: Sequence[BaseOperator | PauliSumOp],
         parameter_values: Sequence[Sequence[float]],
-        parameter_sets: Sequence[set[Parameter]],
+        parameters: Sequence[Sequence[Parameter]],
         **options,
     ) -> EstimatorGradientResult:
         """Compute the estimator gradients on the given circuits."""
         job_circuits, job_observables, job_param_values, metadata = [], [], [], []
         all_n = []
 
-        for circuit, observable, parameter_values_, parameter_set in zip(
-            circuits, observables, parameter_values, parameter_sets
+        for circuit, observable, parameter_values_, parameters_ in zip(
+            circuits, observables, parameter_values, parameters
         ):
             # Indices of parameters to be differentiated
-            indices = [
-                circuit.parameters.data.index(p) for p in circuit.parameters if p in parameter_set
-            ]
-            metadata.append({"parameters": [circuit.parameters[idx] for idx in indices]})
+            indices = [circuit.parameters.data.index(p) for p in parameters_]
+            metadata.append({"parameters": parameters_})
 
             # Combine inputs into a single job to reduce overhead.
             offset = np.identity(circuit.num_parameters)[indices, :]
diff --git a/qiskit_algorithms/gradients/finite_diff_sampler_gradient.py b/qiskit_algorithms/gradients/finite_diff_sampler_gradient.py
index 0bec48fb..c9331230 100644
--- a/qiskit_algorithms/gradients/finite_diff_sampler_gradient.py
+++ b/qiskit_algorithms/gradients/finite_diff_sampler_gradient.py
@@ -86,20 +86,16 @@ def _run(
         self,
         circuits: Sequence[QuantumCircuit],
         parameter_values: Sequence[Sequence[float]],
-        parameter_sets: Sequence[set[Parameter]],
+        parameters: Sequence[Sequence[Parameter]],
         **options,
     ) -> SamplerGradientResult:
         """Compute the sampler gradients on the given circuits."""
         job_circuits, job_param_values, metadata = [], [], []
         all_n = []
-        for circuit, parameter_values_, parameter_set in zip(
-            circuits, parameter_values, parameter_sets
-        ):
+        for circuit, parameter_values_, parameters_ in zip(circuits, parameter_values, parameters):
             # Indices of parameters to be differentiated
-            indices = [
-                circuit.parameters.data.index(p) for p in circuit.parameters if p in parameter_set
-            ]
-            metadata.append({"parameters": [circuit.parameters[idx] for idx in indices]})
+            indices = [circuit.parameters.data.index(p) for p in parameters_]
+            metadata.append({"parameters": parameters_})
             # Combine inputs into a single job to reduce overhead.
             offset = np.identity(circuit.num_parameters)[indices, :]
             if self._method == "central":
diff --git a/qiskit_algorithms/gradients/lin_comb_estimator_gradient.py b/qiskit_algorithms/gradients/lin_comb_estimator_gradient.py
index f7c51af8..e558b9c0 100644
--- a/qiskit_algorithms/gradients/lin_comb_estimator_gradient.py
+++ b/qiskit_algorithms/gradients/lin_comb_estimator_gradient.py
@@ -100,44 +100,44 @@ def _run(
         circuits: Sequence[QuantumCircuit],
         observables: Sequence[BaseOperator | PauliSumOp],
         parameter_values: Sequence[Sequence[float]],
-        parameter_sets: Sequence[set[Parameter]],
+        parameters: Sequence[Sequence[Parameter]],
         **options,
     ) -> EstimatorGradientResult:
         """Compute the estimator gradients on the given circuits."""
-        g_circuits, g_parameter_values, g_parameter_sets = self._preprocess(
-            circuits, parameter_values, parameter_sets, self.SUPPORTED_GATES
+        g_circuits, g_parameter_values, g_parameters = self._preprocess(
+            circuits, parameter_values, parameters, self.SUPPORTED_GATES
         )
         results = self._run_unique(
-            g_circuits, observables, g_parameter_values, g_parameter_sets, **options
+            g_circuits, observables, g_parameter_values, g_parameters, **options
         )
-        return self._postprocess(results, circuits, parameter_values, parameter_sets)
+        return self._postprocess(results, circuits, parameter_values, parameters)
 
     def _run_unique(
         self,
         circuits: Sequence[QuantumCircuit],
         observables: Sequence[BaseOperator | PauliSumOp],
         parameter_values: Sequence[Sequence[float]],
-        parameter_sets: Sequence[set[Parameter]],
+        parameters: Sequence[Sequence[Parameter]],
         **options,
     ) -> EstimatorGradientResult:
         """Compute the estimator gradients on the given circuits."""
         job_circuits, job_observables, job_param_values, metadata = [], [], [], []
         all_n = []
-        for circuit, observable, parameter_values_, parameter_set in zip(
-            circuits, observables, parameter_values, parameter_sets
+        for circuit, observable, parameter_values_, parameters_ in zip(
+            circuits, observables, parameter_values, parameters
         ):
             # Prepare circuits for the gradient of the specified parameters.
-            meta = {"parameters": [p for p in circuit.parameters if p in parameter_set]}
+            meta = {"parameters": parameters_}
             circuit_key = _circuit_key(circuit)
             if circuit_key not in self._lin_comb_cache:
+                # Cache the circuits for the linear combination of unitaries.
+                # We only cache the circuits for the specified parameters in the future.
                 self._lin_comb_cache[circuit_key] = _make_lin_comb_gradient_circuit(
                     circuit, add_measurement=False
                 )
             lin_comb_circuits = self._lin_comb_cache[circuit_key]
             gradient_circuits = []
-            for param in circuit.parameters:
-                if param not in parameter_set:
-                    continue
+            for param in parameters_:
                 gradient_circuits.append(lin_comb_circuits[param])
             n = len(gradient_circuits)
             # Make the observable as :class:`~qiskit.quantum_info.SparsePauliOp` and
diff --git a/qiskit_algorithms/gradients/lin_comb_qgt.py b/qiskit_algorithms/gradients/lin_comb_qgt.py
index 2efbc4a2..dddd0860 100644
--- a/qiskit_algorithms/gradients/lin_comb_qgt.py
+++ b/qiskit_algorithms/gradients/lin_comb_qgt.py
@@ -119,32 +119,31 @@ def _run(
         self,
         circuits: Sequence[QuantumCircuit],
         parameter_values: Sequence[Sequence[float]],
-        parameter_sets: Sequence[set[Parameter]],
+        parameters: Sequence[Sequence[Parameter]],
         **options,
     ) -> QGTResult:
         """Compute the QGT on the given circuits."""
-        g_circuits, g_parameter_values, g_parameter_sets = self._preprocess(
-            circuits, parameter_values, parameter_sets, self.SUPPORTED_GATES
+        g_circuits, g_parameter_values, g_parameters = self._preprocess(
+            circuits, parameter_values, parameters, self.SUPPORTED_GATES
         )
-        results = self._run_unique(g_circuits, g_parameter_values, g_parameter_sets, **options)
-        return self._postprocess(results, circuits, parameter_values, parameter_sets)
+        results = self._run_unique(g_circuits, g_parameter_values, g_parameters, **options)
+        return self._postprocess(results, circuits, parameter_values, parameters)
 
     def _run_unique(
         self,
         circuits: Sequence[QuantumCircuit],
         parameter_values: Sequence[Sequence[float]],
-        parameter_sets: Sequence[set[Parameter]],
+        parameters: Sequence[Sequence[Parameter]],
         **options,
     ) -> QGTResult:
         """Compute the QGTs on the given circuits."""
         job_circuits, job_observables, job_param_values, metadata = [], [], [], []
         all_n, all_m, phase_fixes = [], [], []
 
-        for circuit, parameter_values_, parameter_set in zip(
-            circuits, parameter_values, parameter_sets
-        ):
+        for circuit, parameter_values_, parameters_ in zip(circuits, parameter_values, parameters):
             # Prepare circuits for the gradient of the specified parameters.
-            meta = {"parameters": [p for p in circuit.parameters if p in parameter_set]}
+            parameters_ = [p for p in circuit.parameters if p in parameters_]
+            meta = {"parameters": parameters_}
             metadata.append(meta)
 
             # Compute the first term in the QGT
@@ -157,11 +156,10 @@ def _run_unique(
             lin_comb_qgt_circuits = self._lin_comb_qgt_circuit_cache[circuit_key]
 
             qgt_circuits = []
-            parameters = [p for p in circuit.parameters if p in parameter_set]
-            rows, cols = np.triu_indices(len(parameters))
+            rows, cols = np.triu_indices(len(parameters_))
             for row, col in zip(rows, cols):
-                param_i = parameters[row]
-                param_j = parameters[col]
+                param_i = parameters_[row]
+                param_j = parameters_[col]
                 qgt_circuits.append(lin_comb_qgt_circuits[(param_i, param_j)])
 
             observable = SparsePauliOp.from_list([("I" * circuit.num_qubits, 1)])
@@ -174,13 +172,13 @@ def _run_unique(
                 job_circuits.extend(qgt_circuits * 2)
                 job_observables.extend([observable_1] * n + [observable_2] * n)
                 job_param_values.extend([parameter_values_] * 2 * n)
-                all_m.append(len(parameter_set))
+                all_m.append(len(parameters_))
                 all_n.append(2 * n)
             else:
                 job_circuits.extend(qgt_circuits)
                 job_observables.extend([observable_1] * n)
                 job_param_values.extend([parameter_values_] * n)
-                all_m.append(len(parameter_set))
+                all_m.append(len(parameters_))
                 all_n.append(n)
 
         # Run the single job with all circuits.
@@ -200,7 +198,7 @@ def _run_unique(
                 circuits=circuits,
                 observables=phase_fix_obs,
                 parameter_values=parameter_values,
-                parameters=parameter_sets,
+                parameters=parameters,
                 **options,
             )
 
diff --git a/qiskit_algorithms/gradients/lin_comb_sampler_gradient.py b/qiskit_algorithms/gradients/lin_comb_sampler_gradient.py
index 6e1e45b8..45cccb67 100644
--- a/qiskit_algorithms/gradients/lin_comb_sampler_gradient.py
+++ b/qiskit_algorithms/gradients/lin_comb_sampler_gradient.py
@@ -78,42 +78,39 @@ def _run(
         self,
         circuits: Sequence[QuantumCircuit],
         parameter_values: Sequence[Sequence[float]],
-        parameter_sets: Sequence[set[Parameter]],
+        parameters: Sequence[Sequence[Parameter]],
         **options,
     ) -> SamplerGradientResult:
         """Compute the estimator gradients on the given circuits."""
-        g_circuits, g_parameter_values, g_parameter_sets = self._preprocess(
-            circuits, parameter_values, parameter_sets, self.SUPPORTED_GATES
+        g_circuits, g_parameter_values, g_parameters = self._preprocess(
+            circuits, parameter_values, parameters, self.SUPPORTED_GATES
         )
-        results = self._run_unique(g_circuits, g_parameter_values, g_parameter_sets, **options)
-        return self._postprocess(results, circuits, parameter_values, parameter_sets)
+        results = self._run_unique(g_circuits, g_parameter_values, g_parameters, **options)
+        return self._postprocess(results, circuits, parameter_values, parameters)
 
     def _run_unique(
         self,
         circuits: Sequence[QuantumCircuit],
         parameter_values: Sequence[Sequence[float]],
-        parameter_sets: Sequence[set[Parameter] | None],
+        parameters: Sequence[Sequence[Parameter]],
         **options,
     ) -> SamplerGradientResult:
         """Compute the sampler gradients on the given circuits."""
         job_circuits, job_param_values, metadata = [], [], []
         all_n = []
-        for circuit, parameter_values_, parameter_set in zip(
-            circuits, parameter_values, parameter_sets
-        ):
+        for circuit, parameter_values_, parameters_ in zip(circuits, parameter_values, parameters):
             # Prepare circuits for the gradient of the specified parameters.
-            meta = {"parameters": [p for p in circuit.parameters if p in parameter_set]}
-            metadata.append(meta)
+            metadata.append({"parameters": parameters_})
             circuit_key = _circuit_key(circuit)
             if circuit_key not in self._lin_comb_cache:
+                # Cache the circuits for the linear combination of unitaries.
+                # We only cache the circuits for the specified parameters in the future.
                 self._lin_comb_cache[circuit_key] = _make_lin_comb_gradient_circuit(
                     circuit, add_measurement=True
                 )
             lin_comb_circuits = self._lin_comb_cache[circuit_key]
             gradient_circuits = []
-            for param in circuit.parameters:
-                if param not in parameter_set:
-                    continue
+            for param in parameters_:
                 gradient_circuits.append(lin_comb_circuits[param])
             # Combine inputs into a single job to reduce overhead.
             n = len(gradient_circuits)
diff --git a/qiskit_algorithms/gradients/param_shift_estimator_gradient.py b/qiskit_algorithms/gradients/param_shift_estimator_gradient.py
index 7ab0c8d8..89244149 100644
--- a/qiskit_algorithms/gradients/param_shift_estimator_gradient.py
+++ b/qiskit_algorithms/gradients/param_shift_estimator_gradient.py
@@ -15,7 +15,7 @@
 
 from __future__ import annotations
 
-from typing import Sequence
+from collections.abc import Sequence
 
 from qiskit.circuit import Parameter, QuantumCircuit
 from qiskit.opflow import PauliSumOp
@@ -60,36 +60,36 @@ def _run(
         circuits: Sequence[QuantumCircuit],
         observables: Sequence[BaseOperator | PauliSumOp],
         parameter_values: Sequence[Sequence[float]],
-        parameter_sets: Sequence[set[Parameter]],
+        parameters: Sequence[Sequence[Parameter]],
         **options,
     ) -> EstimatorGradientResult:
         """Compute the gradients of the expectation values by the parameter shift rule."""
-        g_circuits, g_parameter_values, g_parameter_sets = self._preprocess(
-            circuits, parameter_values, parameter_sets, self.SUPPORTED_GATES
+        g_circuits, g_parameter_values, g_parameters = self._preprocess(
+            circuits, parameter_values, parameters, self.SUPPORTED_GATES
         )
         results = self._run_unique(
-            g_circuits, observables, g_parameter_values, g_parameter_sets, **options
+            g_circuits, observables, g_parameter_values, g_parameters, **options
         )
-        return self._postprocess(results, circuits, parameter_values, parameter_sets)
+        return self._postprocess(results, circuits, parameter_values, parameters)
 
     def _run_unique(
         self,
         circuits: Sequence[QuantumCircuit],
         observables: Sequence[BaseOperator | PauliSumOp],
         parameter_values: Sequence[Sequence[float]],
-        parameter_sets: Sequence[set[Parameter]],
+        parameters: Sequence[Sequence[Parameter]],
         **options,
     ) -> EstimatorGradientResult:
         """Compute the estimator gradients on the given circuits."""
         job_circuits, job_observables, job_param_values, metadata = [], [], [], []
         all_n = []
-        for circuit, observable, parameter_values_, parameter_set in zip(
-            circuits, observables, parameter_values, parameter_sets
+        for circuit, observable, parameter_values_, parameters_ in zip(
+            circuits, observables, parameter_values, parameters
         ):
-            metadata.append({"parameters": [p for p in circuit.parameters if p in parameter_set]})
+            metadata.append({"parameters": parameters_})
             # Make parameter values for the parameter shift rule.
             param_shift_parameter_values = _make_param_shift_parameter_values(
-                circuit, parameter_values_, parameter_set
+                circuit, parameter_values_, parameters_
             )
             # Combine inputs into a single job to reduce overhead.
             n = len(param_shift_parameter_values)
diff --git a/qiskit_algorithms/gradients/param_shift_sampler_gradient.py b/qiskit_algorithms/gradients/param_shift_sampler_gradient.py
index ff8301c2..af40e5f0 100644
--- a/qiskit_algorithms/gradients/param_shift_sampler_gradient.py
+++ b/qiskit_algorithms/gradients/param_shift_sampler_gradient.py
@@ -58,33 +58,31 @@ def _run(
         self,
         circuits: Sequence[QuantumCircuit],
         parameter_values: Sequence[Sequence[float]],
-        parameter_sets: Sequence[set[Parameter]],
+        parameters: Sequence[Sequence[Parameter]],
         **options,
     ) -> SamplerGradientResult:
         """Compute the estimator gradients on the given circuits."""
-        g_circuits, g_parameter_values, g_parameter_sets = self._preprocess(
-            circuits, parameter_values, parameter_sets, self.SUPPORTED_GATES
+        g_circuits, g_parameter_values, g_parameters = self._preprocess(
+            circuits, parameter_values, parameters, self.SUPPORTED_GATES
         )
-        results = self._run_unique(g_circuits, g_parameter_values, g_parameter_sets, **options)
-        return self._postprocess(results, circuits, parameter_values, parameter_sets)
+        results = self._run_unique(g_circuits, g_parameter_values, g_parameters, **options)
+        return self._postprocess(results, circuits, parameter_values, parameters)
 
     def _run_unique(
         self,
         circuits: Sequence[QuantumCircuit],
         parameter_values: Sequence[Sequence[float]],
-        parameter_sets: Sequence[set[Parameter]],
+        parameters: Sequence[Sequence[Parameter]],
         **options,
     ) -> SamplerGradientResult:
         """Compute the sampler gradients on the given circuits."""
         job_circuits, job_param_values, metadata = [], [], []
         all_n = []
-        for circuit, parameter_values_, parameter_set in zip(
-            circuits, parameter_values, parameter_sets
-        ):
-            metadata.append({"parameters": [p for p in circuit.parameters if p in parameter_set]})
+        for circuit, parameter_values_, parameters_ in zip(circuits, parameter_values, parameters):
+            metadata.append({"parameters": parameters_})
             # Make parameter values for the parameter shift rule.
             param_shift_parameter_values = _make_param_shift_parameter_values(
-                circuit, parameter_values_, parameter_set
+                circuit, parameter_values_, parameters_
             )
             # Combine inputs into a single job to reduce overhead.
             n = len(param_shift_parameter_values)
diff --git a/qiskit_algorithms/gradients/qfi.py b/qiskit_algorithms/gradients/qfi.py
index ee0adf32..5a0e083a 100644
--- a/qiskit_algorithms/gradients/qfi.py
+++ b/qiskit_algorithms/gradients/qfi.py
@@ -89,20 +89,20 @@ def run(
 
         if parameters is None:
             # If parameters is None, we calculate the gradients of all parameters in each circuit.
-            parameter_sets = [set(circuit.parameters) for circuit in circuits]
+            parameters = [circuit.parameters for circuit in circuits]
         else:
             # If parameters is not None, we calculate the gradients of the specified parameters.
             # None in parameters means that the gradients of all parameters in the corresponding
             # circuit are calculated.
-            parameter_sets = [
-                set(parameters_) if parameters_ is not None else set(circuits[i].parameters)
-                for i, parameters_ in enumerate(parameters)
+            parameters = [
+                params if params is not None else circuits[i].parameters
+                for i, params in enumerate(parameters)
             ]
         # The priority of run option is as follows:
         # options in ``run`` method > QFI's default options > QGT's default setting.
         opts = copy(self._default_options)
         opts.update_options(**options)
-        job = AlgorithmJob(self._run, circuits, parameter_values, parameter_sets, **opts.__dict__)
+        job = AlgorithmJob(self._run, circuits, parameter_values, parameters, **opts.__dict__)
         job.submit()
         return job
 
@@ -110,7 +110,7 @@ def _run(
         self,
         circuits: Sequence[QuantumCircuit],
         parameter_values: Sequence[Sequence[float]],
-        parameter_sets: Sequence[Sequence[Parameter] | None] | None = None,
+        parameters: Sequence[Sequence[Parameter]],
         **options,
     ) -> QFIResult:
         """Compute the QFI on the given circuits."""
@@ -119,7 +119,7 @@ def _run(
             self._qgt.derivative_type,
             DerivativeType.REAL,
         )
-        job = self._qgt.run(circuits, parameter_values, parameter_sets, **options)
+        job = self._qgt.run(circuits, parameter_values, parameters, **options)
         self._qgt.derivative_type = temp_derivative_type
 
         try:
diff --git a/qiskit_algorithms/gradients/reverse_gradient/reverse_gradient.py b/qiskit_algorithms/gradients/reverse_gradient/reverse_gradient.py
index 004d09f1..e8d5f1f3 100644
--- a/qiskit_algorithms/gradients/reverse_gradient/reverse_gradient.py
+++ b/qiskit_algorithms/gradients/reverse_gradient/reverse_gradient.py
@@ -77,24 +77,24 @@ def _run(
         circuits: Sequence[QuantumCircuit],
         observables: Sequence[BaseOperator | PauliSumOp],
         parameter_values: Sequence[Sequence[float]],
-        parameter_sets: Sequence[set[Parameter]],
+        parameters: Sequence[Sequence[Parameter]],
         **options,
     ) -> EstimatorGradientResult:
         """Compute the gradients of the expectation values by the parameter shift rule."""
-        g_circuits, g_parameter_values, g_parameter_sets = self._preprocess(
-            circuits, parameter_values, parameter_sets, self.SUPPORTED_GATES
+        g_circuits, g_parameter_values, g_parameters = self._preprocess(
+            circuits, parameter_values, parameters, self.SUPPORTED_GATES
         )
         results = self._run_unique(
-            g_circuits, observables, g_parameter_values, g_parameter_sets, **options
+            g_circuits, observables, g_parameter_values, g_parameters, **options
         )
-        return self._postprocess(results, circuits, parameter_values, parameter_sets)
+        return self._postprocess(results, circuits, parameter_values, parameters)
 
     def _run_unique(
         self,
         circuits: Sequence[QuantumCircuit],
         observables: Sequence[BaseOperator | PauliSumOp],
         parameter_values: Sequence[Sequence[float]],
-        parameter_sets: Sequence[set[Parameter]],
+        parameters: Sequence[Sequence[Parameter]],
         **options,  # pylint: disable=unused-argument
     ) -> EstimatorGradientResult:
         num_gradients = len(circuits)
@@ -104,25 +104,25 @@ def _run_unique(
         for i in range(num_gradients):
             # temporary variables for easier access
             circuit = circuits[i]
-            parameters = parameter_sets[i]
+            parameters_ = parameters[i]
             observable = observables[i]
             values = parameter_values[i]
 
             # the metadata only contains the parameters as there are no run configs here
             metadata.append(
                 {
-                    "parameters": [p for p in circuits[i].parameters if p in parameter_sets[i]],
+                    "parameters": parameters_,
                     "derivative_type": self.derivative_type,
                 }
             )
 
             # keep track of the parameter order of the circuit, as the circuit splitting might
             # produce a list of unitaries in a different order
-            original_parameter_order = [p for p in circuit.parameters if p in parameters]
+            # original_parameter_order = [p for p in circuit.parameters if p in parameters_]
 
             # split the circuit and generate lists of unitaries [U_1, U_2, ...] and
             # parameters [p_1, p_2, ...] in these unitaries
-            unitaries, paramlist = split(circuit, parameters=parameters)
+            unitaries, paramlist = split(circuit, parameters=parameters_)
 
             parameter_binds = dict(zip(circuit.parameters, values))
             bound_circuit = bind(circuit, parameter_binds)
@@ -132,7 +132,7 @@ def _run_unique(
             lam = _evolve_by_operator(observable, phi)
 
             # store gradients in a dictionary to return them in the correct order
-            grads = {param: 0j for param in original_parameter_order}
+            grads = {param: 0j for param in parameters_}
 
             num_parameters = len(unitaries)
             for j in reversed(range(num_parameters)):
diff --git a/qiskit_algorithms/gradients/spsa_estimator_gradient.py b/qiskit_algorithms/gradients/spsa_estimator_gradient.py
index 96fff0d7..5b17713f 100644
--- a/qiskit_algorithms/gradients/spsa_estimator_gradient.py
+++ b/qiskit_algorithms/gradients/spsa_estimator_gradient.py
@@ -76,20 +76,18 @@ def _run(
         circuits: Sequence[QuantumCircuit],
         observables: Sequence[BaseOperator | PauliSumOp],
         parameter_values: Sequence[Sequence[float]],
-        parameter_sets: Sequence[set[Parameter] | None],
+        parameters: Sequence[Sequence[Parameter]],
         **options,
     ) -> EstimatorGradientResult:
         """Compute the estimator gradients on the given circuits."""
         job_circuits, job_observables, job_param_values, metadata, offsets = [], [], [], [], []
         all_n = []
-        for circuit, observable, parameter_values_, parameter_set in zip(
-            circuits, observables, parameter_values, parameter_sets
+        for circuit, observable, parameter_values_, parameters_ in zip(
+            circuits, observables, parameter_values, parameters
         ):
             # Indices of parameters to be differentiated.
-            indices = [
-                circuit.parameters.data.index(p) for p in circuit.parameters if p in parameter_set
-            ]
-            metadata.append({"parameters": [circuit.parameters[idx] for idx in indices]})
+            indices = [circuit.parameters.data.index(p) for p in parameters_]
+            metadata.append({"parameters": parameters_})
             # Make random perturbation vectors.
             offset = [
                 (-1) ** (self._seed.integers(0, 2, len(circuit.parameters)))
diff --git a/qiskit_algorithms/gradients/spsa_sampler_gradient.py b/qiskit_algorithms/gradients/spsa_sampler_gradient.py
index e76bf1c4..622da468 100644
--- a/qiskit_algorithms/gradients/spsa_sampler_gradient.py
+++ b/qiskit_algorithms/gradients/spsa_sampler_gradient.py
@@ -74,20 +74,16 @@ def _run(
         self,
         circuits: Sequence[QuantumCircuit],
         parameter_values: Sequence[Sequence[float]],
-        parameter_sets: Sequence[Sequence[Parameter] | None],
+        parameters: Sequence[Sequence[Parameter]],
         **options,
     ) -> SamplerGradientResult:
         """Compute the sampler gradients on the given circuits."""
         job_circuits, job_param_values, metadata, offsets = [], [], [], []
         all_n = []
-        for circuit, parameter_values_, parameter_set in zip(
-            circuits, parameter_values, parameter_sets
-        ):
+        for circuit, parameter_values_, parameters_ in zip(circuits, parameter_values, parameters):
             # Indices of parameters to be differentiated.
-            indices = [
-                circuit.parameters.data.index(p) for p in circuit.parameters if p in parameter_set
-            ]
-            metadata.append({"parameters": [circuit.parameters[idx] for idx in indices]})
+            indices = [circuit.parameters.data.index(p) for p in parameters_]
+            metadata.append({"parameters": parameters_})
             offset = np.array(
                 [
                     (-1) ** (self._seed.integers(0, 2, len(circuit.parameters)))
@@ -126,11 +122,7 @@ def _run(
                 dist_diffs[j] = dist_diff
             gradient = []
             indices = [circuits[i].parameters.data.index(p) for p in metadata[i]["parameters"]]
-            for j in range(circuits[i].num_parameters):
-                if not j in indices:
-                    continue
-                # the gradient for jth parameter is the average of the gradients of the jth parameter
-                # for each batch.
+            for j in indices:
                 gradient_j = defaultdict(float)
                 for k in range(self._batch_size):
                     for key, value in dist_diffs[k].items():
diff --git a/qiskit_algorithms/gradients/utils.py b/qiskit_algorithms/gradients/utils.py
index 6f5bc156..f3235f2c 100644
--- a/qiskit_algorithms/gradients/utils.py
+++ b/qiskit_algorithms/gradients/utils.py
@@ -17,6 +17,7 @@
 from __future__ import annotations
 
 from collections import defaultdict
+from collections.abc import Sequence
 from dataclasses import dataclass
 from enum import Enum
 
@@ -88,20 +89,20 @@ class LinearCombGradientCircuit:
 def _make_param_shift_parameter_values(
     circuit: QuantumCircuit,
     parameter_values: np.ndarray,
-    parameter_set: set[Parameter],
+    parameters: Sequence[Parameter],
 ) -> list[np.ndarray]:
     """Returns a list of parameter values with offsets for parameter shift rule.
 
     Args:
         circuit: The original quantum circuit
         parameter_values: parameter values to be added to the base parameter values.
-        param_set: set of parameters to be differentiated
+        parameters: The parameters to be shifted.
 
     Returns:
         A list of parameter values with offsets for parameter shift rule.
     """
     plus_offsets, minus_offsets = [], []
-    indices = [idx for idx, param in enumerate(circuit.parameters) if param in parameter_set]
+    indices = [circuit.parameters.data.index(p) for p in parameters]
     offset = np.identity(circuit.num_parameters)[indices, :]
     plus_offsets = parameter_values + offset * np.pi / 2
     minus_offsets = parameter_values - offset * np.pi / 2
@@ -315,7 +316,6 @@ def _assign_unique_parameters(
                 parameter_map[parameter].append(new_parameter, 1)
                 num_gradient_parameters += 1
         gradient_circuit.global_phase = gradient_circuit.global_phase.subs(substitution_map)
-
     return GradientCircuit(gradient_circuit, parameter_map, gradient_parameter_map)
 
 
@@ -346,21 +346,23 @@ def _make_gradient_parameter_values(
     return g_parameter_values
 
 
-def _make_gradient_parameter_set(
+def _make_gradient_parameters(
     gradient_circuit: GradientCircuit,
-    parameter_set: set[Parameter],
-) -> set[Parameter]:
+    parameters: Sequence[Parameter],
+) -> Sequence[Parameter]:
     """Makes parameter set for the gradient circuit.
 
     Args:
         gradient_circuit: The gradient circuit
-        parameters: The parameters for the original circuit
+        parameters: The parameters in the original circuit to calculate gradients
 
     Returns:
-        The parameters for the gradient circuit.
+        The parameters in the gradient circuit to calculate gradients.
     """
-    return {
+    g_parameters = [
         g_parameter
-        for parameter in parameter_set
+        for parameter in parameters
         for g_parameter, _ in gradient_circuit.parameter_map[parameter]
-    }
+    ]
+    # make g_parameters unique and return it.
+    return list(dict.fromkeys(g_parameters))
diff --git a/test/gradients/test_estimator_gradient.py b/test/gradients/test_estimator_gradient.py
index 00b9c807..c8445031 100644
--- a/test/gradients/test_estimator_gradient.py
+++ b/test/gradients/test_estimator_gradient.py
@@ -221,6 +221,30 @@ def test_gradient_parameters(self, grad):
             gradients = gradient.run([qc], [op], [param], parameters=[[a]]).result().gradients[0]
             np.testing.assert_allclose(gradients, correct_results[i], atol=1e-3)
 
+        # parameter order
+        with self.subTest(msg="The order of gradients"):
+            c = Parameter("c")
+            qc = QuantumCircuit(1)
+            qc.rx(a, 0)
+            qc.rz(b, 0)
+            qc.rx(c, 0)
+
+            param_list = [[np.pi / 4, np.pi / 2, np.pi / 3]]
+            correct_results = [
+                [-0.35355339, 0.61237244, -0.61237244],
+                [-0.61237244, 0.61237244, -0.35355339],
+                [-0.35355339, -0.61237244],
+                [-0.61237244, -0.35355339],
+            ]
+            param = [[a, b, c], [c, b, a], [a, c], [c, a]]
+            op = SparsePauliOp.from_list([("Z", 1)])
+            for i, p in enumerate(param):
+                gradient = grad(estimator)
+                gradients = (
+                    gradient.run([qc], [op], param_list, parameters=[p]).result().gradients[0]
+                )
+                np.testing.assert_allclose(gradients, correct_results[i], atol=1e-3)
+
     @data(*gradient_factories)
     def test_gradient_multi_arguments(self, grad):
         """Test the estimator gradient for multiple arguments"""
@@ -260,19 +284,14 @@ def test_gradient_multi_arguments(self, grad):
         np.testing.assert_allclose(gradients2[1], correct_results2[1], atol=1e-3)
         np.testing.assert_allclose(gradients2[2], correct_results2[2], atol=1e-3)
 
-    @data(FiniteDiffEstimatorGradient, ParamShiftEstimatorGradient, LinCombEstimatorGradient)
+    @data(*gradient_factories)
     def test_gradient_validation(self, grad):
         """Test estimator gradient's validation"""
         estimator = Estimator()
         a = Parameter("a")
         qc = QuantumCircuit(1)
         qc.rx(a, 0)
-        if grad is FiniteDiffEstimatorGradient:
-            gradient = grad(estimator, epsilon=1e-6)
-            with self.assertRaises(ValueError):
-                _ = grad(estimator, epsilon=-0.1)
-        else:
-            gradient = grad(estimator)
+        gradient = grad(estimator)
         param_list = [[np.pi / 4], [np.pi / 2]]
         op = SparsePauliOp.from_list([("Z", 1)])
         with self.assertRaises(ValueError):
@@ -302,22 +321,48 @@ def test_spsa_gradient(self):
         np.testing.assert_allclose(gradients, correct_results, atol=1e-3)
 
         # multi parameters
-        gradient = SPSAEstimatorGradient(estimator, epsilon=1e-6, seed=123)
-        param_list2 = [[1, 1], [1, 1], [3, 3]]
-        gradients2 = (
-            gradient.run([qc] * 3, [op] * 3, param_list2, parameters=[None, [b], None])
-            .result()
-            .gradients
-        )
-        correct_results2 = [[-0.84147098, 0.84147098], [0.84147098], [-0.14112001, 0.14112001]]
-        for grad, correct in zip(gradients2, correct_results2):
-            np.testing.assert_allclose(grad, correct, atol=1e-3)
+        with self.subTest(msg="Multiple parameters"):
+            gradient = SPSAEstimatorGradient(estimator, epsilon=1e-6, seed=123)
+            param_list2 = [[1, 1], [1, 1], [3, 3]]
+            gradients2 = (
+                gradient.run([qc] * 3, [op] * 3, param_list2, parameters=[None, [b], None])
+                .result()
+                .gradients
+            )
+            correct_results2 = [[-0.84147098, 0.84147098], [0.84147098], [-0.14112001, 0.14112001]]
+            for grad, correct in zip(gradients2, correct_results2):
+                np.testing.assert_allclose(grad, correct, atol=1e-3)
 
         # batch size
-        correct_results = [[-0.84147098, 0.1682942]]
-        gradient = SPSAEstimatorGradient(estimator, epsilon=1e-6, batch_size=5, seed=123)
-        gradients = gradient.run([qc], [op], param_list).result().gradients
-        np.testing.assert_allclose(gradients, correct_results, atol=1e-3)
+        with self.subTest(msg="Batch size"):
+            correct_results = [[-0.84147098, 0.1682942]]
+            gradient = SPSAEstimatorGradient(estimator, epsilon=1e-6, batch_size=5, seed=123)
+            gradients = gradient.run([qc], [op], param_list).result().gradients
+            np.testing.assert_allclose(gradients, correct_results, atol=1e-3)
+
+        # parameter order
+        with self.subTest(msg="The order of gradients"):
+            gradient = SPSAEstimatorGradient(estimator, epsilon=1e-6, seed=123)
+            c = Parameter("c")
+            qc = QuantumCircuit(1)
+            qc.rx(a, 0)
+            qc.rz(b, 0)
+            qc.rx(c, 0)
+            op = SparsePauliOp.from_list([("Z", 1)])
+            param_list3 = [[np.pi / 4, np.pi / 2, np.pi / 3]]
+            param = [[a, b, c], [c, b, a], [a, c], [c, a]]
+            expected = [
+                [-0.3535525, 0.3535525, 0.3535525],
+                [0.3535525, 0.3535525, -0.3535525],
+                [-0.3535525, 0.3535525],
+                [0.3535525, -0.3535525],
+            ]
+            for i, p in enumerate(param):
+                gradient = SPSAEstimatorGradient(estimator, epsilon=1e-6, seed=123)
+                gradients = (
+                    gradient.run([qc], [op], param_list3, parameters=[p]).result().gradients[0]
+                )
+                np.testing.assert_allclose(gradients, expected[i], atol=1e-3)
 
     @data(ParamShiftEstimatorGradient, LinCombEstimatorGradient)
     def test_gradient_random_parameters(self, grad):
diff --git a/test/gradients/test_qgt.py b/test/gradients/test_qgt.py
index 044fabdb..74f22d46 100644
--- a/test/gradients/test_qgt.py
+++ b/test/gradients/test_qgt.py
@@ -154,7 +154,7 @@ def test_qgt_coefficients(self, qgt_type):
             np.testing.assert_allclose(qgt_result[0], correct_values[i], atol=1e-3)
 
     @data(LinCombQGT, ReverseQGT)
-    def test_qgt_specify_parameters(self, qgt_type):
+    def test_qgt_parameters(self, qgt_type):
         """Test the QGT with specified parameters"""
         args = () if qgt_type == ReverseQGT else (self.estimator,)
         qgt = qgt_type(*args, derivative_type=DerivativeType.REAL)
@@ -164,10 +164,40 @@ def test_qgt_specify_parameters(self, qgt_type):
         qc = QuantumCircuit(1)
         qc.rx(a, 0)
         qc.ry(b, 0)
-        param_list = [np.pi / 4, np.pi / 4]
-        qgt_result = qgt.run([qc], [param_list], [[a]]).result().qgts
+        param_values = [np.pi / 4, np.pi / 4]
+        qgt_result = qgt.run([qc], [param_values], [[a]]).result().qgts
         np.testing.assert_allclose(qgt_result[0], [[1 / 4]], atol=1e-3)
 
+        with self.subTest("Test with different parameter orders"):
+            c = Parameter("c")
+            qc2 = QuantumCircuit(1)
+            qc2.rx(a, 0)
+            qc2.rz(b, 0)
+            qc2.rx(c, 0)
+            param_values = [np.pi / 4, np.pi / 4, np.pi / 4]
+            params = [[a, b, c], [c, b, a], [a, c], [b, a]]
+            expected = [
+                np.array(
+                    [
+                        [0.25, 0.0, 0.1767767],
+                        [0.0, 0.125, -0.08838835],
+                        [0.1767767, -0.08838835, 0.1875],
+                    ]
+                ),
+                np.array(
+                    [
+                        [0.1875, -0.08838835, 0.1767767],
+                        [-0.08838835, 0.125, 0.0],
+                        [0.1767767, 0.0, 0.25],
+                    ]
+                ),
+                np.array([[0.25, 0.1767767], [0.1767767, 0.1875]]),
+                np.array([[0.125, 0.0], [0.0, 0.25]]),
+            ]
+            for i, param in enumerate(params):
+                qgt_result = qgt.run([qc2], [param_values], [param]).result().qgts
+                np.testing.assert_allclose(qgt_result[0], expected[i], atol=1e-3)
+
     @data(LinCombQGT, ReverseQGT)
     def test_qgt_multi_arguments(self, qgt_type):
         """Test the QGT for multiple arguments"""
diff --git a/test/gradients/test_sampler_gradient.py b/test/gradients/test_sampler_gradient.py
index 7cff9c54..b4caac07 100644
--- a/test/gradients/test_sampler_gradient.py
+++ b/test/gradients/test_sampler_gradient.py
@@ -28,7 +28,7 @@
 )
 from qiskit.circuit import Parameter
 from qiskit.circuit.library import EfficientSU2, RealAmplitudes
-from qiskit.circuit.library.standard_gates import RXXGate, RYYGate, RZXGate, RZZGate
+from qiskit.circuit.library.standard_gates import RXXGate
 from qiskit.primitives import Sampler
 from qiskit.result import QuasiDistribution
 from qiskit.test import QiskitTestCase
@@ -60,16 +60,16 @@ def test_single_circuit(self, grad):
         qc.measure_all()
         gradient = grad(sampler)
         param_list = [[np.pi / 4], [0], [np.pi / 2]]
-        correct_results = [
+        expected = [
             [{0: -0.5 / np.sqrt(2), 1: 0.5 / np.sqrt(2)}],
             [{0: 0, 1: 0}],
             [{0: -0.499999, 1: 0.499999}],
         ]
         for i, param in enumerate(param_list):
-            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)
+            gradients = gradient.run([qc], [param]).result().gradients[0]
+            array1 = _quasi2array(gradients, num_qubits=1)
+            array2 = _quasi2array(expected[i], num_qubits=1)
+            np.testing.assert_allclose(array1, array2, atol=1e-3)
 
     @data(*gradient_factories)
     def test_gradient_p(self, grad):
@@ -83,16 +83,16 @@ def test_gradient_p(self, grad):
         qc.measure_all()
         gradient = grad(sampler)
         param_list = [[np.pi / 4], [0], [np.pi / 2]]
-        correct_results = [
+        expected = [
             [{0: -0.5 / np.sqrt(2), 1: 0.5 / np.sqrt(2)}],
             [{0: 0, 1: 0}],
             [{0: -0.499999, 1: 0.499999}],
         ]
         for i, param in enumerate(param_list):
             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)
+            array1 = _quasi2array(gradients, num_qubits=1)
+            array2 = _quasi2array(expected[i], num_qubits=1)
+            np.testing.assert_allclose(array1, array2, atol=1e-3)
 
     @data(*gradient_factories)
     def test_gradient_u(self, grad):
@@ -108,15 +108,15 @@ def test_gradient_u(self, grad):
         qc.measure_all()
         gradient = grad(sampler)
         param_list = [[np.pi / 4, 0, 0], [np.pi / 4, np.pi / 4, np.pi / 4]]
-        correct_results = [
+        expected = [
             [{0: -0.5 / np.sqrt(2), 1: 0.5 / np.sqrt(2)}, {0: 0, 1: 0}, {0: 0, 1: 0}],
             [{0: -0.176777, 1: 0.176777}, {0: -0.426777, 1: 0.426777}, {0: -0.426777, 1: 0.426777}],
         ]
         for i, param in enumerate(param_list):
             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)
+            array1 = _quasi2array(gradients, num_qubits=1)
+            array2 = _quasi2array(expected[i], num_qubits=1)
+            np.testing.assert_allclose(array1, array2, atol=1e-3)
 
     @data(*gradient_factories)
     def test_gradient_efficient_su2(self, grad):
@@ -129,7 +129,7 @@ def test_gradient_efficient_su2(self, grad):
             [np.pi / 4 for param in qc.parameters],
             [np.pi / 2 for param in qc.parameters],
         ]
-        correct_results = [
+        expected = [
             [
                 {
                     0: -0.11963834764831836,
@@ -208,37 +208,30 @@ def test_gradient_efficient_su2(self, grad):
         ]
         for i, param in enumerate(param_list):
             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)
+            array1 = _quasi2array(gradients, num_qubits=2)
+            array2 = _quasi2array(expected[i], num_qubits=2)
+            np.testing.assert_allclose(array1, array2, atol=1e-3)
 
     @data(*gradient_factories)
     def test_gradient_2qubit_gate(self, grad):
         """Test the sampler gradient for 2 qubit gates"""
         sampler = Sampler()
-        for gate in [RXXGate, RYYGate, RZZGate, RZXGate]:
+        for gate in [RXXGate]:
             param_list = [[np.pi / 4], [np.pi / 2]]
-
-            if gate is RZXGate:
-                correct_results = [
-                    [{0: -0.5 / np.sqrt(2), 1: 0, 2: 0.5 / np.sqrt(2), 3: 0}],
-                    [{0: -0.5, 1: 0, 2: 0.5, 3: 0}],
-                ]
-            else:
-                correct_results = [
-                    [{0: -0.5 / np.sqrt(2), 1: 0, 2: 0, 3: 0.5 / np.sqrt(2)}],
-                    [{0: -0.5, 1: 0, 2: 0, 3: 0.5}],
-                ]
-        for i, param in enumerate(param_list):
-            a = Parameter("a")
-            qc = QuantumCircuit(2)
-            qc.append(gate(a), [qc.qubits[0], qc.qubits[1]], [])
-            qc.measure_all()
-            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)
+            correct_results = [
+                [{0: -0.5 / np.sqrt(2), 1: 0, 2: 0, 3: 0.5 / np.sqrt(2)}],
+                [{0: -0.5, 1: 0, 2: 0, 3: 0.5}],
+            ]
+            for i, param in enumerate(param_list):
+                a = Parameter("a")
+                qc = QuantumCircuit(2)
+                qc.append(gate(a), [qc.qubits[0], qc.qubits[1]], [])
+                qc.measure_all()
+                gradient = grad(sampler)
+                gradients = gradient.run([qc], [param]).result().gradients[0]
+                array1 = _quasi2array(gradients, num_qubits=2)
+                array2 = _quasi2array(correct_results[i], num_qubits=2)
+                np.testing.assert_allclose(array1, array2, atol=1e-3)
 
     @data(*gradient_factories)
     def test_gradient_parameter_coefficient(self, grad):
@@ -310,9 +303,9 @@ def test_gradient_parameter_coefficient(self, grad):
 
         for i, param in enumerate(param_list):
             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], 2)
+            array1 = _quasi2array(gradients, num_qubits=2)
+            array2 = _quasi2array(correct_results[i], num_qubits=2)
+            np.testing.assert_allclose(array1, array2, atol=1e-3)
 
     @data(*gradient_factories)
     def test_gradient_parameters(self, grad):
@@ -326,14 +319,50 @@ def test_gradient_parameters(self, grad):
         qc.measure_all()
         gradient = grad(sampler)
         param_list = [[np.pi / 4, np.pi / 2]]
-        correct_results = [
+        expected = [
             [{0: -0.5 / np.sqrt(2), 1: 0.5 / np.sqrt(2)}],
         ]
         for i, param in enumerate(param_list):
             gradients = gradient.run([qc], [param], parameters=[[a]]).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)
+            array1 = _quasi2array(gradients, num_qubits=1)
+            array2 = _quasi2array(expected[i], num_qubits=1)
+            np.testing.assert_allclose(array1, array2, atol=1e-3)
+
+        # parameter order
+        with self.subTest(msg="The order of gradients"):
+            c = Parameter("c")
+            qc = QuantumCircuit(1)
+            qc.rx(a, 0)
+            qc.rz(b, 0)
+            qc.rx(c, 0)
+            qc.measure_all()
+            param_values = [[np.pi / 4, np.pi / 2, np.pi / 3]]
+            params = [[a, b, c], [c, b, a], [a, c], [c, a]]
+            expected = [
+                [
+                    {0: -0.17677666583387008, 1: 0.17677666583378482},
+                    {0: 0.3061861668168149, 1: -0.3061861668167012},
+                    {0: -0.3061861668168149, 1: 0.30618616681678645},
+                ],
+                [
+                    {0: -0.3061861668168149, 1: 0.30618616681678645},
+                    {0: 0.3061861668168149, 1: -0.3061861668167012},
+                    {0: -0.17677666583387008, 1: 0.17677666583378482},
+                ],
+                [
+                    {0: -0.17677666583387008, 1: 0.17677666583378482},
+                    {0: -0.3061861668168149, 1: 0.30618616681678645},
+                ],
+                [
+                    {0: -0.3061861668168149, 1: 0.30618616681678645},
+                    {0: -0.17677666583387008, 1: 0.17677666583378482},
+                ],
+            ]
+            for i, p in enumerate(params):
+                gradients = gradient.run([qc], param_values, parameters=[p]).result().gradients[0]
+                array1 = _quasi2array(gradients, num_qubits=1)
+                array2 = _quasi2array(expected[i], num_qubits=1)
+                np.testing.assert_allclose(array1, array2, atol=1e-3)
 
     @data(*gradient_factories)
     def test_gradient_multi_arguments(self, grad):
@@ -354,33 +383,35 @@ def test_gradient_multi_arguments(self, grad):
             [{0: -0.499999, 1: 0.499999}],
         ]
         gradients = gradient.run([qc, qc2], param_list).result().gradients
-        for j, q_dists in enumerate(gradients):
-            quasi_dist = q_dists[0]
-            for k in quasi_dist:
-                self.assertAlmostEqual(quasi_dist[k], correct_results[j][0][k], 3)
-
-        c = Parameter("c")
-        qc3 = QuantumCircuit(1)
-        qc3.rx(c, 0)
-        qc3.ry(a, 0)
-        qc3.measure_all()
-        param_list2 = [[np.pi / 4], [np.pi / 4, np.pi / 4], [np.pi / 4, np.pi / 4]]
-        gradients = (
-            gradient.run([qc, qc3, qc3], param_list2, parameters=[[a], [c], None])
-            .result()
-            .gradients
-        )
-        correct_results = [
-            [{0: -0.5 / np.sqrt(2), 1: 0.5 / np.sqrt(2)}],
-            [{0: -0.25, 1: 0.25}],
-            [{0: -0.25, 1: 0.25}, {0: -0.25, 1: 0.25}],
-        ]
-        for i, result in enumerate(gradients):
-            for j, q_dists in enumerate(result):
-                for k in q_dists:
-                    self.assertAlmostEqual(q_dists[k], correct_results[i][j][k], 3)
+        for i, q_dists in enumerate(gradients):
+            array1 = _quasi2array(q_dists, num_qubits=1)
+            array2 = _quasi2array(correct_results[i], num_qubits=1)
+            np.testing.assert_allclose(array1, array2, atol=1e-3)
+
+        # parameters
+        with self.subTest(msg="Different parameters"):
+            c = Parameter("c")
+            qc3 = QuantumCircuit(1)
+            qc3.rx(c, 0)
+            qc3.ry(a, 0)
+            qc3.measure_all()
+            param_list2 = [[np.pi / 4], [np.pi / 4, np.pi / 4], [np.pi / 4, np.pi / 4]]
+            gradients = (
+                gradient.run([qc, qc3, qc3], param_list2, parameters=[[a], [c], None])
+                .result()
+                .gradients
+            )
+            correct_results = [
+                [{0: -0.5 / np.sqrt(2), 1: 0.5 / np.sqrt(2)}],
+                [{0: -0.25, 1: 0.25}],
+                [{0: -0.25, 1: 0.25}, {0: -0.25, 1: 0.25}],
+            ]
+            for i, q_dists in enumerate(gradients):
+                array1 = _quasi2array(q_dists, num_qubits=1)
+                array2 = _quasi2array(correct_results[i], num_qubits=1)
+                np.testing.assert_allclose(array1, array2, atol=1e-3)
 
-    @data(FiniteDiffSamplerGradient, ParamShiftSamplerGradient, LinCombSamplerGradient)
+    @data(*gradient_factories)
     def test_gradient_validation(self, grad):
         """Test sampler gradient's validation"""
         sampler = Sampler()
@@ -388,12 +419,7 @@ def test_gradient_validation(self, grad):
         qc = QuantumCircuit(1)
         qc.rx(a, 0)
         qc.measure_all()
-        if grad is FiniteDiffSamplerGradient:
-            gradient = grad(sampler, epsilon=1e-6)
-            with self.assertRaises(ValueError):
-                _ = grad(sampler, epsilon=-0.1)
-        else:
-            gradient = grad(sampler)
+        gradient = grad(sampler)
         param_list = [[np.pi / 4], [np.pi / 2]]
         with self.assertRaises(ValueError):
             gradient.run([qc], param_list)
@@ -410,6 +436,7 @@ def test_spsa_gradient(self):
 
         a = Parameter("a")
         b = Parameter("b")
+        c = Parameter("c")
         qc = QuantumCircuit(2)
         qc.rx(b, 0)
         qc.rx(a, 1)
@@ -424,58 +451,96 @@ def test_spsa_gradient(self):
         gradient = SPSASamplerGradient(sampler, epsilon=1e-6, seed=123)
         for i, param in enumerate(param_list):
             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)
-        # multi parameters
-        param_list2 = [[1, 2], [1, 2], [3, 4]]
-        correct_results2 = [
-            [
-                {0: 0.2273244, 1: -0.6480598, 2: 0.2273244, 3: 0.1934111},
-                {0: -0.2273244, 1: 0.6480598, 2: -0.2273244, 3: -0.1934111},
-            ],
-            [
-                {0: -0.2273244, 1: 0.6480598, 2: -0.2273244, 3: -0.1934111},
-            ],
-            [
-                {0: -0.0141129, 1: -0.0564471, 2: -0.3642884, 3: 0.4348484},
-                {0: 0.0141129, 1: 0.0564471, 2: 0.3642884, 3: -0.4348484},
-            ],
-        ]
-        gradient = SPSASamplerGradient(sampler, epsilon=1e-6, seed=123)
-        gradients = (
-            gradient.run([qc] * 3, param_list2, parameters=[None, [b], None]).result().gradients
-        )
+            array1 = _quasi2array(gradients, num_qubits=2)
+            array2 = _quasi2array(correct_results[i], num_qubits=2)
+            np.testing.assert_allclose(array1, array2, atol=1e-3)
 
-        for i, result in enumerate(gradients):
-            for j, q_dists in enumerate(result):
-                for k in q_dists:
-                    self.assertAlmostEqual(q_dists[k], correct_results2[i][j][k], 3)
+        # multi parameters
+        with self.subTest(msg="Multiple parameters"):
+            param_list2 = [[1, 2], [1, 2], [3, 4]]
+            correct_results2 = [
+                [
+                    {0: 0.2273244, 1: -0.6480598, 2: 0.2273244, 3: 0.1934111},
+                    {0: -0.2273244, 1: 0.6480598, 2: -0.2273244, 3: -0.1934111},
+                ],
+                [
+                    {0: -0.2273244, 1: 0.6480598, 2: -0.2273244, 3: -0.1934111},
+                ],
+                [
+                    {0: -0.0141129, 1: -0.0564471, 2: -0.3642884, 3: 0.4348484},
+                    {0: 0.0141129, 1: 0.0564471, 2: 0.3642884, 3: -0.4348484},
+                ],
+            ]
+            gradient = SPSASamplerGradient(sampler, epsilon=1e-6, seed=123)
+            gradients = (
+                gradient.run([qc] * 3, param_list2, parameters=[None, [b], None]).result().gradients
+            )
+            for i, result in enumerate(gradients):
+                array1 = _quasi2array(result, num_qubits=2)
+                array2 = _quasi2array(correct_results2[i], num_qubits=2)
+                np.testing.assert_allclose(array1, array2, atol=1e-3)
 
         # batch size
-        param_list = [[1, 1]]
-        gradient = SPSASamplerGradient(sampler, epsilon=1e-6, batch_size=4, seed=123)
-        gradients = gradient.run([qc], param_list).result().gradients
-        correct_results3 = [
-            [
-                {
-                    0: -0.1620149622932887,
-                    1: -0.25872053011771756,
-                    2: 0.3723827084675668,
-                    3: 0.04835278392088804,
-                },
-                {
-                    0: -0.1620149622932887,
-                    1: 0.3723827084675668,
-                    2: -0.25872053011771756,
-                    3: 0.04835278392088804,
-                },
+        with self.subTest(msg="Batch size"):
+            param_list = [[1, 1]]
+            gradient = SPSASamplerGradient(sampler, epsilon=1e-6, batch_size=4, seed=123)
+            gradients = gradient.run([qc], param_list).result().gradients
+            correct_results3 = [
+                [
+                    {
+                        0: -0.1620149622932887,
+                        1: -0.25872053011771756,
+                        2: 0.3723827084675668,
+                        3: 0.04835278392088804,
+                    },
+                    {
+                        0: -0.1620149622932887,
+                        1: 0.3723827084675668,
+                        2: -0.25872053011771756,
+                        3: 0.04835278392088804,
+                    },
+                ]
             ]
-        ]
-        for i, result in enumerate(gradients):
-            for j, q_dists in enumerate(result):
-                for k in q_dists:
-                    self.assertAlmostEqual(q_dists[k], correct_results3[i][j][k], 3)
+            for i, q_dists in enumerate(gradients):
+                array1 = _quasi2array(q_dists, num_qubits=2)
+                array2 = _quasi2array(correct_results3[i], num_qubits=2)
+                np.testing.assert_allclose(array1, array2, atol=1e-3)
+
+        # parameter order
+        with self.subTest(msg="The order of gradients"):
+            qc = QuantumCircuit(1)
+            qc.rx(a, 0)
+            qc.rz(b, 0)
+            qc.rx(c, 0)
+            qc.measure_all()
+            param_list = [[np.pi / 4, np.pi / 2, np.pi / 3]]
+            param = [[a, b, c], [c, b, a], [a, c], [c, a]]
+            correct_results = [
+                [
+                    {0: -0.17677624757590138, 1: 0.17677624757590138},
+                    {0: 0.17677624757590138, 1: -0.17677624757590138},
+                    {0: 0.17677624757590138, 1: -0.17677624757590138},
+                ],
+                [
+                    {0: 0.17677624757590138, 1: -0.17677624757590138},
+                    {0: 0.17677624757590138, 1: -0.17677624757590138},
+                    {0: -0.17677624757590138, 1: 0.17677624757590138},
+                ],
+                [
+                    {0: -0.17677624757590138, 1: 0.17677624757590138},
+                    {0: 0.17677624757590138, 1: -0.17677624757590138},
+                ],
+                [
+                    {0: 0.17677624757590138, 1: -0.17677624757590138},
+                    {0: -0.17677624757590138, 1: 0.17677624757590138},
+                ],
+            ]
+            for i, p in enumerate(param):
+                gradient = SPSASamplerGradient(sampler, epsilon=1e-6, seed=123)
+                gradients = gradient.run([qc], param_list, parameters=[p]).result().gradients[0]
+                array1 = _quasi2array(gradients, num_qubits=1)
+                array2 = _quasi2array(correct_results[i], num_qubits=1)
+                np.testing.assert_allclose(array1, array2, atol=1e-3)
 
     @data(ParamShiftSamplerGradient, LinCombSamplerGradient)
     def test_gradient_random_parameters(self, grad):
diff --git a/test/time_evolvers/variational/test_var_qite.py b/test/time_evolvers/variational/test_var_qite.py
index 067d9d32..9c79a002 100644
--- a/test/time_evolvers/variational/test_var_qite.py
+++ b/test/time_evolvers/variational/test_var_qite.py
@@ -79,14 +79,14 @@ def test_run_d_1_with_aux_ops(self):
         ]
 
         thetas_expected_shots = [
-            0.962539564367932,
-            1.87923124299411,
-            2.721576354971394,
-            2.6710294109655432,
-            2.1574398860715926,
-            1.65846084870719,
-            1.872970698074877,
-            1.9438954933172055,
+            0.9392668013702317,
+            1.8756706968454864,
+            2.6915067128662398,
+            2.655420131540562,
+            2.174687086978046,
+            1.6997059390911056,
+            1.8056912289547045,
+            1.939353810908912,
         ]
 
         with self.subTest(msg="Test exact backend."):
@@ -133,7 +133,7 @@ def test_run_d_1_with_aux_ops(self):
 
             parameter_values = evolution_result.parameter_values[-1]
 
-            expected_aux_ops = (-0.220271, 0.266539)
+            expected_aux_ops = (-0.24629853310903974, 0.2518122871921184)
 
             for i, parameter_value in enumerate(parameter_values):
                 np.testing.assert_almost_equal(