Repository: pyRiemann/pyRiemann-qiskit

Issue:
[Child] Investigate QAOA Qiskit implementation

Discussion:

Solution patch:
diff --git a/benchmarks/light_benchmark.py b/benchmarks/light_benchmark.py
index e66aea1..41e7cc7 100644
--- a/benchmarks/light_benchmark.py
+++ b/benchmarks/light_benchmark.py
@@ -11,25 +11,29 @@ to run on Ci with each PRs.
 # Modified from plot_classify_P300_bi.py of pyRiemann
 # License: BSD (3-clause)
 
-from pyriemann.estimation import XdawnCovariances, Shrinkage
-from pyriemann.tangentspace import TangentSpace
-from sklearn.pipeline import make_pipeline
-from sklearn.model_selection import train_test_split
-from sklearn.metrics import balanced_accuracy_score
-from sklearn.preprocessing import LabelEncoder
-from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
-from sklearn.decomposition import PCA
+import sys
+import warnings
+
 from moabb import set_log_level
 from moabb.datasets import bi2012
 from moabb.paradigms import P300
-from pyriemann_qiskit.utils import distance, mean  # noqa
+from pyriemann.estimation import XdawnCovariances, Shrinkage
+from pyriemann.tangentspace import TangentSpace
+from pyriemann_qiskit.classification import QuanticNCH
 from pyriemann_qiskit.pipelines import (
     QuantumClassifierWithDefaultRiemannianPipeline,
     QuantumMDMWithRiemannianPipeline,
 )
-from pyriemann_qiskit.classification import QuanticNCH
-import warnings
-import sys
+from pyriemann_qiskit.utils import distance, mean  # noqa
+from pyriemann_qiskit.utils.hyper_params_factory import create_mixer_rotational_X_gates
+from qiskit_algorithms.optimizers import SPSA
+from sklearn.decomposition import PCA
+from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
+from sklearn.metrics import balanced_accuracy_score
+from sklearn.model_selection import train_test_split
+from sklearn.pipeline import make_pipeline
+from sklearn.preprocessing import LabelEncoder
+
 
 print(__doc__)
 
