From 1c053897d79a3f8550553a3407c6930707edce63 Mon Sep 17 00:00:00 2001
From: Antoine Collas <22830806+antoinecollas@users.noreply.github.com>
Date: Thu, 7 Mar 2024 10:58:22 +0100
Subject: [PATCH 1/9] [WIP] add GFK

---
 skada/__init__.py           |  2 +
 skada/_mapping.py           | 90 +++++++++++++++++++++++++++++++++++++
 skada/tests/test_mapping.py |  5 +++
 3 files changed, 97 insertions(+)

diff --git a/skada/__init__.py b/skada/__init__.py
index 06a26eeb..4cc0d88e 100644
--- a/skada/__init__.py
+++ b/skada/__init__.py
@@ -23,6 +23,7 @@
     MMDLSConSMapping,
     OTMappingAdapter,
     OTMapping,
+    GFKAdapter
 )
 from ._reweight import (
     DiscriminatorReweightDensityAdapter,
@@ -78,6 +79,7 @@
     "MMDLSConSMapping",
     "OTMappingAdapter",
     "OTMapping",
+    "GFKAdapter",
 
     "DiscriminatorReweightDensityAdapter",
     "DiscriminatorReweightDensity",
diff --git a/skada/_mapping.py b/skada/_mapping.py
index 735e03c0..808f1753 100644
--- a/skada/_mapping.py
+++ b/skada/_mapping.py
@@ -925,3 +925,93 @@ def MMDLSConSMapping(
         ),
         base_estimator,
     )
+
+
+class GFKAdapter(BaseAdapter):
+    """Domain Adaptation using an infinite number of subspaces between domains.
+
+    Geodesic Flow Kernel (GFK) maps the source and target domains to an infinite
+    dimensional space with an infinite number of subspaces. The subspaces are
+    taken from the geodesic flow of the Grassmann manifold between the source
+    and target domains.
+
+    See [5]_ for details.
+
+    Parameters
+    ----------
+    n_components : int, default=None
+        Number of components to keep. If None, all components are kept:
+        ``n_components == min(n_samples, n_features)``.
+
+    Attributes
+    ----------
+    `G_sqrtm_` : array-like, shape (n_features, n_features)
+        The square root of kernel matrix.
+
+    References
+    ----------
+    .. [5] Boqing Gong et. al. Geodesic Flow Kernel for Unsupervised Domain Adaptation
+           In CVPR, 2012.
+    """
+
+    def __init__(self, n_components=None):
+        super().__init__()
+        self.n_components = n_components
+
+    def _set_n_components(self, X):
+        if self.n_components is None:
+            self.n_components = min(X.shape[0], X.shape[1])
+
+    def _kernel_computation(self, X_source, X_target):
+        """Kernel computation for the GFK algorithm."""
+        return np.eye(X_source.shape[1])
+
+    def fit(self, X, y=None, sample_domain=None, **kwargs):
+        """Fit adaptation parameters.
+
+        Parameters
+        ----------
+        X : array-like, shape (n_samples, n_features)
+            The source data.
+        y : Ignored
+            Ignored.
+        sample_domain : array-like, shape (n_samples,)
+            The domain labels (same as sample_domain).
+
+        Returns
+        -------
+        self : object
+            Returns self.
+        """
+        X, sample_domain = check_X_domain(X, sample_domain)
+        self._set_n_components(X)
+
+        X_source, X_target = source_target_split(X, sample_domain=sample_domain)
+
+        G = self._kernel_computation(X_source, X_target)
+        self.G_sqrtm_ = _sqrtm(G)
+
+        return self
+
+    def adapt(self, X, y=None, sample_domain=None, **kwargs):
+        """Predict adaptation (weights, sample or labels).
+
+        Parameters
+        ----------
+        X : array-like, shape (n_samples, n_features)
+            The source data.
+        y : Ignored
+            Ignored.
+        sample_domain : array-like, shape (n_samples,)
+            The domain labels (same as sample_domain).
+
+        Returns
+        -------
+        X_t : array-like, shape (n_samples, n_components)
+            The data (same as X).
+        """
+        X, sample_domain = check_X_domain(X, sample_domain)
+
+        X_adapt = X @ self.G_sqrtm_.T
+
+        return X_adapt
diff --git a/skada/tests/test_mapping.py b/skada/tests/test_mapping.py
index 969a9ec6..8299de75 100644
--- a/skada/tests/test_mapping.py
+++ b/skada/tests/test_mapping.py
@@ -23,6 +23,7 @@
     CORALAdapter,
     EntropicOTMapping,
     EntropicOTMappingAdapter,
