diff --git a/README.md b/README.md
index 830e82c..b558de5 100644
--- a/README.md
+++ b/README.md
@@ -20,6 +20,7 @@ The package is part of the [PyKale](https://github.com/pykale/pykale) ecosystem
- Multilinear Principal Component Analysis (`MPCA`) [[1](#references)]
- Transferable / generalizable feature extraction across domains or groups:
- Correlation Alignment (`CORAL`) [[13](#references)]
+ - Geodesic Flow Kernel (`GFK`) [[14](#references)]
- Transfer Component Analysis (`TCA`) [[2](#references)]
- Joint Distribution Adaptation (`JDA`) [[3](#references)]
- Balanced Distribution Adaptation (`BDA`) [[4](#references)]
@@ -38,6 +39,10 @@ The package is part of the [PyKale](https://github.com/pykale/pykale) ecosystem
- scikit-learn style `fit`, `transform`, `predict`, `fit_transform`, and `fit_predict` workflows where applicable.
- Most methods accept additional `covariates` — e.g., domain or group labels — alongside `X` and `y`, with optional one-hot encoding for categorical values; multiblock transformers (CIFE, AJIVE) take `groups` to specify block membership.
+
+
+
+
KaleLinear requires Python 3.10 or later. Core dependencies include:
- [NumPy](http://www.numpy.org/)
@@ -140,6 +145,8 @@ Worked examples for the main transformers and estimators are collected in
[13] Sun, B., Feng, J. and Saenko, K., 2016. [Return of frustratingly easy domain adaptation](https://ojs.aaai.org/index.php/AAAI/article/view/10306). In _Proceedings of the AAAI Conference on Artificial Intelligence_ (Vol. 30, No. 1, pp. 2058-2065).
+[14] Gong, B., Shi, Y., Sha, F. and Grauman, K., 2012. [Geodesic flow kernel for unsupervised domain adaptation](https://ieeexplore.ieee.org/document/6247911). In _IEEE Conference on Computer Vision and Pattern Recognition_ (pp. 2066-2073).
+
## Other open domain adaptation repositories
- [POT: Python Optimal Transport](https://github.com/rflamary/POT)
diff --git a/docs/images/api_map.jpg b/docs/images/api_map.jpg
new file mode 100644
index 0000000..8b6d075
Binary files /dev/null and b/docs/images/api_map.jpg differ
diff --git a/docs/images/api_map.svg b/docs/images/api_map.svg
new file mode 100644
index 0000000..eba5f13
--- /dev/null
+++ b/docs/images/api_map.svg
@@ -0,0 +1,104 @@
+
diff --git a/docs/source/api_transformers.rst b/docs/source/api_transformers.rst
index 1f11563..dab4429 100644
--- a/docs/source/api_transformers.rst
+++ b/docs/source/api_transformers.rst
@@ -31,6 +31,11 @@ The transformer classes are also available through the PyKale-compatible alias
:undoc-members:
:show-inheritance:
+.. autoclass:: GFK
+ :members:
+ :undoc-members:
+ :show-inheritance:
+
.. autoclass:: MIDA
:members:
:undoc-members:
diff --git a/docs/source/introduction.rst b/docs/source/introduction.rst
index afcbd33..0381ff0 100644
--- a/docs/source/introduction.rst
+++ b/docs/source/introduction.rst
@@ -14,7 +14,7 @@ covariates, side information, or unlabeled target samples.
Main Features
-------------
-* Transformer models for learning feature embeddings: CORAL, MPCA, TCA, JDA,
+* Transformer models for learning feature embeddings: CORAL, GFK, MPCA, TCA, JDA,
BDA, MIDA, CIFE, and AJIVE.
* Estimator models for classification and adaptation: LapSVM, LapRLS, ARSVM,
ARRLS, CoIRSVM, CoIRLS, and GSDA.
diff --git a/docs/source/usage.rst b/docs/source/usage.rst
index aa68fe9..d61c423 100644
--- a/docs/source/usage.rst
+++ b/docs/source/usage.rst
@@ -6,7 +6,7 @@ Usage
Domain Adaptation Transformers
------------------------------
-TCA, JDA, BDA, and CORAL take all samples in a single input array and receive
+TCA, JDA, BDA, CORAL, and GFK take all samples in a single input array and receive
domain labels through ``covariates``. Use ``target_covariate`` to identify
which domain label is the target domain. CORAL is asymmetric: after fitting, it
aligns source samples to the target covariance and centers target samples.
diff --git a/kalelinear/embed.py b/kalelinear/embed.py
index 716f480..43e0202 100644
--- a/kalelinear/embed.py
+++ b/kalelinear/embed.py
@@ -1,5 +1,5 @@
"""Embedding models exposed with a PyKale-style API."""
-from kalelinear.transformer import AJIVE, BDA, CIFE, CORAL, JDA, MIDA, MPCA, TCA
+from kalelinear.transformer import AJIVE, BDA, CIFE, CORAL, GFK, JDA, MIDA, MPCA, TCA
-__all__ = ["TCA", "JDA", "BDA", "CORAL", "MIDA", "MPCA", "CIFE", "AJIVE"]
+__all__ = ["TCA", "JDA", "BDA", "CORAL", "GFK", "MIDA", "MPCA", "CIFE", "AJIVE"]
diff --git a/kalelinear/transformer/__init__.py b/kalelinear/transformer/__init__.py
index 6adf575..6719831 100644
--- a/kalelinear/transformer/__init__.py
+++ b/kalelinear/transformer/__init__.py
@@ -1,6 +1,7 @@
from kalelinear.transformer._ajive import AJIVE
from kalelinear.transformer._cife import CIFE
from kalelinear.transformer._coral import CORAL
+from kalelinear.transformer._gfk import GFK
from kalelinear.transformer._jda import BDA, JDA
from kalelinear.transformer._mida import MIDA
from kalelinear.transformer._mpca import MPCA
@@ -11,6 +12,7 @@
"JDA",
"BDA",
"CORAL",
+ "GFK",
"MIDA",
"MPCA",
"CIFE",
diff --git a/kalelinear/transformer/_gfk.py b/kalelinear/transformer/_gfk.py
new file mode 100644
index 0000000..6408e23
--- /dev/null
+++ b/kalelinear/transformer/_gfk.py
@@ -0,0 +1,224 @@
+# =============================================================================
+# @author: Shuo Zhou, The University of Sheffield, shuo.zhou@sheffield.ac.uk
+# =============================================================================
+"""Geodesic Flow Kernel (GFK)."""
+
+from numbers import Integral
+
+import numpy as np
+from sklearn.base import _fit_context, BaseEstimator, ClassNamePrefixFeaturesOutMixin, TransformerMixin
+from sklearn.utils._param_validation import Interval
+from sklearn.utils.validation import check_is_fitted, validate_data
+
+from kalelinear._domain import check_binary_domain_covariates, split_domain_indices
+from kalelinear.transformer._coral import _symmetric_sqrt
+
+# Below this principal angle (radians) the closed-form weights are replaced by
+# their Taylor expansions to avoid division by (near-)zero angles.
+_SMALL_ANGLE = 1e-2
+
+
+def _pca_basis(X, n_components):
+ """Return the top ``n_components`` principal directions of ``X`` as columns."""
+ _, _, vt = np.linalg.svd(X - X.mean(axis=0), full_matrices=False)
+ return vt[:n_components].T
+
+
+def _geodesic_flow_kernel(source_basis, target_basis):
+ """Return the geodesic flow kernel matrix between two subspaces.
+
+ Both inputs are orthonormal bases of shape ``(n_features, n_components)``.
+ The kernel is computed in closed form from the principal angles between
+ the subspaces, without forming the orthogonal complement of the source
+ subspace explicitly. Principal angles close to zero are handled with
+ Taylor expansions, so identical or overlapping subspaces are supported.
+ """
+ # Principal angles: Ps^T Pt = V1 cos(Theta) V^T
+ v1, cos_theta, vt = np.linalg.svd(source_basis.T @ target_basis, full_matrices=False)
+ v = vt.T
+ theta = np.arccos(np.clip(cos_theta, -1.0, 1.0))
+
+ # Component of the target basis orthogonal to the source subspace.
+ u1 = source_basis @ v1
+ target_rotated = target_basis @ v
+ residual = target_rotated - source_basis @ (source_basis.T @ target_rotated)
+
+ small = theta < _SMALL_ANGLE
+ safe_theta = np.where(small, 1.0, theta)
+
+ # sin(2θ) / (2θ)
+ sinc_2theta = np.sinc(2 * theta / np.pi)
+ lambda1 = 1.0 + sinc_2theta
+ # sin(θ) / θ
+ ratio = np.sinc(theta / np.pi)
+ # (2θ - sin(2θ)) / (2 sin²θ θ) -> 2/3 as θ -> 0
+ w_direct = (2 * safe_theta - np.sin(2 * safe_theta)) / (2 * safe_theta * np.sin(safe_theta) ** 2)
+ w_series = (2 / 3 - 2 * theta**2 / 15 + 4 * theta**4 / 315) / (1 - theta**2 / 3 + 2 * theta**4 / 45)
+ weight = np.where(small, w_series, w_direct)
+
+ cross = residual * ratio
+ kernel = (u1 * lambda1) @ u1.T + u1 @ cross.T + cross @ u1.T + (residual * weight) @ residual.T
+ return (kernel + kernel.T) / 2
+
+
+class GFK(ClassNamePrefixFeaturesOutMixin, TransformerMixin, BaseEstimator):
+ """Geodesic Flow Kernel (GFK) for unsupervised domain adaptation.
+
+ GFK models the source and target domains as points (subspaces spanned by
+ their leading principal components) on a Grassmann manifold and
+ integrates over all intermediate subspaces along the geodesic flow
+ connecting them. The resulting positive semi-definite kernel has a closed
+ form in terms of the principal angles :math:`\\theta_i` between the two
+ subspaces:
+
+ .. math::
+
+ G = Q \\begin{bmatrix} V_1 & 0 \\\\ 0 & V_2 \\end{bmatrix}
+ \\begin{bmatrix} \\Lambda_1 & \\Lambda_2 \\\\ \\Lambda_2 & \\Lambda_3 \\end{bmatrix}
+ \\begin{bmatrix} V_1 & 0 \\\\ 0 & V_2 \\end{bmatrix}^\\top Q^\\top,
+
+ with :math:`\\lambda_{1i} = 1 + \\sin(2\\theta_i) / (2\\theta_i)`,
+ :math:`\\lambda_{2i} = (\\cos(2\\theta_i) - 1) / (2\\theta_i)` and
+ :math:`\\lambda_{3i} = 1 - \\sin(2\\theta_i) / (2\\theta_i)`. The kernel
+ is evaluated stably for principal angles close to zero.
+
+ :meth:`transform` embeds samples as :math:`G^{1/2} x`, so that inner
+ products in the embedded space equal the geodesic flow kernel. The same
+ mapping is applied to source and target samples.
+
+ ``covariates`` represent binary domain labels of length ``n_samples``.
+ They must contain both source and target domains during :meth:`fit`.
+ ``target_covariate`` selects which label is treated as the target domain.
+
+ Parameters
+ ----------
+ n_components : int, default=None
+ Dimensionality of the source and target subspaces. It must not exceed
+ ``min(n_features, n_source - 1, n_target - 1)``. If ``None``, it is
+ set to ``max(1, n_features // 2)`` capped by the above limit.
+
+ Attributes
+ ----------
+ source_basis_ : ndarray of shape (n_features, n_components_)
+ Orthonormal basis of the source PCA subspace.
+ target_basis_ : ndarray of shape (n_features, n_components_)
+ Orthonormal basis of the target PCA subspace.
+ principal_angles_ : ndarray of shape (n_components_,)
+ Principal angles (radians) between the source and target subspaces.
+ kernel_ : ndarray of shape (n_features, n_features)
+ Geodesic flow kernel matrix :math:`G`.
+ transform_matrix_ : ndarray of shape (n_features, n_features)
+ Symmetric square root :math:`G^{1/2}` used by :meth:`transform`.
+ n_components_ : int
+ Effective subspace dimensionality.
+ target_covariate_ : object
+ Domain label treated as the target domain.
+
+ References
+ ----------
+ Gong, B., Shi, Y., Sha, F. and Grauman, K., 2012. Geodesic flow kernel
+ for unsupervised domain adaptation. In *IEEE Conference on Computer
+ Vision and Pattern Recognition* (pp. 2066-2073).
+ """
+
+ _parameter_constraints: dict = {
+ "n_components": [Interval(Integral, 1, None, closed="left"), None],
+ }
+
+ def __init__(self, n_components=None):
+ self.n_components = n_components
+
+ @_fit_context(prefer_skip_nested_validation=True)
+ def fit(self, X, y=None, covariates=None, target_covariate=None):
+ """Fit GFK on source and target samples.
+
+ Parameters
+ ----------
+ X : array-like of shape (n_samples, n_features)
+ Source and target samples stacked along the first axis.
+ y : array-like, default=None
+ Ignored. Present for scikit-learn API consistency. GFK is an
+ unsupervised domain adaptation method.
+ covariates : array-like of shape (n_samples,), default=None
+ Binary domain label for every sample. Must contain both source
+ and target domain values.
+ target_covariate : object, default=None
+ Domain label treated as the target domain. Defaults to the last
+ unique covariate value, as in the other KaleLinear adapters.
+
+ Returns
+ -------
+ self : GFK
+ Fitted transformer.
+ """
+ X = validate_data(self, X, dtype=[np.float64, np.float32])
+ n_samples, n_features = X.shape
+
+ if covariates is None:
+ raise ValueError("Covariates must be provided for GFK during `fit`.")
+
+ covariates, unique_covariates = check_binary_domain_covariates(
+ covariates,
+ n_samples,
+ require_numeric=True,
+ error_prefix=f"Covariates for {self.__class__.__name__}",
+ both_domains_message=(
+ f"Covariates for {self.__class__.__name__} must contain both source and target domain values."
+ ),
+ )
+ split = split_domain_indices(covariates, target_covariate)
+
+ X_source = X[split.source_idx]
+ X_target = X[split.target_idx]
+ if X_source.shape[0] < 2:
+ raise ValueError("At least two source samples are required to estimate the source subspace.")
+ if X_target.shape[0] < 2:
+ raise ValueError("At least two target samples are required to estimate the target subspace.")
+
+ max_components = min(n_features, X_source.shape[0] - 1, X_target.shape[0] - 1)
+ if self.n_components is None:
+ n_components = min(max(1, n_features // 2), max_components)
+ elif self.n_components > max_components:
+ raise ValueError(
+ f"n_components={self.n_components} must be less than or equal to "
+ f"min(n_features, n_source - 1, n_target - 1) = {max_components}."
+ )
+ else:
+ n_components = self.n_components
+
+ self.source_basis_ = _pca_basis(X_source, n_components)
+ self.target_basis_ = _pca_basis(X_target, n_components)
+ cos_theta = np.linalg.svd(self.source_basis_.T @ self.target_basis_, compute_uv=False)
+ self.principal_angles_ = np.arccos(np.clip(cos_theta, -1.0, 1.0))
+ self.kernel_ = _geodesic_flow_kernel(self.source_basis_, self.target_basis_)
+ self.transform_matrix_ = _symmetric_sqrt(self.kernel_)
+ self.n_components_ = n_components
+ self.target_covariate_ = split.target_covariate
+ self.domain_values_ = unique_covariates
+ self._n_features_out = n_features
+ return self
+
+ def transform(self, X):
+ """Embed samples with the square root of the geodesic flow kernel.
+
+ Parameters
+ ----------
+ X : array-like of shape (n_samples, n_features)
+ Samples to transform. The same mapping applies to source and
+ target samples.
+
+ Returns
+ -------
+ X_new : ndarray of shape (n_samples, n_features)
+ Transformed samples. Feature dimensionality is preserved.
+ """
+ check_is_fitted(self, "kernel_")
+ X = validate_data(self, X, dtype=[np.float64, np.float32], reset=False)
+ return X @ self.transform_matrix_
+
+ def fit_transform(self, X, y=None, covariates=None, target_covariate=None):
+ """Fit GFK and transform the source and target samples.
+
+ Parameters are the same as for :meth:`fit`.
+ """
+ return self.fit(X, y=y, covariates=covariates, target_covariate=target_covariate).transform(X)
diff --git a/tests/test_public_api.py b/tests/test_public_api.py
index fbe8d0c..f3d79e7 100644
--- a/tests/test_public_api.py
+++ b/tests/test_public_api.py
@@ -9,6 +9,7 @@ def test_embed_module_exposes_transformers():
assert embed.JDA is transformer.JDA
assert embed.BDA is transformer.BDA
assert embed.CORAL is transformer.CORAL
+ assert embed.GFK is transformer.GFK
assert embed.MIDA is transformer.MIDA
assert embed.MPCA is transformer.MPCA
assert embed.CIFE is transformer.CIFE
diff --git a/tests/transformer/test_gfk.py b/tests/transformer/test_gfk.py
new file mode 100644
index 0000000..4d291dc
--- /dev/null
+++ b/tests/transformer/test_gfk.py
@@ -0,0 +1,129 @@
+import numpy as np
+import pytest
+from numpy import testing
+from sklearn.exceptions import NotFittedError
+
+from kalelinear.transformer import GFK
+from kalelinear.transformer._gfk import _geodesic_flow_kernel, _pca_basis
+
+
+@pytest.fixture(scope="module")
+def gfk_data():
+ rng = np.random.default_rng(0)
+ n_source, n_target, n_features = 60, 80, 6
+
+ X_source = rng.normal(size=(n_source, n_features)) * np.array([3.0, 2.0, 1.0, 0.3, 0.2, 0.1])
+ X_target = rng.normal(size=(n_target, n_features)) * np.array([0.1, 0.3, 1.0, 2.0, 3.0, 0.2])
+ domains = np.concatenate((np.zeros(n_source, dtype=int), np.ones(n_target, dtype=int)))
+ X = np.vstack((X_source, X_target))
+ return X, domains, X_source, X_target
+
+
+def _numerical_kernel(source_basis, target_basis, n_steps=4000):
+ """Integrate the geodesic flow numerically (twice the integral, as in Gong et al.)."""
+ n_features = source_basis.shape[0]
+ v1, cos_theta, vt = np.linalg.svd(source_basis.T @ target_basis, full_matrices=False)
+ theta = np.arccos(np.clip(cos_theta, -1, 1))
+ u1 = source_basis @ v1
+ target_rotated = target_basis @ vt.T
+ residual = target_rotated - source_basis @ (source_basis.T @ target_rotated)
+ t = (np.arange(n_steps) + 0.5) / n_steps
+ kernel = np.zeros((n_features, n_features))
+ for ti in t:
+ phi = u1 * np.cos(ti * theta) + residual * (np.sin(ti * theta) / np.sin(theta))
+ kernel += phi @ phi.T
+ return 2 * kernel / n_steps
+
+
+def test_gfk_kernel_matches_numerical_integration(gfk_data):
+ _, _, X_source, X_target = gfk_data
+ ps, pt = _pca_basis(X_source, 3), _pca_basis(X_target, 3)
+
+ testing.assert_allclose(_geodesic_flow_kernel(ps, pt), _numerical_kernel(ps, pt), atol=1e-5)
+
+
+def test_gfk_transform_reproduces_kernel(gfk_data):
+ X, domains, _, _ = gfk_data
+ gfk = GFK(n_components=3).fit(X, covariates=domains, target_covariate=1)
+ z = gfk.transform(X)
+
+ assert z.shape == X.shape
+ testing.assert_allclose(z @ z.T, X @ gfk.kernel_ @ X.T, atol=1e-8)
+ testing.assert_allclose(gfk.kernel_, gfk.kernel_.T)
+ assert np.linalg.eigvalsh(gfk.kernel_).min() > -1e-10
+ testing.assert_allclose(gfk.fit_transform(X, covariates=domains, target_covariate=1), z)
+
+
+def test_gfk_identical_subspaces_gives_scaled_projection(gfk_data):
+ _, _, X_source, _ = gfk_data
+ ps = _pca_basis(X_source, 3)
+ kernel = _geodesic_flow_kernel(ps, ps)
+
+ assert np.isfinite(kernel).all()
+ testing.assert_allclose(kernel, 2 * ps @ ps.T, atol=1e-10)
+
+
+def test_gfk_near_zero_angles_are_continuous(gfk_data):
+ _, _, X_source, _ = gfk_data
+ ps = _pca_basis(X_source, 3)
+ rng = np.random.default_rng(2)
+ direction = rng.normal(size=ps.shape)
+ kernels = []
+ for eps in (1e-1, 2e-2, 1e-2, 5e-3, 1e-6):
+ pt, _ = np.linalg.qr(ps + eps * direction)
+ kernels.append(_geodesic_flow_kernel(ps, pt))
+ for previous, current in zip(kernels, kernels[1:]):
+ assert np.linalg.norm(current - previous) < 0.5
+ testing.assert_allclose(kernels[-1], 2 * ps @ ps.T, atol=1e-4)
+
+
+def test_gfk_attributes_and_feature_names(gfk_data):
+ X, domains, _, _ = gfk_data
+ gfk = GFK().fit(X, covariates=domains)
+
+ assert gfk.n_components_ == X.shape[1] // 2
+ assert gfk.principal_angles_.shape == (gfk.n_components_,)
+ assert gfk.target_covariate_ == 1
+ assert gfk.kernel_.shape == (X.shape[1], X.shape[1])
+ testing.assert_array_equal(gfk.get_feature_names_out(), np.array([f"gfk{i}" for i in range(X.shape[1])]))
+
+
+def test_gfk_transform_requires_fit(gfk_data):
+ X, _, _, _ = gfk_data
+ with pytest.raises(NotFittedError, match="not fitted"):
+ GFK().transform(X)
+
+
+def test_gfk_fit_requires_covariates(gfk_data):
+ X, _, _, _ = gfk_data
+ with pytest.raises(ValueError, match="Covariates must be provided"):
+ GFK().fit(X)
+
+
+def test_gfk_fit_requires_both_domains(gfk_data):
+ X, _, _, _ = gfk_data
+ with pytest.raises(ValueError, match="both source and target"):
+ GFK().fit(X, covariates=np.zeros(X.shape[0], dtype=int))
+
+
+def test_gfk_requires_two_samples_per_domain(gfk_data):
+ X, _, _, _ = gfk_data
+ with pytest.raises(ValueError, match="At least two target samples"):
+ GFK().fit(X[:3], covariates=np.array([0, 0, 1]), target_covariate=1)
+
+
+def test_gfk_validates_parameters(gfk_data):
+ X, domains, _, _ = gfk_data
+ with pytest.raises(ValueError, match="n_components"):
+ GFK(n_components=0).fit(X, covariates=domains)
+ with pytest.raises(ValueError, match="n_components"):
+ GFK(n_components=X.shape[1] + 1).fit(X, covariates=domains)
+ with pytest.raises(ValueError, match="target_covariate"):
+ GFK().fit(X, covariates=domains, target_covariate=2)
+
+
+def test_gfk_transform_validates_features(gfk_data):
+ X, domains, _, _ = gfk_data
+ gfk = GFK().fit(X, covariates=domains)
+ with pytest.raises(ValueError, match="features"):
+ gfk.transform(X[:, :2])
+