diff --git a/README.md b/README.md
index f87a65c..830e82c 100644
--- a/README.md
+++ b/README.md
@@ -19,6 +19,7 @@ The package is part of the [PyKale](https://github.com/pykale/pykale) ecosystem
   - Dimension reduction for multiview tensor data:
     - Multilinear Principal Component Analysis (`MPCA`) [[1](#references)]
   - Transferable / generalizable feature extraction across domains or groups:
+    - Correlation Alignment (`CORAL`) [[13](#references)]
     - Transfer Component Analysis (`TCA`) [[2](#references)]
     - Joint Distribution Adaptation (`JDA`) [[3](#references)]
     - Balanced Distribution Adaptation (`BDA`) [[4](#references)]
@@ -97,7 +98,7 @@ See [CONTRIBUTING.md](CONTRIBUTING.md) for detailed contribution guidelines.
 ### Public API
 
 ```python
-from kalelinear.transformer import BDA, JDA, MIDA, MPCA, TCA
+from kalelinear.transformer import BDA, CORAL, JDA, MIDA, MPCA, TCA
 from kalelinear.estimator import ARRLS, ARSVM, CoIRLS, CoIRSVM, GSDA, LapRLS, LapSVM
 ```
 
@@ -105,6 +106,7 @@ Worked examples for the main transformers and estimators are collected in
 [Tutorials](TUTORIALS.md):
 
 - Learn a domain-invariant embedding with TCA
+- Align source and target features with CORAL
 - Use MIDA with categorical covariates
 - Extract common and individual features across groups with CIFE or AJIVE
 - Train a domain adaptation classifier (ARSVM, ARRLS)
@@ -136,6 +138,8 @@ Worked examples for the main transformers and estimators are collected in
 
 [12] Feng, Q., Jiang, M., Hannig, J. and Marron, J.S., 2018. [Angle-based joint and individual variation explained](https://www.sciencedirect.com/science/article/pii/S0047259X1730204X). _Journal of Multivariate Analysis_, 166, pp.241-265.
 
+[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).
+
 ## Other open domain adaptation repositories
 
 - [POT: Python Optimal Transport](https://github.com/rflamary/POT)
diff --git a/TUTORIALS.md b/TUTORIALS.md
index 8941c36..aa79361 100644
--- a/TUTORIALS.md
+++ b/TUTORIALS.md
@@ -37,6 +37,39 @@ TCA, JDA, and BDA take domain labels through `covariates`. They do not accept
 separate source and target arrays; stack samples into one array and use
 `target_covariate` to identify the target domain.
 
+### Align Source and Target Features with CORAL
+
+CORAL is an asymmetric domain alignment method: it whitens the source
+covariance and recolors it with the target covariance. As with TCA, pass all
+samples in one array with binary domain labels. `fit_transform` centers each
+domain by its own mean, aligns the source samples, and leaves the target
+samples in the same centered space.
+
+```python
+import numpy as np
+from kalelinear.transformer import CORAL
+
+X = np.array(
+    [
+        [-2.0, -1.8],
+        [-1.8, -2.1],
+        [1.9, 1.7],
+        [2.1, 2.0],
+        [-1.4, -1.2],
+        [-1.2, -1.1],
+        [1.2, 1.1],
+        [1.4, 1.3],
+    ]
+)
+domain_labels = np.array([0, 0, 0, 0, 1, 1, 1, 1])
+
+transformer = CORAL()
+z = transformer.fit_transform(X, covariates=domain_labels, target_covariate=1)
+
+z_source = z[domain_labels == 0]
+z_target = z[domain_labels == 1]
+```
+
 ### Use MIDA with Categorical Domain Covariates
 
 ```python
diff --git a/docs/source/api_transformers.rst b/docs/source/api_transformers.rst
index d20ac72..1f11563 100644
--- a/docs/source/api_transformers.rst
+++ b/docs/source/api_transformers.rst
@@ -26,6 +26,11 @@ The transformer classes are also available through the PyKale-compatible alias
    :undoc-members:
    :show-inheritance:
 
+.. autoclass:: CORAL
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
 .. autoclass:: MIDA
    :members:
    :undoc-members:
diff --git a/docs/source/introduction.rst b/docs/source/introduction.rst
index a21aa52..afcbd33 100644
--- a/docs/source/introduction.rst
+++ b/docs/source/introduction.rst
@@ -14,8 +14,8 @@ covariates, side information, or unlabeled target samples.
 Main Features
 -------------
 
-* Transformer models for learning feature embeddings: MPCA, TCA, JDA, BDA,
-  MIDA, CIFE, and AJIVE.
+* Transformer models for learning feature embeddings: CORAL, MPCA, TCA, JDA,
+  BDA, MIDA, CIFE, and AJIVE.
 * Estimator models for classification and adaptation: LapSVM, LapRLS, ARSVM,
   ARRLS, CoIRSVM, CoIRLS, and GSDA.
 * NumPy-compatible inputs and outputs.
diff --git a/docs/source/tutorial.rst b/docs/source/tutorial.rst
index a83442d..898d560 100644
--- a/docs/source/tutorial.rst
+++ b/docs/source/tutorial.rst
@@ -33,6 +33,40 @@ label is the target domain.
    z_source = z[domain_labels == 0]
    z_target = z[domain_labels == 1]
 
+Align Source and Target Features with CORAL
+-------------------------------------------
+
+CORAL is an asymmetric domain alignment method: it whitens the source
+covariance and recolors it with the target covariance. As with TCA, pass all
+samples in one array with binary domain labels. ``fit_transform`` centers each
+domain by its own mean, aligns the source samples, and leaves the target
+samples in the same centered space.
+
+.. code-block:: python
+
+   import numpy as np
+   from kalelinear.transformer import CORAL
+
+   X = np.array(
+       [
+           [-2.0, -1.8],
+           [-1.8, -2.1],
+           [1.9, 1.7],
+           [2.1, 2.0],
+           [-1.4, -1.2],
+           [-1.2, -1.1],
+           [1.2, 1.1],
+           [1.4, 1.3],
+       ]
+   )
+   domain_labels = np.array([0, 0, 0, 0, 1, 1, 1, 1])
+
+   transformer = CORAL()
+   z = transformer.fit_transform(X, covariates=domain_labels, target_covariate=1)
+
+   z_source = z[domain_labels == 0]
+   z_target = z[domain_labels == 1]
+
 Use MIDA with Categorical Covariates
 ------------------------------------
 
diff --git a/docs/source/usage.rst b/docs/source/usage.rst
index 29f1116..aa68fe9 100644
--- a/docs/source/usage.rst
+++ b/docs/source/usage.rst
@@ -6,9 +6,10 @@ Usage
 Domain Adaptation Transformers
 ------------------------------
 
-TCA, JDA, and BDA 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.
+TCA, JDA, BDA, and CORAL 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.
 
 .. code-block:: python
 
diff --git a/kalelinear/embed.py b/kalelinear/embed.py
index 3cb98b5..716f480 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, JDA, MIDA, MPCA, TCA
+from kalelinear.transformer import AJIVE, BDA, CIFE, CORAL, JDA, MIDA, MPCA, TCA
 
-__all__ = ["TCA", "JDA", "BDA", "MIDA", "MPCA", "CIFE", "AJIVE"]
+__all__ = ["TCA", "JDA", "BDA", "CORAL", "MIDA", "MPCA", "CIFE", "AJIVE"]
diff --git a/kalelinear/transformer/__init__.py b/kalelinear/transformer/__init__.py
index b0c44bb..6adf575 100644
--- a/kalelinear/transformer/__init__.py
+++ b/kalelinear/transformer/__init__.py
@@ -1,5 +1,6 @@
 from kalelinear.transformer._ajive import AJIVE
 from kalelinear.transformer._cife import CIFE
+from kalelinear.transformer._coral import CORAL
 from kalelinear.transformer._jda import BDA, JDA
 from kalelinear.transformer._mida import MIDA
 from kalelinear.transformer._mpca import MPCA
@@ -9,6 +10,7 @@
     "TCA",
     "JDA",
     "BDA",
+    "CORAL",
     "MIDA",
     "MPCA",
     "CIFE",
diff --git a/kalelinear/transformer/_coral.py b/kalelinear/transformer/_coral.py
new file mode 100644
index 0000000..992f8af
--- /dev/null
+++ b/kalelinear/transformer/_coral.py
@@ -0,0 +1,224 @@
+# =============================================================================
+# @author: Shuo Zhou, The University of Sheffield, shuo.zhou@sheffield.ac.uk
+# =============================================================================
+"""Correlation Alignment (CORAL)."""
+
+from numbers import Real
+
+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
+
+
+def _regularized_covariance(X, ridge):
+    """Return the sample covariance of ``X`` with ridge regularization.
+
+    The sample covariance is centered and normalized by ``n_samples - 1``,
+    matching the classical CORAL implementation. ``ridge`` is added to the
+    diagonal to keep the whitening step stable when the covariance is
+    (close to) singular.
+    """
+    _, n_features = X.shape
+    covariance = np.cov(X, rowvar=False)
+    return covariance + ridge * np.eye(n_features, dtype=covariance.dtype)
+
+
+def _symmetric_sqrt(matrix):
+    """Return the symmetric positive square root of a symmetric matrix."""
+    eigenvalues, eigenvectors = np.linalg.eigh(matrix)
+    eigenvalues = np.clip(eigenvalues, 0, None)
+    return (eigenvectors * np.sqrt(eigenvalues)) @ eigenvectors.T
+
+
+def _symmetric_inv_sqrt(matrix):
+    """Return the symmetric inverse square root via a PSD pseudo-inverse."""
+    eigenvalues, eigenvectors = np.linalg.eigh(matrix)
+    eigenvalues = np.clip(eigenvalues, 0, None)
+    max_eigenvalue = eigenvalues.max()
+    if max_eigenvalue <= 0:
+        raise ValueError(
+            "The source covariance matrix has no positive eigenvalue. "
+            "CORAL cannot compute a whitening transform; try a positive `lambda_`."
+        )
+
+    tolerance = 10 * np.finfo(eigenvalues.dtype).eps * max(1.0, max_eigenvalue) * eigenvalues.shape[0]
+    keep = eigenvalues > tolerance
+    inverse_sqrt = np.zeros_like(eigenvalues)
+    inverse_sqrt[keep] = 1.0 / np.sqrt(eigenvalues[keep])
+    return (eigenvectors * inverse_sqrt) @ eigenvectors.T
+
+
+class CORAL(ClassNamePrefixFeaturesOutMixin, TransformerMixin, BaseEstimator):
+    """Correlation Alignment (CORAL) for unsupervised domain adaptation.
+
+    CORAL aligns the second-order statistics of the source and target
+    feature distributions by whitening the (mean-centered) source features
+    and recoloring them with the target covariance:
+
+    .. math::
+
+        A = (C_S + \\lambda I)^{-1/2} (C_T + \\lambda I)^{1/2},
+
+    where :math:`C_S` and :math:`C_T` are the source and target sample
+    covariance matrices. Samples are centered by their domain mean during
+    :meth:`fit_transform`, so both domains are embedded in a common
+    mean-centered feature space. Source samples are transformed with the
+    learned alignment :math:`A`; target samples are only centered.
+
+    ``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
+    ----------
+    lambda_ : float, default=1e-5
+        Regularization added to the diagonal of the source and target
+        covariance matrices before computing the alignment. Use a larger
+        value when the covariance estimates are unstable (e.g. very few
+        samples or near-constant features).
+
+    Attributes
+    ----------
+    source_mean_ : ndarray of shape (n_features,)
+        Mean of the source training samples.
+    target_mean_ : ndarray of shape (n_features,)
+        Mean of the target training samples.
+    source_covariance_ : ndarray of shape (n_features, n_features)
+        Regularized source covariance matrix.
+    target_covariance_ : ndarray of shape (n_features, n_features)
+        Regularized target covariance matrix.
+    alignment_ : ndarray of shape (n_features, n_features)
+        Linear whitening/recoloring transform applied to source features.
+    target_covariate_ : object
+        Domain label treated as the target domain.
+
+    References
+    ----------
+    Sun, B., Feng, J. and Saenko, K., 2016. Return of Frustratingly Easy
+    Domain Adaptation. In *Proceedings of the AAAI Conference on Artificial
+    Intelligence*.
+    """
+
+    _parameter_constraints: dict = {
+        "lambda_": [Interval(Real, 0, None, closed="left")],
+    }
+
+    def __init__(self, lambda_=1e-5):
+        self.lambda_ = lambda_
+
+    @_fit_context(prefer_skip_nested_validation=True)
+    def fit(self, X, y=None, covariates=None, target_covariate=None):
+        """Fit CORAL 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. CORAL 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 : CORAL
+            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 CORAL 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 covariance.")
+        if X_target.shape[0] < 2:
+            raise ValueError("At least two target samples are required to estimate the target covariance.")
+
+        self.source_mean_ = X_source.mean(axis=0)
+        self.target_mean_ = X_target.mean(axis=0)
+        self.source_covariance_ = _regularized_covariance(X_source, self.lambda_)
+        self.target_covariance_ = _regularized_covariance(X_target, self.lambda_)
+        self.alignment_ = _symmetric_inv_sqrt(self.source_covariance_) @ _symmetric_sqrt(self.target_covariance_)
+        self.target_covariate_ = split.target_covariate
+        self.domain_values_ = unique_covariates
+        self._n_features_out = n_features
+        return self
+
+    def transform(self, X, covariates=None):
+        """Align new samples to the learned target feature space.
+
+        Parameters
+        ----------
+        X : array-like of shape (n_samples, n_features)
+            Samples to transform.
+        covariates : array-like of shape (n_samples,), default=None
+            Binary domain labels for the new samples. Source rows are
+            whitened with the source covariance and recolored with the target
+            covariance; target rows are centered by the fitted target mean.
+            When ``None``, all rows are treated as source samples.
+
+        Returns
+        -------
+        X_new : ndarray of shape (n_samples, n_features)
+            Transformed samples. Feature dimensionality is preserved.
+        """
+        check_is_fitted(self, "alignment_")
+        X = validate_data(self, X, dtype=[np.float64, np.float32], reset=False)
+
+        if covariates is None:
+            return (X - self.source_mean_) @ self.alignment_
+
+        covariates = np.asarray(covariates)
+        if covariates.ndim == 2 and covariates.shape[1] == 1:
+            covariates = covariates.reshape(-1)
+        if covariates.ndim != 1:
+            raise ValueError(f"Covariates for {self.__class__.__name__} must be a 1D array of domain labels.")
+        if covariates.shape[0] != X.shape[0]:
+            raise ValueError("Covariates and X must have the same number of samples.")
+        if not (np.issubdtype(covariates.dtype, np.number) or np.issubdtype(covariates.dtype, np.bool_)):
+            raise ValueError(f"Covariates for {self.__class__.__name__} should be numeric or boolean domain labels.")
+
+        unknown_values = np.setdiff1d(np.unique(covariates), self.domain_values_)
+        if unknown_values.size:
+            raise ValueError(
+                f"Covariates for {self.__class__.__name__} contain domain values "
+                f"not seen at fit time: {unknown_values.tolist()}."
+            )
+
+        source_mask = covariates != self.target_covariate_
+        X_new = (X - self.target_mean_).copy()
+        if np.any(source_mask):
+            X_new[source_mask] = (X[source_mask] - self.source_mean_) @ self.alignment_
+        return X_new
+
+    def fit_transform(self, X, y=None, covariates=None, target_covariate=None):
+        """Fit CORAL and transform the source and target samples.
+
+        Parameters are the same as for :meth:`fit`. Domain labels are
+        forwarded to :meth:`transform`, so source samples are aligned and
+        target samples are centered in the returned array.
+        """
+        self.fit(X, y=y, covariates=covariates, target_covariate=target_covariate)
+        return self.transform(X, covariates=covariates)
diff --git a/tests/test_public_api.py b/tests/test_public_api.py
index 708c32a..fbe8d0c 100644
--- a/tests/test_public_api.py
+++ b/tests/test_public_api.py
@@ -8,6 +8,7 @@ def test_embed_module_exposes_transformers():
     assert embed.TCA is transformer.TCA
     assert embed.JDA is transformer.JDA
     assert embed.BDA is transformer.BDA
+    assert embed.CORAL is transformer.CORAL
     assert embed.MIDA is transformer.MIDA
     assert embed.MPCA is transformer.MPCA
     assert embed.CIFE is transformer.CIFE
diff --git a/tests/transformer/test_coral.py b/tests/transformer/test_coral.py
new file mode 100644
index 0000000..8bab28d
--- /dev/null
+++ b/tests/transformer/test_coral.py
@@ -0,0 +1,149 @@
+import numpy as np
+import pytest
+from numpy import testing
+from sklearn.exceptions import NotFittedError
+
+from kalelinear.transformer import CORAL
+
+
+@pytest.fixture(scope="module")
+def coral_data():
+    rng = np.random.default_rng(0)
+    n_source, n_target, n_features = 80, 120, 3
+
+    X_source = rng.normal(size=(n_source, n_features)) * np.array([0.5, 1.0, 2.0])
+    X_target = rng.normal(size=(n_target, n_features)) * np.array([2.0, 1.0, 0.5])
+    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 test_coral_aligns_source_covariance_to_target(coral_data):
+    X, domains, X_source, X_target = coral_data
+    coral = CORAL(lambda_=1e-6)
+
+    z = coral.fit_transform(X, covariates=domains, target_covariate=1)
+    z_source = z[domains == 0]
+    z_target = z[domains == 1]
+
+    covariance_before = np.linalg.norm(np.cov(X_source, rowvar=False) - np.cov(X_target, rowvar=False))
+    covariance_after = np.linalg.norm(np.cov(z_source, rowvar=False) - np.cov(z_target, rowvar=False))
+
+    assert z_source.shape == (X_source.shape[0], X.shape[1])
+    assert z_target.shape == (X_target.shape[0], X.shape[1])
+    assert covariance_after < covariance_before / 100
+    assert np.isfinite(z).all()
+
+
+def test_coral_centers_and_aligns_rows_independently(coral_data):
+    X, domains, X_source, X_target = coral_data
+    coral = CORAL(lambda_=1e-6).fit(X, covariates=domains, target_covariate=1)
+
+    expected_source = (X_source - coral.source_mean_) @ coral.alignment_
+    expected_target = X_target - coral.target_mean_
+
+    testing.assert_allclose(
+        coral.transform(X_source),
+        expected_source,
+    )
+    testing.assert_allclose(
+        coral.transform(X_target, covariates=np.ones(X_target.shape[0], dtype=int)),
+        expected_target,
+    )
+
+    # Row-wise domain labels must reproduce the same rows as `fit_transform`.
+    permutation = np.random.default_rng(1).permutation(X.shape[0])
+    z_permuted = coral.transform(X[permutation], covariates=domains[permutation])
+    z_expected = coral.fit_transform(X, covariates=domains, target_covariate=1)
+    testing.assert_allclose(z_permuted, z_expected[permutation])
+
+
+def test_coral_fit_transform_default_lambda(coral_data):
+    X, domains, _, _ = coral_data
+    coral = CORAL()
+
+    z = coral.fit_transform(X, covariates=domains, target_covariate=1)
+
+    assert z.shape == X.shape
+    assert np.isfinite(z).all()
+    assert coral.source_covariance_.shape == (X.shape[1], X.shape[1])
+    assert coral.target_covariance_.shape == (X.shape[1], X.shape[1])
+    assert coral.alignment_.shape == (X.shape[1], X.shape[1])
+    assert coral._n_features_out == X.shape[1]
+    testing.assert_array_equal(
+        coral.get_feature_names_out(),
+        np.array([f"coral{i}" for i in range(X.shape[1])]),
+    )
+
+
+def test_coral_transform_requires_fit(coral_data):
+    X, _, _, _ = coral_data
+    with pytest.raises(NotFittedError, match="not fitted"):
+        CORAL().transform(X)
+
+
+def test_coral_fit_requires_covariates(coral_data):
+    X, _, _, _ = coral_data
+    with pytest.raises(ValueError, match="Covariates must be provided"):
+        CORAL().fit(X)
+
+
+def test_coral_fit_requires_both_domains(coral_data):
+    X, _, X_source, _ = coral_data
+    with pytest.raises(ValueError, match="both source and target"):
+        CORAL().fit(X, covariates=np.zeros(X.shape[0], dtype=int))
+    with pytest.raises(ValueError, match="both source and target"):
+        CORAL().fit(X_source, covariates=np.zeros(X_source.shape[0], dtype=int))
+
+
+def test_coral_requires_two_samples_per_domain(coral_data):
+    X, _, _, _ = coral_data
+    domains = np.array([0, 0, 1])
+    with pytest.raises(ValueError, match="At least two target samples"):
+        CORAL().fit(X[:3], covariates=domains, target_covariate=1)
+
+
+def test_coral_validates_target_covariate(coral_data):
+    X, domains, _, _ = coral_data
+    with pytest.raises(ValueError, match="target_covariate"):
+        CORAL().fit(X, covariates=domains, target_covariate=2)
+
+
+def test_coral_validates_lambda(coral_data):
+    X, domains, _, _ = coral_data
+    with pytest.raises(ValueError, match="lambda_"):
+        CORAL(lambda_=-1).fit(X, covariates=domains)
+
+
+def test_coral_transform_validates_inputs(coral_data):
+    X, domains, X_source, X_target = coral_data
+    coral = CORAL().fit(X, covariates=domains, target_covariate=1)
+
+    with pytest.raises(ValueError, match="features"):
+        coral.transform(X[:, :2])
+
+    with pytest.raises(ValueError, match="numeric or boolean"):
+        coral.transform(X, covariates=domains.astype(str))
+
+    with pytest.raises(ValueError, match="not seen at fit time"):
+        coral.transform(X, covariates=np.full(X.shape[0], 2, dtype=int))
+
+    with pytest.raises(ValueError, match="same number of samples"):
+        coral.transform(X, covariates=np.zeros(X_target.shape[0] + 1, dtype=int))
+
+    with pytest.raises(ValueError, match="1D array"):
+        coral.transform(X_source, covariates=np.zeros((X_source.shape[0], 1, 1)))
+
+
+def test_coral_accepts_boolean_domain_labels(coral_data):
+    X, _, _, _ = coral_data
+    domains = np.concatenate((np.zeros(X.shape[0] - 2, dtype=bool), np.ones(2, dtype=bool)))
+    # Keep at least two samples per domain.
+    domains[:2] = False
+
+    coral = CORAL(lambda_=1e-6)
+    z = coral.fit_transform(X, covariates=domains, target_covariate=True)
+
+    assert coral.target_covariate_ is True
+    assert z.shape == X.shape
+    assert np.isfinite(z).all()