+    GFKAdapter,
     LinearOTMapping,
     LinearOTMappingAdapter,
     MMDLSConSMapping,
@@ -67,6 +68,7 @@
             MMDLSConSMapping(),
             marks=pytest.mark.skipif(not torch, reason="PyTorch not installed"),
         ),
+        make_da_pipeline(GFKAdapter(), LogisticRegression()),
     ],
 )
 def test_mapping_estimator(estimator, da_blobs_dataset):
@@ -117,6 +119,7 @@ def test_mapping_estimator(estimator, da_blobs_dataset):
             MMDLSConSMapping(Ridge()),
             marks=pytest.mark.skipif(not torch, reason="PyTorch not installed"),
         ),
+        make_da_pipeline(GFKAdapter(), Ridge()),
     ],
 )
 def test_reg_mapping_estimator(estimator, da_reg_dataset):
@@ -169,6 +172,7 @@ def _base_test_new_X_adapt(estimator, da_dataset):
             MMDLSConSMappingAdapter(gamma=1e-3),
             marks=pytest.mark.skipif(not torch, reason="PyTorch not installed"),
         ),
+        GFKAdapter(),
     ],
 )
 def test_new_X_adapt(estimator, da_dataset):
@@ -188,6 +192,7 @@ def test_new_X_adapt(estimator, da_dataset):
             MMDLSConSMappingAdapter(gamma=1e-3),
             marks=pytest.mark.skipif(not torch, reason="PyTorch not installed"),
         ),
+        GFKAdapter(),
     ],
 )
 def test_reg_new_X_adapt(estimator, da_reg_dataset):

From 76654312e4c2ec1fdcece727357a651409360f20 Mon Sep 17 00:00:00 2001
From: Antoine Collas <22830806+antoinecollas@users.noreply.github.com>
Date: Thu, 7 Mar 2024 11:55:20 +0100
Subject: [PATCH 2/9] Compute kernel

---
 skada/_mapping.py           | 65 ++++++++++++++++++++++++++++++++++++-
 skada/tests/test_mapping.py |  2 +-
 2 files changed, 65 insertions(+), 2 deletions(-)

diff --git a/skada/_mapping.py b/skada/_mapping.py
index 808f1753..4d4ffb57 100644
--- a/skada/_mapping.py
+++ b/skada/_mapping.py
@@ -9,6 +9,7 @@
 
 import numpy as np
 from ot import da
+from sklearn.decomposition import PCA
 from sklearn.metrics.pairwise import pairwise_distances
 from sklearn.svm import SVC
 
@@ -927,6 +928,41 @@ def MMDLSConSMapping(
     )
 
 
