Repository: pyRiemann/pyRiemann

Issue:
Add TSA to transfer learning module
In the transfer learning module, add Tangent Space Alignement (TSA), see
https://www.frontiersin.org/articles/10.3389/fnhum.2022.1049985/pdf

Hints from issue comments:
@qbarthelemy: Code shared by first author
```
import numpy as np
from sklearn.base import BaseEstimator, TransformerMixin
from sklearn.decomposition import PCA
from pyriemann.estimation import Covariances
from pyriemann.transfer import decode_domains
from sklearn.preprocessing import LabelEncoder

class TSAxPCA(BaseEstimator, TransformerMixin):
    def __init__(self, target_domain, eVar=0.999, n_clusters=3, n_components=1):
        self.target_domain = target_domain
        self.eVar = eVar
        self.source_means = None
        self.target_means = None
        self.u = None
        self.vh = None

        self.incomplete_source = None
        self.incomplete_target = None

        self.n_clusters = n_clusters
        self.n_components = n_components

        self.name = "TSAxPCA"

    def fit(self, X, y_enc):
        X, y, domains = decode_domains(X, y_enc)
        X_source = X[domains != self.target_domain]
        y_source = y[domains != self.target_domain]
        X_target = X[domains == self.target_domain]
        y_target = y[domains == self.target_domain]
        self.fit_source(X_source, y_source)
        self.fit_target(X_target, y_target)

    def fit_source(self, X, y):
        y = LabelEncoder().fit_transform(y=y)
        self.incomplete_source = False
        self.source_pcas = [PCA(n_components=(self.n_components if self.n_components != "max" else X.shape[1])) for i in range(len(np.unique(y)))]
        self.source_means = np.array([np.ndarray.mean(X[y == i, :], axis=0) for i in np.unique(y)])
        for i in np.unique(y):
            Xi = X[y == i]
            ni = len(Xi)
            if ni < self.n_clusters:
                self.incomplete_source = True
                break
            self.source_pcas[i].fit(Xi)
            X_pca = self.source_pcas[i].transform(Xi)
            for j in range(self.n_components):
                args = np.argsort(X_pca[:, j], axis=0)
                self.source_means = np.vstack([self.source_means,
                                               [np.ndarray.mean(Xi[args[round(k*ni/self.n_clusters):round((k+1)*ni/self.n_clusters)]], axis=0) for k in range(self.n_clusters)]])
        if self.target_means is not None:
            if self.incomplete_source or self.incomplete_target:
                print("There is not enough point for either target or source. Compute basic TSA rotation...")
                self.u, self.vh = self.compute_incomplete_rotation(len(np.unique(y)))
            else:
                self.u, self.vh = self.compute_rotations()
        return self

    def compute_rotations(self):
        c12 = np.dot(self.source_means.T, self.target_means)
        u, s, vh = np.linalg.svd(c12)
        if self.eVar < 2:
            n_features = np.ndarray.sum(np.ndarray.cumsum(s) < self.eVar * np.ndarray.sum(s)) + 1
        else:
            n_features = self.eVar
        u = u[:, :n_features]
        vh = vh[:n_features, :]
        return u, vh

    def compute_incomplete_rotation(self, n_classes):
        c12 = np.dot(self.source_means[:n_classes].T, self.target_means[:n_classes])
        u, s, vh = np.linalg.svd(c12)
        if self.eVar < 2:
            n_features = np.ndarray.sum(np.ndarray.cumsum(s) < self.eVar * np.ndarray.sum(s)) + 1
        else:
            n_features = self.eVar
        u = u[:, :n_features]
        vh = vh[:n_features, :]
        return u, vh

    def fit_target(self, X, y):
        y = LabelEncoder().fit_transform(y=y)
        self.incomplete_target = False
        self.target_means = np.array([np.ndarray.mean(X[y == i, :], axis=0) for i in np.unique(y)])
        for i in np.unique(y):
            Xi = X[y == i]
            ni = len(Xi)
            if ni < self.n_clusters:
                self.incomplete_target = True
                break
            X_pca = self.source_pcas[i].transform(Xi)
            for j in range(self.n_components):
                args = np.argsort(X_pca[:, j], axis=0)
                self.target_means = np.vstack([self.target_means,
                                               [np.ndarray.mean(Xi[args[round(k*ni/self.n_clusters):round((k+1)*ni/self.n_clusters)]], axis=0) for k in range(self.n_clusters)]])
        if self.source_means is not None:
            if self.incomplete_source or self.incomplete_target:
                print("There is not enough point for either target or source. Compute basic TSA rotation...")
                self.u, self.vh = self.compute_incomplete_rotation(len(np.unique(y)))
            else:
                self.u, self.vh = self.compute_rotations()
        return self

    def transform(self, X, y=None):
        return np.dot(np.dot(X, self.vh.T), self.u.T)

    def transform_source(self, X, y=None):
        return np.dot(np.dot(X, self.u), self.u.T)


class TSAxPCA_generalized(BaseEstimator, TransformerMixin):
    def __init__(self, target_domain, eVar=0.999, n_clusters=3, n_components=1, equal_dims=False):
        self.target_domain = target_domain
        self.eVar = eVar
        self.n_clusters = n_clusters
        self.n_components = n_components
        self.equal_dims = equal_dims
        self.source_means = None
        self.target_means = None
        self.u = None
        self.vh = None
        self.source_eig = []
        self.source_vals = []
        self.target_eig = []
        self.target_vals = []
        self.source_pcas = []

        self.incomplete_source = None
        self.incomplete_target = None

        self.name = "TSAxPCA_generalized"

    def fit(self, X, y_enc):
        X, y, domains = decode_domains(X, y_enc)
        X_source = X[domains != self.target_domain]
        y_source = y[domains != self.target_domain]
        X_target = X[domains == self.target_domain]
        y_target = y[domains == self.target_domain]
        self.fit_source(X_source, y_source)
        self.fit_target(X_target, y_target)

    def fit_source(self, X, y):
        self.incomplete_source = False
        self.source_pcas = [PCA(n_components=(self.n_components if self.n_components != "max" else X.shape[1])) for i in range(len(np.unique(y)))]
        self.source_means = np.array([np.ndarray.mean(X[y == yi, :], axis=0) for yi in np.unique(y)])
        for i_class, yi in enumerate(np.unique(y)):
            Xi = X[y == yi]
            ni = len(Xi)
            cov = Covariances(estimator='oas').fit_transform(Xi.T.reshape(1, Xi.shape[1], Xi.shape[0]))
            vals, eig = np.linalg.eigh(cov[0])
            arg = np.argsort(np.abs(vals))[::-1]
            self.source_vals.append(vals[arg])
            self.source_eig.append(eig[:, arg])
            if not self.equal_dims:
                Xi_pca = np.dot(Xi, self.source_eig[i_class])
                for i, vect in enumerate(self.source_eig[i_class].T[:self.n_components]):
                    sign = np.sign(np.dot(vect, np.ones(len(vect))))
                    args = np.argsort(sign*Xi_pca[:, i], axis=0)
                    self.source_means = np.vstack([self.source_means,
                                                   [np.ndarray.mean(Xi[args[round(k * ni / self.n_clusters):round(
                                                       (k + 1) * ni / self.n_clusters)]], axis=0) for k in
                                                    range(self.n_clusters)]])
            else:
                Xi = X[y == yi]
                ni = len(Xi)
                if ni < self.n_clusters:
                    self.incomplete_target = True
                    break
                self.source_pcas[i_class].fit(Xi)
                X_pca = self.source_pcas[i_class].transform(Xi)
                for j in range(self.n_components):
                    args = np.argsort(X_pca[:, j], axis=0)
                    self.source_means = np.vstack([self.source_means,
                                                   [np.ndarray.mean(Xi[args[round(k * ni / self.n_clusters):round(
                                                       (k + 1) * ni / self.n_clusters)]], axis=0) for k in
                                                    range(self.n_clusters)]])
        if self.target_means is not None:
            self.u, self.vh = self.compute_rotations()
        return self

    def compute_rotations(self):
        c12 = np.dot(self.source_means.T, self.target_means)
        u, s, vh = np.linalg.svd(c12)
        if self.eVar < 2:
            n_features = np.min([np.ndarray.sum(np.ndarray.cumsum(s) < self.eVar * np.ndarray.sum(s))+1, len(u), len(vh)])
        else:
            n_features = np.min([self.eVar, len(u), len(vh)])
        u = u[:, :n_features]
        vh = vh[:n_features, :]
        return u, vh

    def compute_incomplete_rotation(self, n_classes):
        c12 = np.dot(self.source_means[:n_classes].T, self.target_means[:n_classes])
        u, s, vh = np.linalg.svd(c12)
        if self.eVar < 2:
            n_features = np.min([np.ndarray.sum(np.ndarray.cumsum(s) < self.eVar * np.ndarray.sum(s))+1, len(u), len(vh)])
        else:
            n_features = self.eVar
        u = u[:, :n_features]
        vh = vh[:n_features, :]
        return u, vh

    def fit_target(self, X, y):
        self.incomplete_target = False
        self.target_means = np.array([np.ndarray.mean(X[y == i, :], axis=0) for i in np.unique(y)])
        for i_class, yi in enumerate(np.unique(y)):
            Xi = X[y == yi]
            ni = len(Xi)
            cov = Covariances(estimator='oas').fit_transform(Xi.T.reshape(1, Xi.shape[1], Xi.shape[0]))
            vals, eig = np.linalg.eigh(cov[0])
            arg = np.argsort(np.abs(vals))[::-1]
            self.target_vals.append(vals[arg])
            self.target_eig.append(eig[:, arg])
            if not self.equal_dims:
                Xi_pca = np.dot(Xi, self.target_eig[i_class])
                for i, vect in enumerate(self.target_eig[i_class].T[:self.n_components]):
                    sign = np.sign(np.dot(vect, np.ones(len(vect))))
                    args = np.argsort(sign * Xi_pca[:, i], axis=0)
                    self.target_means = np.vstack([self.target_means,
                                                   [np.ndarray.mean(Xi[args[round(k * ni / self.n_clusters):round(
                                                       (k + 1) * ni / self.n_clusters)]], axis=0) for k in
                                                    range(self.n_clusters)]])
            elif len(self.source_eig) == 0:
                raise(ValueError("The source is not fitted yet. Please use the 'fit' function on Source data before fitting target data."))
            else:
                Xi = X[y == yi]
                ni = len(Xi)
                if ni < self.n_clusters:
                    self.incomplete_target = True
                    break
                X_pca = self.source_pcas[i_class].transform(Xi)
                for j in range(self.n_components):
                    args = np.argsort(X_pca[:, j], axis=0)
                    self.target_means = np.vstack([self.target_means,
                                                   [np.ndarray.mean(Xi[args[round(k * ni / self.n_clusters):round(
                                                       (k + 1) * ni / self.n_clusters)]], axis=0) for k in
                                                    range(self.n_clusters)]])
        if self.source_means is not None:
            if self.incomplete_source or self.incomplete_target:
                print("There is not enough point for either target or source. Compute basic TSA rotation.")
                self.u, self.vh = self.compute_incomplete_rotation(len(np.unique(y)))
            else:
                self.u, self.vh = self.compute_rotations()
        return self

    def transform(self, X, y=None):
        return np.dot(np.dot(X, self.vh.T), self.u.T)

    def transform_source(self, X, y=None):
        return np.dot(np.dot(X, self.u), self.u.T)

def rot_data(X, y_enc, **kwargs):
    X, y, domains = decode_domains(X, y_enc)
    Xs = X[domains != self.target_domain]
    ys = y[domains != self.target_domain]
    Xt = X[domains == self.target_domain]
    yt = y[domains == self.target_domain]
    tsa = TSAxPCA(**kwargs)
    tsa.fit(Xs, ys)
    tsa.fit_target(Xt, yt)
    Xt_rot, _, _ = tsa.transform(Xt)
    return Xt_rot


def align_data(X, y_enc, **kwargs):
    X, y, domains = decode_domains(X, y_enc)
    Xs = X[domains != self.target_domain]
    ys = y[domains != self.target_domain]
    Xt = X[domains == self.target_domain]
    yt = y[domains == self.target_domain]
    Xt_rec = Xt - np.mean(Xt, axis=0) + np.mean(Xs, axis=0)
    Xt_rot = rot_data(Xs, ys, Xt_rec, yt, **kwargs)
    Xt_sca = Xt_rot / np.mean(np.linalg.norm(Xt_rot, axis=1), axis=0) * np.mean(np.linalg.norm(Xs, axis=1), axis=0)
    return Xt_sca


def align_test(Xs, ys, Xt_train, yt_train, Xt_test, yt_test=None, **kwargs):
    Xt_rec_train = Xt_train - np.mean(Xt_train, axis=0) + np.mean(Xs, axis=0)
    Xt_rec_test = Xt_test - np.mean(Xt_train, axis=0) + np.mean(Xs, axis=0)
    tsa = TSAxPCA(**kwargs)
    tsa.fit(Xs, ys)
    tsa.fit_target(Xt_rec_train, yt_train)
    Xt_rot_train, _, _ = tsa.transform(Xt_rec_train)
    Xt_rot_test, _, _ = tsa.transform(Xt_rec_test)

    Xt_sca_test = Xt_rot_test / np.mean(np.linalg.norm(Xt_rot_train, axis=1), axis=0) * np.mean(np.linalg.norm(Xs, axis=1), axis=0)
    return Xt_sca_test

```@qbarthelemy: Example

```
import matplotlib.pyplot as plt
import numpy as np
from pyriemann.tangentspace import tangent_space, untangent_space
from pyriemann.utils.mean import mean_riemann, mean_logeuclid
from pyriemann.estimation import Covariances
from pyriemann.transfer import encode_domains

from _tsa import TSAxPCA


n_classes = 2
n_clusters = 500
n_dims = 2

source_label = "source"
target_label = "target"

n_trials_source = n_classes*n_clusters
n_trials_target = n_trials_source
fake_mean = 0.5
fake_std = 0.5
theta = np.pi/5
fake_rot = np.array([[np.cos(theta), -np.sin(theta)],
                     [np.sin(theta),  np.cos(theta)]])

############# ADDED FALSE SOURCE ##################

Xs_means = 5 * np.array([[-1 / np.sqrt(n_classes * (n_classes - 1)) for _ in range(0, i)] +
                     [np.sqrt((n_classes - 1) / n_classes)] +
                     [-1 / np.sqrt(n_classes * (n_classes - 1)) for _ in range(i + 1, n_classes)] +
                     [0 for _ in range(n_dims - n_classes)] for i in range(n_classes)])
Xs_ts_list = [Xs_means[0] + (-(n_clusters - 1) + 2 * i) / (n_clusters - 1) * np.array(
    [0, 5] + [0 for _ in range(n_dims - n_classes)]) for i in
              range(n_clusters)] + \
             [Xs_means[1] + (-(n_clusters - 1) + 2 * i) / (n_clusters - 1) * np.array(
                 [5, 0] + [0 for _ in range(n_dims - n_classes)]) for i in
              range(n_clusters)]
ys = np.array([0 for _ in range(n_clusters)] + [1 for _ in range(n_clusters)])
Xs_ts = np.array(Xs_ts_list)
Xs_ts = Xs_ts + np.random.randn(n_clusters*n_classes, n_dims)

mean_source = np.mean(Xs_ts, axis=0)
Xs_ts = Xs_ts - mean_source


std_source = np.mean(np.linalg.norm(Xs_ts, axis=1), axis=0)
Xs_ts = Xs_ts / std_source

###### Create the target as a modification of the source data #############
yt = ys
# Xt_ts = np.multiply(1+5*np.array([[np.cos((i+1)*np.pi/2), np.sin((i+1)*np.pi/2), 0] for i in yt]), (np.random.normal(scale=0.15, size=(n_trials_target, 3)) + 0.4 * np.array([[i - 0.5, i - 0.5, 0] for i in yt])))
Xt_ts = Xs_ts
Xt_ts = np.dot(Xt_ts, fake_rot)

X = np.append(Xs_ts, Xt_ts, axis=0)
y = np.append(ys, yt)
yDom = np.append(np.repeat(source_label, n_trials_source), np.repeat(target_label, n_trials_target))

X_enc, y_enc = encode_domains(X, y, yDom)


tsa = TSAxPCA(target_label, n_clusters=5, eVar=1.0)

tsa.fit(X_enc, y_enc=y_enc)
Xt_ts_aligned = tsa.transform(X[yDom == target_label])

plt.figure()
plt.title("Source")
plt.scatter(Xs_ts[:, 0], Xs_ts[:, 1])
plt.xlim([-2, 2])
plt.ylim([-2, 2])

plt.figure()
plt.title("Target")
plt.scatter(Xt_ts[:, 0], Xt_ts[:, 1])
plt.xlim([-2, 2])
plt.ylim([-2, 2])

plt.figure()
plt.title("Target aligned")
plt.scatter(Xt_ts_aligned[:, 0], Xt_ts_aligned[:, 1])
plt.xlim([-2, 2])
plt.ylim([-2, 2])

plt.show()
```@agramfort: cc @apmellot @antoinecollas@qbarthelemy: Closed with merge of #320

