diff --git a/Dockerfile b/Dockerfile
index a2501d96..464f0abf 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -24,16 +24,16 @@ RUN mkdir /root/mne_data
 RUN mkdir /home/mne_data
 
 ## Workaround for firestore
-RUN pip install protobuf==4.25.2
-RUN pip install google_cloud_firestore==2.14.0
+RUN pip install protobuf==4.25.3
+RUN pip install google_cloud_firestore==2.15.0
 ### Missing __init__ file in protobuf
-RUN touch /usr/local/lib/python3.9/site-packages/protobuf-4.25.2-py3.9.egg/google/__init__.py
+RUN touch /usr/local/lib/python3.9/site-packages/protobuf-4.25.3-py3.9.egg/google/__init__.py
 ## google.cloud.location is never used in these files, and is missing in path.
-RUN sed -i 's/from google.cloud.location import locations_pb2//g' '/usr/local/lib/python3.9/site-packages/google_cloud_firestore-2.14.0-py3.9.egg/google/cloud/firestore_v1/services/firestore/client.py'
-RUN sed -i 's/from google.cloud.location import locations_pb2//g' '/usr/local/lib/python3.9/site-packages/google_cloud_firestore-2.14.0-py3.9.egg/google/cloud/firestore_v1/services/firestore/transports/base.py'
-RUN sed -i 's/from google.cloud.location import locations_pb2//g' '/usr/local/lib/python3.9/site-packages/google_cloud_firestore-2.14.0-py3.9.egg/google/cloud/firestore_v1/services/firestore/transports/grpc.py'
-RUN sed -i 's/from google.cloud.location import locations_pb2//g' '/usr/local/lib/python3.9/site-packages/google_cloud_firestore-2.14.0-py3.9.egg/google/cloud/firestore_v1/services/firestore/transports/grpc_asyncio.py'
-RUN sed -i 's/from google.cloud.location import locations_pb2//g' '/usr/local/lib/python3.9/site-packages/google_cloud_firestore-2.14.0-py3.9.egg/google/cloud/firestore_v1/services/firestore/transports/rest.py'
-RUN sed -i 's/from google.cloud.location import locations_pb2//g' '/usr/local/lib/python3.9/site-packages/google_cloud_firestore-2.14.0-py3.9.egg/google/cloud/firestore_v1/services/firestore/async_client.py'
+RUN sed -i 's/from google.cloud.location import locations_pb2//g' '/usr/local/lib/python3.9/site-packages/google_cloud_firestore-2.15.0-py3.9.egg/google/cloud/firestore_v1/services/firestore/client.py'
+RUN sed -i 's/from google.cloud.location import locations_pb2//g' '/usr/local/lib/python3.9/site-packages/google_cloud_firestore-2.15.0-py3.9.egg/google/cloud/firestore_v1/services/firestore/transports/base.py'
+RUN sed -i 's/from google.cloud.location import locations_pb2//g' '/usr/local/lib/python3.9/site-packages/google_cloud_firestore-2.15.0-py3.9.egg/google/cloud/firestore_v1/services/firestore/transports/grpc.py'
+RUN sed -i 's/from google.cloud.location import locations_pb2//g' '/usr/local/lib/python3.9/site-packages/google_cloud_firestore-2.15.0-py3.9.egg/google/cloud/firestore_v1/services/firestore/transports/grpc_asyncio.py'
+RUN sed -i 's/from google.cloud.location import locations_pb2//g' '/usr/local/lib/python3.9/site-packages/google_cloud_firestore-2.15.0-py3.9.egg/google/cloud/firestore_v1/services/firestore/transports/rest.py'
+RUN sed -i 's/from google.cloud.location import locations_pb2//g' '/usr/local/lib/python3.9/site-packages/google_cloud_firestore-2.15.0-py3.9.egg/google/cloud/firestore_v1/services/firestore/async_client.py'
 
 ENTRYPOINT [ "python", "/examples/ERP/classify_P300_bi.py" ]
diff --git a/benchmarks/light_benchmark.py b/benchmarks/light_benchmark.py
index 9f2a136b..9ad973e8 100644
--- a/benchmarks/light_benchmark.py
+++ b/benchmarks/light_benchmark.py
@@ -11,7 +11,7 @@
 # Modified from plot_classify_P300_bi.py of pyRiemann
 # License: BSD (3-clause)
 