+def _gsvd(A, B):
+    """Generalized singular value decomposition.
+
+    Compute the generalized singular value decomposition of two matrices A and B.
+
+    Parameters
+    ----------
+    A : array-like, shape (m, n)
+        The first matrix.
+    B : array-like, shape (m, n)
+        The second matrix.
+
+    Returns
+    -------
+    U_A : array-like, shape (m, min(m, n))
+        The left singular vectors of A.
+    U_B : array-like, shape (m, min(m, n))
+        The left singular vectors of B.
+    S_A : array-like, shape (min(m, n),)
+        The singular values of A.
+    S_B : array-like, shape (min(m, n),)
+        The singular values of B.
+    Vt : array-like, shape (min(m, n), n)
+        The right singular vectors of A and B.
+    """
+    # TODO: implement the true gsvd
+    U_A, S_A, V_A = np.linalg.svd(A)
+    U_B, S_B, V_B = np.linalg.svd(B)
+
+    S_A = np.diag(S_A)
+    S_B = np.diag(S_B)
+
+    return U_A, U_B, S_A, S_B, V_A
+
+
 class GFKAdapter(BaseAdapter):
     """Domain Adaptation using an infinite number of subspaces between domains.
 
@@ -962,9 +998,36 @@ def _set_n_components(self, X):
         if self.n_components is None:
             self.n_components = min(X.shape[0], X.shape[1])
 
+    def _compute_pca_subspace(self, X):
+        pca = PCA(n_components=self.n_components)
+        pca.fit(X)
+        return pca.components_.T
+
     def _kernel_computation(self, X_source, X_target):
         """Kernel computation for the GFK algorithm."""
-        return np.eye(X_source.shape[1])
+        P_S = self._compute_pca_subspace(X_source)
+        P_S_fat, _, _ = np.linalg.svd(P_S, full_matrices=True)
+        R_S = P_S_fat[:, self.n_components :]
+        P_T = self._compute_pca_subspace(X_target)
+
+        U_1, U_2, Gamma, Sigma, Vt = _gsvd(P_S.T @ P_T, R_S.T @ P_T)
+
+        theta = np.arccos(Gamma)
+        Lambda_1 = 1 + (np.sin(2 * theta) / (2 * theta))
+        Lambda_2 = (np.cos(2 * theta) - 1) / (2 * theta)
+        Lambda_3 = 1 - (np.sin(2 * theta) / (2 * theta))
+
+        A = np.block([P_S @ U_1, R_S @ U_2])
+        B = np.block(
+            [
+                [np.diag(Lambda_1), np.diag(Lambda_2)],
+                [np.diag(Lambda_2), np.diag(Lambda_3)],
+            ]
+        )
+        C = np.block([[U_1.T @ P_S.T], [U_2.T @ R_S.T]])
+        G = A @ B @ C
+
+        return G
 
     def fit(self, X, y=None, sample_domain=None, **kwargs):
         """Fit adaptation parameters.
diff --git a/skada/tests/test_mapping.py b/skada/tests/test_mapping.py
index 8299de75..d92ca45f 100644
--- a/skada/tests/test_mapping.py
+++ b/skada/tests/test_mapping.py
@@ -68,7 +68,7 @@
             MMDLSConSMapping(),
             marks=pytest.mark.skipif(not torch, reason="PyTorch not installed"),
         ),
-        make_da_pipeline(GFKAdapter(), LogisticRegression()),
+        make_da_pipeline(GFKAdapter(n_components=1), LogisticRegression()),
     ],
 )
 def test_mapping_estimator(estimator, da_blobs_dataset):

From 3e42dec43b859e3be1166fdb46015c272beb7cc8 Mon Sep 17 00:00:00 2001
From: Antoine Collas <22830806+antoinecollas@users.noreply.github.com>
Date: Thu, 7 Mar 2024 14:46:04 +0100
Subject: [PATCH 3/9] check minimum eigenvalue in _sqrtm

---
 skada/_mapping.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/skada/_mapping.py b/skada/_mapping.py