Solution patch:
diff --git a/doc/api.rst b/doc/api.rst
index 25c0da1..59cd386 100644
--- a/doc/api.rst
+++ b/doc/api.rst
@@ -51,7 +51,7 @@ Classification
 
     MDM
     FgMDM
-    TSclassifier
+    TSClassifier
     KNearestNeighbor
     SVC
     MeanField
@@ -157,7 +157,7 @@ Transfer Learning
     TLRegressor
     TLDummy
     TLCenter
-    TLStretch
+    TLScale
     TLRotate
     MDWM
 
diff --git a/doc/whatsnew.rst b/doc/whatsnew.rst
index a450a7a..28f2e63 100644
--- a/doc/whatsnew.rst
+++ b/doc/whatsnew.rst
@@ -15,6 +15,10 @@ v0.8.dev
 
 - Enhance :func:`pyriemann.utils.mean.mean_covariance` to support "power" and "poweuclid" metrics. :pr:`329` by :user:`qbarthelemy`
 
+- Add tangent space alignment (TSA) in transfer learning module.
+  ``TLStretch`` is deprecated and renamed into :class:`pyriemann.transfer._estimators.TLScale`,
+  and ``TSclassifier`` into :class:`pyriemann.classification.TSClassifier`. :pr:`320` by :user:`qbarthelemy`
+
 v0.7 (October 2024)
 -------------------
 
diff --git a/examples/biosignal-mi/plot_single.py b/examples/biosignal-mi/plot_single.py
index c418291..922abae 100644
--- a/examples/biosignal-mi/plot_single.py
+++ b/examples/biosignal-mi/plot_single.py
@@ -19,7 +19,7 @@ from sklearn.linear_model import LogisticRegression
 from sklearn.model_selection import cross_val_score, KFold
 from sklearn.pipeline import Pipeline
 
-from pyriemann.classification import MDM, TSclassifier
+from pyriemann.classification import MDM, TSClassifier
 from pyriemann.estimation import Covariances
 
 
@@ -93,7 +93,7 @@ print("MDM Classification accuracy: %f / Chance level: %f" % (np.mean(scores),
 # Classification with Tangent Space Logistic Regression
 # -----------------------------------------------------
 
-clf = TSclassifier()
+clf = TSClassifier()
 # Use scikit-learn Pipeline with cross_val_score function
 scores = cross_val_score(clf, cov_data_train, labels, cv=cv, n_jobs=1)
 
diff --git a/examples/transfer/plot_bci_example.py b/examples/transfer/plot_bci_example.py
index 4f2820c..e295463 100644
--- a/examples/transfer/plot_bci_example.py
+++ b/examples/transfer/plot_bci_example.py
@@ -10,7 +10,7 @@ from several source subjects to train a single classifier using all of their
 data points pooled together. We compare the results of simply mixing all
 covariances from all source subjects without any care (dummy) versus
 transforming the covariances of all subjects so that they are centered around
-the Identity matrix (recenter) [1]_. We use data from the Physionet BCI
+the identity matrix (recenter) [1]_. We use data from the Physionet BCI
 database and compare the classification performance of MDM with each strategy.
 """
 
@@ -114,7 +114,11 @@ n_splits = 5  # How many times to split the target domain into train/test
 tl_cv = TLSplitter(
     target_domain="",
     cv=StratifiedShuffleSplit(
-        n_splits=n_splits, train_size=0.10, random_state=42))
+        n_splits=n_splits,
+        train_size=0.10,
+        random_state=42,
+    ),
+)
 
 # We consider two types of pipelines for transfer learning
 # dct : no transformation of dataset between the domains
@@ -140,7 +144,7 @@ for i_subject in tqdm(range(n_subjects)):
         X_enc_train, X_enc_test = X_enc[train_idx], X_enc[test_idx]
         y_enc_train, y_enc_test = y_enc[train_idx], y_enc[test_idx]
 
-        # (1) Dummy pipeline: no transfer learning at all.
+        # Dummy pipeline: no transfer learning at all.
         # Classifier is trained only with points from the source domain.
         domain_weight_dummy = {}
         for d in np.unique(domains):
@@ -160,7 +164,7 @@ for i_subject in tqdm(range(n_subjects)):
         pipeline.fit(X_enc_train, y_enc_train)
         scores_cv["dummy"].append(pipeline.score(X_enc_test, y_enc_test))
 
-        # (2) Recentering pipeline: recenter the data from each domain to
+        # Recentering pipeline: recenter the data from each domain to
         # identity [1]_.
         # Classifier is trained only with points from the source domain.
         domain_weight_rct = {}
diff --git a/examples/transfer/plot_comparison_pipelines.py b/examples/transfer/plot_comparison_pipelines.py
index 0149b1e..42db67e 100644
--- a/examples/transfer/plot_comparison_pipelines.py
+++ b/examples/transfer/plot_comparison_pipelines.py
@@ -3,50 +3,68 @@
 Comparison of pipelines for transfer learning
 ====================================================================
 
-We compare the classification performance of MDM on five different strategies
-for transfer learning. These include re-centering the datasets as done in [1]_,
-matching the statistical distributions in a semi-supervised way with Riemannian
-Procrustes Analysis [2]_, and improving the MDM classifier with a weighting
-scheme (MDWM) [3]_. All data points are simulated from a toy model based on the
+We compare the classification performance of MDM on different strategies
+for transfer learning.
+
+Matrices are simulated from a toy model based on the
 Riemannian Gaussian distribution and the differences in statistics between
 source and target distributions are determined by a set of parameters that have
 control over the distance between the centers of each dataset, the angle of
 rotation between the means of each class, and the differences in dispersion
-of the data points from each dataset.
+of the matrices from each dataset.
 """
 
 from tqdm import tqdm
 
-import numpy as np
 import matplotlib.pyplot as plt
-from sklearn.pipeline import make_pipeline
+import numpy as np
+from sklearn.linear_model import LogisticRegression
 from sklearn.model_selection import StratifiedShuffleSplit
+from sklearn.pipeline import make_pipeline
 
 from pyriemann.classification import MDM
 from pyriemann.datasets.simulated import make_classification_transfer
+from pyriemann.tangentspace import TangentSpace
 from pyriemann.transfer import (
     TLSplitter,
     TLDummy,
     TLCenter,
-    TLStretch,
+    TLScale,
     TLRotate,
     TLClassifier,
-    MDWM
+    MDWM,
 )
 
 
 ###############################################################################
+# Pipelines
+# ---------
+#
+# We consider several pipelines for transfer learning:
+#
+# - calib: use only data from target-train partition,
+#   classifier is trained only with matrices from the target domain.
+# - dummy: no transfer learning at all, ie no transformation of data between
+#   domains, classifier is trained only with matrices from the source domain.
+# - rct: recenter data from each domain to the identity matrix [1]_,
+#   classifier is trained only with matrices from the source domain.
+# - rpa: match the statistical distributions in a semi-supervised way with
+#   Riemannian Procrustes Analysis (RPA) [2]_: center, stretch and rotate
+#   matrices in manifold.
+#   Classifier is trained with matrices from source and target.
+# - mdwm: improve the MDM classifier with a weighting strategy, giving the
+#   minimum distance to weighted mean (MDWM) [3]_.
+# - tsa: align tangent vectors by Procrustes analysis [4]_: center,
+#   normalize and rotate vectors in tangent space.
+methods = ["calib", "dummy", "rct", "rpa", "mdwm", "tsa"]
+scores = {meth: [] for meth in methods}
+
+# Base classifier to consider in manifold
+clf_base = MDM()
+
 # Choose seed for reproducible results
 seed = 100
 
-# We consider several types of pipeline for transfer learning
-# dummy : no transformation of dataset between the domains
-# rct : re-center the data points from each domain to the Identity
-# rpa : re-center, stretch and rotate (Riemannian Procrustes Analysis)
-# mdwm : transfer learning with minimum distance to weighted mean
-# calibration : use only data from target-train partition for the classifier
-scores = {meth: [] for meth in ['dummy', 'rct', 'rpa', 'mdwm', 'calibration']}
-
 # Create a dataset with two domains, each with two classes both datasets
 # are generated by the same generative procedure with the SPD Gaussian
 # and one of them is transformed by a matrix A, i.e. X <- A @ X @ A.T
@@ -60,17 +78,16 @@ X_enc, y_enc = make_classification_transfer(
 )
 
 # Object for splitting the datasets into training and validation partitions
-# the training set is composed of all data points from the source domain
+# the training set is composed of all matrices from the source domain
 # plus a partition of the target domain whose size we can control
-target_domain = 'target_domain'
+target_domain = "target_domain"
 n_splits = 5  # how many times to split the target domain into train/test
 tl_cv = TLSplitter(
     target_domain=target_domain,
     cv=StratifiedShuffleSplit(n_splits=n_splits, random_state=seed),
 )
 
-# Which base classifier to consider
-clf_base = MDM()
+###############################################################################
 
 # Vary the proportion of the target domain for training
 target_train_frac_array = np.linspace(0.01, 0.20, 10)
@@ -89,74 +106,89 @@ for target_train_frac in tqdm(target_train_frac_array):
         X_enc_train, X_enc_test = X_enc[train_idx], X_enc[test_idx]
         y_enc_train, y_enc_test = y_enc[train_idx], y_enc[test_idx]
 
-        # (1) Dummy pipeline: no transfer learning at all.
-        # Classifier is trained only with samples from the source dataset.
+        # Calibration
+        pipeline = make_pipeline(
+            TLClassifier(
+                target_domain=target_domain,
+                estimator=clf_base,
+                domain_weight={"source_domain": 0.0, "target_domain": 1.0},
+            ),
+        )
+
+        pipeline.fit(X_enc_train, y_enc_train)
+        scores_cv["calib"].append(pipeline.score(X_enc_test, y_enc_test))
+
+        # Dummy
         pipeline = make_pipeline(
             TLDummy(),
             TLClassifier(
                 target_domain=target_domain,
                 estimator=clf_base,
-                domain_weight={'source_domain': 1.0, 'target_domain': 0.0},
+                domain_weight={"source_domain": 1.0, "target_domain": 0.0},
             ),
         )
 
-        # Fit and get accuracy score
         pipeline.fit(X_enc_train, y_enc_train)
-        scores_cv['dummy'].append(pipeline.score(X_enc_test, y_enc_test))
+        scores_cv["dummy"].append(pipeline.score(X_enc_test, y_enc_test))
 
-        # (2) RCT pipeline: recenter data from each domain to identity [1]_.
-        # Classifier is trained only with points from the source domain.
+        # Recentering pipeline
         pipeline = make_pipeline(
             TLCenter(target_domain=target_domain),
             TLClassifier(
                 target_domain=target_domain,
                 estimator=clf_base,
-                domain_weight={'source_domain': 1.0, 'target_domain': 0.0},
+                domain_weight={"source_domain": 1.0, "target_domain": 0.0},
             ),
         )
 
         pipeline.fit(X_enc_train, y_enc_train)
-        scores_cv['rct'].append(pipeline.score(X_enc_test, y_enc_test))
+        scores_cv["rct"].append(pipeline.score(X_enc_test, y_enc_test))
 
-        # (3) RPA pipeline: recenter, stretch, and rotate [2]_.
-        # Classifier is trained with points from source and target.
+        # RPA pipeline
         pipeline = make_pipeline(
             TLCenter(target_domain=target_domain),
-            TLStretch(
+            TLScale(
                 target_domain=target_domain,
                 final_dispersion=1,
                 centered_data=True,
             ),
-            TLRotate(target_domain=target_domain, metric='euclid'),
+            TLRotate(target_domain=target_domain, metric="euclid"),
             TLClassifier(
                 target_domain=target_domain,
                 estimator=clf_base,
-                domain_weight={'source_domain': 0.5, 'target_domain': 0.5},
+                domain_weight={"source_domain": 0.5, "target_domain": 0.5},
             ),
         )
 
         pipeline.fit(X_enc_train, y_enc_train)
-        scores_cv['rpa'].append(pipeline.score(X_enc_test, y_enc_test))
+        scores_cv["rpa"].append(pipeline.score(X_enc_test, y_enc_test))
+
+        # MDWM pipeline
+        domain_tradeoff = 1 - np.exp(-100 * target_train_frac)
+        pipeline = MDWM(
+            domain_tradeoff=domain_tradeoff,
+            target_domain=target_domain,
+            metric="riemann",
+        )
 
-        # (4) MDWM pipeline
-        domain_tradeoff = 1 - np.exp(-25*target_train_frac)
-        pipeline = MDWM(domain_tradeoff=domain_tradeoff,
-                        target_domain=target_domain)
         pipeline.fit(X_enc_train, y_enc_train)
-        scores_cv['mdwm'].append(pipeline.score(X_enc_test, y_enc_test))
+        scores_cv["mdwm"].append(pipeline.score(X_enc_test, y_enc_test))
 
-        # (5) Calibration: use only data from target-train partition.
-        # Classifier is trained only with points from the target domain.
+        # TSA pipeline
         pipeline = make_pipeline(
+            TangentSpace(metric="riemann"),
+            TLCenter(target_domain=target_domain),
+            TLScale(target_domain=target_domain),
+            TLRotate(target_domain=target_domain),
             TLClassifier(
                 target_domain=target_domain,
-                estimator=clf_base,
-                domain_weight={'source_domain': 0.0, 'target_domain': 1.0},
+                estimator=LogisticRegression(),
+                domain_weight={"source_domain": 0.5, "target_domain": 0.5},
             ),
         )
 
         pipeline.fit(X_enc_train, y_enc_train)
-        scores_cv['calibration'].append(pipeline.score(X_enc_test, y_enc_test))
+        scores_cv["tsa"].append(pipeline.score(X_enc_test, y_enc_test))
 
     # Get the average score of each pipeline
     for meth in scores.keys():
@@ -166,8 +198,10 @@ for target_train_frac in tqdm(target_train_frac_array):
 for meth in scores.keys():
     scores[meth] = np.array(scores[meth])
 
-
 ###############################################################################
+# Results
+# -------
+#
 # Plot the results, reproducing Figure 2 of [2]_.
 
 fig, ax = plt.subplots(figsize=(6.7, 5.7))
@@ -176,17 +210,17 @@ for meth in scores.keys():
         target_train_frac_array,
         scores[meth],
         label=meth,
-        lw=3.0)
-ax.legend(loc='lower right')
-ax.set_ylim(0.45, 0.75)
+        lw=3.0 if meth == "calib" else 2.0,
+    )
+ax.legend(loc="lower right")
+ax.set_ylim(0.5, 0.75)
 ax.set_yticks([0.5, 0.6, 0.7])
 ax.set_xlim(0.00, 0.21)
 ax.set_xticks([0.01, 0.05, 0.10, 0.15, 0.20])
 ax.set_xticklabels([1, 5, 10, 15, 20])
-ax.set_xlabel('Percentage of training partition in target domain')
-ax.set_ylabel('Classification accuracy')
-ax.set_title('Comparison of transfer learning pipelines')
-
+ax.set_xlabel("Percentage of training partition in target domain")
+ax.set_ylabel("Classification accuracy")
+ax.set_title("Comparison of transfer learning pipelines")
 plt.show()
 
 
@@ -208,3 +242,8 @@ plt.show()
 #    <https://hal.uvsq.fr/hal-01911092>`_
 #    E Kalunga et al, 26th European Signal Processing Conference (EUSIPCO 2018)
 #    Sep 2018, Rome, Italy, pp.1685-1689
+# .. [4] `Tangent space alignment: Transfer learning for brain-computer
+#    interface
+#    <https://www.frontiersin.org/articles/10.3389/fnhum.2022.1049985/pdf>`_
+#    A. Bleuzé, J. Mattout and M. Congedo, Frontiers in Human Neuroscience,
+#    2022
diff --git a/examples/transfer/plot_rpa_steps.py b/examples/transfer/plot_rpa_steps.py
index 8877d21..d0629ad 100644
--- a/examples/transfer/plot_rpa_steps.py
+++ b/examples/transfer/plot_rpa_steps.py
@@ -1,10 +1,10 @@
 """
 ====================================================================
-Plot the data transformations in the Riemannian Procrustes Analysis
+Data transformations in the Riemannian Procrustes Analysis
 ====================================================================
 
-Use the SpectralEmbedding module to plot in 2D the transformations on the data
-points from source and target domains when applying the Riemannian Procrustes
+Use the SpectralEmbedding module to plot in 2D the transformations on the
+matrices from source and target domains when applying the Riemannian Procrustes
 Analysis [1]_ to match their statistics.
 """
 
@@ -79,19 +79,23 @@ for axi, step, title in zip(ax, steps, titles):
     axi.scatter(
         S_source[y_source == "1"][:, 0],
         S_source[y_source == "1"][:, 1],
-        c="C0", s=50, alpha=0.50)
+        c="C0", s=50, alpha=0.50,
+    )
     axi.scatter(
         S_source[y_source == "2"][:, 0],
         S_source[y_source == "2"][:, 1],
-        c="C1", s=50, alpha=0.50)
+        c="C1", s=50, alpha=0.50,
+    )
     axi.scatter(
         S_target[y_target == "1"][:, 0],
         S_target[y_target == "1"][:, 1],
-        c="C0", s=50, alpha=0.30, marker="^")
+        c="C0", s=50, alpha=0.30, marker="^",
+    )
     axi.scatter(
         S_target[y_target == "2"][:, 0],
         S_target[y_target == "2"][:, 1],
-        c="C1", s=50, alpha=0.30, marker="^")
+        c="C1", s=50, alpha=0.30, marker="^",
+    )
     axi.scatter(S[-1, 0], S[-1, 1], c="k", s=80, marker="*")
     axi.set_xlim(-0.60, +1.60)
     axi.set_ylim(-1.10, +1.25)