-from pyriemann.estimation import XdawnCovariances
+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
@@ -78,15 +78,17 @@
 )
 
 pipelines["RG_VQC"] = QuantumClassifierWithDefaultRiemannianPipeline(
-    shots=100, spsa_trials=5, two_local_reps=2, params={"seed": 42}
+    shots=100, spsa_trials=1, two_local_reps=2, params={"seed": 42}
 )
 
 pipelines["QMDM_mean"] = QuantumMDMWithRiemannianPipeline(
-    convex_metric="mean", quantum=True
+    metric={"mean": "euclid_cpm", "distance": "euclid"},
+    quantum=True,
+    regularization=Shrinkage(shrinkage=0.9),
 )
 
 pipelines["QMDM_dist"] = QuantumMDMWithRiemannianPipeline(
-    convex_metric="distance", quantum=True
+    metric={"mean": "logeuclid", "distance": "logeuclid_cpm"}, quantum=True
 )
 
 pipelines["RG_LDA"] = make_pipeline(
diff --git a/doc/api.rst b/doc/api.rst
index 71a05d71..06e1d570 100644
--- a/doc/api.rst
+++ b/doc/api.rst
@@ -96,7 +96,8 @@ Mean
 .. autosummary::
     :toctree: generated/
 
-    fro_mean_convex
+    mean_euclid_cpm
+    mean_logeuclid_cpm
 
 Distance
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -106,7 +107,7 @@ Distance
 .. autosummary::
     :toctree: generated/
 
-    logeucl_dist_convex
+    distance_logeuclid_cpm
 
 Docplex
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~
diff --git a/doc/requirements.txt b/doc/requirements.txt
index 5639a962..7545fe2c 100644
--- a/doc/requirements.txt
+++ b/doc/requirements.txt
@@ -10,6 +10,7 @@ imbalanced-learn==0.11.0
 joblib
 pandas
 cvxpy==1.4.1
+qiskit==0.45.0
 qiskit_machine_learning==0.6.1
 qiskit-ibm-provider==0.7.3
 qiskit-optimization==0.5.0
diff --git a/examples/ERP/classify_P300_bi_illiteracy.py b/examples/ERP/classify_P300_bi_illiteracy.py
index 368ad6a6..5bde4378 100644
--- a/examples/ERP/classify_P300_bi_illiteracy.py
+++ b/examples/ERP/classify_P300_bi_illiteracy.py
@@ -172,19 +172,19 @@ def placeholder(key):
 placeholder(PIP.xDAWN_LDA.value)
 
 pipelines[PIP.ERPCov_CvxMDM_Dist.value] = QuantumMDMWithRiemannianPipeline(
-    convex_metric="distance", quantum=False
+    metric="distance", quantum=False
 )
 placeholder(PIP.ERPCov_CvxMDM_Dist.value)
 
 # Quantum Pipelines
 
 pipelines[PIP.ERPCov_QMDM_Dist.value] = QuantumMDMWithRiemannianPipeline(
-    convex_metric="distance", quantum=True
+    metric="distance", quantum=True
 )
 placeholder(PIP.ERPCov_QMDM_Dist.value)
 
 pipelines[PIP.ERPCov_QMDM_Dist.value] = QuantumMDMWithRiemannianPipeline(
-    convex_metric="distance", quantum=True
+    metric="distance", quantum=True
 )
 placeholder(PIP.ERPCov_QMDM_Dist.value)
 
diff --git a/examples/ERP/classify_P300_bi_quantum_mdm.py b/examples/ERP/classify_P300_bi_quantum_mdm.py
index 60a50260..8ecba8cc 100644
--- a/examples/ERP/classify_P300_bi_quantum_mdm.py
+++ b/examples/ERP/classify_P300_bi_quantum_mdm.py
@@ -6,8 +6,8 @@
 The mean and the distance in MDM algorithm are formulated as
 optimization problems. These optimization problems are translated
 to Qiskit using Docplex and additional glue code. These optimizations
-are enabled when we use convex mean or convex distance. This is set
-using the 'convex_metric' parameter of the QuantumMDMWithRiemannianPipeline.
+are enabled when we use cpm mean or cpm distance. This is set
+using the 'metric' parameter of the QuantumMDMWithRiemannianPipeline.
 
 Classification can be run either on emulation or real quantum computer.
 
@@ -43,7 +43,7 @@
 from moabb.evaluations import WithinSessionEvaluation
 from moabb.paradigms import P300
 
-# inject convex distance and mean to pyriemann (if not done already)
+# inject cpm distance and mean to pyriemann (if not done already)
 from pyriemann_qiskit.utils import distance, mean  # noqa
 from pyriemann_qiskit.pipelines import (
     QuantumMDMVotingClassifier,
@@ -107,15 +107,15 @@
 
 pipelines = {}
 
-pipelines["mean=convex/distance=euclid"] = QuantumMDMWithRiemannianPipeline(
-    convex_metric="mean", quantum=quantum
+pipelines["mean=logeuclid_cpm/distance=logeuclid"] = QuantumMDMWithRiemannianPipeline(
+    metric="mean", quantum=quantum
 )
 
-pipelines["mean=logeuclid/distance=convex"] = QuantumMDMWithRiemannianPipeline(
-    convex_metric="distance", quantum=quantum
+pipelines["mean=logeuclid/distance=logeuclid_cpm"] = QuantumMDMWithRiemannianPipeline(
+    metric="distance", quantum=quantum
 )
 
-pipelines["Voting convex"] = QuantumMDMVotingClassifier(quantum=quantum)
+pipelines["Voting logeuclid_cpm"] = QuantumMDMVotingClassifier(quantum=quantum)
 
 ##############################################################################
 # Run evaluation
diff --git a/examples/MI/classify_alexmi_with_quantum_pipeline.py b/examples/MI/classify_alexmi_with_quantum_pipeline.py
index dbcb60c9..ac482c49 100644
--- a/examples/MI/classify_alexmi_with_quantum_pipeline.py
+++ b/examples/MI/classify_alexmi_with_quantum_pipeline.py
@@ -20,7 +20,7 @@
 from moabb.evaluations import WithinSessionEvaluation
 from moabb.paradigms import MotorImagery
 
-# inject convex distance and mean to pyriemann (if not done already)
+# inject cpm distance and mean to pyriemann (if not done already)
 from pyriemann_qiskit.utils import distance, mean  # noqa
 from pyriemann_qiskit.pipelines import (
     QuantumMDMWithRiemannianPipeline,
@@ -68,8 +68,8 @@
 pipelines = {}
 
 # Will run QAOA under the hood
-pipelines["mean=logeuclid/distance=convex"] = QuantumMDMWithRiemannianPipeline(
-    convex_metric="distance", quantum=True
+pipelines["mean=logeuclid/distance=cpm"] = QuantumMDMWithRiemannianPipeline(
+    metric="distance", quantum=True
 )
 
 # Classical baseline for evaluation
diff --git a/pyriemann_qiskit/classification.py b/pyriemann_qiskit/classification.py
index 8a9b6f23..98cf4df6 100644
--- a/pyriemann_qiskit/classification.py
+++ b/pyriemann_qiskit/classification.py
@@ -21,6 +21,7 @@
 from qiskit_ibm_provider import IBMProvider, least_busy
 from qiskit_machine_learning.algorithms import QSVC, VQC, PegasosQSVC
 from qiskit_machine_learning.kernels.quantum_kernel import QuantumKernel
+from qiskit_optimization.algorithms import CobylaOptimizer
 from sklearn.base import BaseEstimator, ClassifierMixin
 from sklearn.svm import SVC
 
@@ -582,7 +583,7 @@ class QuanticMDM(QuanticClassifierBase):
 
     """Quantum-enhanced MDM classifier
 
-    This class is a convex implementation of the Minimum Distance to Mean (MDM)
+    This class is a quantic implementation of the Minimum Distance to Mean (MDM)
     [1]_, which can run with quantum optimization.
     Only log-Euclidean distance between trial and class prototypes is supported
     at the moment, but any type of metric can be used for centroid estimation.
@@ -593,11 +594,13 @@ class QuanticMDM(QuanticClassifierBase):
     .. versionchanged:: 0.1.0
         Fix: copy estimator not keeping base class parameters.
     .. versionchanged:: 0.2.0
-        Add seed parameter
+        Add seed parameter.
+        Add regularization parameter.
+        Add classical_optimizer parameter.
 
     Parameters
     ----------
-    metric : string | dict, default={"mean": 'logeuclid', "distance": 'convex'}
+    metric : string | dict, default={"mean": 'logeuclid', "distance": 'cpm'}
         The type of metric used for centroid and distance estimation.
         see `mean_covariance` for the list of supported metric.
         the metric could be a dict with two keys, `mean` and `distance` in
@@ -606,7 +609,7 @@ class QuanticMDM(QuanticClassifierBase):
         the mean in order to boost the computional speed and 'riemann' for the
         distance in order to keep the good sensitivity for the classification.
     quantum : bool (default: True)
-        Only applies if `metric` contains a convex distance or mean.
+        Only applies if `metric` contains a cpm distance or mean.
 
         - If true will run on local or remote backend
           (depending on q_account_token value),
@@ -624,6 +627,10 @@ class QuanticMDM(QuanticClassifierBase):
         Random seed for the simulation
     upper_bound : int (default: 7)
         The maximum integer value for matrix normalization.
+    regularization: MixinTransformer (defulat: None)
+        Additional post-processing to regularize means.
+    classical_optimizer : OptimizationAlgorithm
+        An instance of OptimizationAlgorithm [3]_
 
     See Also
     --------
@@ -642,26 +649,32 @@ class QuanticMDM(QuanticClassifierBase):
         A. Barachant, S. Bonnet, M. Congedo and C. Jutten. 9th International
         Conference Latent Variable Analysis and Signal Separation
         (LVA/ICA 2010), LNCS vol. 6365, 2010, p. 629-636.
+    .. [3] \
+        https://qiskit-community.github.io/qiskit-optimization/stubs/qiskit_optimization.algorithms.OptimizationAlgorithm.html#optimizationalgorithm
     """
 
     def __init__(
         self,
-        metric={"mean": "logeuclid", "distance": "convex"},
+        metric={"mean": "logeuclid", "distance": "logeuclid_cpm"},
         quantum=True,
         q_account_token=None,
         verbose=True,
         shots=1024,
         seed=None,
         upper_bound=7,
+        regularization=None,
+        classical_optimizer=CobylaOptimizer(rhobeg=2.1, rhoend=0.000001),
     ):
         QuanticClassifierBase.__init__(
             self, quantum, q_account_token, verbose, shots, None, seed
         )
         self.metric = metric
         self.upper_bound = upper_bound
+        self.regularization = regularization
+        self.classical_optimizer = classical_optimizer
 
     def _init_algo(self, n_features):
-        self._log("Convex MDM initiating algorithm")
+        self._log("Quantic MDM initiating algorithm")
         classifier = MDM(metric=self.metric)
         if self.quantum:
             self._log("Using NaiveQAOAOptimizer")
@@ -670,10 +683,17 @@ def _init_algo(self, n_features):
             )
         else:
             self._log("Using ClassicalOptimizer (COBYLA)")
-            self._optimizer = ClassicalOptimizer()
+            self._optimizer = ClassicalOptimizer(self.classical_optimizer)
         set_global_optimizer(self._optimizer)
         return classifier
 
+    def _train(self, X, y):
+        QuanticClassifierBase._train(self, X, y)
+        if self.regularization is not None:
+            self._classifier.covmeans_ = self.regularization.fit_transform(
+                self._classifier.covmeans_
+            )
+
     def predict(self, X):
         """Calculates the predictions.
 
diff --git a/pyriemann_qiskit/pipelines.py b/pyriemann_qiskit/pipelines.py
index 56c1c5d0..fb4ac338 100644
--- a/pyriemann_qiskit/pipelines.py
+++ b/pyriemann_qiskit/pipelines.py
@@ -4,9 +4,11 @@
 from sklearn.decomposition import PCA
 from sklearn.pipeline import make_pipeline
 from sklearn.ensemble import VotingClassifier
+from qiskit_optimization.algorithms import CobylaOptimizer
 from pyriemann.estimation import XdawnCovariances, ERPCovariances
 from pyriemann.tangentspace import TangentSpace
 from pyriemann.preprocessing import Whitening
+from pyriemann_qiskit.utils.mean import is_cpm_mean
 from pyriemann_qiskit.utils.filtering import NoDimRed
 from pyriemann_qiskit.utils.hyper_params_factory import (
     # gen_zz_feature_map,
@@ -304,20 +306,14 @@ def _create_pipe(self):
 
 class QuantumMDMWithRiemannianPipeline(BasePipeline):
 
-    """MDM with Riemannian pipeline adapted for convex metrics.
+    """MDM with Riemannian pipeline adapted for cpm metrics.
 
     It can run on classical or quantum optimizer.
 
     Parameters
     ----------
-    convex_metric : string (default: "distance")
-        `metric` passed to the inner QuanticMDM depends on the
-        `convex_metric` as follows (convex_metric => metric):
-
-        - "distance" => {mean=logeuclid, distance=convex},
-        - "mean" => {mean=convex, distance=euclid},
-        - "both" => {mean=convex, distance=convex},
-        - other => same as "distance".
+    metric : string | dict, default={"mean": 'logeuclid', "distance": 'logeuclid_cpm'}
+        The type of metric used for centroid and distance estimation.
     quantum : bool (default: True)
         - If true will run on local or remote backend
           (depending on q_account_token value),
@@ -333,6 +329,10 @@ class QuantumMDMWithRiemannianPipeline(BasePipeline):
         Number of repetitions of each circuit, for sampling.
     upper_bound : int (default: 7)
         The maximum integer value for matrix normalization.
+    regularization: MixinTransformer (defulat: None)
+        Additional post-processing to regularize means.
+    classical_optimizer : OptimizationAlgorithm
+        An instance of OptimizationAlgorithm [1]_
 
     Attributes
     ----------
@@ -342,40 +342,49 @@ class QuantumMDMWithRiemannianPipeline(BasePipeline):
     Notes
     -----
     .. versionadded:: 0.1.0
+    .. versionchanged:: 0.2.0
+        Add regularization parameter.
+        Add classical_optimizer parameter.
+        Change metric, so you can pass the kernel of your choice\
+            as when using MDM.
 
     See Also
     --------
     QuanticMDM
 
+    References
+    ----------
+    .. [1] \
+        https://qiskit-community.github.io/qiskit-optimization/stubs/qiskit_optimization.algorithms.OptimizationAlgorithm.html#optimizationalgorithm
+
     """
 
     def __init__(
         self,
-        convex_metric="distance",
+        metric={"mean": "logeuclid", "distance": "logeuclid_cpm"},
         quantum=True,
         q_account_token=None,
         verbose=True,
         shots=1024,
         upper_bound=7,
+        regularization=None,
+        classical_optimizer=CobylaOptimizer(rhobeg=2.1, rhoend=0.000001),
     ):
-        self.convex_metric = convex_metric
+        self.metric = metric
         self.quantum = quantum
         self.q_account_token = q_account_token
         self.verbose = verbose
         self.shots = shots
         self.upper_bound = upper_bound
+        self.regularization = regularization
+        self.classical_optimizer = classical_optimizer
 
         BasePipeline.__init__(self, "QuantumMDMWithRiemannianPipeline")
 
     def _create_pipe(self):
-        if self.convex_metric == "both":
-            metric = {"mean": "convex", "distance": "convex"}
-        elif self.convex_metric == "mean":
-            metric = {"mean": "convex", "distance": "euclid"}
-        else:
-            metric = {"mean": "logeuclid", "distance": "convex"}
-
-        if metric["mean"] == "convex":
+        print(self.metric)
+        print(self.metric["mean"])
+        if is_cpm_mean(self.metric["mean"]):
             if self.quantum:
                 covariances = XdawnCovariances(
                     nfilter=1, estimator="scm", xdawn_estimator="lwf"
@@ -389,12 +398,14 @@ def _create_pipe(self):
             filtering = NoDimRed()
 
         clf = QuanticMDM(
-            metric=metric,
+            metric=self.metric,
             quantum=self.quantum,
             q_account_token=self.q_account_token,
             verbose=self.verbose,
             shots=self.shots,
             upper_bound=self.upper_bound,
+            regularization=self.regularization,
+            classical_optimizer=self.classical_optimizer,
         )
 
         return make_pipeline(covariances, filtering, clf)
@@ -407,8 +418,8 @@ class QuantumMDMVotingClassifier(BasePipeline):
     Voting classifier with two configurations of
     QuantumMDMWithRiemannianPipeline:
 
-    - with mean = convex and distance = euclid,
-    - with mean = logeuclid and distance = convex.
+    - with mean = euclid_cpm and distance = euclid,
+    - with mean = logeuclid and distance = logeuclid_cpm.
 
     Parameters
     ----------
@@ -460,16 +471,16 @@ def __init__(
         BasePipeline.__init__(self, "QuantumMDMVotingClassifier")
 
     def _create_pipe(self):
-        clf_mean_logeuclid_dist_convex = QuantumMDMWithRiemannianPipeline(
-            "distance",
+        clf_mean_logeuclid_dist_cpm = QuantumMDMWithRiemannianPipeline(
+            {"mean": "logeuclid", "distance": "logeuclid_cpm"},
             self.quantum,
             self.q_account_token,
             self.verbose,
             self.shots,
             self.upper_bound,
         )
-        clf_mean_convex_dist_euclid = QuantumMDMWithRiemannianPipeline(
-            "mean",
+        clf_mean_cpm_dist_euclid = QuantumMDMWithRiemannianPipeline(
+            {"mean": "euclid_cpm", "distance": "euclid"},
             self.quantum,
             self.q_account_token,
             self.verbose,
@@ -480,8 +491,8 @@ def _create_pipe(self):
         return make_pipeline(
             VotingClassifier(
                 [
-                    ("mean_logeuclid_dist_convex", clf_mean_logeuclid_dist_convex),
-                    ("mean_convex_dist_euclid ", clf_mean_convex_dist_euclid),
+                    ("mean_logeuclid_dist_cpm", clf_mean_logeuclid_dist_cpm),
+                    ("mean_cpm_dist_euclid ", clf_mean_cpm_dist_euclid),
                 ],
                 voting="soft",
             )
diff --git a/pyriemann_qiskit/utils/__init__.py b/pyriemann_qiskit/utils/__init__.py
index c44c297b..b4704ce1 100644
--- a/pyriemann_qiskit/utils/__init__.py
+++ b/pyriemann_qiskit/utils/__init__.py
@@ -18,7 +18,8 @@
     add_moabb_dataframe_results_to_caches,
     convert_caches_to_dataframes,
 )
-from .distance import logeucl_dist_convex
+from . import distance
+from . import mean
 
 __all__ = [
     "hyper_params_factory",
@@ -36,7 +37,8 @@
     "NaiveQAOAOptimizer",
     "set_global_optimizer",
     "get_global_optimizer",
-    "logeucl_dist_convex",
+    "distance",
+    "mean",
     "FirebaseConnector",
     "Cache",
     "generate_caches",
diff --git a/pyriemann_qiskit/utils/distance.py b/pyriemann_qiskit/utils/distance.py
index cc1970dc..eb2b8df9 100644
--- a/pyriemann_qiskit/utils/distance.py
+++ b/pyriemann_qiskit/utils/distance.py
@@ -2,28 +2,49 @@
 from docplex.mp.model import Model
 from pyriemann_qiskit.utils.docplex import ClassicalOptimizer, get_global_optimizer
 from pyriemann.classification import MDM
-from pyriemann.utils.distance import distance_functions
+from pyriemann.utils.distance import distance_functions, distance_logeuclid
 from pyriemann.utils.base import logm
+from pyriemann.utils.mean import mean_logeuclid
+from typing_extensions import deprecated
 
 
-def logeucl_dist_convex(X, y, optimizer=ClassicalOptimizer()):
-    """Convex formulation of the MDM algorithm with log-Euclidean metric.
+@deprecated(
+    "logeucl_dist_convex is deprecated and will be removed in 0.3.0; "
+    "please use distance_logeuclid_cpm."
+)
+def logeucl_dist_convex():
+    pass
+
+
+def distance_logeuclid_cpm(A, B, optimizer=ClassicalOptimizer(), return_weights=False):
+    """Log-Euclidean distance to a convex hull of SPD matrices.
+
+    Log-Euclidean distance between a SPD matrix B and the convex hull of a set
+    of SPD matrices A [1]_, formulated as a Constraint Programming Model (CPM)
+    [2]_.
 
     Parameters
     ----------
-    X : ndarray, shape (n_classes, n_channels, n_channels)
+    A : ndarray, shape (n_matrices, n_channels, n_channels)
         Set of SPD matrices.
-    y : ndarray, shape (n_channels, n_channels)
-        A trial
+    B : ndarray, shape (n_channels, n_channels)
+        SPD matrix.
     optimizer: pyQiskitOptimizer
       An instance of pyQiskitOptimizer.
+    return_weights : bool, default=False
+        Whether to return optimized weights.
 
     Returns
     -------
-    weights : ndarray, shape (n_classes,)
-        The weights associated with each class.
-        Higher the weight, closer it is to the class prototype.
-        Weights are not normalized.
+    distance : float
+        Log-Euclidean distance between the SPD matrix B and the convex hull of
+        the set of SPD matrices A, defined as the distance between B and the
+        matrix of the convex hull closest to matrix B.
+    weights : ndarray, shape (n_matrices,)
+        If return_weights is True,
+        it returns the optimized weights for the set of SPD matrices A.
+        Using these weights, the weighted Log-Euclidean mean of set A
+        provides the matrix of the convex hull closest to matrix B.
 
     Notes
     -----
@@ -32,53 +53,92 @@ def logeucl_dist_convex(X, y, optimizer=ClassicalOptimizer()):
     References
     ----------
     .. [1] \
-        http://ibmdecisionoptimization.github.io/docplex-doc/mp/_modules/docplex/mp/model.html#Model
-    """
+        K. Zhao, A. Wiliem, S. Chen, and B. C. Lovell,
+        ‘Convex Class Model on Symmetric Positive Definite Manifolds’,
+        Image and Vision Computing, 2019.
+    .. [2] \
+        http://ibmdecisionoptimization.github.io/docplex-doc/cp/creating_model.html
 
-    optimizer = get_global_optimizer(optimizer)
-
-    n_classes, _, _ = X.shape
-    classes = range(n_classes)
+    """
+    n_matrices, _, _ = A.shape
+    matrices = range(n_matrices)
 
     def log_prod(m1, m2):
         return np.nansum(logm(m1).flatten() * logm(m2).flatten())
 
     prob = Model()
-
+    optimizer = get_global_optimizer(optimizer)
     # should be part of the optimizer
-    w = optimizer.get_weights(prob, classes)
-
-    _2VecLogYD = 2 * prob.sum(w[i] * log_prod(y, X[i]) for i in classes)
+    w = optimizer.get_weights(prob, matrices)
 
-    wtDw = prob.sum(
-        w[i] * w[j] * log_prod(X[i], X[j]) for i in classes for j in classes
+    wtLogAtLogAw = prob.sum(
+        w[i] * w[j] * log_prod(A[i], A[j]) for i in matrices for j in matrices
     )
-
-    objectives = wtDw - _2VecLogYD
+    wLogBLogA = prob.sum(w[i] * log_prod(B, A[i]) for i in matrices)
+    objectives = wtLogAtLogAw - 2 * wLogBLogA
 
     prob.set_objective("min", objectives)
+    prob.add_constraint(prob.sum(w) == 1)
+    weights = optimizer.solve(prob, reshape=False)
 
-    result = optimizer.solve(prob, reshape=False)
+    # compute nearest matrix and distance
+    C = mean_logeuclid(A, weights)
+    distance = distance_logeuclid(C, B)
 
-    return 1 - result
+    if return_weights:
+        return distance, weights
+    return distance
 
 
 _mdm_predict_distances_original = MDM._predict_distances
 
 
 def predict_distances(mdm, X):
-    if mdm.metric_dist == "convex":
+    if mdm.metric_dist == "logeuclid_cpm":
         centroids = np.array(mdm.covmeans_)
-        return np.array([logeucl_dist_convex(centroids, x) for x in X])
+
+        weights = [
+            distance_logeuclid_cpm(centroids, x, return_weights=True)[1] for x in X
+        ]
+        return 1 - np.array(weights)
     else:
         return _mdm_predict_distances_original(mdm, X)
 
 
+def is_cpm_dist(string):
+    """Indicates if the distance is a CPM distance.
+
+    Return True is "string" represents a Constraint Programming Model (CPM) [1]_
+    distance available in the library.
+
+    Parameters
+    ----------
+    string: str
+        A string representation of the distance.
+
+    Returns
+    -------
+    is_cpm_dist : boolean
+        True if "string" represents a CPM distance available in the library.
+
+    Notes
+    -----
+    .. versionadded:: 0.2.0
+
+    References
+    ----------
+    .. [1] \
+        http://ibmdecisionoptimization.github.io/docplex-doc/cp/creating_model.html
+
+    """
+    return "_cpm" in string and string in distance_functions
+
+
 MDM._predict_distances = predict_distances
 
 # This is only for validation inside the MDM.
 # In fact, we override the _predict_distances method
-# inside MDM to directly use logeucl_dist_convex when the metric is "convex"
+# inside MDM to directly use distance_logeuclid_cpm when the metric is "logeuclid_cpm"
 # This is due to the fact the the signature of this method is different from
 # the usual distance functions.
-distance_functions["convex"] = logeucl_dist_convex
+distance_functions["logeuclid_cpm"] = distance_logeuclid_cpm
diff --git a/pyriemann_qiskit/utils/docplex.py b/pyriemann_qiskit/utils/docplex.py
index 32397b44..5e94b2d6 100644
--- a/pyriemann_qiskit/utils/docplex.py
+++ b/pyriemann_qiskit/utils/docplex.py
@@ -307,19 +307,32 @@ class ClassicalOptimizer(pyQiskitOptimizer):
 
     """Wrapper for the classical Cobyla optimizer.
 
+    Attributes
+    ----------
+    optimizer : OptimizationAlgorithm
+        An instance of OptimizationAlgorithm [1]_
+
     Notes
     -----
     .. versionadded:: 0.0.2
     .. versionchanged:: 0.0.4
+    .. versionchanged:: 0.2.0
+        Add attribute `optimizer`.
 
     See Also
     --------
     pyQiskitOptimizer
 
+    References
+    ----------
+    .. [1] \
+        https://qiskit-community.github.io/qiskit-optimization/stubs/qiskit_optimization.algorithms.OptimizationAlgorithm.html#optimizationalgorithm
+
     """
 
-    def __init__(self):
+    def __init__(self, optimizer=CobylaOptimizer(rhobeg=2.1, rhoend=0.000001)):
         pyQiskitOptimizer.__init__(self)
+        self.optimizer = optimizer
 
     """Helper to create a docplex representation of a
     covariance matrix variable.
@@ -360,7 +373,7 @@ def covmat_var(self, prob, channels, name):
         return square_cont_mat_var(prob, channels, name)
 
     def _solve_qp(self, qp, reshape=True):
-        result = CobylaOptimizer(rhobeg=2.1, rhoend=0.000001).solve(qp).x
+        result = self.optimizer.solve(qp).x
         if reshape:
             n_channels = int(math.sqrt(result.shape[0]))
             return np.reshape(result, (n_channels, n_channels))
diff --git a/pyriemann_qiskit/utils/mean.py b/pyriemann_qiskit/utils/mean.py
index 15cdfb3c..e169efd0 100644
--- a/pyriemann_qiskit/utils/mean.py
+++ b/pyriemann_qiskit/utils/mean.py
@@ -1,69 +1,149 @@
+from typing_extensions import deprecated
 from docplex.mp.model import Model
 from pyriemann.utils.mean import mean_functions
 from pyriemann_qiskit.utils.docplex import ClassicalOptimizer, get_global_optimizer
-from pyriemann.estimation import Shrinkage
+from pyriemann.utils.base import logm, expm
+from qiskit_optimization.algorithms import ADMMOptimizer
 
 
-def fro_mean_convex(
-    covmats, sample_weight=None, optimizer=ClassicalOptimizer(), shrink=True
-):
-    """Convex formulation of the mean with Frobenius distance.
+@deprecated(
+    "fro_mean_convex is deprecated and will be removed in 0.3.0; "
+    "please use mean_euclid_cpm."
+)
+def fro_mean_convex():
+    pass
+
+
+def mean_euclid_cpm(X, sample_weight=None, optimizer=ClassicalOptimizer()):
+    """Euclidean mean with Constraint Programming Model.
+
+    Constraint Programming Model (CPM) [1]_ formulation of the mean
+    with Euclidean distance.
 
     Parameters
     ----------
-    covmats: ndarray, shape (n_matrices, n_channels, n_channels)
+    X : ndarray, shape (n_matrices, n_channels, n_channels)
         Set of SPD matrices.
-    sample_weights:  None | ndarray, shape (n_matrices,), default=None
+    sample_weights : None | ndarray, shape (n_matrices,), default=None
         Weights for each matrix. Never used in practice.
         It is kept only for standardization with pyRiemann.
-    optimizer: pyQiskitOptimizer
+    optimizer : pyQiskitOptimizer
         An instance of pyQiskitOptimizer.
-    shrink: boolean (default: true)
-        If True, it applies shrinkage regularization [2]_
-        of the resulting covariance matrix.
 
     Returns
     -------
     mean : ndarray, shape (n_channels, n_channels)
-        Convex-optimized Frobenius mean.
+        CPM-optimized Euclidean mean.
 
     Notes
     -----
     .. versionadded:: 0.0.3
     .. versionchanged:: 0.0.4
         Add regularization of the results.
+    .. versionchanged:: 0.2.0
+        Rename from `fro_mean_convex` to `mean_euclid_cpm`
+        Remove shrinkage
 
     References
     ----------
     .. [1] \
-        http://ibmdecisionoptimization.github.io/docplex-doc/mp/_modules/docplex/mp/model.html#Model
-    .. [2] \
-        https://pyriemann.readthedocs.io/en/v0.4/generated/pyriemann.estimation.Shrinkage.html
+        http://ibmdecisionoptimization.github.io/docplex-doc/cp/creating_model.html
     """
 
     optimizer = get_global_optimizer(optimizer)
 
-    n_trials, n_channels, _ = covmats.shape
+    n_matrices, n_channels, _ = X.shape
     channels = range(n_channels)
-    trials = range(n_trials)
+    matrices = range(n_matrices)
 
     prob = Model()
 
     X_mean = optimizer.covmat_var(prob, channels, "fro_mean")
 
-    def _fro_dist(A, B):
+    def _dist_euclid(A, B):
         A = optimizer.convert_covmat(A)
         return prob.sum_squares(A[r, c] - B[r, c] for r in channels for c in channels)
 
-    objectives = prob.sum(_fro_dist(covmats[i], X_mean) for i in trials)
+    objectives = prob.sum(_dist_euclid(X[i], X_mean) for i in matrices)
 
     prob.set_objective("min", objectives)
 
     result = optimizer.solve(prob)
 
-    if shrink:
-        return Shrinkage(shrinkage=0.9).transform([result])[0]
     return result
 
 
-mean_functions["convex"] = fro_mean_convex
+def mean_logeuclid_cpm(
+    X, sample_weight=None, optimizer=ClassicalOptimizer(optimizer=ADMMOptimizer())
+):
+    """Log-Euclidean mean with Constraint Programming Model.
+
+    Constraint Programming Model (CPM) [2]_ formulation of the mean
+    with Log-Euclidean distance [1]_.
+
+    Parameters
+    ----------
+    X : ndarray, shape (n_matrices, n_channels, n_channels)
+        Set of SPD matrices.
+    sample_weights : None | ndarray, shape (n_matrices,), default=None
+        Weights for each matrix. Never used in practice.
+        It is kept only for standardization with pyRiemann.
+    optimizer : pyQiskitOptimizer
+        An instance of pyQiskitOptimizer.
+
+    Returns
+    -------
+    mean : ndarray, shape (n_channels, n_channels)
+        CPM-optimized Log-Euclidean mean.
+
+    Notes
+    -----
+    .. versionadded:: 0.2.0
+
+    References
+    ----------
+    .. [1] \
+        Geometric means in a novel vector space structure on
+        symmetric positive-definite matrices
+        V. Arsigny, P. Fillard, X. Pennec, and N. Ayache.
+        SIAM Journal on Matrix Analysis and Applications. Volume 29, Issue 1 (2007).
+    .. [2] \
+        http://ibmdecisionoptimization.github.io/docplex-doc/cp/creating_model.html
+    """
+
+    log_X = logm(X)
+    result = mean_euclid_cpm(log_X, sample_weight, optimizer)
+    return expm(result)
+
+
+def is_cpm_mean(string):
+    """Indicates if the mean is a CPM mean.
+
+    Return True is "string" represents a Constraint Programming Model (CPM) [1]_
+    mean available in the library.
+
+    Parameters
+    ----------
+    string: str
+        A string representation of the mean.
+
+    Returns
+    -------
+    is_cpm_mean : boolean
+        True if "string" represents a CPM mean aailable in the library.
+
+    Notes
+    -----
+    .. versionadded:: 0.2.0
+
+    References
+    ----------
+    .. [1] \
+        http://ibmdecisionoptimization.github.io/docplex-doc/cp/creating_model.html
+
+    """
+    return "_cpm" in string and string in mean_functions
+
+
+mean_functions["euclid_cpm"] = mean_euclid_cpm
+mean_functions["logeuclid_cpm"] = mean_logeuclid_cpm
diff --git a/requirements.txt b/requirements.txt
index bc45cb22..02c17745 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -2,6 +2,7 @@ numpy<1.27
 cython
 scikit-learn==1.3.2
 git+https://github.com/pyRiemann/pyRiemann#egg=pyriemann
+qiskit==0.45.0
 qiskit_machine_learning==0.6.1
 qiskit-ibm-provider==0.7.3
 qiskit-optimization==0.5.0
diff --git a/tests/conftest.py b/tests/conftest.py
index 333f39a2..ea97f52b 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -4,7 +4,7 @@
 import pytest
 import numpy as np
 from functools import partial
-from pyriemann.datasets import make_covariances
+from pyriemann.datasets import make_matrices
 from pyriemann_qiskit.datasets import get_mne_sample
 from operator import itemgetter
 
@@ -36,7 +36,9 @@ def rndstate():
 @pytest.fixture
 def get_covmats(rndstate):
     def _gen_cov(n_matrices, n_channels):
-        return make_covariances(n_matrices, n_channels, rndstate, return_params=False)
+        return make_matrices(
+            n_matrices, n_channels, "spd", rndstate, return_params=False
+        )
 
     return _gen_cov
 
@@ -44,7 +46,9 @@ def _gen_cov(n_matrices, n_channels):
 @pytest.fixture
 def get_covmats_params(rndstate):
     def _gen_cov_params(n_matrices, n_channels):
-        return make_covariances(n_matrices, n_channels, rndstate, return_params=True)
+        return make_matrices(
+            n_matrices, n_channels, "spd", rndstate, return_params=True
+        )
 
     return _gen_cov_params
 
@@ -94,11 +98,13 @@ def _get_dataset(n_samples, n_features, n_classes, type="bin"):
             samples = _get_separable_feats(n_samples, n_features, n_classes)
             labels = _get_labels(n_samples, n_classes)
         elif type == "rand_cov":
-            samples = make_covariances(n_samples, n_features, 0, return_params=False)
+            samples = make_matrices(
+                n_samples, n_features, "spd", 0, return_params=False
+            )
             labels = _get_labels(n_samples, n_classes)
         elif type == "bin_cov":
-            samples_0 = make_covariances(
-                n_samples // n_classes, n_features, 0, return_params=False
+            samples_0 = make_matrices(
+                n_samples // n_classes, n_features, "spd", 0, return_params=False
             )
             samples = np.concatenate(
                 [samples_0 * (i + 1) for i in range(n_classes)], axis=0
diff --git a/tests/test_utils_distance.py b/tests/test_utils_distance.py
index e59f4ddc..757781dc 100644
--- a/tests/test_utils_distance.py
+++ b/tests/test_utils_distance.py
@@ -3,8 +3,8 @@
 from pyriemann_qiskit.utils import (
     ClassicalOptimizer,
     NaiveQAOAOptimizer,
-    logeucl_dist_convex,
 )
+from pyriemann_qiskit.utils.distance import distance_logeuclid_cpm
 from pyriemann_qiskit.datasets import get_mne_sample
 from pyriemann.classification import MDM
 from pyriemann.estimation import XdawnCovariances
@@ -13,7 +13,7 @@
 
 
 def test_performance():
-    metric = {"mean": "logeuclid", "distance": "convex"}
+    metric = {"mean": "logeuclid", "distance": "logeuclid_cpm"}
 
     clf = make_pipeline(XdawnCovariances(), MDM(metric=metric))
     skf = StratifiedKFold(n_splits=3)
@@ -23,10 +23,11 @@ def test_performance():
 
 
 @pytest.mark.parametrize("optimizer", [ClassicalOptimizer(), NaiveQAOAOptimizer()])
-def test_logeucl_dist_convex(optimizer):
+def test_distance_logeuclid_cpm(optimizer):
     X_0 = np.array([[0.9, 1.1], [0.9, 1.1]])
     X_1 = X_0 + 1
     X = np.stack((X_0, X_1))
     y = (X_0 + X_1) / 3
-    distances = logeucl_dist_convex(X, y, optimizer=optimizer)
+    _, weights = distance_logeuclid_cpm(X, y, optimizer=optimizer, return_weights=True)
+    distances = 1 - weights
     assert distances.argmin() == 0
diff --git a/tests/test_utils_mean.py b/tests/test_utils_mean.py
index 4fe1ac52..2854ad5d 100644
--- a/tests/test_utils_mean.py
+++ b/tests/test_utils_mean.py
@@ -1,78 +1,112 @@
 import pytest
 import numpy as np
-from pyriemann.utils.mean import mean_euclid
-from pyriemann.classification import MDM
-from pyriemann.estimation import XdawnCovariances
+from pyriemann.utils.mean import mean_euclid, mean_logeuclid
+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 fro_mean_convex
+from pyriemann_qiskit.utils.mean import mean_euclid_cpm, mean_logeuclid_cpm
 from pyriemann_qiskit.utils import ClassicalOptimizer, NaiveQAOAOptimizer
-
-
-def test_performance(get_covmats, get_labels):
-    metric = {"mean": "convex", "distance": "euclid"}
-
-    clf = make_pipeline(XdawnCovariances(), MDM(metric=metric))
-    skf = StratifiedKFold(n_splits=5)
-    n_matrices, n_channels, n_classes = 100, 3, 2
-    covset = get_covmats(n_matrices, n_channels)
-    labels = get_labels(n_matrices, n_classes)
+from pyriemann_qiskit.classification import QuanticMDM
+from pyriemann_qiskit.datasets import get_mne_sample
+from qiskit_optimization.algorithms import ADMMOptimizer
+
+
+@pytest.mark.parametrize(
+    "kernel",
+    [
+        ({"mean": "euclid_cpm", "distance": "euclid"}, Shrinkage(shrinkage=0.9)),
+        ({"mean": "logeuclid_cpm", "distance": "logeuclid"}, Shrinkage(shrinkage=0.9)),
+    ],
+)
+def test_performance(kernel):
+    metric, regularization = kernel
+    clf = make_pipeline(
+        XdawnCovariances(),
+        QuanticMDM(metric=metric, regularization=regularization, quantum=False),
+    )
+    skf = StratifiedKFold(n_splits=3)
+
+    covset, labels = get_mne_sample()
 
     score = cross_val_score(clf, covset, labels, cv=skf, scoring="roc_auc")
     assert score.mean() > 0
 
 
-def test_mean_convex_vs_euclid(get_covmats):
-    """Test that euclidian and convex mean returns close results"""
+@pytest.mark.parametrize(
+    "means", [(mean_euclid, mean_euclid_cpm), (mean_logeuclid, mean_logeuclid_cpm)]
+)
+def test_analytic_vs_cpm_mean(get_covmats, means):
+    """Test that analytic and cpm mean returns close results"""
+    analytic_mean, cpm_mean = means
     n_trials, n_channels = 5, 3
     covmats = get_covmats(n_trials, n_channels)
-    C = fro_mean_convex(covmats, shrink=False)
-    C_euclid = mean_euclid(covmats)
-    assert np.allclose(C, C_euclid, atol=0.001)
+    C = cpm_mean(covmats)
+    C_analytic = analytic_mean(covmats)
+    assert np.allclose(C, C_analytic, atol=0.00001)
 
 
-def test_mean_convex_shape(get_covmats):
+@pytest.mark.parametrize("mean", [mean_euclid_cpm, mean_logeuclid_cpm])
+def test_mean_cpm_shape(get_covmats, mean):
     """Test the shape of mean"""
     n_trials, n_channels = 5, 3
     covmats = get_covmats(n_trials, n_channels)
-    C = fro_mean_convex(covmats)
+    C = mean(covmats)
     assert C.shape == (n_channels, n_channels)
 
 
-@pytest.mark.parametrize("optimizer", [ClassicalOptimizer(), NaiveQAOAOptimizer()])
-def test_mean_convex_all_zeros(optimizer):
+@pytest.mark.parametrize(
+    "optimizer",
+    [ClassicalOptimizer(optimizer=ADMMOptimizer()), NaiveQAOAOptimizer()],
+)
+@pytest.mark.parametrize("mean", [mean_euclid_cpm])
+def test_mean_cpm_all_zeros(optimizer, mean):
     """Test that the mean of covariance matrices containing zeros
     is a matrix filled with zeros"""
     n_trials, n_channels = 5, 2
     covmats = np.zeros((n_trials, n_channels, n_channels))
-    C = fro_mean_convex(covmats, optimizer=optimizer, shrink=False)
+    C = mean(covmats, optimizer=optimizer)
     assert np.allclose(covmats[0], C, atol=0.001)
 
 
-def test_mean_convex_all_ones():
+@pytest.mark.parametrize(
+    "optimizer",
+    [ClassicalOptimizer(optimizer=ADMMOptimizer()), NaiveQAOAOptimizer()],
+)
+@pytest.mark.parametrize("mean", [mean_euclid_cpm])
+def test_mean_cpm_all_ones(optimizer, mean):
     """Test that the mean of covariance matrices containing ones
     is a matrix filled with ones"""
     n_trials, n_channels = 5, 2
     covmats = np.ones((n_trials, n_channels, n_channels))
-    C = fro_mean_convex(covmats, shrink=False)
+    C = mean(covmats, optimizer=optimizer)
     assert np.allclose(covmats[0], C, atol=0.001)
 
 
-def test_mean_convex_all_equals():
+@pytest.mark.parametrize(
+    "optimizer",
+    [ClassicalOptimizer(optimizer=ADMMOptimizer()), NaiveQAOAOptimizer()],
+)
+@pytest.mark.parametrize("mean", [mean_euclid_cpm])
+def test_mean_cpm_all_equals(optimizer, mean):
     """Test that the mean of covariance matrices filled with the same value
     is a matrix identical to the input"""
     n_trials, n_channels, value = 5, 2, 2.5
     covmats = np.full((n_trials, n_channels, n_channels), value)
-    C = fro_mean_convex(covmats, shrink=False)
+    C = mean(covmats, optimizer=optimizer)
     assert np.allclose(covmats[0], C, atol=0.001)
 
 
-def test_mean_convex_mixed():
+@pytest.mark.parametrize(
+    "optimizer",
+    [ClassicalOptimizer(optimizer=ADMMOptimizer()), NaiveQAOAOptimizer()],
+)
+@pytest.mark.parametrize("mean", [mean_euclid_cpm])
+def test_mean_cpm_mixed(optimizer, mean):
     """Test that the mean of covariances matrices with zero and ones
     is a matrix filled with 0.5"""
     n_trials, n_channels = 5, 2
     covmats_0 = np.zeros((n_trials, n_channels, n_channels))
     covmats_1 = np.ones((n_trials, n_channels, n_channels))
     expected_mean = np.full((n_channels, n_channels), 0.5)
-    C = fro_mean_convex(np.concatenate((covmats_0, covmats_1), axis=0), shrink=False)
+    C = mean(np.concatenate((covmats_0, covmats_1), axis=0), optimizer=optimizer)
     assert np.allclose(expected_mean, C, atol=0.001)