def METHOD_NAME(
    self,
    circuits: Sequence[QuantumCircuit],
    parameter_values: Sequence[Sequence[float]],
    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_, parameters_ in zip(circuits, parameter_values, parameters):
        # Indices of parameters to be differentiated.
        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)))
                for _ in range(self._batch_size)
            ]
        )
        plus = [parameter_values_ + self._epsilon * offset_ for offset_ in offset]
        minus = [parameter_values_ - self._epsilon * offset_ for offset_ in offset]
        offsets.append(offset)
        # Combine inputs into a single job to reduce overhead.
        n = 2 * self._batch_size
        job_circuits.extend([circuit] * n)
        job_param_values.extend(plus + minus)
        all_n.append(n)
    # Run the single job with all circuits.
    job = self._sampler.run(job_circuits, job_param_values, **options)
    try:
        results = job.result()
    except Exception as exc:
        raise AlgorithmError("Sampler job failed.") from exc
    # Compute the gradients.
    gradients = []
    partial_sum_n = 0
    for i, n in enumerate(all_n):
        dist_diffs = {}
        result = results.quasi_dists[partial_sum_n : partial_sum_n + n]
        for j, (dist_plus, dist_minus) in enumerate(zip(result[: n // 2], result[n // 2 :])):
            dist_diff: dict[int, float] = defaultdict(float)
            for key, value in dist_plus.items():
                dist_diff[key] += value / (2 * self._epsilon)
            for key, value in dist_minus.items():
                dist_diff[key] -= value / (2 * self._epsilon)
            dist_diffs[j] = dist_diff
        gradient = []
        indices = [circuits[i].parameters.data.index(p) for p in metadata[i]["parameters"]]
        for j in indices:
            gradient_j: dict[int, float] = defaultdict(float)
            for k in range(self._batch_size):
                for key, value in dist_diffs[k].items():
                    gradient_j[key] += value * offsets[i][k][j]
            gradient_j = {key: value / self._batch_size for key, value in gradient_j.items()}
            gradient.append(gradient_j)
        gradients.append(gradient)
        partial_sum_n += n
    opt = self._get_local_options(options)
    return SamplerGradientResult(gradients=gradients, metadata=metadata, options=opt)