@@ -114,7 +118,6 @@ pipelines["RG_LDA"] = make_pipeline(
 pipelines["NCH_MIN_HULL"] = make_pipeline(
     XdawnCovariances(
         nfilter=3,
-        # classes=[labels_dict["Target"]],
         estimator="lwf",
         xdawn_estimator="scm",
     ),
@@ -128,6 +131,26 @@ pipelines["NCH_MIN_HULL"] = make_pipeline(
     ),
 )
 
+
+pipelines["NCH_MIN_HULL_QAOACV"] = make_pipeline(
+    XdawnCovariances(
+        nfilter=3,
+        estimator="lwf",
+        xdawn_estimator="scm",
+    ),
+    QuanticNCH(
+        n_hulls_per_class=1,
+        n_samples_per_hull=3,
+        n_jobs=12,
+        subsampling="min",
+        quantum=True,
+        # Provide create_mixer to force QAOA-CV optimization
+        create_mixer=create_mixer_rotational_X_gates(0),
+        shots=100,
+        qaoa_optimizer=SPSA(maxiter=500),
+    ),
+)
+
 ##############################################################################
 # Compute score
 # --------------
diff --git a/doc/api.rst b/doc/api.rst
index 0cdf65f..e89a092 100644
--- a/doc/api.rst
+++ b/doc/api.rst
@@ -93,6 +93,7 @@ Hyper-parameters generation
     create_mixer_rotational_X_gates
     create_mixer_rotational_XY_gates
     create_mixer_rotational_XZ_gates
+    create_mixer_qiskit_default
 
 Preprocessing
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -169,6 +170,7 @@ Math
     cov_to_corr_matrix
     union_of_diff
     to_xyz
+    is_pauli_identity
 
 Firebase
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~
diff --git a/examples/toys_dataset/plot_qaoa_cv.py b/examples/toys_dataset/plot_qaoa_cv.py
index 2dfee07..ed4d80b 100644
--- a/examples/toys_dataset/plot_qaoa_cv.py
+++ b/examples/toys_dataset/plot_qaoa_cv.py
@@ -18,12 +18,11 @@ import math
 
 from docplex.mp.model import Model
 import matplotlib.pyplot as plt
-from qiskit.primitives import BackendSampler
-from qiskit_aer import AerSimulator
 from qiskit_algorithms.optimizers import COBYLA, SPSA
 
 from pyriemann_qiskit.utils.docplex import QAOACVOptimizer
 from pyriemann_qiskit.utils.hyper_params_factory import (
+    create_mixer_qiskit_default,
     create_mixer_rotational_X_gates,
     create_mixer_rotational_XY_gates,
     create_mixer_rotational_XZ_gates,
@@ -46,16 +45,16 @@ def run_qaoa_cv(n_reps, optimizer, create_mixer):
     # objective function to minimize
     mdl.minimize((x - 0.83 + y + 2 * z) ** 2)
 
-    # Define the BackendSampler (previously QuantumInstance)
-    backend = AerSimulator(method="statevector", cuStateVec_enable=True)
-    quantum_instance = BackendSampler(
-        backend, options={"shots": 200, "seed_simulator": 42}
+    # Instanciate the QAOA-CV
+    # Note: if quantum_instance is None, it will be created inside the optimizer.
+    qaoa_cv = QAOACVOptimizer(
+        create_mixer, n_reps, quantum_instance=None, optimizer=optimizer
     )
-    quantum_instance.transpile_options["seed_transpiler"] = 42
 
-    # Instanciate the QAOA-CV
-    qaoa_cv = QAOACVOptimizer(create_mixer, n_reps, quantum_instance, optimizer)
-    solution = qaoa_cv.solve(mdl)
+    # reshape is when working with covariance matrices
+    # So the vector of solution is reshaped into a matrix
+    # (this is not the case here)
+    solution = qaoa_cv.solve(mdl, reshape=False)
 
     # print the time, the solution (that it the value for our three variable)
     # and the minimum of the objective function
@@ -96,6 +95,7 @@ ret = {}
 for angle in range(n_angles):
     angle = math.pi * angle / n_angles
     mixers = [
+        create_mixer_qiskit_default(angle),
         create_mixer_rotational_X_gates(angle),
         create_mixer_rotational_XY_gates(angle),
         create_mixer_rotational_XZ_gates(angle),
diff --git a/pyriemann_qiskit/classification.py b/pyriemann_qiskit/classification.py
index 683d036..42da0e1 100644
--- a/pyriemann_qiskit/classification.py
+++ b/pyriemann_qiskit/classification.py
@@ -28,8 +28,10 @@ from sklearn.svm import SVC
 from .datasets import get_feature_dimension
 from .utils.docplex import (
     set_global_optimizer,
+    get_global_optimizer,
     ClassicalOptimizer,
     NaiveQAOAOptimizer,
+    QAOACVOptimizer,
 )
 from .utils.distance import (
     distance_functions,
@@ -619,6 +621,37 @@ class QuanticVQC(QuanticClassifierBase):
         return 0
 
 
+def _get_docplex_optimizer_from_params_bag(
+    logger,
+    quantum,
+    quantum_instance,
+    upper_bound,
+    qaoa_optimizer,
+    classical_optimizer,
+    create_mixer,
+    n_reps,
+):
+    if quantum:
+        if create_mixer:
+            logger._log("Using QAOACVOptimizer")
+            return QAOACVOptimizer(
+                create_mixer=create_mixer,
+                n_reps=n_reps,
+                quantum_instance=quantum_instance,
+                optimizer=qaoa_optimizer,
+            )
+        else:
+            logger._log("Using NaiveQAOAOptimizer")
+            return NaiveQAOAOptimizer(
+                quantum_instance=quantum_instance,
+                upper_bound=upper_bound,
+                optimizer=qaoa_optimizer,
+            )
+    else:
+        logger._log("Using ClassicalOptimizer (COBYLA)")
+        return ClassicalOptimizer(classical_optimizer)
+
+
 class QuanticMDM(QuanticClassifierBase):
 
     """Quantum-enhanced MDM classifier
@@ -639,6 +672,8 @@ class QuanticMDM(QuanticClassifierBase):
         Add classical_optimizer parameter.
     .. versionchanged:: 0.3.0
         Add qaoa_optimizer parameter.
+    .. versionchanged:: 0.4.0
+        Add QAOA-CV optimization.
 
     Parameters
     ----------
@@ -675,6 +710,13 @@ class QuanticMDM(QuanticClassifierBase):
     qaoa_optimizer : SciPyOptimizer, default=SLSQP()
         An instance of a scipy optimizer to find the optimal weights for the
         parametric circuit (ansatz).
+    create_mixer : None | Callable[int, QuantumCircuit], default=None
+        A delegate that takes into input an angle and returns a QuantumCircuit.
+        This circuit is the mixer operatior for the QAOA-CV algorithm.
+        If None and quantum, the NaiveQAOAOptimizer will be used instead.
+    n_reps : int, default=3
+        The number of time the mixer and cost operator are repeated in the QAOA-CV
+        circuit.
 
     See Also
     --------
@@ -709,6 +751,8 @@ class QuanticMDM(QuanticClassifierBase):
         regularization=None,
         classical_optimizer=CobylaOptimizer(rhobeg=2.1, rhoend=0.000001),
         qaoa_optimizer=SLSQP(),
+        create_mixer=None,
+        n_reps=3,
     ):
         QuanticClassifierBase.__init__(
             self, quantum, q_account_token, verbose, shots, None, seed
@@ -718,6 +762,8 @@ class QuanticMDM(QuanticClassifierBase):
         self.regularization = regularization
         self.classical_optimizer = classical_optimizer
         self.qaoa_optimizer = qaoa_optimizer
+        self.create_mixer = create_mixer
+        self.n_reps = n_reps
 
     @staticmethod
     def _override_predict_distance(mdm):
@@ -753,16 +799,16 @@ class QuanticMDM(QuanticClassifierBase):
         classifier._predict_distances = QuanticMDM._override_predict_distance(
             classifier
         )
-        if self.quantum:
-            self._log("Using NaiveQAOAOptimizer")
-            self._optimizer = NaiveQAOAOptimizer(
-                quantum_instance=self._quantum_instance,
-                upper_bound=self.upper_bound,
-                optimizer=self.qaoa_optimizer,
-            )
-        else:
-            self._log("Using ClassicalOptimizer (COBYLA)")
-            self._optimizer = ClassicalOptimizer(self.classical_optimizer)
+        self._optimizer = _get_docplex_optimizer_from_params_bag(
+            self,
+            self.quantum,
+            self._quantum_instance if self.quantum else None,
+            self.upper_bound,
+            self.qaoa_optimizer,
+            self.classical_optimizer,
+            self.create_mixer,
+            self.n_reps,
+        )
         set_global_optimizer(self._optimizer)
         return classifier
 
@@ -953,9 +999,15 @@ class NearestConvexHull(BaseEstimator, ClassifierMixin, TransformerMixin):
                 print("Not running in parallel")
 
         if parallel:
-            dists = Parallel(n_jobs=self.n_jobs)(
-                delayed(self._process_sample)(x) for x in X
-            )
+            # Get global optimizer in this process
+            optimizer = get_global_optimizer(default=None)
+
+            def job(x):
+                # Set the global optimizer inside the new process
+                set_global_optimizer(optimizer)
+                return self._process_sample(x)
+
+            dists = Parallel(n_jobs=self.n_jobs)(delayed(job)(x) for x in X)
 
         else:
             for x in X:
@@ -1018,6 +1070,8 @@ class QuanticNCH(QuanticClassifierBase):
     .. versionadded:: 0.2.0
     .. versionchanged:: 0.3.0
         Add qaoa_optimizer parameter.
+    .. versionchanged:: 0.4.0
+        Add QAOA-CV optimization.
 
     Parameters
     ----------
@@ -1059,6 +1113,13 @@ class QuanticNCH(QuanticClassifierBase):
     qaoa_optimizer : SciPyOptimizer, default=SLSQP()
         An instance of a scipy optimizer to find the optimal weights for the
         parametric circuit (ansatz).
+    create_mixer : None | Callable[int, QuantumCircuit], default=None
+        A delegate that takes into input an angle and returns a QuantumCircuit.
+        This circuit is the mixer operatior for the QAOA-CV algorithm.
+        If None and quantum, the NaiveQAOAOptimizer will be used instead.
+    n_reps : int, default=3
+        The number of time the mixer and cost operator are repeated in the QAOA-CV
+        circuit.
 
     References
     ----------
@@ -1081,6 +1142,8 @@ class QuanticNCH(QuanticClassifierBase):
         n_samples_per_hull=10,
         subsampling="min",
         qaoa_optimizer=SLSQP(),
+        create_mixer=None,
+        n_reps=3,
     ):
         QuanticClassifierBase.__init__(
             self, quantum, q_account_token, verbose, shots, None, seed
@@ -1093,6 +1156,8 @@ class QuanticNCH(QuanticClassifierBase):
         self.n_jobs = n_jobs
         self.subsampling = subsampling
         self.qaoa_optimizer = qaoa_optimizer
+        self.create_mixer = create_mixer
+        self.n_reps = n_reps
 
     def _init_algo(self, n_features):
         self._log("Nearest Convex Hull Classifier initiating algorithm")
@@ -1104,16 +1169,16 @@ class QuanticNCH(QuanticClassifierBase):
             subsampling=self.subsampling,
         )
 
-        if self.quantum:
-            self._log("Using NaiveQAOAOptimizer")
-            self._optimizer = NaiveQAOAOptimizer(
-                quantum_instance=self._quantum_instance,
-                upper_bound=self.upper_bound,
-                optimizer=self.qaoa_optimizer,
-            )
-        else:
-            self._log("Using ClassicalOptimizer")
-            self._optimizer = ClassicalOptimizer(self.classical_optimizer)
+        self._optimizer = _get_docplex_optimizer_from_params_bag(
+            self,
+            self.quantum,
+            self._quantum_instance if self.quantum else None,
+            self.upper_bound,
+            self.qaoa_optimizer,
+            self.classical_optimizer,
+            self.create_mixer,
+            self.n_reps,
+        )
 
         # sets the optimizer for the distance functions
         # used in NearestConvexHull class
diff --git a/pyriemann_qiskit/utils/docplex.py b/pyriemann_qiskit/utils/docplex.py
index 10a016a..da1baa8 100644
--- a/pyriemann_qiskit/utils/docplex.py
+++ b/pyriemann_qiskit/utils/docplex.py
@@ -15,13 +15,15 @@ from qiskit.circuit.library import QAOAAnsatz
 from qiskit.primitives import BackendSampler
 from qiskit.quantum_info import Statevector
 from qiskit_algorithms import QAOA
-from qiskit_algorithms.optimizers import SLSQP
+from qiskit_algorithms.optimizers import SLSQP, SPSA
 from qiskit_optimization.algorithms import CobylaOptimizer, MinimumEigenOptimizer
-from qiskit_optimization.converters import IntegerToBinary
+from qiskit_optimization.converters import IntegerToBinary, LinearEqualityToPenalty
 from qiskit_optimization.translators import from_docplex_mp
 from qiskit_optimization.problems import VarType
 from sklearn.preprocessing import MinMaxScaler
 
+from .hyper_params_factory import create_mixer_rotational_X_gates
+from .math import is_pauli_identity
 from .quantum_provider import get_simulator
 
 
@@ -64,8 +66,10 @@ def get_global_optimizer(default):
     return _global_optimizer[0] if _global_optimizer[0] is not None else default
 
 
-def square_cont_mat_var(prob, channels, name="cont_covmat"):
-    """Creates a 2-dimensional dictionary of continuous decision variables,
+def square_cont_mat_var(prob, channels, name="cont_spdmat"):
+    """ Docplex square continous matrix
+
+    Creates a 2-dimensional dictionary of continuous decision variables,
     indexed by pairs of key objects.
     The dictionary represents a square matrix of size
     len(channels) x len(channels).
@@ -104,8 +108,10 @@ def square_cont_mat_var(prob, channels, name="cont_covmat"):
     )
 
 
-def square_int_mat_var(prob, channels, upper_bound=7, name="int_covmat"):
-    """Creates a 2-dimensional dictionary of integer decision variables,
+def square_int_mat_var(prob, channels, upper_bound=7, name="int_spdmat"):
+    """ Docplex square integer matrix
+
+    Creates a 2-dimensional dictionary of integer decision variables,
     indexed by pairs of key objects.
     The dictionary represents a square matrix of size
     len(channels) x len(channels).
@@ -146,8 +152,10 @@ def square_int_mat_var(prob, channels, upper_bound=7, name="int_covmat"):
     )
 
 
-def square_bin_mat_var(prob, channels, name="bin_covmat"):
-    """Creates a 2-dimensional dictionary of binary decision variables,
+def square_bin_mat_var(prob, channels, name="bin_spdmat"):
+    """ Docplex square binary matrix
+
+    Creates a 2-dimensional dictionary of binary decision variables,
     indexed by pairs of key objects.
     The dictionary represents a square matrix of size
     len(channels) x len(channels).
@@ -189,7 +197,7 @@ class pyQiskitOptimizer:
     """Wrapper for Qiskit optimizer.
 
     This class is an abstract class which provides an interface
-    for running our docplex model independently of the optimizer type
+    for running a docplex model independently of the optimizer type
     (such as classical or quantum optimizer).
 
     Notes
@@ -201,114 +209,119 @@ class pyQiskitOptimizer:
     def __init__(self):
         pass
 
-    """Hook to apply some transformation on a covariance matrix.
-
-    Parameters
-    ----------
-    covmat : ndarray, shape (n_features, n_features)
-        The covariance matrix.
-
-    Returns
-    -------
-    transformed_covmat : ndarray, shape (n_features, n_features)
-        A transformation of the covariance matrix.
-
-    Notes
-    -----
-    .. versionadded:: 0.0.2
-    """
-
-    def convert_covmat(self, covmat):
-        return covmat
-
-    """Helper to create a docplex representation of a
-    covariance matrix variable.
-
-    Parameters
-    ----------
-    prob : Model
-        An instance of the docplex model [1]_
-    channels : list
-        The list of channels. A channel can be any Python object,
-        such as channels'name or number but None.
-    name : string
-        An custom name for the variable. The name is used internally by docplex
-        and may appear if your print the model to a file for example.
-
-    Returns
-    -------
-    docplex_covmat : dict
-        A square matrix of decision variables representing
-        our covariance matrix.
-
-    Notes
-    -----
-    .. versionadded:: 0.0.2
-
-    References
-    ----------
-    .. [1] \
-        http://ibmdecisionoptimization.github.io/docplex-doc/mp/_modules/docplex/mp/model.html#Model
-
-    """
-
-    def covmat_var(self, prob, channels, name):
+    def convert_spdmat(self, X):
+        """Convert a SPD matrix
+
+        Hook to apply some transformation on a SPD matrix.
+
+        Parameters
+        ----------
+        X : ndarray, shape (n_features, n_features)
+            A SPD matrix.
+
+        Returns
+        -------
+        X_new : ndarray, shape (n_features, n_features)
+            A transformation of the SPD matrix.
+
+        Notes
+        -----
+        .. versionadded:: 0.0.2
+        .. versionchanged:: 0.4.0
+            rename convert_covmat to convert_spdmat
+        """
+        return X
+
+    def spdmat_var(self, prob, channels, name):
+        """ Create docplex matrix variable
+
+        Helper to create a docplex representation of a
+        SPD matrix variable.
+
+        Parameters
+        ----------
+        prob : Model
+            An instance of the docplex model [1]_.
+        channels : list
+            The list of channels. A channel can be any Python object,
+            such as channels'name or number but None.
+        name : string
+            An custom name for the variable. The name is used internally by docplex
+            and may appear if your print the model to a file for example.
+
+        Returns
+        -------
+        docplex_spdmat : dict
+            A docplex representation of a SPD matrix.
+
+        Notes
+        -----
+        .. versionadded:: 0.0.2
+        .. versionchanged:: 0.4.0
+            rename covmat_var to spdmat_var
+
+        References
+        ----------
+        .. [1] \
+            http://ibmdecisionoptimization.github.io/docplex-doc/mp/_modules/docplex/mp/model.html#Model
+
+        """
         raise NotImplementedError()
 
     def _solve_qp(self, qp, reshape=True):
         raise NotImplementedError()
 
-    """Solve the docplex problem.
-
-    Parameters
-    ----------
-    prob : Model
-        An instance of the docplex model [1]_
+    def solve(self, prob, reshape=True):
+        """Solve the docplex problem.
 
-    Returns
-    -------
-    result : OptimizationResult
-        The result of the optimization.
+        Parameters
+        ----------
+        prob : Model
+            An instance of the docplex model [1]_
 
-    Notes
-    -----
-    .. versionadded:: 0.0.2
-    .. versionchanged:: 0.0.4
+        Returns
+        -------
+        result : OptimizationResult
+            The result of the optimization.
 
-    References
-    ----------
-    .. [1] \
-        http://ibmdecisionoptimization.github.io/docplex-doc/mp/_modules/docplex/mp/model.html#Model
+        Notes
+        -----
+        .. versionadded:: 0.0.2
+        .. versionchanged:: 0.0.4
 
-    """
+        References
+        ----------
+        .. [1] \
+            http://ibmdecisionoptimization.github.io/docplex-doc/mp/_modules/docplex/mp/model.html#Model
 
-    def solve(self, prob, reshape=True):
+        """
         qp = from_docplex_mp(prob)
         return self._solve_qp(qp, reshape)
 
-    """Helper to create a docplex representation of a
-    weight vector.
+    def get_weights(self, prob, classes):
+        """Weights variable
 
-    Parameters
-    ----------
-    prob : Model
-        An instance of the docplex model [1]_
-    classes : list
-        The classes.
+        Helper to create a docplex representation of a
+        weight vector.
 
-    Returns
-    -------
-    docplex_weights : dict
-        A vector of decision variables representing
-        our weights.
+        Parameters
+        ----------
+        prob : Model
+            An instance of the docplex model [1]_
+        classes : list
+            The classes.
 
-    Notes
-    -----
-    .. versionadded:: 0.0.4
+        Returns
+        -------
+        docplex_weights : dict
+            A vector of decision variables representing
+            weights.
 
-    """
+        Notes
+        -----
+        .. versionadded:: 0.0.4
 
-    def get_weights(self, prob, classes):
+        """
         raise NotImplementedError()
 
 
@@ -343,42 +356,43 @@ class ClassicalOptimizer(pyQiskitOptimizer):
         pyQiskitOptimizer.__init__(self)
         self.optimizer = optimizer
 
-    """Helper to create a docplex representation of a
-    covariance matrix variable.
-
-    Parameters
-    ----------
-    prob : Model
-        An instance of the docplex model [1]_
-    channels : list
-        The list of channels. A channel can be any Python object,
-        such as channels'name or number but None.
-    name : string
-        An custom name for the variable. The name is used internally by docplex
-        and may appear if your print the model to a file for example.
-
-    Returns
-    -------
-    docplex_covmat : dict
-        A square matrix of continuous decision variables representing
-        our covariance matrix.
-
-    See Also
-    -----
-    square_cont_mat_var
-
-    Notes
-    -----
-    .. versionadded:: 0.0.2
-
-    References
-    ----------
-    .. [1] \
-        http://ibmdecisionoptimization.github.io/docplex-doc/mp/_modules/docplex/mp/model.html#Model
-
-    """
-
-    def covmat_var(self, prob, channels, name):
+    def spdmat_var(self, prob, channels, name):
+        """ Create docplex matrix variable
+
+        Helper to create a docplex representation of a
+        SPD matrix variable.
+
+        Parameters
+        ----------
+        prob : Model
+            An instance of the docplex model [1]_
+        channels : list
+            The list of channels. A channel can be any Python object,
+            such as channels'name or number but None.
+        name : string
+            An custom name for the variable. The name is used internally by docplex
+            and may appear if your print the model to a file for example.
+
+        Returns
+        -------
+        docplex_spdmat : dict
+            A docplex representation of a SPD matrix with continuous variables.
+
+        See Also
+        -----
+        square_cont_mat_var
+
+        Notes
+        -----
+        .. versionadded:: 0.0.2
+        .. versionchanged:: 0.4.0
+            rename covmat_var to spdmat_var
+
+        References
+        ----------
+        .. [1] \
+            http://ibmdecisionoptimization.github.io/docplex-doc/mp/_modules/docplex/mp/model.html#Model
+        """
         return square_cont_mat_var(prob, channels, name)
 
     def _solve_qp(self, qp, reshape=True):
@@ -388,29 +402,29 @@ class ClassicalOptimizer(pyQiskitOptimizer):
             return np.reshape(result, (n_channels, n_channels))
         return result
 
-    """Helper to create a docplex representation of a
-    weight vector.
+    def get_weights(self, prob, classes):
+        """Weights variabpe
 
-    Parameters
-    ----------
-    prob : Model
-        An instance of the docplex model [1]_
-    classes : list
-        The classes.
+        Helper to create a docplex representation of a
+        weight vector.
 
-    Returns
-    -------
-    docplex_weights : dict
-        A vector of continuous decision variables representing
-        our weights.
+        Parameters
+        ----------
+        prob : Model
+            An instance of the docplex model [1]_
+        classes : list
+            The classes.
 
-    Notes
-    -----
-    .. versionadded:: 0.0.4
+        Returns
+        -------
+        docplex_weights : dict
+            A vector of continuous decision variables representing weights.
 
-    """
+        Notes
+        -----
+        .. versionadded:: 0.0.4
 
-    def get_weights(self, prob, classes):
+        """
         w = prob.continuous_var_matrix(
             keys1=[1], keys2=classes, name="weight", lb=0, ub=1
         )
@@ -418,6 +432,20 @@ class ClassicalOptimizer(pyQiskitOptimizer):
         return w
 
 
+def _get_quantum_instance(self):
+    if self.quantum_instance is None:
+        backend = get_simulator()
+        seed = 42
+        shots = 1024
+        quantum_instance = BackendSampler(
+            backend, options={"shots": shots, "seed_simulator": seed}
+        )
+        quantum_instance.transpile_options["seed_transpiler"] = seed
+    else:
+        quantum_instance = self.quantum_instance
+    return quantum_instance
+
+
 class NaiveQAOAOptimizer(pyQiskitOptimizer):
 
     """Wrapper for the quantum optimizer QAOA.
@@ -458,83 +486,78 @@ class NaiveQAOAOptimizer(pyQiskitOptimizer):
         self.quantum_instance = quantum_instance
         self.optimizer = optimizer
 
-    """Transform all values in the covariance matrix
-    to integers.
+    def convert_spdmat(self, X):
+        """Convert a SPD matrix
 
-    Example:
-    0.123 -> 1230
+        Transform all values in the SPD matrix to integers.
 
-    Parameters
-    ----------
-    covmat : ndarray, shape (n_features, n_features)
-        The covariance matrix.
+        Example:
+        0.123 -> 1230
 
-    Returns
-    -------
-    transformed_covmat : ndarray, shape (n_features, n_features)
-        A transformation of the covariance matrix.
+        Parameters
+        ----------
+        X : ndarray, shape (n_features, n_features)
+            The SPD matrix.
 
-    Notes
-    -----
-    .. versionadded:: 0.0.2
+        Returns
+        -------
+        transformed_X : ndarray, shape (n_features, n_features)
+            A transformation of the SPD matrix.
 
-    """
+        Notes
+        -----
+        .. versionadded:: 0.0.2
+        .. versionchanged:: 0.4.0
+            rename convert_covmat to convert_spdmat
 
-    def convert_covmat(self, covmat):
-        corr = normalize(covmat, "corr")
+        """
+        corr = normalize(X, "corr")
         return np.round(corr * self.upper_bound, 0)
 
-    """Helper to create a docplex representation of a
-    covariance matrix variable.
-
-    Parameters
-    ----------
-    prob : Model
-        An instance of the docplex model [1]_
-    channels : list
-        The list of channels. A channel can be any Python object,
-        such as channels'name or number but None.
-    name : string
-        An custom name for the variable. The name is used internally by docplex
-        and may appear if your print the model to a file for example.
-
-    Returns
-    -------
-    docplex_covmat : dict
-        A square matrix of integer decision variables representing
-        our covariance matrix.
-
-    See Also
-    -----
-    square_int_mat_var
-
-    Notes
-    -----
-    .. versionadded:: 0.0.2
-
-    References
-    ----------
-    .. [1] \
-        http://ibmdecisionoptimization.github.io/docplex-doc/mp/_modules/docplex/mp/model.html#Model
-
-    """
-
-    def covmat_var(self, prob, channels, name):
+    def spdmat_var(self, prob, channels, name):
+        """ Create docplex matrix variable
+
+        Helper to create a docplex representation of a
+        SPD matrix variable.
+
+        Parameters
+        ----------
+        prob : Model
+            An instance of the docplex model [1]_
+        channels : list
+            The list of channels. A channel can be any Python object,
+            such as channels'name or number but None.
+        name : string
+            An custom name for the variable. The name is used internally by docplex
+            and may appear if your print the model to a file for example.
+
+        Returns
+        -------
+        docplex_spdmat : dict
+            A docplex representation of a SPD matrix with integer variables.
+
+        See Also
+        -----
+        square_int_mat_var
+
+        Notes
+        -----
+        .. versionadded:: 0.0.2
+        .. versionchanged:: 0.4.0
+            rename covmat_var to spdmat_var
+
+        References
+        ----------
+        .. [1] \
+            http://ibmdecisionoptimization.github.io/docplex-doc/mp/_modules/docplex/mp/model.html#Model
+
+        """
         return square_int_mat_var(prob, channels, self.upper_bound, name)
 
     def _solve_qp(self, qp, reshape=True):
         conv = IntegerToBinary()
         qubo = conv.convert(qp)
-        if self.quantum_instance is None:
-            backend = get_simulator()
-            seed = 42
-            shots = 1024
-            quantum_instance = BackendSampler(
-                backend, options={"shots": shots, "seed_simulator": seed}
-            )
-            quantum_instance.transpile_options["seed_transpiler"] = seed
-        else:
-            quantum_instance = self.quantum_instance
+        quantum_instance = _get_quantum_instance(self)
 
         self.evaluated_values_ = []
 
@@ -554,29 +577,28 @@ class NaiveQAOAOptimizer(pyQiskitOptimizer):
             return np.reshape(result, (n_channels, n_channels))
         return result
 
-    """Helper to create a docplex representation of a
-    weight vector.
+    def get_weights(self, prob, classes):
+        """Get weights variable
 
-    Parameters
-    ----------
-    prob : Model
-        An instance of the docplex model [1]_
-    classes : list
-        The classes.
+        Helper to create a docplex representation of a weight vector.
 
-    Returns
-    -------
-    docplex_weights : dict
-        A vector of integer decision variables representing
-        our weights.
+        Parameters
+        ----------
+        prob : Model
+            An instance of the docplex model [1]_
+        classes : list
+            The classes.
 
-    Notes
-    -----
-    .. versionadded:: 0.0.4
+        Returns
+        -------
+        docplex_weights : dict
+            A vector of integer decision variables representing weights.
 
-    """
+        Notes
+        -----
+        .. versionadded:: 0.0.4
 
-    def get_weights(self, prob, classes):
+        """
         w = prob.integer_var_matrix(
             keys1=[1], keys2=classes, name="weight", lb=0, ub=self.upper_bound
         )
@@ -590,16 +612,17 @@ class QAOACVOptimizer(pyQiskitOptimizer):
 
     Parameters
     ----------
-    create_mixer : Callable[int,QuantumCircuit]
-        A delegate that takes into input an angle and returns a QuantumCircuit
-    n_reps : int
+    create_mixer : Callable[int, QuantumCircuit], \
+            default=create_mixer_rotational_X_gates(0)
+        A delegate that takes into input an angle and returns a QuantumCircuit.
+    n_reps : int, default=3
         The number of repetitions for the QAOA ansatz.
         It defines how many time Mixer and Cost operators will be repeated
         in the circuit.
     quantum_instance : QuantumInstance, default=None
         A quantum backend instance.
         If None, AerSimulator will be used.
-    optimizer : SciPyOptimizer, default=SLSQP()
+    optimizer : SciPyOptimizer, default=SPSA()
         An instance of a scipy optimizer to find the optimal weights for the
         parametric circuit (ansatz).
 
@@ -626,9 +649,16 @@ class QAOACVOptimizer(pyQiskitOptimizer):
     See Also
     --------
     pyQiskitOptimizer
+    create_mixer_rotational_X_gates
     """
 
-    def __init__(self, create_mixer, n_reps, quantum_instance=None, optimizer=SLSQP()):
+    def __init__(
+        self,
+        create_mixer=create_mixer_rotational_X_gates(0),
+        n_reps=3,
+        quantum_instance=None,
+        optimizer=SPSA(),
+    ):
         self.n_reps = n_reps
         self.create_mixer = create_mixer
         self.quantum_instance = quantum_instance
@@ -645,9 +675,77 @@ class QAOACVOptimizer(pyQiskitOptimizer):
                 # print(scaler.data_min_, scaler.data_max_)
                 scalers.append(scaler)
                 v.vartype = VarType.BINARY
-        return scalers
+                v.lowerbound = 0
+                v.upperbound = 1
+        conv = LinearEqualityToPenalty()
+        qp = conv.convert(qp)
+
+        return qp, scalers
+
+    def spdmat_var(self, prob, channels, name):
+        """ Create docplex matrix variable
+
+        Parameters
+        ----------
+        prob : Model
+            An instance of the docplex model [1]_.
+        channels : list
+            The list of channels. A channel can be any Python object,
+            such as channels'name or number but None.
+        name : string
+            An custom name for the variable. The name is used internally by docplex
+            and may appear if your print the model to a file for example.
+
+        Returns
+        -------
+        docplex_spdmat : dict
+            A docplex representation of a SPD matrix with continuous variables.
+
+        See Also
+        -----
+        square_cont_mat_var
+
+        Notes
+        -----
+        .. versionadded:: 0.4.0
+        .. versionchanged:: 0.4.0
+            rename covmat_var to spdmat_var
+
+        References
+        ----------
+        .. [1] \
+            http://ibmdecisionoptimization.github.io/docplex-doc/mp/_modules/docplex/mp/model.html#Model
+
+        """
+        return ClassicalOptimizer.spdmat_var(self, prob, channels, name)
+
+    def get_weights(self, prob, classes):
+        """Get weights variable
+
+        Helper to create a docplex representation of a weight vector.
+
+        Parameters
+        ----------
+        prob : Model
+            An instance of the docplex model [1]_
+        classes : list
+            The classes.
+
+        Returns
+        -------
+        docplex_weights : dict
+            A vector of integer decision variables representing weights.
+
+        Notes
+        -----
+        .. versionadded:: 0.4.0
+
+        """
+        return ClassicalOptimizer.get_weights(self, prob, classes)
 
     def _solve_qp(self, qp, reshape=True):
+        quantum_instance = _get_quantum_instance(self)
+
         n_var = qp.get_num_vars()
         # Extract the objective function from the docplex model
         # We want the object expression with continuous variable
@@ -655,7 +753,7 @@ class QAOACVOptimizer(pyQiskitOptimizer):
 
         # Convert continous variable to binary ones
         # Get scalers corresponding to the definition range of each variables
-        scalers = QAOACVOptimizer.prepare_model(qp)
+        qp, scalers = QAOACVOptimizer.prepare_model(qp)
 
         # Check all variables are converted to binary, and scalers are registered
         # print(qp.prettyprint(), scalers)
@@ -664,7 +762,14 @@ class QAOACVOptimizer(pyQiskitOptimizer):
         # Get operator associated with model
         cost, _offset = qp.to_ising()
 
-        mixer = self.create_mixer(cost.num_qubits)
+        # If the cost operator is a Pauli identity
+        # or the cost operator has no parameters
+        # the number of parameters in the QAOAAnsatz will be 0.
+        # We will then create a mixer with parameters
+        # So we get some parameters in the circuit to optimize
+        cost_op_has_no_parameter = is_pauli_identity(cost) or len(cost.parameters) == 0
+
+        mixer = self.create_mixer(cost.num_qubits, use_params=cost_op_has_no_parameter)
 
         # QAOA cirtcuit
         ansatz_0 = QAOAAnsatz(
@@ -690,7 +795,7 @@ class QAOACVOptimizer(pyQiskitOptimizer):
 
             # p is in the range [0, 1].
             # We now need to scale it in the definition
-            # range of our continuous variables
+            # range of the continuous variables
             p = scalers[i].inverse_transform([[p]])[0][0]
             return p
 
@@ -699,14 +804,14 @@ class QAOACVOptimizer(pyQiskitOptimizer):
         self.y_ = []
 
         def loss(params):
-            job = self.quantum_instance.run(ansatz, params)
+            job = quantum_instance.run(ansatz, params)
             var_hat = [prob(job, i) for i in range(n_var)]
             cost = objective_expr.evaluate(var_hat)
             self.x_.append(len(self.x_))
             self.y_.append(cost)
             return cost
 
-        # Initial guess for the parameters
+        # Initial guess for the parameters.
         initial_guess = np.array([1, 1] * self.n_reps)
 
         # minimize function to search for the optimal parameters
@@ -718,11 +823,15 @@ class QAOACVOptimizer(pyQiskitOptimizer):
         self.optim_params_ = result.x
 
         # running QAOA circuit with optimal parameters
-        job = self.quantum_instance.run(ansatz, self.optim_params_)
-        solution = [prob(job, i) for i in range(n_var)]
+        job = quantum_instance.run(ansatz, self.optim_params_)
+        solution = np.array([prob(job, i) for i in range(n_var)])
         self.minimum_ = objective_expr.evaluate(solution)
 
         optimized_circuit = ansatz_0.assign_parameters(self.optim_params_)
         self.state_vector_ = Statevector(optimized_circuit)
 
+        if reshape:
+            n_channels = int(math.sqrt(solution.shape[0]))
+            return np.reshape(solution, (n_channels, n_channels))
+
         return solution
diff --git a/pyriemann_qiskit/utils/hyper_params_factory.py b/pyriemann_qiskit/utils/hyper_params_factory.py
index a07f09d..ec5140a 100644
--- a/pyriemann_qiskit/utils/hyper_params_factory.py
+++ b/pyriemann_qiskit/utils/hyper_params_factory.py
@@ -1,7 +1,6 @@
 import inspect
-import math
 
-from qiskit.circuit import QuantumCircuit, QuantumRegister
+from qiskit.circuit import QuantumCircuit, QuantumRegister, Parameter
 from qiskit.circuit.library import (
     ZZFeatureMap,
     ZFeatureMap,
@@ -344,6 +343,44 @@ def get_spsa_parameters(spsa):
     )
 
 
+def create_mixer_qiskit_default(_angle):
+    r"""Return the default mixing operator with QAOA.
+    (qiskit implementation)
+
+    .. math::
+        H_X = \sum_{i}^{N} X_i
+
+    See [1]_ for details.
+
+    Parameters
+    ----------
+    angle : float
+        The angle of the gates' rotation.
+        Not used. Just kept for compatibility with other
+        `create_mixer`
+
+    Returns
+    -------
+    mixer : Callable[[int, boolean], QuantumCircuit]
+        A method that returns None,
+        forcing qiskit to use its own implementation of the mixer.
+
+    Notes
+    -----
+    .. versionadded:: 0.4.0
+
+    References
+    ----------
+    .. [1] \
+        https://dice.cyfronet.pl/papers/JPlewa_JSienko_msc_v2.pdf
+    """
+
+    def mixer_qiskit_default(n_qubits, use_params=False):
+        return None
+
+    return mixer_qiskit_default
+
+
 def create_mixer_rotational_X_gates(angle):
     r"""Return the default mixing operator with QAOA.
 
@@ -359,8 +396,11 @@ def create_mixer_rotational_X_gates(angle):
 
     Returns
     -------
-    mixer : QuantumCircuit
-        The mixer.
+    mixer : Callable[[int, boolean], QuantumCircuit]
+        A method that create the mixer with the correct angle.
+        This method takes into parameters the number of qubits in the circuit,
+        and a parameter `use_params`.
+        if True, `use_params` create a Parameter in the circuit (same for all gates)
 
     Notes
     -----
@@ -372,12 +412,16 @@ def create_mixer_rotational_X_gates(angle):
         https://dice.cyfronet.pl/papers/JPlewa_JSienko_msc_v2.pdf
     """
 
-    def mixer_X(n_qubits):
+    def mixer_X(n_qubits, use_params=False):
         qr = QuantumRegister(n_qubits)
         mixer = QuantumCircuit(qr)
-
+        p = Parameter("beta")
         for qr_ in qr:
-            mixer.rx(angle, qr_)
+            if use_params:
+                mixer.rx(p + angle, qr_)
+            else:
+                mixer.rx(angle, qr_)
+        return mixer
 
     return mixer_X
 
@@ -397,8 +441,11 @@ def create_mixer_rotational_XY_gates(angle):
 
     Returns
     -------
-    mixer : QuantumCircuit
-        The mixer.
+    mixer : Callable[[int, boolean], QuantumCircuit]
+        A method that create the mixer with the correct angle.
+        This method takes into parameters the number of qubits in the circuit,
+        and a parameter `use_params`.
+        if True, `use_params` create a Parameter in the circuit (same for all gates)
 
     Notes
     -----
@@ -410,15 +457,17 @@ def create_mixer_rotational_XY_gates(angle):
         https://dice.cyfronet.pl/papers/JPlewa_JSienko_msc_v2.pdf
     """
 
-    def mixer_XY(n_qubits):
+    def mixer_XY(n_qubits, use_params):
         qr = QuantumRegister(n_qubits)
         mixer = QuantumCircuit(qr)
-
+        p = Parameter("beta")
         for i in range(n_qubits - 1):
-            mixer.rx(angle, qr[i])
-            mixer.rx(angle, qr[i + 1])
-            mixer.ry(angle, qr[i])
-            mixer.ry(angle, qr[i + 1])
+            if use_params:
+                mixer.rxx(p + angle, qr[i], qr[i + 1])
+            else:
+                mixer.rxx(angle, qr[i], qr[i + 1])
+            mixer.ryy(0 + angle, qr[i], qr[i + 1])
+        return mixer
 
     return mixer_XY
 
@@ -438,8 +487,11 @@ def create_mixer_rotational_XZ_gates(angle):
 
     Returns
     -------
-    mixer : QuantumCircuit
-        The mixer.
+    mixer : Callable[[int, boolean], QuantumCircuit]
+        A method that create the mixer with the correct angle.
+        This method takes into parameters the number of qubits in the circuit,
+        and a parameter `use_params`.
+        if True, `use_params` create a Parameter in the circuit (same for all gates)
 
     Notes
     -----
@@ -451,14 +503,17 @@ def create_mixer_rotational_XZ_gates(angle):
         https://dice.cyfronet.pl/papers/JPlewa_JSienko_msc_v2.pdf
     """
 
-    def mixer_XZ(n_qubits):
+    def mixer_XZ(n_qubits, use_params):
         qr = QuantumRegister(n_qubits)
         mixer = QuantumCircuit(qr)
 
-        for i in range(1, n_qubits - 1):
-            mixer.rz(angle, qr[i - 1])
-            mixer.rx(angle, qr[i])
-            mixer.rx(angle + math.pi / 2, qr[i])
-            mixer.rz(angle, qr[i + 1])
+        p = Parameter("beta")
+        for i in range(n_qubits - 1):
+            if use_params:
+                mixer.rzx(p + angle, qr[i], qr[i + 1])
+            else:
+                mixer.rzx(angle, qr[i], qr[i + 1])
+            mixer.rzx(0 + angle, qr[i], qr[i + 1]).inverse(annotated=True)
+        return mixer
 
     return mixer_XZ
diff --git a/pyriemann_qiskit/utils/math.py b/pyriemann_qiskit/utils/math.py
index 3929676..5a63de0 100644
--- a/pyriemann_qiskit/utils/math.py
+++ b/pyriemann_qiskit/utils/math.py
@@ -3,6 +3,7 @@ from typing_extensions import deprecated
 
 import numpy as np
 from pyriemann.utils.covariance import normalize
+from qiskit.quantum_info import Pauli, SparsePauliOp
 
 
 @deprecated(
@@ -94,3 +95,37 @@ def to_xyz(X):
     if X.shape[1] != 2 and X.shape[2] != 2:
         raise ValueError("SPD matrices must have size 2 x 2")
     return np.array([[spd[0, 0], spd[0, 1], spd[1, 1]] for spd in X])
+
+
+def is_pauli_identity(operator):
+    """Check if operator is pauli identity
+
+    Implementation from qiskit [1]_.
+
+    Parameters
+    ----------
+    operator : Union[QuantumCircuit,Pauli,SparsePauliOp]
+        An operator.
+
+    Returns
+    -------
+    is_pauli : boolean
+        True if the operator is a pauli identity
+
+    Notes
+    -----
+    .. versionadded:: 0.4.0
+
+    References
+    ----------
+    .. [1] \
+        https://github.com/Qiskit/qiskit/blob/stable/1.2/qiskit/circuit/library/n_local/evolved_operator_ansatz.py#L248
+    """
+    if isinstance(operator, SparsePauliOp):
+        if len(operator.paulis) == 1:
+            operator = operator.paulis[0]  # check if the single Pauli is identity below
+        else:
+            return False
+    if isinstance(operator, Pauli):
+        return not np.any(np.logical_or(operator.x, operator.z))
+    return False
diff --git a/pyriemann_qiskit/utils/mean.py b/pyriemann_qiskit/utils/mean.py
index a91e056..207269b 100644
--- a/pyriemann_qiskit/utils/mean.py
+++ b/pyriemann_qiskit/utils/mean.py
@@ -60,10 +60,10 @@ def qmean_euclid(X, sample_weight=None, optimizer=ClassicalOptimizer()):
 
     prob = Model()
 
-    X_mean = optimizer.covmat_var(prob, channels, "fro_mean")
+    X_mean = optimizer.spdmat_var(prob, channels, "fro_mean")
 
     def _dist_euclid(A, B):
-        A = optimizer.convert_covmat(A)
+        A = optimizer.convert_spdmat(A)
         return prob.sum_squares(A[r, c] - B[r, c] for r in channels for c in channels)
 
     objectives = prob.sum(_dist_euclid(X[i], X_mean) for i in matrices)


Test patch:
diff --git a/tests/test_docplex.py b/tests/test_docplex.py
index c75ae37..d7f1e54 100644
--- a/tests/test_docplex.py
+++ b/tests/test_docplex.py
@@ -7,6 +7,7 @@ from pyriemann_qiskit.utils.docplex import (
     square_bin_mat_var,
     ClassicalOptimizer,
     NaiveQAOAOptimizer,
+    QAOACVOptimizer,
 )
 
 
@@ -38,6 +39,8 @@ def test_get_square_cont_var(square_mat_var):
     assert type(first_element.vartype) is expected_result_type
 
 
-@pytest.mark.parametrize("optimizer", [ClassicalOptimizer, NaiveQAOAOptimizer])
+@pytest.mark.parametrize(
+    "optimizer", [ClassicalOptimizer, NaiveQAOAOptimizer, QAOACVOptimizer]
+)
 def test_optimizer_creation(optimizer):
     assert optimizer()
diff --git a/tests/test_utils_distance.py b/tests/test_utils_distance.py
index 530d6d4..13bfc6c 100644
--- a/tests/test_utils_distance.py
+++ b/tests/test_utils_distance.py
@@ -3,6 +3,7 @@ import numpy as np
 from pyriemann_qiskit.utils.docplex import (
     ClassicalOptimizer,
     NaiveQAOAOptimizer,
+    QAOACVOptimizer,
 )
 from pyriemann_qiskit.utils.distance import (
     qdistance_logeuclid_to_convex_hull,
@@ -32,7 +33,9 @@ def test_performance(metric):
     assert score.mean() > 0
 
 
-@pytest.mark.parametrize("optimizer", [ClassicalOptimizer(), NaiveQAOAOptimizer()])
+@pytest.mark.parametrize(
+    "optimizer", [ClassicalOptimizer(), NaiveQAOAOptimizer(), QAOACVOptimizer()]
+)
 def test_qdistance_logeuclid_to_convex_hull(optimizer, get_covmats):
     n_trials, n_channels = 5, 3
     covmats = get_covmats(n_trials, n_channels)
@@ -45,7 +48,9 @@ def test_qdistance_logeuclid_to_convex_hull(optimizer, get_covmats):
     assert dist == pytest.approx(0, rel=1e-5, abs=1e-5)
 
 
-@pytest.mark.parametrize("optimizer", [ClassicalOptimizer(), NaiveQAOAOptimizer()])
+@pytest.mark.parametrize(
+    "optimizer", [ClassicalOptimizer(), NaiveQAOAOptimizer(), QAOACVOptimizer()]
+)
 def test_weight_logeuclid_to_convex_hull(optimizer):
     X_0 = np.array([[0.9, 1.1], [0.9, 1.1]])
     X_1 = X_0 + 1
diff --git a/tests/test_utils_mean.py b/tests/test_utils_mean.py
index 00c9b99..b717f89 100644
--- a/tests/test_utils_mean.py
+++ b/tests/test_utils_mean.py
@@ -5,7 +5,11 @@ from pyriemann.estimation import XdawnCovariances, Shrinkage
 from sklearn.pipeline import make_pipeline
 from sklearn.model_selection import StratifiedKFold, cross_val_score
 from pyriemann_qiskit.utils.mean import qmean_euclid, qmean_logeuclid
-from pyriemann_qiskit.utils.docplex import ClassicalOptimizer, NaiveQAOAOptimizer
+from pyriemann_qiskit.utils.docplex import (
+    ClassicalOptimizer,
+    NaiveQAOAOptimizer,
+    QAOACVOptimizer,
+)
 from pyriemann_qiskit.classification import QuanticMDM
 from pyriemann_qiskit.datasets import get_mne_sample
 from qiskit_optimization.algorithms import ADMMOptimizer
@@ -56,7 +60,11 @@ def test_mean_cpm_shape(get_covmats, mean):
 
 @pytest.mark.parametrize(
     "optimizer",
-    [ClassicalOptimizer(optimizer=ADMMOptimizer()), NaiveQAOAOptimizer()],
+    [
+        ClassicalOptimizer(optimizer=ADMMOptimizer()),
+        NaiveQAOAOptimizer(),
+        QAOACVOptimizer(),
+    ],
 )
 @pytest.mark.parametrize("mean", [qmean_euclid])
 def test_mean_cpm_all_zeros(optimizer, mean):
@@ -70,7 +78,11 @@ def test_mean_cpm_all_zeros(optimizer, mean):
 
 @pytest.mark.parametrize(
     "optimizer",
-    [ClassicalOptimizer(optimizer=ADMMOptimizer()), NaiveQAOAOptimizer()],
+    [
+        ClassicalOptimizer(optimizer=ADMMOptimizer()),
+        NaiveQAOAOptimizer(),
+        QAOACVOptimizer(),
+    ],
 )
 @pytest.mark.parametrize("mean", [qmean_euclid])
 def test_mean_cpm_all_ones(optimizer, mean):
@@ -84,7 +96,11 @@ def test_mean_cpm_all_ones(optimizer, mean):
 
 @pytest.mark.parametrize(
     "optimizer",
-    [ClassicalOptimizer(optimizer=ADMMOptimizer()), NaiveQAOAOptimizer()],
+    [
+        ClassicalOptimizer(optimizer=ADMMOptimizer()),
+        NaiveQAOAOptimizer(),
+        QAOACVOptimizer(),
+    ],
 )
 @pytest.mark.parametrize("mean", [qmean_euclid])
 def test_mean_cpm_all_equals(optimizer, mean):
@@ -98,7 +114,11 @@ def test_mean_cpm_all_equals(optimizer, mean):
 
 @pytest.mark.parametrize(
     "optimizer",
-    [ClassicalOptimizer(optimizer=ADMMOptimizer()), NaiveQAOAOptimizer()],
+    [
+        ClassicalOptimizer(optimizer=ADMMOptimizer()),
+        NaiveQAOAOptimizer(),
+        QAOACVOptimizer(),
+    ],
 )
 @pytest.mark.parametrize("mean", [qmean_euclid])
 def test_mean_cpm_mixed(optimizer, mean):