@@ -103,7 +107,6 @@ ax[0].scatter([], [], c="C1", label="source - class 1")
 ax[0].scatter([], [], marker="^", c="C0", label="target - class 0")
 ax[0].scatter([], [], marker="^", c="C1", label="target - class 1")
 ax[0].legend(loc="upper right")
-
 plt.show()
 
 
diff --git a/pyriemann/classification.py b/pyriemann/classification.py
index 7cc1574..c9de2f6 100644
--- a/pyriemann/classification.py
+++ b/pyriemann/classification.py
@@ -9,6 +9,7 @@ from sklearn.utils.extmath import softmax
 from sklearn.linear_model import LogisticRegression
 from sklearn.pipeline import make_pipeline
 
+from .utils import deprecated
 from .utils.kernel import kernel
 from .utils.mean import mean_covariance
 from .utils.distance import distance
@@ -341,7 +342,7 @@ class FgMDM(BaseEstimator, ClassifierMixin, TransformerMixin):
         return self._mdm.transform(cov)
 
 
-class TSclassifier(BaseEstimator, ClassifierMixin):
+class TSClassifier(BaseEstimator, ClassifierMixin):
     """Classification in the tangent space.
 
     Project data in the tangent space and apply a classifier on the projected
@@ -396,7 +397,7 @@ class TSclassifier(BaseEstimator, ClassifierMixin):
         self.clf = clf
 
     def fit(self, X, y, sample_weight=None):
-        """Fit TSclassifier.
+        """Fit TsClassifier.
 
         Parameters
         ----------
@@ -409,8 +410,8 @@ class TSclassifier(BaseEstimator, ClassifierMixin):
 
         Returns
         -------
-        self : TSclassifier instance
-            The TSclassifier instance.
+        self : TSClassifier instance
+            The TSClassifier instance.
         """
         if not isinstance(self.clf, ClassifierMixin):
             raise TypeError("clf must be a ClassifierMixin")
@@ -456,6 +457,14 @@ class TSclassifier(BaseEstimator, ClassifierMixin):
         return self._pipe.predict_proba(X)
 
 
+@deprecated(
+    "TSclassifier is deprecated and will be removed in 0.10.0; "
+    "please use TSClassifier."
+)
+class TSclassifier(TSClassifier):
+    pass
+
+
 class KNearestNeighbor(MDM):
     """Classification by k-nearest neighbors.
 
diff --git a/pyriemann/transfer/__init__.py b/pyriemann/transfer/__init__.py
index 01d8f1f..1ee71ea 100644
--- a/pyriemann/transfer/__init__.py
+++ b/pyriemann/transfer/__init__.py
@@ -7,7 +7,7 @@ from ._tools import (
 from ._estimators import (
     TLDummy,
     TLCenter,
-    TLStretch,
+    TLScale,
     TLRotate,
     TLEstimator,
     TLClassifier,
@@ -18,10 +18,10 @@ from ._estimators import (
 __all__ = [
     "encode_domains",
     "decode_domains",
-    "TLDummy",
     "TLSplitter",
+    "TLDummy",
     "TLCenter",
-    "TLStretch",
+    "TLScale",
     "TLRotate",
     "TLEstimator",
     "TLClassifier",
diff --git a/pyriemann/transfer/_estimators.py b/pyriemann/transfer/_estimators.py
index 594f9d4..8310c20 100644
--- a/pyriemann/transfer/_estimators.py
+++ b/pyriemann/transfer/_estimators.py
@@ -1,48 +1,57 @@
-import numpy as np
+import warnings
+
 from joblib import Parallel, delayed
+import numpy as np
 from sklearn.base import (
     BaseEstimator,
     TransformerMixin,
     is_classifier,
     is_regressor
 )
-from sklearn.pipeline import Pipeline
+from sklearn.decomposition import PCA
 from sklearn.metrics import accuracy_score, r2_score
+from sklearn.pipeline import Pipeline
 
-from ..utils.mean import mean_covariance, mean_riemann
-from ..utils.distance import distance
+from ._rotate import _get_rotation_manifold, _get_rotation_tangentspace
+from ._tools import decode_domains
+from ..classification import MDM
+from ..preprocessing import Whitening
+from ..utils import deprecated
 from ..utils.base import invsqrtm, powm, sqrtm
+from ..utils.distance import distance
 from ..utils.geodesic import geodesic
+from ..utils.mean import mean_covariance, mean_riemann
 from ..utils.utils import check_weights, check_metric
-from ._rotate import _get_rotation_matrix
-from ..classification import MDM
-from ..preprocessing import Whitening
-from ._tools import decode_domains
+
+
+###############################################################################
+
+
+def _check_inputs(X):
+    if X.ndim not in [2, 3]:
+        raise ValueError(f"Input must be a 2d or a 3d array (Got {X.ndim}).")
 
 
 class TLDummy(BaseEstimator, TransformerMixin):
-    """No transformation on data for transfer learning.
+    """No transformation for transfer learning.
 
-    No transformation of the data points between the domains.
-    This is what we call the Direct Center Transfer (DCT) method.
+    No transformation of data between the domains.
 
     Notes
     -----
     .. versionadded:: 0.4
     """
 
-    def __init__(self):
-        pass
-
-    def fit(self, X, y_enc):
+    def fit(self, X, y_enc=None):
         """Do nothing.
 
         Parameters
         ----------
-        X : ndarray, shape (n_matrices, n_channels, n_channels)
-            Set of SPD matrices.
-        y_enc : ndarray, shape (n_matrices,)
-            Extended labels for each matrix.
+        X : ndarray, shape (n_matrices, n_channels, n_channels) or \
+                shape (n_vectors, n_ts)
+            Set of SPD matrices or tangent vectors.
+        y_enc : None
+            Not used, here for compatibility with sklearn API.
 
         Returns
         -------
@@ -51,54 +60,61 @@ class TLDummy(BaseEstimator, TransformerMixin):
         """
         return self
 
-    def transform(self, X, y_enc=None):
+    def transform(self, X):
         """Do nothing.
 
         Parameters
         ----------
-        X : ndarray, shape (n_matrices, n_channels, n_channels)
-            Set of SPD matrices.
-        y_enc : None
-            Not used, here for compatibility with sklearn API.
+        X : ndarray, shape (n_matrices, n_channels, n_channels) or \
+                shape (n_vectors, n_ts)
+            Set of SPD matrices or tangent vectors.
 
         Returns
         -------
-        X : ndarray, shape (n_matrices, n_classes)
-            Same set of SPD matrices as in the input.
+        X_new : ndarray, shape (n_matrices, n_channels, n_channels) or \
+                shape (n_vectors, n_ts)
+            Same data as in the input.
         """
         return X
 
-    def fit_transform(self, X, y_enc):
+    def fit_transform(self, X, y_enc=None):
         """Do nothing.
 
         Parameters
         ----------
-        X : ndarray, shape (n_matrices, n_channels, n_channels)
-            Set of SPD matrices.
-        y_enc : ndarray, shape (n_matrices,)
-            Extended labels for each matrix.
+        X : ndarray, shape (n_matrices, n_channels, n_channels) or \
+                shape (n_vectors, n_ts)
+            Set of SPD matrices or tangent vectors.
+        y_enc : None
+            Not used, here for compatibility with sklearn API.
 
         Returns
         -------
-        X : ndarray, shape (n_matrices, n_classes)
-            Set of SPD matrices with mean in the Identity.
+        X_new : ndarray, shape (n_matrices, n_channels, n_channels) or \
+                shape (n_vectors, n_ts)
+            Same data as in the input.
         """
-        return self.fit(X, y_enc).transform(X, y_enc)
+        return self.fit(X, y_enc).transform(X)
 
 
 class TLCenter(BaseEstimator, TransformerMixin):
-    """Recenter data for transfer learning.
+    """Centering for transfer learning.
+
+    For inputs in matrix manifold, it recenters the matrices from each domain
+    to the identity matrix on manifold, ie it makes the mean of the matrices of
+    each domain become the identity [1]_.
+    This operation corresponds to a whitening when the matrices represent the
+    spatial covariance matrices of multivariate signals.
 
-    Recenter the data points from each domain to the Identity on manifold, ie
-    make the mean of the datasets become the identity. This operation
-    corresponds to a whitening step if the SPD matrices represent the spatial
-    covariance matrices of multivariate signals.
+    For inputs in tangent space, it recenters the tangent vectors from each
+    domain to the origin of tangent space, ie it makes the mean of the vectors
+    of each domain become zero.
 
     .. note::
        Using .fit() and then .transform() will give different results than
        .fit_transform(). In fact, .fit_transform() should be applied on the
-       training dataset (target and source) and .transform() on the test
-       partition of the target dataset.
+       training set (target and source domains),
+       and .transform() on the target domain of the test set.
 
     Parameters
     ----------
@@ -109,15 +125,22 @@ class TLCenter(BaseEstimator, TransformerMixin):
           target domain;
         * else, ``transform()`` recenters matrices to the last fitted domain.
     metric : str, default="riemann"
-        Metric used for mean estimation. For the list of supported metrics,
+        For inputs in manifold,
+        metric used for mean estimation. For the list of supported metrics,
         see :func:`pyriemann.utils.mean.mean_covariance`.
         Note, however, that only when using the "riemann" metric that we are
-        ensured to re-center the matrices precisely to the Identity.
+        ensured to re-center the matrices precisely to the identity.
 
     Attributes
     ----------
-    recenter_ : dict
-        If fit, dictionary with key=domain_name and value=domain_mean.
+    centers_ : dict
+        Dictionary with key=domain_name and value=domain_center.
+
+    Notes
+    -----
+    .. versionadded:: 0.4
+    .. versionchanged:: 0.8
+        Added support for tangent space centering.
 
     References
     ----------
@@ -130,10 +153,6 @@ class TLCenter(BaseEstimator, TransformerMixin):
         A Euclidean Space Data Alignment Approach
         <https://arxiv.org/abs/1808.05464>`_
         He He and Dongrui Wu, IEEE Transactions on Biomedical Engineering, 2019
-
-    Notes
-    -----
-    .. versionadded:: 0.4
     """
 
     def __init__(self, target_domain, metric="riemann"):
@@ -141,130 +160,182 @@ class TLCenter(BaseEstimator, TransformerMixin):
         self.target_domain = target_domain
         self.metric = metric
 
+    @property
+    @deprecated(
+        "Attribute `recenter_` is deprecated and will be removed in 0.10.0; "
+        "please use `centers_`."
+    )
+    def recenter_(self):
+        return self.centers_
+
     def fit(self, X, y_enc, sample_weight=None):
         """Fit TLCenter.
 
-        For each domain, calculates the mean of matrices of this domain.
+        For each domain, it calculates the mean of matrices or vectors of this
+        domain.
 
         Parameters
         ----------
-        X : ndarray, shape (n_matrices, n_channels, n_channels)
-            Set of SPD matrices.
-        y_enc : ndarray, shape (n_matrices,)
-            Extended labels for each matrix.
-        sample_weight : None | ndarray, shape (n_matrices,), default=None
-            Weights for each matrix. If None, it uses equal weights.
+        X : ndarray, shape (n_matrices, n_channels, n_channels) or \
+                shape (n_vectors, n_ts)
+            Set of SPD matrices or tangent vectors.
+        y_enc : ndarray, shape (n_matrices,) or shape (n_vectors,)
+            Extended labels for each matrix or vector.
+        sample_weight : None | ndarray, shape (n_matrices,) or \
+                shape (n_vectors,), default=None
+            Weights for each matrix or vector. If None, it uses equal weights.
 
         Returns
         -------
         self : TLCenter instance
             The TLCenter instance.
         """
+        _check_inputs(X)
         _, _, domains = decode_domains(X, y_enc)
-        n_matrices, _, _ = X.shape
-        sample_weight = check_weights(sample_weight, n_matrices)
+        sample_weight = check_weights(sample_weight, X.shape[0])
+
+        self.centers_ = {}
 
-        self.recenter_ = {}
         for d in np.unique(domains):
             idx = domains == d
-            self.recenter_[d] = Whitening(metric=self.metric).fit(
-                X[idx], sample_weight=sample_weight[idx]
-            )
+
+            if X.ndim == 3:
+                self.centers_[d] = Whitening(metric=self.metric).fit(
+                    X[idx],
+                    sample_weight=sample_weight[idx],
+                )
+
+            else:
+                self.centers_[d] = np.mean(X[idx], axis=0)
+
         return self
 
-    def transform(self, X, y_enc=None):
-        """Re-center matrices in the target domain.
+    def transform(self, X):
+        """Center in the target domain.
 
         .. note::
            This method is designed for using at test time,
-           recentering all matrices in target domain.
+           recentering all inputs in target domain, or in the last fitted
+           domain.
 
         Parameters
         ----------
-        X : ndarray, shape (n_matrices, n_channels, n_channels)
-            Set of SPD matrices.
-        y_enc : None
-            Not used, here for compatibility with sklearn API.
+        X : ndarray, shape (n_matrices, n_channels, n_channels) or \
+                shape (n_vectors, n_ts)
+            Set of SPD matrices or tangent vectors.
 
         Returns
         -------
-        X : ndarray, shape (n_matrices, n_classes)
-            Set of recentered SPD matrices.
+        X_new : ndarray, shape (n_matrices, n_classes) or \
+                shape (n_vectors, n_ts)
+            Set of centered SPD matrices or tangent vectors in target domain.
         """
+        _check_inputs(X)
+
         # if target domain is specified, use it
         if self.target_domain != "":
             target_domain = self.target_domain
         # else, use last calibrated domain as target domain
         else:
-            target_domain = list(self.recenter_.keys())[-1]
+            target_domain = list(self.centers_.keys())[-1]
+
+        if X.ndim == 3:
+            X_new = self.centers_[target_domain].transform(X)
+        else:
+            X_new = X - self.centers_[self.target_domain]
 
-        X_rct = self.recenter_[target_domain].transform(X)
-        return X_rct
+        return X_new
 
     def fit_transform(self, X, y_enc, sample_weight=None):
-        """Fit TLCenter and then transform matrices.
+        """Fit TLCenter and then center each domain.
 
-        For each domain, calculates the mean of matrices of this domain and
-        then recenters them to Identity.
+        For each domain, it calculates the mean of matrices or vectors of this
+        domain, and then recenters them to identity matrix or to null vector.
 
         .. note::
-           This method is designed for using at training time. The output for
-           .fit_transform() will be different than using .fit() and
-           .transform() separately.
+           This method is designed for using at training time.
+           The output for .fit_transform() will be different
+           than using .fit() and .transform() separately.
 
         Parameters
         ----------
-        X : ndarray, shape (n_matrices, n_channels, n_channels)
-            Set of SPD matrices.
-        y_enc : ndarray, shape (n_matrices,)
-            Extended labels for each matrix.
-        sample_weight : None | ndarray, shape (n_matrices,), default=None
-            Weights for each matrix. If None, it uses equal weights.
+        X : ndarray, shape (n_matrices, n_channels, n_channels) or \
+                shape (n_vectors, n_ts)
+            Set of SPD matrices or tangent vectors.
+        y_enc : ndarray, shape (n_matrices,) or shape (n_vectors,)
+            Extended labels for each matrix or vector.
+        sample_weight : None | ndarray, shape (n_matrices,) or \
+                shape (n_vectors,), default=None
+            Weights for each matrix or vector. If None, it uses equal weights.
 
         Returns
         -------