index 4d4ffb57..dce9671c 100644
--- a/skada/_mapping.py
+++ b/skada/_mapping.py
@@ -518,6 +518,8 @@ def _sqrtm(C):
         Matrix inverse square root of C.
     """
     eigvals, eigvecs = np.linalg.eigh(C)
+    if np.any(eigvals < 0):
+        eigvals += 1e-14
     return (eigvecs * np.sqrt(eigvals)) @ eigvecs.T
 
 

From 2ca01c8c9dc91bf5acd4392975f19fe9598a7485 Mon Sep 17 00:00:00 2001
From: Antoine Collas <22830806+antoinecollas@users.noreply.github.com>
Date: Thu, 7 Mar 2024 14:46:27 +0100
Subject: [PATCH 4/9] rm np.diag in moke _gsvd

---
 skada/_mapping.py | 4 ----
 1 file changed, 4 deletions(-)

diff --git a/skada/_mapping.py b/skada/_mapping.py
index dce9671c..92f65050 100644
--- a/skada/_mapping.py
+++ b/skada/_mapping.py
@@ -958,10 +958,6 @@ def _gsvd(A, B):
     # TODO: implement the true gsvd
     U_A, S_A, V_A = np.linalg.svd(A)
     U_B, S_B, V_B = np.linalg.svd(B)
-
-    S_A = np.diag(S_A)
-    S_B = np.diag(S_B)
-
     return U_A, U_B, S_A, S_B, V_A
 
 

From ea933c04ab7dcc818d7bdda4240218b6a50004f5 Mon Sep 17 00:00:00 2001
From: Antoine Collas <22830806+antoinecollas@users.noreply.github.com>
Date: Thu, 7 Mar 2024 14:47:03 +0100
Subject: [PATCH 5/9] maximum n_components is min(X.shape[0], X.shape[1]) - 1

---
 skada/_mapping.py | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/skada/_mapping.py b/skada/_mapping.py
index 92f65050..ee5843fa 100644
--- a/skada/_mapping.py
+++ b/skada/_mapping.py
@@ -975,7 +975,7 @@ class GFKAdapter(BaseAdapter):
     ----------
     n_components : int, default=None
         Number of components to keep. If None, all components are kept:
-        ``n_components == min(n_samples, n_features)``.
+        ``n_components == min(n_samples, n_features) - 1``.
 
     Attributes
     ----------
@@ -994,7 +994,11 @@ def __init__(self, n_components=None):
 
     def _set_n_components(self, X):
         if self.n_components is None:
-            self.n_components = min(X.shape[0], X.shape[1])
+            self.n_components = min(X.shape[0], X.shape[1]) - 1
+        elif self.n_components > min(X.shape[0], X.shape[1]) - 1:
+            raise ValueError(
+                "n_components must be less than min(n_samples, n_features) - 1"
+            )
 
     def _compute_pca_subspace(self, X):
         pca = PCA(n_components=self.n_components)

From 3767224a567c308a854bedfe7bdac660f49a7f9b Mon Sep 17 00:00:00 2001
From: Antoine Collas <22830806+antoinecollas@users.noreply.github.com>
Date: Thu, 7 Mar 2024 14:47:53 +0100
Subject: [PATCH 6/9] fix lambda computations

---
 skada/_mapping.py | 13 +++++++++----
 1 file changed, 9 insertions(+), 4 deletions(-)

diff --git a/skada/_mapping.py b/skada/_mapping.py
index ee5843fa..efe546bf 100644
--- a/skada/_mapping.py
+++ b/skada/_mapping.py
@@ -1012,12 +1012,17 @@ def _kernel_computation(self, X_source, X_target):
         R_S = P_S_fat[:, self.n_components :]
         P_T = self._compute_pca_subspace(X_target)
 
-        U_1, U_2, Gamma, Sigma, Vt = _gsvd(P_S.T @ P_T, R_S.T @ P_T)
+        U_1, U_2, Gamma, _, Vt = _gsvd(P_S.T @ P_T, R_S.T @ P_T)
 
         theta = np.arccos(Gamma)
-        Lambda_1 = 1 + (np.sin(2 * theta) / (2 * theta))
-        Lambda_2 = (np.cos(2 * theta) - 1) / (2 * theta)
-        Lambda_3 = 1 - (np.sin(2 * theta) / (2 * theta))
+
+        Lambda_1 = 1 + np.sinc(2 * theta / np.pi)
+        Lambda_3 = 1 - np.sinc(2 * theta / np.pi)
+
+        EPS_TOL = 1e-14
+        Lambda_2 = np.zeros_like(theta)
+        Lambda_2[theta > EPS_TOL] = (np.cos(2 * theta) - 1) / (2 * theta)
+        Lambda_2[theta <= EPS_TOL] = 0
 
         A = np.block([P_S @ U_1, R_S @ U_2])
         B = np.block(

From c71038b25f919aa3c534ca7c2ff49061d8665fb4 Mon Sep 17 00:00:00 2001
From: Antoine Collas <22830806+antoinecollas@users.noreply.github.com>
Date: Fri, 8 Mar 2024 11:56:13 +0100
Subject: [PATCH 7/9] add gsvd implementation

---
 skada/_mapping.py           | 39 ++++++++++++++++++++++++++-----------
 skada/tests/test_mapping.py | 26 +++++++++++++++++++++++++
 2 files changed, 54 insertions(+), 11 deletions(-)

diff --git a/skada/_mapping.py b/skada/_mapping.py
index efe546bf..57d8c129 100644
--- a/skada/_mapping.py
+++ b/skada/_mapping.py
@@ -937,28 +937,45 @@ def _gsvd(A, B):
 
     Parameters
     ----------
-    A : array-like, shape (m, n)
+    A : array-like, shape (d, d)
         The first matrix.
-    B : array-like, shape (m, n)
+    B : array-like, shape (d, d)
         The second matrix.
 
     Returns
     -------
-    U_A : array-like, shape (m, min(m, n))
+    U_A : array-like, shape (d, d)
         The left singular vectors of A.
-    U_B : array-like, shape (m, min(m, n))
+    U_B : array-like, shape (d, d)
         The left singular vectors of B.
-    S_A : array-like, shape (min(m, n),)
+    S_A : array-like, shape (d, d)
         The singular values of A.
-    S_B : array-like, shape (min(m, n),)
+    S_B : array-like, shape (d, d)
         The singular values of B.
-    Vt : array-like, shape (min(m, n), n)
+    Vt : array-like, shape (d, d)
         The right singular vectors of A and B.
     """