-        X : ndarray, shape (n_matrices, n_classes)
-            Set of recentered SPD matrices.
+        X_new : ndarray, shape (n_matrices, n_channels, n_channels) or \
+                shape (n_vectors, n_ts)
+            Set of centered SPD matrices or tangent vectors in each domain.
         """
-        self.fit(X, y_enc, sample_weight)
+        self.fit(X, y_enc, sample_weight=sample_weight)
         _, _, domains = decode_domains(X, y_enc)
 
-        X_rct = np.zeros_like(X)
+        X_new = np.zeros_like(X)
         for d in np.unique(domains):
             idx = domains == d
-            X_rct[idx] = self.recenter_[d].transform(X[idx])
-        return X_rct
 
+            if X.ndim == 3:
+                X_new[idx] = self.centers_[d].transform(X[idx])
+            else:
+                X_new[idx] = X[idx] - self.centers_[d]
+
+        return X_new
 
-class TLStretch(BaseEstimator, TransformerMixin):
-    """Stretch data for transfer learning.
 
-    Change the dispersion of the datapoints around their geometric mean
-    for each dataset so that they all have the same desired value.
+class TLScale(BaseEstimator, TransformerMixin):
+    """Scaling for transfer learning.
+
+    For inputs in matrix manifold, it stretches the matrices from each domain
+    around their mean so that the dispersion of the matrices of each domain is
+    equal to one [1]_.
+
+    For inputs in tangent space, it scales the tangent vectors from each domain
+    so that the mean of norms of vectors of each domain is equal to one.
 
     .. note::
        Using .fit() and then .transform() will give different results than
        .fit_transform(). In fact, .fit_transform() should be applied on the
-       training dataset (target and source) and .transform() on the test
-       partition of the target dataset.
+       training set (target and source domains),
+       and .transform() on the target domain of the test set.
 
     Parameters
     ----------
     target_domain : str
         Domain to consider as target.
     dispersion : float, default=1.0
-        Target value for the dispersion of the data points.
+        For inputs in manifold, target value for the dispersion of the
+        matrices.
     centered_data : bool, default=False
-        Whether the data has been re-centered to the Identity beforehand.
+        For inputs in manifold, whether the matrices have been re-centered to
+        the identity matrix beforehand.
     metric : str, default="riemann"
-        Metric used for calculating the dispersion.
+        For inputs in manifold, metric used for calculating the dispersion.
         For the list of supported metrics,
         see :func:`pyriemann.utils.distance.distance`.
+        The stretching operation in manifold is properly defined only for the
+        "riemann" metric.
 
     Attributes
     ----------
-    dispersions_ : dict
-        Dictionary with key=domain_name and value=domain_dispersion.
+    scales_ : dict
+        Dictionary with key=domain_name and value=domain_scale.
+
+    See Also
+    --------
+    TLCenter
+
+    Notes
+    -----
+    .. versionadded:: 0.4
+    .. versionchanged:: 0.8
+        Added support for tangent space scaling.
 
     References
     ----------
@@ -273,10 +344,6 @@ class TLStretch(BaseEstimator, TransformerMixin):
         <https://hal.archives-ouvertes.fr/hal-01971856>`_
         PLC Rodrigues et al, IEEE Transactions on Biomedical Engineering,
         vol. 66, no. 8, pp. 2390-2401, December, 2018
-
-    Notes
-    -----
-    .. versionadded:: 0.4
     """
 
     def __init__(
@@ -292,46 +359,63 @@ class TLStretch(BaseEstimator, TransformerMixin):
         self.centered_data = centered_data
         self.metric = metric
 
+    @property
+    @deprecated(
+        "Attribute `dispersions_` is deprecated and will be removed in 0.10.0;"
+        " please use `scales_`."
+    )
+    def dispersions_(self):
+        return self.scales_
+
     def fit(self, X, y_enc, sample_weight=None):
-        """Fit TLStretch.
+        """Fit TLScale.
 
-        Calculate the dispersion around the mean for each domain.
+        For each domain, it calculates the scaling of this domain,
+        ie the dispersion around the mean of matrices,
+        or the mean of the norm of vectors.
 
         Parameters
         ----------
-        X : ndarray, shape (n_matrices, n_channels, n_channels)
-            Set of SPD matrices.
-        y_enc : ndarray, shape (n_matrices,)
-            Extended labels for each matrix.
-        sample_weight : None | ndarray, shape (n_matrices,), default=None
-            Weights for each matrix. If None, it uses equal weights.
+        X : ndarray, shape (n_matrices, n_channels, n_channels) or \
+                shape (n_vectors, n_ts)
+            Set of SPD matrices or tangent vectors.
+        y_enc : ndarray, shape (n_matrices,) or shape (n_vectors,)
+            Extended labels for each matrix or vector.
+        sample_weight : None | ndarray, shape (n_matrices,) or \
+                shape (n_vectors,), default=None
+            Weights for each matrix or vector. If None, it uses equal weights.
 
         Returns
         -------
-        self : TLStretch instance
-            The TLStretch instance.
+        self : TLScale instance
+            The TLScale instance.
         """
+        _check_inputs(X)
         _, _, domains = decode_domains(X, y_enc)
-        n_matrices, n_channels, _ = X.shape
-        sample_weight = check_weights(sample_weight, n_matrices)
+        sample_weight = check_weights(sample_weight, X.shape[0])
 
-        self._means, self.dispersions_ = {}, {}
+        self._means, self.scales_ = {}, {}
         for d in np.unique(domains):
             idx = domains == d
             sample_weight_d = check_weights(sample_weight[idx], np.sum(idx))
-            if self.centered_data:
-                self._means[d] = np.eye(n_channels)
-            else:
-                self._means[d] = mean_riemann(
-                    X[idx], sample_weight=sample_weight_d
+
+            if X.ndim == 3:
+                if self.centered_data:
+                    self._means[d] = np.eye(X.shape[-1])
+                else:
+                    self._means[d] = mean_riemann(
+                        X[idx], sample_weight=sample_weight_d
+                    )
+                dist = distance(
+                    X[idx],
+                    self._means[d],
+                    metric=self.metric,
+                    squared=True,
                 )
-            dist = distance(
-                X[idx],
-                self._means[d],
-                metric=self.metric,
-                squared=True,
-            )
-            self.dispersions_[d] = np.sum(sample_weight_d * np.squeeze(dist))
+                self.scales_[d] = np.sum(sample_weight_d * np.squeeze(dist))
+
+            else:
+                self.scales_[d] = np.mean(np.linalg.norm(X[idx], axis=1))
 
         return self
 
@@ -346,128 +430,169 @@ class TLStretch(BaseEstimator, TransformerMixin):
     def _strech(self, X, dispersion_in, dispersion_out):
         return powm(X, np.sqrt(dispersion_out / dispersion_in))
 
-    def transform(self, X, y_enc=None):
-        """Stretch the data points in the target domain.
+    def transform(self, X):
+        """Scale in the target domain.
 
         .. note::
-           The stretching operation is properly defined only for the riemann
-           metric.
+           This method is designed for using at test time,
+           scaling all inputs in target domain.
 
         Parameters
         ----------
-        X : ndarray, shape (n_matrices, n_channels, n_channels)
-            Set of SPD matrices.
-        y_enc : None
-            Not used, here for compatibility with sklearn API.
+        X : ndarray, shape (n_matrices, n_channels, n_channels) or \
+                shape (n_vectors, n_ts)
+            Set of SPD matrices or tangent vectors.
 
         Returns
         -------
-        X : ndarray, shape (n_matrices, n_classes)
-            Set of SPD matrices with desired final dispersion.
+        X_new : ndarray, shape (n_matrices, n_classes) or \
+                shape (n_vectors, n_ts)
+            Set of scaled SPD matrices or tangent vectors in target domain.
         """
+        _check_inputs(X)
 
-        if not self.centered_data:
-            # center matrices to Identity
-            X = self._center(X, self._means[self.target_domain])
+        if X.ndim == 3:
+            if not self.centered_data:  # center matrices to identity
+                X = self._center(X, self._means[self.target_domain])
 
-        # stretch
-        X_str = self._strech(
-            X, self.dispersions_[self.target_domain], self.final_dispersion
-        )
+            # stretch
+            X_new = self._strech(
+                X, self.scales_[self.target_domain], self.final_dispersion
+            )
+
+            if not self.centered_data:  # re-center back to previous mean
+                X_new = self._uncenter(X_new, self._means[self.target_domain])
 
-        if not self.centered_data:
-            # re-center back to previous mean
-            X_str = self._uncenter(X_str, self._means[self.target_domain])
+        else:
+            X_new = X / self.scales_[self.target_domain]
 
-        return X_str
+        return X_new
 
     def fit_transform(self, X, y_enc, sample_weight=None):
-        """Fit TLStretch and then transform data points.
+        """Fit TLScale and then scale each domain.
 
-        Calculate the dispersion around the mean for each domain and then
-        stretch the data points to the desired final dispersion.
+        For each domain, it calculates the dispersion around the mean of this
+        domain, and then stretches them to the desired final dispersion.
+        For vectors, it scales them so that the mean of norms of vectors of
+        each domain is equal to one.
 
         .. note::
-           This method is designed for using at training time. The output for
-           .fit_transform() will be different than using .fit() and
-           .transform() separately.
+           This method is designed for using at training time.
+           The output for .fit_transform() will be different
+           than using .fit() and .transform() separately.
 
         Parameters
         ----------
-        X : ndarray, shape (n_matrices, n_channels, n_channels)
-            Set of SPD matrices.
-        y_enc : ndarray, shape (n_matrices,)
-            Extended labels for each matrix.
-        sample_weight : None | ndarray, shape (n_matrices,), default=None
-            Weights for each matrix. If None, it uses equal weights.
+        X : ndarray, shape (n_matrices, n_channels, n_channels) or \
+                shape (n_vectors, n_ts)
+            Set of SPD matrices or tangent vectors.
+        y_enc : ndarray, shape (n_matrices,) or shape (n_vectors,)
+            Extended labels for each matrix or vector.
+        sample_weight : None | ndarray, shape (n_matrices,) or \
+                shape (n_vectors,), default=None
+            Weights for each matrix or vector. If None, it uses equal weights.
 
         Returns
         -------
-        X : ndarray, shape (n_matrices, n_classes)
-            Set of SPD matrices with desired final dispersion.
+        X_new : ndarray, shape (n_matrices, n_channels, n_channels) or \
+                shape (n_vectors, n_ts)
+            Set of scaled SPD matrices or tangent vectors in each domain.
         """
-
-        # used during fit, in pipeline
-        self.fit(X, y_enc, sample_weight)
+        self.fit(X, y_enc, sample_weight=sample_weight)
         _, _, domains = decode_domains(X, y_enc)
 
-        X_str = np.zeros_like(X)
+        X_new = np.zeros_like(X)
         for d in np.unique(domains):
             idx = domains == d
 
-            if not self.centered_data:
-                # re-center matrices to Identity
-                X[idx] = self._center(X[idx], self._means[d])
+            if X.ndim == 3:
+                if not self.centered_data:  # re-center matrices to identity
+                    X[idx] = self._center(X[idx], self._means[d])
 
-            # stretch
-            X_str[idx] = self._strech(
-                X[idx], self.dispersions_[d], self.final_dispersion
-            )
+                # stretch
+                X_new[idx] = self._strech(
+                    X[idx], self.scales_[d], self.final_dispersion
+                )
+
+                if not self.centered_data:  # re-center back to previous mean
+                    X_new[idx] = self._uncenter(X_new[idx], self._means[d])
+
+            else:
+                X_new[idx] = X[idx] / self.scales_[d]
+
+        return X_new
 
-            if not self.centered_data:
-                # re-center back to previous mean
-                X_str[idx] = self._uncenter(X_str[idx], self._means[d])
 
-        return X_str
+@deprecated(
+    "TLStretch is deprecated and will be removed in 0.10.0; "
+    "please use TLScale."
+)
+class TLStretch(TLScale):
+    pass
 
 
 class TLRotate(BaseEstimator, TransformerMixin):
-    """Rotate data for transfer learning.
+    """Rotation for transfer learning.
+
+    For inputs in matrix manifold, it rotates the matrices from each source
+    domain so to match its class means with those from the target domain.
+    The loss function for this matching is described in [1]_ and the
+    optimization procedure for minimizing it in [2]_.
 
-    Rotate the data points from each source domain so to match its class means
-    with those from the target domain. The loss function for this matching was
-    first proposed in [1]_ and the optimization procedure for mininimizing it
-    follows the presentation from [2]_.
+    For inputs in tangent space, it rotates the tangent vectors from source
+    domain so to match its class means with those from the target domain [3]_.
+    Current implementation supports only one source domain.
 
     .. note::
-       The data points from each domain must have been re-centered to the
-       identity before calculating the rotation.
+       The inputs from each domain must have been centered to the before
+       calculating the rotation.
 
     .. note::
        Using .fit() and then .transform() will give different results than
        .fit_transform(). In fact, .fit_transform() should be applied on the
-       training dataset (target and source) and .transform() on the test
-       partition of the target dataset.
+       training set (target and source domains),
+       and .transform() on the target domain of the test set.
 
     Parameters
     ----------
     target_domain : str
         Domain to consider as target.
     weights : None | array, shape (n_classes,), default=None
-        Weights to assign for each class. If None, then give the same weight
-        for each class.
+        Weights to assign for each class. If None, it uses equal weights.
     metric : {"euclid", "riemann"}, default="euclid"
-        Metric for the distance to minimize between class means.
+        For inputs in manifold, distance to minimize between class means.
     n_jobs : int, default=1
-        The number of jobs to use for the computation. This works by computing
-        the rotation matrix for each source domain in parallel. If -1 all CPUs
-        are used.
+        For inputs in manifold, the number of jobs to use for the computation.
+        This works by computing the rotation matrix for each source domain in
+        parallel. If -1 all CPUs are used.
+    expl_var : float, default=0.999
+        For inputs in tangent space, dimension reduction applied to the cross
+        product matrix during Procrustes analysis.
+        If float in (0,1], percentage of variance that needs to be explained.
+        Else, number of components.
+    n_components : int | "max", default=1
+        For inputs in tangent space,
+        parameter `n_components` used in `sklearn.decomposition.PCA`.
+        If int, number of components to keep in PCA.
+        If "max", all components are kept.
+    n_clusters : int, default=3
+        For inputs in tangent space, number of clusters used to split data.
 
     Attributes
     ----------
     rotations_ : dict
         Dictionary with key=domain_name and value=domain_rotation_matrix.
 
+    See Also
+    --------
+    TLCenter
+
+    Notes
+    -----
+    .. versionadded:: 0.4
+    .. versionchanged:: 0.8
+        Added support for tangent space rotation.
+
     References
     ----------
     .. [1] `Riemannian Procrustes analysis: transfer learning for
@@ -478,56 +603,89 @@ class TLRotate(BaseEstimator, TransformerMixin):
     .. [2] `An introduction to optimization on smooth manifolds
         <https://www.nicolasboumal.net/book/>`_
         N. Boumal. To appear with Cambridge University Press. June, 2022
-
-    Notes
-    -----
-    .. versionadded:: 0.4
+    .. [3] `Tangent space alignment: Transfer learning for brain-computer
+        interface
+        <https://www.frontiersin.org/articles/10.3389/fnhum.2022.1049985/pdf>`_
+        A. Bleuzé, J. Mattout and M. Congedo, Frontiers in Human Neuroscience,
+        2022
     """
 
-    def __init__(self, target_domain, weights=None, metric="euclid", n_jobs=1):
+    def __init__(
+        self,
+        target_domain,
+        weights=None,
+        metric="euclid",
+        n_jobs=1,
+        expl_var=0.999,
+        n_components=1,
+        n_clusters=3,
+    ):
         """Init"""
         self.target_domain = target_domain
         self.weights = weights
         self.metric = metric
         self.n_jobs = n_jobs
+        self.expl_var = expl_var
+        self.n_components = n_components
+        self.n_clusters = n_clusters
 
     def fit(self, X, y_enc, sample_weight=None):
         """Fit TLRotate.
 
-        Calculate the rotations matrices to transform each source domain into
-        the target domain.
+        It calculates the rotations matrices to transform each source domain
+        into the target domain.
 
         Parameters
         ----------
-        X : ndarray, shape (n_matrices, n_channels, n_channels)
-            Set of SPD matrices.
-        y_enc : ndarray, shape (n_matrices,)
-            Extended labels for each matrix.
-        sample_weight : None | ndarray, shape (n_matrices,), default=None
-            Weights for each matrix. If None, it uses equal weights.
+        X : ndarray, shape (n_matrices, n_channels, n_channels) or \
+                shape (n_vectors, n_ts)
+            Set of SPD matrices or tangent vectors.
+        y_enc : ndarray, shape (n_matrices,) or shape (n_vectors,)
+            Extended labels for each matrix or vector.
+        sample_weight : None | ndarray, shape (n_matrices,) or \
+                shape (n_vectors,), default=None
+            Weights for each matrix or vector. If None, it uses equal weights.
 
         Returns
         -------
         self : TLRotate instance
             The TLRotate instance.
         """
-
+        _check_inputs(X)
         _, _, domains = decode_domains(X, y_enc)
-        n_matrices, _, _ = X.shape
-        sample_weight = check_weights(sample_weight, n_matrices)
+        sample_weight = check_weights(sample_weight, X.shape[0])
+
+        if X.ndim == 3:
+            self._fit_manifold(X, y_enc, domains, sample_weight)
+        else:
+            if len(np.unique(domains)) > 2:
+                raise ValueError(
+                    "Current implementation supports only one source domain"
+                )
+            self._fit_tangentspace(X, y_enc, domains)
+
+        return self
 
+    def _fit_manifold(self, X, y_enc, domains, sample_weight):
+        """Fit TLRotate on manifold.
+
+        It computes the mean M of SPD matrices X for each class of each domain.
+        It computes the rotation for each source domain so that its class means
+        are close to those from the target domain.
+        """
         idx = domains == self.target_domain
         X_target, y_target = X[idx], y_enc[idx]
         M_target = np.stack([
-            mean_riemann(X_target[y_target == label],
-                         sample_weight=sample_weight[idx][y_target == label])
-            for label in np.unique(y_target)
+            mean_riemann(
+                X_target[y_target == label],
+                sample_weight=sample_weight[idx][y_target == label],
+            ) for label in np.unique(y_target)
         ])
 
         source_domains = np.unique(domains)
         source_domains = source_domains[source_domains != self.target_domain]
         rotations = Parallel(n_jobs=self.n_jobs)(
-            delayed(_get_rotation_matrix)(
+            delayed(_get_rotation_manifold)(
                 np.stack([
                     mean_riemann(
                         X[domains == d][y_enc[domains == d] == label],
@@ -546,75 +704,159 @@ class TLRotate(BaseEstimator, TransformerMixin):
         for d, rot in zip(source_domains, rotations):
             self.rotations_[d] = rot
 
-        return self
+    def _fit_tangentspace(self, X, y_enc, domains):
+        """Fit TLRotate in tangent space.
+
+        It computes anchors for each class of each domain.
+        It computes the rotation for source domain so that its anchors
+        are close to those from the target domain.
+        """
+        _, y, _ = decode_domains(X, y_enc)
+        X_src = X[domains != self.target_domain]
+        y_src = y[domains != self.target_domain]
+        X_tgt = X[domains == self.target_domain]
+        y_tgt = y[domains == self.target_domain]
+
+        if len(np.unique(y_src)) != len(np.unique(y_tgt)):
+            raise ValueError(
+                "The number of classes in each domain don't match"
+            )
+
+        if self.n_components == "max":
+            self.n_components = X.shape[1]
+        self._src_pca = [
+            PCA(n_components=self.n_components) for _ in np.unique(y_src)
+        ]
+
+        src_anchors, src_valid = self._get_anchors(X_src, y_src, fit_pca=True)
+        tgt_anchors, tgt_valid = self._get_anchors(X_tgt, y_tgt, fit_pca=False)
+
+        if not src_valid or not tgt_valid:
+            n_classes = len(np.unique(y_tgt))
+            src_anchors = src_anchors[:n_classes]
+            tgt_anchors = tgt_anchors[:n_classes]
+            if not src_valid:
+                warnings.warn("Not enough vectors for source domain")
+            if not tgt_valid:
+                warnings.warn("Not enough vectors for target domain")
+
+        self.rotations_ = _get_rotation_tangentspace(
+            src_anchors,
+            tgt_anchors,
+            self.expl_var,
+        )
+
+    def _get_anchors(self, X, y, fit_pca):
+        """Get anchors.
+
+        It computes anchors for each class of each domain, ie class means and
+        clusters along principal components.
+        """
+        is_valid = True
+
+        classes = np.unique(y)
+        anchors = np.array([
+            np.mean(X[y == c], axis=0) for c in classes
+        ])
+
+        for i, c in enumerate(classes):
+            Xc = X[y == c]
+            n_vectors_c = len(Xc)
 
-    def transform(self, X, y_enc=None):
-        """Rotate the data points in the target domain.
+            if n_vectors_c < self.n_clusters:
+                is_valid = False
+                break
+            r = n_vectors_c / self.n_clusters
+            if fit_pca:
+                self._src_pca[i].fit(Xc)
 
-        The rotations are done from source to target, so in this step the data
-        points suffer no transformation at all.
+            Xc_pca = self._src_pca[i].transform(Xc)
+            for j in range(self.n_components):
+                inds = np.argsort(Xc_pca[:, j], axis=0)
+                anchor = [
+                    np.mean(Xc[inds[round(k*r):round((k+1)*r)]], axis=0)
+                    for k in range(self.n_clusters)
+                ]
+                anchors = np.vstack([anchors, anchor])
+
+        return anchors, is_valid
+
+    def transform(self, X):
+        """Rotate in the target domain, ie do nothing.
+
+        .. note::
+           This method is designed for using at test time on target data.
+           No transformation is applied since rotations are done from source
+           to target domain.
 
         Parameters
         ----------
-        X : ndarray, shape (n_matrices, n_channels, n_channels)
-            Set of SPD matrices.
-        y_enc : None
-            Not used, here for compatibility with sklearn API.
+        X : ndarray, shape (n_matrices, n_channels, n_channels) or \
+                shape (n_vectors, n_ts)
+            Set of SPD matrices or tangent vectors.
 
         Returns
         -------
-        X : ndarray, shape (n_matrices, n_classes)
-            Same set of SPD matrices as in the input.
+        X : ndarray, shape (n_matrices, n_channels, n_channels) or \
+                shape (n_vectors, n_ts)
+            Same data as in the input.
         """
-
-        # used during inference on target domain
         return X
 
     def fit_transform(self, X, y_enc, sample_weight=None):
-        """Fit TLRotate and then transform data points.
+        """Fit TLRotate and then rotate each source domain to target domain.
 
-        Calculate the rotation matrix for matching each source domain to the
-        target domain.
+        It calculates and applies the rotation matrix for matching each source
+        domain to the target domain.
 
         .. note::
-           This method is designed for using at training time. The output for
-           .fit_transform() will be different than using .fit() and
-           .transform() separately.
+           This method is designed for using at training time.
+           The output for .fit_transform() will be different
+           than using .fit() and .transform() separately.
 
         Parameters
         ----------
-        X : ndarray, shape (n_matrices, n_channels, n_channels)
-            Set of SPD matrices.
-        y_enc : ndarray, shape (n_matrices,)
-            Extended labels for each matrix.
-        sample_weight : None | ndarray, shape (n_matrices,), default=None
-            Weights for each matrix. If None, it uses equal weights.
+        X : ndarray, shape (n_matrices, n_channels, n_channels) or \
+                shape (n_vectors, n_ts)
+            Set of SPD matrices or tangent vectors.
+        y_enc : ndarray, shape (n_matrices,) or shape (n_vectors,)
+            Extended labels for each matrix or vector.
+        sample_weight : None | ndarray, shape (n_matrices,) or \
+                shape (n_vectors,), default=None
+            Weights for each matrix or vector. If None, it uses equal weights.
 
         Returns
         -------
-        X : ndarray, shape (n_matrices, n_classes)
-            Set of SPD matrices after rotation step.
+        X_new : ndarray, shape (n_matrices, n_classes)
+            Set of rotated SPD matrices or tangent vectors to target domain.
         """
-
-        # used during fit in pipeline, rotate each source domain
-        self.fit(X, y_enc, sample_weight)
+        self.fit(X, y_enc, sample_weight=sample_weight)
         _, _, domains = decode_domains(X, y_enc)
 
-        X_rot = np.zeros_like(X)
+        X_new = np.zeros_like(X)
         for d in np.unique(domains):
             idx = domains == d
-            if d != self.target_domain:
-                X_rot[idx] = self.rotations_[d] @ X[idx] @ self.rotations_[d].T
+
+            if d == self.target_domain:
+                X_new[idx] = X[idx]
             else:
-                X_rot[idx] = X[idx]
-        return X_rot
+                if X.ndim == 3:
+                    X_new[idx] = self.rotations_[d] @ X[idx] \
+                        @ self.rotations_[d].T
+                else:
+                    X_new[idx] = X[idx] @ self.rotations_
+
+        return X_new
+
+
+###############################################################################
 
 
 class TLEstimator(BaseEstimator):
     """Transfer learning wrapper for estimators.
 
-    This is a wrapper for any BaseEstimator (i.e. classifier or regressor) that
-    converts extended labels used in Transfer Learning into the usual y array
+    This is a wrapper for any BaseEstimator (classifier or regressor) that
+    converts extended labels used in transfer learning into the usual y array
     to train a classifier/regressor of choice.
 
     Parameters
@@ -622,13 +864,18 @@ class TLEstimator(BaseEstimator):
     target_domain : str
         Domain to consider as target.
     estimator : BaseEstimator
-        The estimator to apply on matrices. It can be any regressor or
-        classifier from pyRiemann.
+        The estimator to apply on data. It can be any regressor or classifier
+        from pyRiemann.
     domain_weight : None | dict, default=None
-        Weights to combine matrices from each domain to train the estimator.
+        Weights to combine data from each domain to train the estimator.
         The dict contains key=domain_name and value=weight_to_assign.
         If None, it uses equal weights.
 
+    See Also
+    --------
+    TLClassifier
+    TLRegressor
+
     Notes
     -----
     .. versionadded:: 0.4
@@ -645,10 +892,11 @@ class TLEstimator(BaseEstimator):
 
         Parameters
         ----------
-        X : ndarray, shape (n_matrices, n_channels, n_channels)
-            Set of SPD matrices.
-        y_enc : ndarray, shape (n_matrices,)
-            Extended labels for each matrix.
+        X : ndarray, shape (n_matrices, n_channels, n_channels) or \
+                shape (n_vectors, n_ts)
+            Set of SPD matrices or tangent vectors.
+        y_enc : ndarray, shape (n_matrices,) or shape (n_vectors,)
+            Extended labels for each matrix or vector.
 
         Returns
         -------
@@ -688,13 +936,14 @@ class TLEstimator(BaseEstimator):
 
         Parameters
         ----------
-        X : ndarray, shape (n_matrices, n_channels, n_channels)
-            Set of SPD matrices.
+        X : ndarray, shape (n_matrices, n_channels, n_channels) or \
+                shape (n_vectors, n_ts)
+            Set of SPD matrices or tangent vectors.
 
         Returns
         -------
-        pred : ndarray, shape (n_matrices,)
-            Predictions for each matrix according to the estimator.
+        pred : ndarray, shape (n_matrices,) or shape (n_vectors,)
+            Predictions according to the estimator.
         """
         return self.estimator.predict(X)
 
@@ -703,7 +952,7 @@ class TLClassifier(TLEstimator):
     """Transfer learning wrapper for classifiers.
 
     This is a wrapper for any classifier that converts extended labels used in
-    Transfer Learning into the usual y array to train a classifier of choice.
+    transfer learning into the usual y array to train a classifier of choice.
 
     Parameters
     ----------
@@ -712,10 +961,14 @@ class TLClassifier(TLEstimator):
     estimator : BaseClassifier
         The classifier to apply on matrices.
     domain_weight : None | dict, default=None
-        Weights to combine matrices from each domain to train the classifier.
+        Weights to combine data from each domain to train the classifier.
         The dict contains key=domain_name and value=weight_to_assign.
         If None, it uses equal weights.
 
+    See Also
+    --------
+    TLRegressor
+
     Notes
     -----
     .. versionadded:: 0.4
@@ -726,10 +979,11 @@ class TLClassifier(TLEstimator):
 
         Parameters
         ----------
-        X : ndarray, shape (n_matrices, n_channels, n_channels)
-            Set of SPD matrices.
-        y_enc : ndarray, shape (n_matrices,)
-            Extended labels for each matrix.
+        X : ndarray, shape (n_matrices, n_channels, n_channels) or \
+                shape (n_vectors, n_ts)
+            Set of SPD matrices or tangent vectors.
+        y_enc : ndarray, shape (n_matrices,) or shape (n_vectors,)
+            Extended labels for each matrix or vector.
 
         Returns
         -------
@@ -746,13 +1000,15 @@ class TLClassifier(TLEstimator):
 
         Parameters
         ----------
-        X : ndarray, shape (n_matrices, n_channels, n_channels)
-            Set of SPD matrices.
+        X : ndarray, shape (n_matrices, n_channels, n_channels) or \
+                shape (n_vectors, n_ts)
+            Set of SPD matrices or tangent vectors.
 
         Returns
         -------
-        pred : ndarray, shape (n_matrices, n_classes)
-            Predictions for each matrix.
+        pred : ndarray, shape (n_matrices, n_classes) or \
+                shape (n_vectors, n_classes)
+            Predictions for each matrix or vector.
         """
         return self.estimator.predict_proba(X)
 
@@ -761,10 +1017,11 @@ class TLClassifier(TLEstimator):
 
         Parameters
         ----------
-        X : ndarray, shape (n_matrices, n_channels, n_channels)
-            Test set of SPD matrices.
-        y_enc : ndarray, shape (n_matrices,)
-            Extended true labels for each matrix.
+        X : ndarray, shape (n_matrices, n_channels, n_channels) or \
+                shape (n_vectors, n_ts)
+            Set of SPD matrices or tangent vectors.
+        y_enc : ndarray, shape (n_matrices,) or shape (n_vectors,)
+            Extended labels for each matrix or vector.
 
         Returns
         -------
@@ -780,7 +1037,7 @@ class TLRegressor(TLEstimator):
     """Transfer learning wrapper for regressors.
 
     This is a wrapper for any regressor that converts extended labels used in
-    Transfer Learning into the usual y array to train a regressor of choice.
+    transfer learning into the usual y array to train a regressor of choice.
 
     Parameters
     ----------
@@ -789,10 +1046,14 @@ class TLRegressor(TLEstimator):
     estimator : BaseRegressor
         The regressor to apply on matrices.
     domain_weight : None | dict, default=None
-        Weights to combine matrices from each domain to train the regressor.
+        Weights to combine data from each domain to train the regressor.
         The dict contains key=domain_name and value=weight_to_assign.
         If None, it uses equal weights.
 
+    See Also
+    --------
+    TLClassifier
+
     Notes
     -----
     .. versionadded:: 0.4
@@ -803,10 +1064,11 @@ class TLRegressor(TLEstimator):
 
         Parameters
         ----------
-        X : ndarray, shape (n_matrices, n_channels, n_channels)
-            Set of SPD matrices.
-        y_enc : ndarray, shape (n_matrices,)
-            Extended labels for each matrix.
+        X : ndarray, shape (n_matrices, n_channels, n_channels) or \
+                shape (n_vectors, n_ts)
+            Set of SPD matrices or tangent vectors.
+        y_enc : ndarray, shape (n_matrices,) or shape (n_vectors,)
+            Extended labels for each matrix or vector.
 
         Returns
         -------
@@ -823,10 +1085,11 @@ class TLRegressor(TLEstimator):
 
         Parameters
         ----------
-        X : ndarray, shape (n_matrices, n_channels, n_channels)
-            Test set of SPD matrices.
-        y_enc : ndarray, shape (n_matrices,)
-            Extended true values for each matrix.
+        X : ndarray, shape (n_matrices, n_channels, n_channels) or \
+                shape (n_vectors, n_ts)
+            Set of SPD matrices or tangent vectors.
+        y_enc : ndarray, shape (n_matrices,) or shape (n_vectors,)
+            Extended labels for each matrix or vector.
 
         Returns
         -------
@@ -838,6 +1101,9 @@ class TLRegressor(TLEstimator):
         return r2_score(y_true.astype(float), y_pred)
 
 
+###############################################################################
+
+
 class MDWM(MDM):
     """Classification by Minimum Distance to Weighted Mean.
 
@@ -940,8 +1206,9 @@ class MDWM(MDM):
 
         if not 0 <= self.domain_tradeoff <= 1:
             raise ValueError(
-                "Value domain_tradeoff must be included in [0, 1] (Got %d)"
-                % self.domain_tradeoff)
+                "Parameter domain_tradeoff must be included in [0, 1] "
+                f"(Got {self.domain_tradeoff})"
+            )
 
         X_dec, y_dec, domains = decode_domains(X, y_enc)
         X_src = X_dec[domains != self.target_domain]
@@ -963,19 +1230,19 @@ class MDWM(MDM):
         self.source_means_ = np.stack(
             Parallel(n_jobs=self.n_jobs)(
                 delayed(mean_covariance)(
-                    X_src[y_src == ll],
+                    X_src[y_src == c],
                     metric=self.metric_mean,
-                    sample_weight=sample_weight[y_src == ll],
-                ) for ll in self.classes_
+                    sample_weight=sample_weight[y_src == c],
+                ) for c in self.classes_
             )
         )
 
         self.target_means_ = np.stack(
             Parallel(n_jobs=self.n_jobs)(
                 delayed(mean_covariance)(
-                    X_tgt[y_tgt == ll],
+                    X_tgt[y_tgt == c],
                     metric=self.metric_mean,
-                ) for ll in self.classes_
+                ) for c in self.classes_
             )
         )
 
@@ -985,6 +1252,7 @@ class MDWM(MDM):
             self.domain_tradeoff,
             metric=self.metric_mean,
         )
+
         return self
 
     def score(self, X, y_enc, sample_weight=None):
@@ -993,9 +1261,11 @@ class MDWM(MDM):
         Parameters
         ----------
         X : ndarray, shape (n_matrices, n_channels, n_channels)
-            Set of SPD matrices.
+            Test set of SPD matrices.
         y_enc : ndarray, shape (n_matrices,)
-            Extended labels for each matrix.
+            Extended true labels for each matrix.
+        sample_weight : None | ndarray, shape (n_matrices,), default=None
+            Weights for each matrix. If None, it uses equal weights.
 
         Returns
         -------
diff --git a/pyriemann/transfer/_rotate.py b/pyriemann/transfer/_rotate.py
index 53400b2..dc0aab6 100644
--- a/pyriemann/transfer/_rotate.py
+++ b/pyriemann/transfer/_rotate.py
@@ -1,11 +1,14 @@
 import warnings
+
 import numpy as np
+
 from ..utils.distance import distance
 from ..utils.utils import check_weights
 
 
 def _project(X, U):
-    """
+    """Projection.
+
     Project U on the tangent space of orthogonal matrices with base point X.
     """
 
@@ -13,7 +16,8 @@ def _project(X, U):
 
 
 def _retract(X, v):
-    """
+    """Retraction.
+
     Retraction taking tangent vector v at base point X back to manifold of
     orthogonal matrices.
     """
@@ -26,13 +30,13 @@ def _retract(X, v):
 
 
 def _loss(Q, X, Y, weights, metric="euclid"):
-    """Loss function for estimating the rotation matrix in RPA."""
+    """Loss function for estimating the rotation matrix."""
 
     return weights @ distance(X, Q @ Y @ Q.T, metric=metric)
 
 
 def _grad(Q, X, Y, weights, metric="euclid"):
-    """Gradient of loss function between class means."""
+    """Gradient of loss function for estimating the rotation matrix."""
 
     if metric == "euclid":
         return np.einsum("a,abc->bc", weights, -4 * (X - Q @ Y @ Q.T) @ Q @ Y)
@@ -92,8 +96,8 @@ def _warm_start(X, Y, weights, metric='euclid'):
 
 def _run_minimization(
     Q_ini,
-    M_source,
-    M_target,
+    X_source,
+    X_target,
     weights=None,
     metric="euclid",
     tol_step=1e-9,
@@ -107,11 +111,12 @@ def _run_minimization(
     for _ in range(maxiter):
 
         # get the current value for the loss function
-        F_1 = _loss(Q_1, M_target, M_source, weights, metric=metric)
+        F_1 = _loss(Q_1, X_target, X_source, weights, metric=metric)
 
         # get the direction of steepest descent
         direction = _project(
-            Q_1, _grad(Q_1, M_target, M_source, weights, metric=metric)
+            Q_1,
+            _grad(Q_1, X_target, X_source, weights, metric=metric),
         )
 
         # backtracking line search
@@ -119,13 +124,13 @@ def _run_minimization(
         tau = 0.50
         r = 1e-4
         Q = _retract(Q_1, -alpha * direction)
-        F = _loss(Q, M_target, M_source, weights, metric=metric)
+        F = _loss(Q, X_target, X_source, weights, metric=metric)
         for _ in range(maxiter_linesearch):
             if F_1 - F > r * alpha * np.linalg.norm(direction)**2:
                 break
             alpha = tau * alpha
             Q = _retract(Q_1, -alpha * direction)
-            F = _loss(Q, M_target, M_source, weights, metric=metric)
+            F = _loss(Q, X_target, X_source, weights, metric=metric)
 
         # test if the step size is small
         crit = np.linalg.norm(-alpha * direction)
@@ -142,41 +147,40 @@ def _run_minimization(
     return Q, F
 
 
-def _get_rotation_matrix(
-    M_source,
-    M_target,
+def _get_rotation_manifold(
+    X_source,
+    X_target,
     weights=None,
     metric="euclid",
     tol_step=1e-9,
     maxiter=10_000,
     maxiter_linesearch=32,
 ):
-    r"""Calculate rotation matrix for the Riemannian Procustes Analysis.
+    r"""Calculate rotation matrix in the manifold.
 
-    Get the rotation matrix Q that minimizes the loss function:
+    Riemannian Procustes analysis in the matrix manifold calculates the
+    rotation matrix :math:`\mathbf{Q}` that minimizes the loss function:
 
     .. math::
-        L(Q) = \sum_i w_i delta^2(M_{target_i}, Q M_{source_i} Q^T)
+        L(Q) = \sum_i w_i \delta^2(X^\text{target}_i, Q X^\text{source}_i Q^T)
 
-    The solution can then be used to transform the data points from the source
-    domain so that their class means are close to the those from the target
-    domain. This manifold optimization problem was first defined in Eq(35) of
-    [1]_. The optimization procedure follows the usual setup for optimization
-    on manifolds as described in [2]_.
+    This rotation transforms SPD matrices :math:`\mathbf{M^\text{source}}` from
+    the source domain close to matrices :math:`\mathbf{M^\text{target}}` from
+    the target domain.
+    This manifold optimization problem was first defined in Eq.(35) of [1]_.
+    The optimization procedure follows the usual setup for optimization on
+    manifolds as described in [2]_.
 
     Parameters
     ----------
-    M_source : ndarray, shape (n_classes, n, n)
-        Set with the means of the n_classes from the source domain.
-    M_target : ndarray, shape (n_classes, n, n)
-        Set with the means of the n_classes from the target domain.
-    weights : None | array, shape (n_classes,), default=None
-        Weights for each class. If None, then give the same weight for each
-        class.
+    X_source : ndarray, shape (n_matrices, n, n)
+        Set of SPD matrices from the source domain.
+    X_target : ndarray, shape (n_matrices, n, n)
+        Set of SPD matrices from the target domain.
+    weights : None | array, shape (n_matrices,), default=None
+        Weights for each pair of matrices. If None, it uses equal weights.
     metric : {"euclid", "riemann"}, default="euclid"
-        Which type of distance to minimize between the class means, either
-        "euclid" for Euclidean distance or "riemann" for the affine-invariant
-        Riemannian distance between SPD matrices.
+        Distance to minimize between SPD matrices.
     tol_step : float, default 1e-9
         Stopping criterion based on the norm of the descent direction.
     maxiter : int, default=10_000
@@ -187,9 +191,13 @@ def _get_rotation_matrix(
     Returns
     -------
     Q : ndarray, shape (n, n)
-        Orthogonal matrix that minimizes the distance between the class means
+        Orthogonal matrix that minimizes the distance between the SPD matrices
         for the source and target domains.
 
+    Notes
+    -----
+    .. versionadded:: 0.4
+
     References
     ----------
     .. [1] `Riemannian Procrustes analysis: transfer
@@ -200,37 +208,104 @@ def _get_rotation_matrix(
     .. [2] `An introduction to optimization on smooth manifolds
         <https://www.nicolasboumal.net/book/>`_
         N. Boumal. To appear with Cambridge University Press. June, 2022
-
-    Notes
-    -----
-    .. versionadded:: 0.4
     """
 
-    if M_source.shape[0] != M_target.shape[0]:
-        raise ValueError("The number of classes in each domain don't match")
-    if M_source.shape[1:] != M_target.shape[1:]:
+    if X_source.shape[0] != X_target.shape[0]:
+        raise ValueError("The number of matrices in each domain don't match")
+    if X_source.shape[1:] != X_target.shape[1:]:
         raise ValueError("The number of channels in each domain don't match")
 
-    weights = check_weights(weights, len(M_source))
+    weights = check_weights(weights, len(X_source))
 
     # initialize the solution with an educated guess
-    Q_ini = _warm_start(M_target, M_source, weights, metric=metric)
+    Q_ini = _warm_start(X_target, X_source, weights, metric=metric)
     if np.linalg.det(Q_ini) < 0:  # make sure it is a rotation
         Q_ini[[0, 1], :] = Q_ini[[1, 0], :]
 
     # run the optimization procedure on rotation matrices
     Q_posdet, F_posdet = _run_minimization(
-        Q_ini, M_source, M_target, weights, metric, tol_step, maxiter,
-        maxiter_linesearch)
+        Q_ini,
+        X_source,
+        X_target,
+        weights,
+        metric,
+        tol_step,
+        maxiter,
+        maxiter_linesearch,
+    )
 
     # run the optimization procedure on reflection matrices
     Q_ini[[0, 1], :] = Q_ini[[1, 0], :]
     Q_negdet, F_negdet = _run_minimization(
-        Q_ini, M_source, M_target, weights, metric, tol_step, maxiter,
-        maxiter_linesearch)
+        Q_ini,
+        X_source,
+        X_target,
+        weights,
+        metric,
+        tol_step,
+        maxiter,
+        maxiter_linesearch,
+    )
 
     # check which of the options yield a smaller loss value
     if F_posdet < F_negdet:
         return Q_posdet
     else:
         return Q_negdet
+
+
+###############################################################################
+
+
+def _get_rotation_tangentspace(X_source, X_target, expl_var):
+    r"""Calculate rotation matrix in the tangent space.
+
+    Procustes analysis in the Euclidean tangent space calculates the rotation
+    matrix :math:`\mathbf{Q}` that minimizes:
+
+    .. math::
+        || X^\text{source} Q -  X^\text{target} ||_F^2
+
+    From code provided by the authors [1]_.
+
+    Parameters
+    ----------
+    X_source : ndarray, shape (n_vectors, n_ts)
+        Set of tangent vectors from the source domain.
+    X_target : ndarray, shape (n_vectors, n_ts)
+        Set of tangent vectors from the target domain.
+    expl_var : float
+        Dimension reduction applied to the cross product matrix during
+        Procrustes analysis.
+        If float in (0,1], percentage of variance that needs to be explained.
+        Else, number of components.
+
+    Returns
+    -------
+    Q : ndarray, shape (n_ts, n_ts)
+        Orthogonal matrix that minimizes the distance between the tangent
+        vectors for the source and target domains.
+
+    Notes
+    -----
+    .. versionadded:: 0.8
+
+    References
+    ----------
+    .. [1] `Tangent space alignment: Transfer learning for brain-computer
+        interface
+        <https://www.frontiersin.org/articles/10.3389/fnhum.2022.1049985/pdf>`_
+        A. Bleuzé, J. Mattout and M. Congedo, Frontiers in Human Neuroscience,
+        2022
+    """
+    C = X_source.T @ X_target
+    u, s, vh = np.linalg.svd(C)
+
+    if expl_var <= 1:
+        n_comps = np.sum(np.cumsum(s) < expl_var * np.sum(s)) + 1
+    else:
+        n_comps = int(expl_var)
+
+    u = u[:, :n_comps]
+    vh = vh[:n_comps, :]
+    return vh.T @ u.T
diff --git a/pyriemann/transfer/_tools.py b/pyriemann/transfer/_tools.py
index 8217d58..e664cd4 100644
--- a/pyriemann/transfer/_tools.py
+++ b/pyriemann/transfer/_tools.py
@@ -2,30 +2,32 @@ import numpy as np
 
 
 def encode_domains(X, y, domain):
-    r"""Encode the domains of the matrices in the labels.
+    r"""Encode the domains of the data in the labels.
 
     We handle the possibility of having different domains for the datasets by
-    extending the labels of the matrices and including this information to
-    them. For instance, if we have a matrix X with class `left_hand` on the
-    `domain_01` then its extended label will be `domain_01/left_hand`. Note
-    that if the classes were integers at first, they will be converted to
+    extending the labels of the data and including this information to them.
+    For instance, if we have a datum X with class `left_hand` on the
+    `domain_01` then its extended label will be `domain_01/left_hand`.
+    Note that if the classes were integers at first, they will be converted to
     strings.
 
     Parameters
     ----------
-    X : ndarray, shape (n_matrices, n_channels, n_channels)
-        Set of SPD matrices.
-    y : ndarray, shape (n_matrices,)
-        Labels for each matrix.
-    domain : ndarray, shape (n_matrices,)
-        Domains for each matrix.
+    X : ndarray, shape (n_matrices, n_channels, n_channels) or \
+            shape (n_vectors, n_ts)
+        Set of SPD matrices or tangent vectors.
+    y : ndarray, shape (n_matrices,) or shape (n_vectors,)
+        Labels for each matrix or vector.
+    domain : ndarray, shape (n_matrices,) or shape (n_vectors,)
+        Domains for each matrix or vector.
 
     Returns
     -------
-    X_enc : ndarray, shape (n_matrices, n_channels, n_channels)
-        The same set of SPD matrices given as input.
-    y_enc : ndarray, shape (n_matrices,)
-        Extended labels for each matrix.
+    X_enc : ndarray, shape (n_matrices, n_channels, n_channels) or \
+            shape (n_vectors, n_ts)
+        The same data given as input.
+    y_enc : ndarray, shape (n_matrices,) or shape (n_vectors,)
+        Extended labels for each matrix or vector.
 
     See Also
     --------
@@ -38,33 +40,35 @@ def encode_domains(X, y, domain):
     if len(y) != len(domain):
         raise ValueError("Input lengths don't match")
 
-    y_enc = [str(d_) + '/' + str(y_) for (d_, y_) in zip(domain, y)]
+    y_enc = [str(d_) + "/" + str(y_) for (d_, y_) in zip(domain, y)]
     return X, np.array(y_enc)
 
 
 def decode_domains(X_enc, y_enc):
-    """Decode the domains of the matrices in the labels.
+    """Decode the domains of the data in the labels.
 
     We handle the possibility of having different domains for the datasets by
-    encoding the domain information into the labels of the matrices. This
-    method converts the data into its original form, with a separate data
+    encoding the domain information into the labels of the data.
+    This method converts the data into its original form, with a separate data
     structure for labels and for domains.
 
     Parameters
     ----------
-    X_enc : ndarray, shape (n_matrices, n_channels, n_channels)
-        Set of SPD matrices.
-    y_enc : ndarray, shape (n_matrices,)
-        Extended labels for each matrix.
+    X_enc : ndarray, shape (n_matrices, n_channels, n_channels) or \
+            shape (n_vectors, n_ts)
+        Set of SPD matrices or tangent vectors.
+    y_enc : ndarray, shape (n_matrices,) or shape (n_vectors,)
+        Extended labels for each matrix or vector.
 
     Returns
     -------
-    X : ndarray, shape (n_matrices, n_channels, n_channels)
-        Set of SPD matrices.
-    y : ndarray, shape (n_matrices,)
-        Labels for each matrix.
-    domain : ndarray, shape (n_matrices,)
-        Domains for each matrix.
+    X : ndarray, shape (n_matrices, n_channels, n_channels) or \
+            shape (n_vectors, n_ts)
+        The same data given as input.
+    y : ndarray, shape (n_matrices,) or shape (n_vectors,)
+        Labels for each matrix or vector.
+    domain : ndarray, shape (n_matrices,) or shape (n_vectors,)
+        Domains for each matrix or vector.
 
     See Also
     --------
@@ -76,7 +80,7 @@ def decode_domains(X_enc, y_enc):
     """
     y, domain = [], []
     for y_enc_ in y_enc:
-        y_dec_ = y_enc_.split('/')
+        y_dec_ = y_enc_.split("/")
         domain.append(y_dec_[-2])
         y.append(y_dec_[-1])
     return X_enc, np.array(y), np.array(domain)
@@ -86,9 +90,9 @@ class TLSplitter():
     """Class for handling the cross-validation splits of multi-domain data.
 
     This is a wrapper to sklearn's cross-validation iterators [1]_ which
-    ensures the handling of domain information with the data points. In fact,
+    ensures the handling of domain information with the data. In fact,
     the data from source domain is always fully available in the training
-    partition whereas the random splits are done on the data points from the
+    partition whereas the random splits are done on the data from the
     target domain.
 
     Parameters
@@ -96,18 +100,17 @@ class TLSplitter():
     target_domain : str
         Domain considered as target.
     cv : None | BaseCrossValidator | BaseShuffleSplit, default=None
-        An instance of a cross validation iterator from sklearn.
-
-    References
-    ----------
-    .. [1] https://scikit-learn.org/stable/modules/cross_validation.html#cross-validation-iterators
+        An instance of a cross-validation iterator from sklearn.
 
     Notes
     -----
     .. versionadded:: 0.4
+
+    References
+    ----------
+    .. [1] https://scikit-learn.org/stable/modules/cross_validation.html#cross-validation-iterators
     """  # noqa
     def __init__(self, target_domain, cv):
-
         self.target_domain = target_domain
         self.cv = cv
 
@@ -116,10 +119,11 @@ class TLSplitter():
 
         Parameters
         ----------
-        X : ndarray, shape (n_matrices, n_channels, n_channels)
-            Set of SPD matrices.
-        y : ndarray, shape (n_matrices,)
-            Extended labels for each matrix.
+        X : ndarray, shape (n_matrices, n_channels, n_channels) or \
+                shape (n_vectors, n_ts)
+            Set of SPD matrices or tangent vectors.
+        y : ndarray, shape (n_matrices,) or shape (n_vectors,)
+            Extended labels for each matrix or vector.
 
         Yields
         ------
@@ -129,24 +133,25 @@ class TLSplitter():
             The testing set indices for that split.
         """
 
-        # decode the domains of the data points
+        # decode the domains of the data
         X, y, domain = decode_domains(X, y)
 
-        # indentify the indices of the target dataset
+        # identify the indices of the target dataset
         idx_source = np.where(domain != self.target_domain)[0]
         idx_target = np.where(domain == self.target_domain)[0]
         y_target = y[idx_target]
 
-        # index of training-split for the target data points
+        # index of training-split for the target data
         ss_target = self.cv.split(idx_target, y_target)
         for train_sub_idx_target, test_sub_idx_target in ss_target:
             train_idx = np.concatenate(
-                [idx_source, idx_target[train_sub_idx_target]])
+                [idx_source, idx_target[train_sub_idx_target]]
+            )
             test_idx = idx_target[test_sub_idx_target]
             yield train_idx, test_idx
 
     def get_n_splits(self, X=None, y=None):
-        """Returns the number of splitting iterations in the cross-validator.
+        """Return the number of splitting iterations in the cross-validator.
 
         Parameters
         ----------
@@ -158,6 +163,6 @@ class TLSplitter():
         Returns
         -------
         n_splits : int
-            Returns the number of splitting iterations in the cross-validator.
+            Number of splitting iterations in the cross-validator.
         """
         return self.cv.n_splits
diff --git a/pyriemann/utils/base.py b/pyriemann/utils/base.py
index 66952bc..9a4419a 100644
--- a/pyriemann/utils/base.py
+++ b/pyriemann/utils/base.py
@@ -1,7 +1,7 @@
 """Base functions for SPD/HPD matrices."""
 
 import numpy as np
-from numpy.core.numerictypes import typecodes
+
 from .test import is_pos_def
 
 
@@ -9,7 +9,7 @@ def _matrix_operator(C, operator):
     """Matrix function."""
     if not isinstance(C, np.ndarray) or C.ndim < 2:
         raise ValueError("Input must be at least a 2D ndarray")
-    if C.dtype.char in typecodes['AllFloat'] and (
+    if C.dtype.char in np.typecodes['AllFloat'] and (
             np.isinf(C).any() or np.isnan(C).any()):
         raise ValueError(
             "Matrices must be positive definite. Add "
diff --git a/requirements.txt b/requirements.txt
index 95cd64d..274c1e6 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,4 +1,4 @@
-numpy!=1.24.0
+numpy>=2.0.0
 scipy
 scikit-learn>=0.24
 joblib
diff --git a/setup.py b/setup.py
index 04f2dfa..44f5eaf 100644
--- a/setup.py
+++ b/setup.py
@@ -38,7 +38,7 @@ setup(
     platforms="any",
     python_requires=">=3.9",
     install_requires=[
-        "numpy!=1.24.0",
+        "numpy>=2.0.0",
         "scipy",
         "scikit-learn>=0.24",
         "joblib",


Test patch:
diff --git a/tests/test_classification.py b/tests/test_classification.py
index 35245d2..7bc836b 100644
--- a/tests/test_classification.py
+++ b/tests/test_classification.py
@@ -17,13 +17,13 @@ from pyriemann.classification import (
     MDM,
     FgMDM,
     KNearestNeighbor,
-    TSclassifier,
+    TSClassifier,
     SVC,
     MeanField,
     class_distinctiveness,
 )
 
-rclf = [MDM, FgMDM, KNearestNeighbor, TSclassifier, SVC, MeanField]
+rclf = [MDM, FgMDM, KNearestNeighbor, TSClassifier, SVC, MeanField]
 
 
 @pytest.mark.parametrize(
@@ -181,7 +181,7 @@ def test_metric_str(classif, metric, get_mats, get_labels):
     mats = get_mats(n_matrices, n_channels, "spd")
     clf = classif(metric=metric)
 
-    if classif in [SVC, FgMDM, TSclassifier] \
+    if classif in [SVC, FgMDM, TSClassifier] \
             and metric not in ["euclid", "logchol", "logeuclid", "riemann"]:
         with pytest.raises((KeyError, ValueError)):
             clf.fit(mats, labels).predict(mats)
@@ -266,7 +266,7 @@ def test_tsclassifier(metric_mean, metric_map, get_mats, get_labels):
     labels = get_labels(n_matrices, n_classes)
     mats = get_mats(n_matrices, n_channels, "spd")
 
-    clf = TSclassifier(metric={"mean": metric_mean, "map": metric_map})
+    clf = TSClassifier(metric={"mean": metric_mean, "map": metric_map})
     clf.fit(mats, labels).predict(mats)
 
 
@@ -276,7 +276,7 @@ def test_tsclassifier_fit(get_mats, get_labels):
     mats = get_mats(n_matrices, n_channels, "spd")
     labels = get_labels(n_matrices, n_classes)
 
-    clf = TSclassifier(clf=DummyClassifier())
+    clf = TSClassifier(clf=DummyClassifier())
     clf.fit(mats, labels).predict(mats)
 
 
@@ -286,7 +286,7 @@ def test_tsclassifier_clf_error(get_mats, get_labels):
     mats = get_mats(n_matrices, n_channels, "spd")
     labels = get_labels(n_matrices, n_classes)
     with pytest.raises(TypeError):
-        TSclassifier(clf=Covariances()).fit(mats, labels)
+        TSClassifier(clf=Covariances()).fit(mats, labels)
 
 
 def test_svc_params():
diff --git a/tests/test_transfer.py b/tests/test_transfer.py
index 77c1e54..9c0f261 100644
--- a/tests/test_transfer.py
+++ b/tests/test_transfer.py
@@ -1,9 +1,11 @@
 import numpy as np
-from numpy.testing import assert_array_equal
+from numpy.testing import assert_array_equal, assert_array_almost_equal
 import pytest
 from pytest import approx
+from sklearn.linear_model import LogisticRegression, LinearRegression
 from sklearn.model_selection import KFold, StratifiedShuffleSplit
 from sklearn.pipeline import make_pipeline, Pipeline
+from sklearn.svm import LinearSVC, LinearSVR
 
 from pyriemann.datasets.simulated import (
     make_classification_transfer,
@@ -13,7 +15,7 @@ from pyriemann.classification import (
     MDM,
     FgMDM,
     KNearestNeighbor,
-    TSclassifier,
+    TSClassifier,
     SVC,
     MeanField,
 )
@@ -21,24 +23,49 @@ from pyriemann.regression import KNearestNeighborRegressor, SVR
 from pyriemann.transfer import (
     decode_domains,
     encode_domains,
+    TLSplitter,
     TLDummy,
     TLCenter,
-    TLStretch,
+    TLScale,
     TLRotate,
-    TLSplitter,
     TLClassifier,
     TLRegressor,
     MDWM,
 )
 from pyriemann.utils.distance import distance, distance_riemann
 from pyriemann.utils.mean import mean_covariance, mean_riemann
+from pyriemann.utils.tangentspace import tangent_space
 from pyriemann.utils.utils import check_weights
 
 rndstate = 1234
 
 
+###############################################################################
+
+
+def make_classification_transfer_tangspace(rndstate, n_vectors_d=50, n_ts=10):
+    n_vectors = 2 * n_vectors_d
+    X = rndstate.randn(n_vectors, n_ts)
+    y = rndstate.randint(0, 3, size=n_vectors)
+    domains = np.array(n_vectors_d * ["src"] + n_vectors_d * ["tgt"])
+    _, y_enc = encode_domains(X, y, domains)
+    return X, y_enc
+
+
+def make_regression_transfer(rndstate, n_matrices, n_channels):
+    X = make_matrices(n_matrices, n_channels, "spd", rs=rndstate)
+    y = rndstate.randn(n_matrices)
+    domain = ["source_domain"] * (n_matrices // 2) \
+        + ["target_domain"] * (n_matrices // 2)
+    _, y_enc = encode_domains(X, y, domain)
+    return X, y_enc
+
+
+###############################################################################
+
+
 def test_encode_decode_domains(rndstate):
-    """Test encoding and decoding of domains for data points"""
+    """Test encoding and decoding of domains"""
     n_matrices, n_channels = 4, 2
     X = make_matrices(n_matrices, n_channels, "spd", rs=rndstate)
     y = np.array(["left_hand", "right_hand", "left_hand", "right_hand"])
@@ -55,23 +82,50 @@ def test_encode_decode_domains(rndstate):
     assert (domain == d_dec).all()
 
 
-def test_tldummy(rndstate):
-    """Test TLDummy"""
-    dum = TLDummy()
+@pytest.mark.parametrize(
+    "cv",
+    [
+        KFold(n_splits=5, shuffle=True),
+        StratifiedShuffleSplit(n_splits=5, train_size=0.80),
+    ]
+)
+def test_tlsplitter(rndstate, cv):
+    """Test wrapper for cross-validation"""
+    X, y_enc = make_classification_transfer(
+        n_matrices=25, class_sep=5, class_disp=1.0, random_state=rndstate
+    )
+
+    cv = TLSplitter(target_domain="target_domain", cv=cv)
+    train_idx, test_idx = next(cv.split(X, y_enc))
+
+    assert len(train_idx) == 90  # 50 from source and 4/5*100 from target
+    assert len(test_idx) == 10  # 1/5*100 from target
+    assert cv.get_n_splits() == 5
+
+
+###############################################################################
+
+
+@pytest.mark.parametrize("space", ["manifold", "tangentspace"])
+def test_tldummy(rndstate, space):
     X, y_enc = make_classification_transfer(
         n_matrices=5, random_state=rndstate
     )
+    if space == "tangentspace":
+        X = tangent_space(X, Cref=mean_riemann(X), metric="riemann")
+
+    dum = TLDummy()
     dum.fit(X, y_enc)
     dum.fit_transform(X, y_enc)
+    dum.transform(X)
 
 
 @pytest.mark.parametrize("metric", ["riemann", "euclid"])
 @pytest.mark.parametrize("use_weight", [True, False])
 @pytest.mark.parametrize("target_domain", ["target_domain", ""])
-def test_tlcenter(rndstate, get_weights, metric, use_weight, target_domain):
-    """Test pipeline for recentering data to Identity"""
-    # check if the global mean of the domains is indeed Identity
-    rct = TLCenter(target_domain=target_domain, metric=metric)
+def test_tlcenter_manifold(rndstate, get_weights,
+                           metric, use_weight, target_domain):
+    """Test centering matrices to identity"""
     X, y_enc = make_classification_transfer(
         n_matrices=25, random_state=rndstate
     )
@@ -80,12 +134,15 @@ def test_tlcenter(rndstate, get_weights, metric, use_weight, target_domain):
     else:
         weights = None
 
+    _, _, domain = decode_domains(X, y_enc)
+
+    rct = TLCenter(target_domain=target_domain, metric=metric)
+
     # Test fit
     rct.fit(X, y_enc, sample_weight=weights)
 
-    # Test fit_transform
+    # Test fit_transform, if the mean of each domain is indeed identity
     X_rct = rct.fit_transform(X, y_enc, sample_weight=weights)
-    _, _, domain = decode_domains(X_rct, y_enc)
     for d in np.unique(domain):
         idx = domain == d
         Xd = X_rct[idx]
@@ -96,15 +153,24 @@ def test_tlcenter(rndstate, get_weights, metric, use_weight, target_domain):
             Md = mean_covariance(Xd, metric=metric)
         assert Md == pytest.approx(np.eye(2))
 
+    # Test transform
+    rct.transform(X)
+
+    # Test deprecated attribute
+    import warnings
+    with warnings.catch_warnings(record=True) as w:
+        warnings.simplefilter("always")
+        rct.recenter_
+        assert len(w) == 1
+        assert issubclass(w[-1].category, DeprecationWarning)
+
 
 @pytest.mark.parametrize("metric", ["riemann", "euclid"])
 @pytest.mark.parametrize("use_weight", [True, False])
 @pytest.mark.parametrize("target_domain", ["target_domain", ""])
-def test_tlcenter_fit_transf(rndstate, get_weights,
-                             metric, use_weight, target_domain):
+def test_tlcenter_manifold_fit_transf(rndstate, get_weights,
+                                      metric, use_weight, target_domain):
     """Test .fit_transform() versus .fit().transform()"""
-    # check if the global mean of the domains is indeed Identity
-    rct = TLCenter(target_domain=target_domain, metric=metric)
     X, y_enc = make_classification_transfer(
         n_matrices=25, random_state=rndstate
     )
@@ -120,10 +186,12 @@ def test_tlcenter_fit_transf(rndstate, get_weights,
     X1, y1 = X[:50], y_enc[:50]
     X2, y2 = X[50:], y_enc[50:]
 
+    rct = TLCenter(target_domain=target_domain, metric=metric)
     X1_rct = rct.fit_transform(X1, y1, sample_weight=weights_1)
     X2_rct = rct.fit(X2, y2, sample_weight=weights_2).transform(X2)
     X_rct = np.concatenate((X1_rct, X2_rct))
 
+    # Test if the mean of each domain is indeed identity
     for d in np.unique(domain):
         idx = domain == d
         Xd = X_rct[idx]
@@ -136,22 +204,40 @@ def test_tlcenter_fit_transf(rndstate, get_weights,
         assert Md == pytest.approx(np.eye(2))
 
 
+def test_tlcenter_tangentspace(rndstate):
+    """Test centering tangent vectors to origin"""
+    n_vectors_d, n_ts = 50, 10
+    X, y_enc = make_classification_transfer_tangspace(
+        rndstate, n_vectors_d, n_ts
+    )
+
+    _, _, domain = decode_domains(X, y_enc)
+
+    tlctr = TLCenter(target_domain="tgt")
+
+    tlctr.fit(X, y_enc)
+
+    # Test if the mean of each domain is indeed zero
+    X_rct = tlctr.fit_transform(X, y_enc)
+    assert X_rct.shape == X.shape
+    for d in np.unique(domain):
+        idx = domain == d
+        md = np.mean(X_rct[idx], axis=0)
+        assert_array_almost_equal(md, np.zeros((n_ts)))
+
+    X_rct = tlctr.transform(X)
+    assert X_rct.shape == X.shape
+
+
 @pytest.mark.parametrize("use_centered_data", [True, False])
 @pytest.mark.parametrize("metric", ["riemann"])
 @pytest.mark.parametrize("use_weight", [True, False])
-def test_tlstretch(rndstate, get_weights,
-                   use_centered_data, metric, use_weight):
-    """Test pipeline for stretching data"""
-    # check if the dispersion of the dataset indeed decreases to 1
-    tlstr = TLStretch(
-        target_domain="target_domain",
-        final_dispersion=1.0,
-        centered_data=use_centered_data,
-        metric=metric
-    )
+def test_tlscale_manifold(rndstate, get_weights,
+                          use_centered_data, metric, use_weight):
+    """Test scaling matrices to a target dispersion"""
     X, y_enc = make_classification_transfer(
-        n_matrices=25, class_disp=2.0, random_state=rndstate)
-
+        n_matrices=25, class_disp=2.0, random_state=rndstate
+    )
     if use_weight:
         weights = get_weights(len(y_enc))
     else:
@@ -160,9 +246,17 @@ def test_tlstretch(rndstate, get_weights,
         tlrct = TLCenter(target_domain="target_domain", metric=metric)
         X = tlrct.fit_transform(X, y_enc, sample_weight=weights)
 
+    _, _, domain = decode_domains(X, y_enc)
+
+    tlstr = TLScale(
+        target_domain="target_domain",
+        final_dispersion=1.0,
+        centered_data=use_centered_data,
+        metric=metric
+    )
     X_str = tlstr.fit_transform(X, y_enc, sample_weight=weights)
 
-    _, _, domain = decode_domains(X_str, y_enc)
+    # Test if the dispersion of each domain is indeed 1
     for d in np.unique(domain):
         idx = domain == d
         Xd = X_str[idx]
@@ -176,67 +270,131 @@ def test_tlstretch(rndstate, get_weights,
             disp = np.mean(distance(Xd, Md, metric=metric, squared=True))
         assert np.isclose(disp, 1.0)
 
+    tlstr.transform(X)
+
+    # Test deprecated attribute
+    import warnings
+    with warnings.catch_warnings(record=True) as w:
+        warnings.simplefilter("always")
+        tlstr.dispersions_
+        assert len(w) == 1
+        assert issubclass(w[-1].category, DeprecationWarning)
+
+
+def test_tlscale_tangentspace(rndstate):
+    """Test scaling vectors to unit norm"""
+    n_vectors_d, n_ts = 50, 10
+    X, y_enc = make_classification_transfer_tangspace(
+        rndstate, n_vectors_d, n_ts
+    )
+
+    _, _, domain = decode_domains(X, y_enc)
+
+    tlscl = TLScale(target_domain="tgt")
+
+    tlscl.fit(X, y_enc)
+
+    # Test if the mean of norms of each domain is indeed 1
+    X_scl = tlscl.fit_transform(X, y_enc)
+    assert X_scl.shape == X.shape
+    for d in np.unique(domain):
+        idx = domain == d
+        assert np.mean(np.linalg.norm(X_scl[idx], axis=1)) == pytest.approx(1)
+
+    X_scl = tlscl.transform(X)
+    assert X_scl.shape == X.shape
+
 
 @pytest.mark.parametrize("metric", ["euclid", "riemann"])
 @pytest.mark.parametrize("use_weight", [True, False])
-def test_tlrotate(rndstate, get_weights, metric, use_weight):
-    """Test fit_transform method for rotating the datasets"""
-    # check if the distance between the classes of each domain is reduced
+def test_tlrotate_manifold(rndstate, get_weights, metric, use_weight):
+    """Test rotating matrices"""
     X, y_enc = make_classification_transfer(
-        n_matrices=50, class_sep=3, class_disp=1.0, theta=np.pi / 2,
-        random_state=rndstate)
-
+        n_matrices=50,
+        class_sep=3,
+        class_disp=1.0,
+        theta=np.pi / 2,
+        random_state=rndstate,
+    )
     if use_weight:
         weights = get_weights(len(y_enc))
     else:
         weights = None
     weights = check_weights(weights, len(y_enc))
 
+    _, y, domain = decode_domains(X, y_enc)
+
     rct = TLCenter(target_domain="target_domain", metric=metric)
     X_rct = rct.fit_transform(X, y_enc, sample_weight=weights)
     rot = TLRotate(target_domain="target_domain", metric=metric)
     X_rot = rot.fit_transform(X_rct, y_enc, sample_weight=weights)
+    for k in rot.rotations_.keys():
+        assert rot.rotations_[k].shape == (2, 2)
+    assert_array_equal(
+        X_rot[domain == "target_domain"],
+        X_rct[domain == "target_domain"],
+    )
 
-    _, y, domain = decode_domains(X_rot, y_enc)
+    # check if the distance between the classes of each domain is reduced
     for label in np.unique(y):
         d = "source_domain"
         M_rct_label_source = mean_riemann(
             X_rct[domain == d][y[domain == d] == label],
-            sample_weight=weights[domain == d][y[domain == d] == label])
+            sample_weight=weights[domain == d][y[domain == d] == label]
+        )
         M_rot_label_source = mean_riemann(
             X_rot[domain == d][y[domain == d] == label],
-            sample_weight=weights[domain == d][y[domain == d] == label])
+            sample_weight=weights[domain == d][y[domain == d] == label]
+        )
         d = "target_domain"
         M_rct_label_target = mean_riemann(
             X_rct[domain == d][y[domain == d] == label],
-            sample_weight=weights[domain == d][y[domain == d] == label])
+            sample_weight=weights[domain == d][y[domain == d] == label]
+        )
         M_rot_label_target = mean_riemann(
             X_rot[domain == d][y[domain == d] == label],
-            sample_weight=weights[domain == d][y[domain == d] == label])
+            sample_weight=weights[domain == d][y[domain == d] == label]
+        )
         d_rct = distance_riemann(M_rct_label_source, M_rct_label_target)
         d_rot = distance_riemann(M_rot_label_source, M_rot_label_target)
         assert d_rot <= d_rct
 
+    X_rot = rot.transform(X_rct)
+    assert_array_equal(X_rot, X_rct)
 
-@pytest.mark.parametrize(
-    "cv",
-    [
-        KFold(n_splits=5, shuffle=True),
-        StratifiedShuffleSplit(n_splits=5, train_size=0.80),
-    ]
-)
-def test_tlsplitter(rndstate, cv):
-    """Test wrapper for cross-validation in transfer learning"""
-    X, y_enc = make_classification_transfer(
-        n_matrices=25, class_sep=5, class_disp=1.0, random_state=rndstate)
-    cv = TLSplitter(
-        target_domain="target_domain",
-        cv=cv,
+
+@pytest.mark.parametrize("n_components", [1, 3, "max"])
+@pytest.mark.parametrize("n_clusters", [1, 2, 5])
+def test_tlrotate_tangentspace(rndstate, n_components, n_clusters):
+    """Test rotating vectors"""
+    n_vectors_d, n_ts = 50, 10
+    X, y_enc = make_classification_transfer_tangspace(
+        rndstate, n_vectors_d, n_ts
     )
-    train_idx, test_idx = next(cv.split(X, y_enc))
-    assert len(train_idx) == 90  # 50 from source and 4/5*100 from target
-    assert len(test_idx) == 10  # 1/5*100 from target
-    assert cv.get_n_splits() == 5
+
+    _, _, domain = decode_domains(X, y_enc)
+
+    tlrot = TLRotate(
+        target_domain="tgt",
+        n_components=n_components,
+        n_clusters=n_clusters,
+    )
+
+    tlrot.fit(X, y_enc)
+    assert tlrot.rotations_.shape == (n_ts, n_ts)
+
+    X_rot = tlrot.fit_transform(X, y_enc)
+    assert X_rot.shape == X.shape
+    assert_array_equal(
+        X_rot[domain == "target_domain"],
+        X[domain == "target_domain"],
+    )
+
+    X_rot = tlrot.transform(X)
+    assert_array_equal(X_rot, X)
+
+
+###############################################################################
 
 
 @pytest.mark.parametrize(
@@ -247,11 +405,11 @@ def test_tlsplitter(rndstate, cv):
     ]
 )
 @pytest.mark.parametrize(
-    "source_domain, target_domain",
+    "source_weight, target_weight",
     [(1.0, 0.0), (0.0, 1.0), (1.0, 1.0)],
 )
-def test_tlclassifier_mdm(rndstate, clf, source_domain, target_domain):
-    """Test wrapper for MDM classifier in transfer learning"""
+def test_tlclassifier_mdm(rndstate, clf, source_weight, target_weight):
+    """Test wrapper for MDM classifier"""
     n_matrices = 40
     X, y_enc = make_classification_transfer(
         n_matrices=n_matrices // 4,
@@ -260,11 +418,14 @@ def test_tlclassifier_mdm(rndstate, clf, source_domain, target_domain):
         random_state=rndstate,
     )
 
-    tlclf = TLClassifier(target_domain="target_domain", estimator=clf)
-    tlclf.domain_weight = {
-        "source_domain": source_domain,
-        "target_domain": target_domain,
-    }
+    tlclf = TLClassifier(
+        target_domain="target_domain",
+        estimator=clf,
+        domain_weight={
+            "source_domain": source_weight,
+            "target_domain": target_weight,
+        },
+    )
     tlclf.fit(X, y_enc)
 
     X, y, domain = decode_domains(X, y_enc)
@@ -273,13 +434,13 @@ def test_tlclassifier_mdm(rndstate, clf, source_domain, target_domain):
     X_target = X[domain == "target_domain"]
     y_target = y[domain == "target_domain"]
 
-    if source_domain == 1.0 and target_domain == 0.0:
+    if source_weight == 1.0 and target_weight == 0.0:
         X_0 = X_source[y_source == tlclf.estimator.classes_[0]]
         X_1 = X_source[y_source == tlclf.estimator.classes_[1]]
-    elif source_domain == 0.0 and target_domain == 1.0:
+    elif source_weight == 0.0 and target_weight == 1.0:
         X_0 = X_target[y_target == tlclf.estimator.classes_[0]]
         X_1 = X_target[y_target == tlclf.estimator.classes_[1]]
-    elif source_domain == 1.0 and target_domain == 1.0:
+    elif source_weight == 1.0 and target_weight == 1.0:
         X_0 = X[y == tlclf.estimator.classes_[0]]
         X_1 = X[y == tlclf.estimator.classes_[1]]
 
@@ -300,32 +461,61 @@ def test_tlclassifier_mdm(rndstate, clf, source_domain, target_domain):
         make_pipeline(MDM(metric="riemann")),
         FgMDM(),
         make_pipeline(KNearestNeighbor(metric="logeuclid")),
-        TSclassifier(),
+        TSClassifier(),
         SVC(metric="riemann"),
         make_pipeline(MeanField(metric="logeuclid")),
     ]
 )
 @pytest.mark.parametrize(
-    "source_domain, target_domain",
+    "source_weight, target_weight",
     [(1.0, 0.0), (0.0, 1.0), (1.0, 1.0)],
 )
-def test_tlclassifiers(rndstate, clf, source_domain, target_domain):
-    """Test wrapper for classifiers in transfer learning"""
-    n_classes, n_matrices = 2, 40
+def test_tlclassifier_manifold(rndstate, clf, source_weight, target_weight):
+    """Test wrapper for classifiers in manifold"""
+    n_matrices = 40
+    X, y_enc = make_classification_transfer(
+        n_matrices=n_matrices // 4,
+        class_sep=5,
+        class_disp=1.0,
+        random_state=rndstate,
+    )
+
+    tlclassifier(clf, X, y_enc, source_weight, target_weight)
+
+
+@pytest.mark.parametrize("clf", [LinearSVC(), LogisticRegression()])
+@pytest.mark.parametrize(
+    "source_weight, target_weight",
+    [(1.0, 0.0), (0.0, 1.0), (1.0, 1.0)],
+)
+def test_tlclassifier_tangentspace(rndstate,
+                                   clf, source_weight, target_weight):
+    """Test wrapper for classifiers in tangent space"""
+    n_matrices = 40
     X, y_enc = make_classification_transfer(
         n_matrices=n_matrices // 4,
         class_sep=5,
         class_disp=1.0,
         random_state=rndstate,
     )
+    X = tangent_space(X, Cref=mean_riemann(X), metric="riemann")
+
+    tlclassifier(clf, X, y_enc, source_weight, target_weight)
+
+
+def tlclassifier(clf, X, y_enc, source_weight, target_weight, n_classes=2):
+    n_inputs = X.shape[0]
 
     if hasattr(clf, "probability"):
         clf.set_params(**{"probability": True})
-    tlclf = TLClassifier(target_domain="target_domain", estimator=clf)
-    tlclf.domain_weight = {
-        "source_domain": source_domain,
-        "target_domain": target_domain,
-    }
+    tlclf = TLClassifier(
+        target_domain="target_domain",
+        estimator=clf,
+        domain_weight={
+            "source_domain": source_weight,
+            "target_domain": target_weight,
+        },
+    )
 
     # test fit
     tlclf.fit(X, y_enc)
@@ -333,12 +523,13 @@ def test_tlclassifiers(rndstate, clf, source_domain, target_domain):
 
     # test predict
     predicted = tlclf.predict(X)
-    assert predicted.shape == (n_matrices,)
+    assert predicted.shape == (n_inputs,)
 
     # test predict_proba
-    probabilities = tlclf.predict_proba(X)
-    assert probabilities.shape == (n_matrices, n_classes)
-    assert probabilities.sum(axis=1) == approx(np.ones(n_matrices))
+    if hasattr(clf, "predict_proba"):
+        probabilities = tlclf.predict_proba(X)
+        assert probabilities.shape == (n_inputs, n_classes)
+        assert probabilities.sum(axis=1) == approx(np.ones(n_inputs))
 
     # test score
     tlclf.score(X, y_enc)
@@ -353,44 +544,61 @@ def test_tlclassifiers(rndstate, clf, source_domain, target_domain):
     ]
 )
 @pytest.mark.parametrize(
-    "source_domain, target_domain",
+    "source_weight, target_weight",
     [(1.0, 0.0), (0.0, 1.0), (1.0, 1.0)],
 )
-def test_tlregressors(rndstate, get_weights,
-                      reg, source_domain, target_domain):
-    """Test wrapper for regressors in transfer learning"""
+def test_tlregressor_manifold(rndstate, reg, source_weight, target_weight):
+    """Test wrapper for regressors in manifold"""
+    n_matrices, n_channels = 40, 2
+    X, y_enc = make_regression_transfer(
+        rndstate, n_matrices=n_matrices, n_channels=n_channels
+    )
+
+    tlregressor(reg, X, y_enc, source_weight, target_weight)
 
-    def make_regression_transfer(n_matrices, n_channels, rs):
-        X = make_matrices(
-            n_matrices, n_channels, "spd", rs=rndstate
-        )
-        y = get_weights(n_matrices)
-        domain = np.array(20 * ["source_domain"] + 20 * ["target_domain"])
-        _, y_enc = encode_domains(X, y, domain)
-        return X, y_enc
 
+@pytest.mark.parametrize("reg", [LinearSVR(), LinearRegression()])
+@pytest.mark.parametrize(
+    "source_weight, target_weight",
+    [(1.0, 0.0), (0.0, 1.0), (1.0, 1.0)],
+)
+def test_tlregressor_tangentspace(rndstate, reg, source_weight, target_weight):
+    """Test wrapper for regressors in tangent space"""
     n_matrices, n_channels = 40, 2
     X, y_enc = make_regression_transfer(
-        n_matrices=n_matrices, n_channels=n_channels, rs=rndstate
+        rndstate, n_matrices=n_matrices, n_channels=n_channels
     )
+    X = tangent_space(X, Cref=mean_riemann(X), metric="riemann")
+
+    tlregressor(reg, X, y_enc, source_weight, target_weight)
+
 
-    tlreg = TLRegressor(target_domain="target_domain", estimator=reg)
-    tlreg.domain_weight = {
-        "source_domain": source_domain,
-        "target_domain": target_domain,
-    }
+def tlregressor(reg, X, y_enc, source_weight, target_weight):
+    n_inputs = X.shape[0]
+
+    tlreg = TLRegressor(
+        target_domain="target_domain",
+        estimator=reg,
+        domain_weight={
+            "source_domain": source_weight,
+            "target_domain": target_weight,
+        },
+    )
 
     # test fit
     tlreg.fit(X, y_enc)
 
     # test predict
     predicted = tlreg.predict(X)
-    assert predicted.shape == (n_matrices,)
+    assert predicted.shape == (n_inputs,)
 
     # test score
     tlreg.score(X, y_enc)
 
 
+###############################################################################
+
+
 @pytest.mark.parametrize("domain_tradeoff", [0, 1])  # 0.5
 @pytest.mark.parametrize("metric", ["euclid", "logeuclid", "riemann"])
 @pytest.mark.parametrize("n_jobs", [-1, 1, 2])
@@ -445,6 +653,10 @@ def test_mdwm(rndstate, domain_tradeoff, metric, n_jobs):
     assert probabilities.shape == (n_matrices, n_classes)
     assert probabilities.sum(axis=1) == approx(np.ones(n_matrices))
 
+    # test transform
+    dists = clf.transform(X)
+    assert dists.shape == (n_matrices, n_classes)
+
     # test score
     clf.score(X, y_enc)
 
@@ -460,7 +672,7 @@ def test_mdwm_weights(rndstate, get_weights):
     weights = get_weights(n_matrices // n_classes)
 
     clf = MDWM(
-        domain_tradeoff=0.5,
+        domain_tradeoff=0.4,
         target_domain="target_domain",
         metric="riemann",
     )