-    # TODO: implement the true gsvd
-    U_A, S_A, V_A = np.linalg.svd(A)
-    U_B, S_B, V_B = np.linalg.svd(B)
-    return U_A, U_B, S_A, S_B, V_A
+    d = A.shape[0]
+    if A.shape != B.shape or A.shape[1] != d:
+        raise ValueError("A and B must have the same shape and be square matrices")
+
+    # computation as described in the construction section of
+    # https://en.wikipedia.org/wiki/Generalized_singular_value_decomposition
+    C = np.block([[A], [B]])
+    k = np.linalg.matrix_rank(C)
+    P, D, Vt = np.linalg.svd(C, full_matrices=False)
+    Vt = Vt[:d]
+    P11 = P[:d, :k]
+    P21 = P[d:, :k]
+    U_A, Sigma1, Wt = np.linalg.svd(P11, full_matrices=False)
+    P21_W = P21 @ Wt.T
+    Sigma2 = np.linalg.norm(P21_W, axis=0)
+    U_B = P21_W / Sigma2
+
+    S_A = np.diag(Sigma1) @ Wt @ np.diag(D)
+    S_B = np.diag(Sigma2) @ Wt @ np.diag(D)
+
+    return U_A, U_B, S_A, S_B, Vt
 
 
 class GFKAdapter(BaseAdapter):
diff --git a/skada/tests/test_mapping.py b/skada/tests/test_mapping.py
index d92ca45f..39b1c508 100644
--- a/skada/tests/test_mapping.py
+++ b/skada/tests/test_mapping.py
@@ -33,6 +33,7 @@
     make_da_pipeline,
     source_target_split,
 )
+from skada._mapping import _gsvd
 from skada.datasets import DomainAwareDataset
 
 
@@ -197,3 +198,28 @@ def test_new_X_adapt(estimator, da_dataset):
 )
 def test_reg_new_X_adapt(estimator, da_reg_dataset):
     _base_test_new_X_adapt(estimator, da_reg_dataset)
+
+
+def test_gsvd():
+    d = 10
+    A = np.random.multivariate_normal(np.zeros(d), np.eye(d), size=d)
+    B = np.random.multivariate_normal(np.zeros(d), np.eye(d), size=d)
+
+    # gsvd
+    U_A, U_B, S_A, S_B, Vt = _gsvd(A, B)
+
+    # shape of the matrices
+    assert U_A.shape == (d, d)
+    assert U_B.shape == (d, d)
+    assert S_A.shape == (d, d)
+    assert S_B.shape == (d, d)
+    assert Vt.shape == (d, d)
+
+    # orthogonality
+    assert np.allclose(U_A.T @ U_A, np.eye(d))
+    assert np.allclose(U_B.T @ U_B, np.eye(d))
+    assert np.allclose(Vt @ Vt.T, np.eye(d))
+
+    # check the reconstruction
+    assert np.allclose(A, U_A @ S_A @ Vt)
+    assert np.allclose(B, U_B @ S_B @ Vt)

From 55656573c4d8953a53c51989ea0c7d7974a4f28f Mon Sep 17 00:00:00 2001
From: Antoine Collas <22830806+antoinecollas@users.noreply.github.com>
Date: Fri, 8 Mar 2024 12:08:32 +0100
Subject: [PATCH 8/9] add test on singular values

---
 skada/tests/test_mapping.py | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/skada/tests/test_mapping.py b/skada/tests/test_mapping.py
index 39b1c508..badcb074 100644
--- a/skada/tests/test_mapping.py
+++ b/skada/tests/test_mapping.py
@@ -220,6 +220,10 @@ def test_gsvd():
     assert np.allclose(U_B.T @ U_B, np.eye(d))
     assert np.allclose(Vt @ Vt.T, np.eye(d))
 
+    # check condition on S_A and S_B
+    cond = S_A.T @ S_A + S_B.T @ S_B
+    assert np.allclose(cond, np.diag(np.diag(cond)))
+
     # check the reconstruction
     assert np.allclose(A, U_A @ S_A @ Vt)
     assert np.allclose(B, U_B @ S_B @ Vt)

From a283554cd9c3cf2e4bc5de892e9e7004750c7ae2 Mon Sep 17 00:00:00 2001
From: Antoine Collas <22830806+antoinecollas@users.noreply.github.com>
Date: Fri, 8 Mar 2024 15:45:45 +0100
Subject: [PATCH 9/9] add gsvd test based on GFK

---
 skada/tests/test_mapping.py | 18 ++++++++++++++++++
 1 file changed, 18 insertions(+)

diff --git a/skada/tests/test_mapping.py b/skada/tests/test_mapping.py
index badcb074..7abfed4d 100644
--- a/skada/tests/test_mapping.py
+++ b/skada/tests/test_mapping.py
@@ -227,3 +227,21 @@ def test_gsvd():
     # check the reconstruction
     assert np.allclose(A, U_A @ S_A @ Vt)
     assert np.allclose(B, U_B @ S_B @ Vt)
+
+    n, d, k = 100, 10, 5
+    Xs = np.random.multivariate_normal(np.zeros(d), np.eye(d), size=n)
+    Xt = np.random.multivariate_normal(np.zeros(d), np.eye(d), size=n)
+
+    Ps, _, _ = np.linalg.svd(Xs.T, full_matrices=False)
+    Ps, Rs = Ps[:, :k], Ps[:, k:]
+    Pt, _, _ = np.linalg.svd(Xt.T, full_matrices=False)
+    Pt = Pt[:, :k]
+
+    U1, U2, Gamma, Sigma, Vt = _gsvd(Ps.T @ Pt, -Rs.T @ Pt)
+
+    # assert Gamma and Sigma are diagonal
+    import ipdb
+
+    ipdb.set_trace()
+    assert np.allclose(Gamma, np.diag(np.diag(Gamma)))
+    assert np.allclose(Sigma, np.diag(np.diag(Sigma)))