diff --git a/README.md b/README.md index 2c50427..f87a65c 100644 --- a/README.md +++ b/README.md @@ -9,28 +9,35 @@ [![PyPI version](https://img.shields.io/pypi/v/kalelinear?color=blue)](https://pypi.org/project/kalelinear/) [![PyPI downloads](https://pepy.tech/badge/kalelinear)](https://pepy.tech/project/kalelinear) -`kalelinear` is a Python library for learning harmonized or individualized models from multi-source/multi-view data in linear or reproducing kernel Hilbert spaces (RKHS). It provides NumPy-based methods for leveraging related data distributions and structural assumptions, including transfer learning, domain adaptation, manifold regularization, and group-aware learning, through a [`scikit-learn`](https://github.com/scikit-learn/scikit-learn) style API. +KaleLinear is a Python library for non-deep machine learning that learns transferable, shared, or group-specific models from data across multiple sources, groups, blocks, or views. It provides NumPy-based methods in linear or reproducing kernel Hilbert spaces (RKHS), including transfer learning, domain adaptation, manifold regularization, and group-aware learning, through a [`scikit-learn`](https://github.com/scikit-learn/scikit-learn) style API. The package is part of the [PyKale](https://github.com/pykale/pykale) ecosystem and focuses on linear and kernel methods for data characterized by covariates (e.g., domain labels, group labels, side information), unlabeled target samples, or tensor structures. -## What's included - -- Transformer models for learning feature embeddings: - - Multilinear Principal Component Analysis (MPCA) [[1](#references)] - - Transfer Component Analysis (TCA) [[2](#references)] - - Joint Distribution Adaptation (JDA) [[3](#references)] - - Balanced Distribution Adaptation (BDA) [[4](#references)] - - Maximum Independence Domain Adaptation (MIDA) [[5](#references)] -- Estimator models for classification: - - Manifold Regularization Learning Framework (LapSVM, LapRLS) [[6](#references)] - - Adaptation Regularization Learning Framework (ARSVM, ARRLS) [[7](#references)] - - Covariate Independence Regularized Learning Framework (CoIRSVM, CoIRLS) [[8](#references)][[9](#references)] - - Group-specific Discriminant Analysis (GSDA) [[9](#references)][[10](#references)] -- NumPy-compatible inputs and outputs. +## Key features + +- Feature transformation models for data embedding via `kalelinear.transformer` (PyKale-style alias: `kalelinear.embed`): + - Dimension reduction for multiview tensor data: + - Multilinear Principal Component Analysis (`MPCA`) [[1](#references)] + - Transferable / generalizable feature extraction across domains or groups: + - Transfer Component Analysis (`TCA`) [[2](#references)] + - Joint Distribution Adaptation (`JDA`) [[3](#references)] + - Balanced Distribution Adaptation (`BDA`) [[4](#references)] + - Maximum Independence Domain Adaptation (`MIDA`) [[5](#references)] + - Common (or shared or joint) and individual feature separation / extraction across groups or blocks: + - Common and Individual Feature Extraction (`CIFE`) [[11](#references)] + - Angle-based Joint and Individual Variation Explained (`AJIVE`) [[12](#references)] +- Estimator models for prediction via `kalelinear.estimator` (PyKale-style alias: `kalelinear.predict`): + - Predictive models that generalize across domains or groups: + - Manifold Regularization Learning Framework (`LapSVM`, `LapRLS`) [[6](#references)] + - Adaptation Regularization Learning Framework (`ARSVM`, `ARRLS`) [[7](#references)] + - Covariate Independence Regularized Learning Framework (`CoIRSVM`, `CoIRLS`) [[8](#references)][[9](#references)] + - Group-specific predictive models: + - Group-specific Discriminant Analysis (`GSDA`) [[9](#references)][[10](#references)] +- Lightweight: plain NumPy array inputs and outputs — no deep-learning framework or GPU required. - scikit-learn style `fit`, `transform`, `predict`, `fit_transform`, and `fit_predict` workflows where applicable. -- Optional covariate encoding for categorical domain or group labels. +- 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: +KaleLinear requires Python 3.10 or later. Core dependencies include: - [NumPy](http://www.numpy.org/) - [SciPy](https://www.scipy.org/) @@ -99,6 +106,7 @@ Worked examples for the main transformers and estimators are collected in - Learn a domain-invariant embedding with TCA - Use MIDA with categorical covariates +- Extract common and individual features across groups with CIFE or AJIVE - Train a domain adaptation classifier (ARSVM, ARRLS) - Train a manifold-regularized classifier (LapSVM, LapRLS) @@ -124,6 +132,10 @@ Worked examples for the main transformers and estimators are collected in [10] Zhou, S., Luo, J., Jiang, Y., Wang, H., Lu, H. and Gong, G., 2025. [Group-specific discriminant analysis enhances detection of sex differences in brain functional network lateralization](https://academic.oup.com/gigascience/article/doi/10.1093/gigascience/giaf082/8244707). _GigaScience_, 14, p.giaf082. +[11] Zhou, G., Cichocki, A., Zhang, Y. and Mandic, D., 2016. [Group component analysis for multiblock data: Common and individual feature extraction](https://ieeexplore.ieee.org/abstract/document/7310871). _IEEE Transactions on Neural Networks and Learning Systems_, 27(11), pp.2426-2439. + +[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. + ## Other open domain adaptation repositories - [POT: Python Optimal Transport](https://github.com/rflamary/POT) @@ -136,4 +148,4 @@ Worked examples for the main transformers and estimators are collected in ## License -`kalelinear` is released under the MIT License. See [LICENSE](LICENSE) for details. +KaleLinear is released under the MIT License. See [LICENSE](LICENSE) for details. diff --git a/TUTORIALS.md b/TUTORIALS.md index 320e125..075561a 100644 --- a/TUTORIALS.md +++ b/TUTORIALS.md @@ -101,3 +101,30 @@ clf = LapSVM(kernel="linear") clf.fit(X_train, ys) y_pred = clf.predict(X_target) ``` + +## Learn Common and Individual Features Across Groups + +CIFE and AJIVE decompose multiblock data into a common feature subspace shared +by all blocks plus block-specific individual subspaces. Blocks can be passed +either as a stacked matrix with a `groups` array of block ids or as a list of +block matrices sharing the same feature space. + +```python +import numpy as np +from kalelinear.transformer import CIFE, AJIVE + +rng = np.random.default_rng(0) +X = np.vstack([rng.normal(size=(50, 10)) for _ in range(3)]) +groups = np.repeat([0, 1, 2], 50) + +cife = CIFE(random_state=0) +z_common = cife.fit_transform(X, groups=groups) +z_individual = cife.transform_individual(X, groups=groups) + +ajive = AJIVE(n_resamples=50, random_state=0) +z_common_ajive = ajive.fit_transform(X, groups=groups) +``` + +`transform` returns the common feature scores of the samples, and +`transform_individual` returns a list with the block-specific individual +scores. diff --git a/docs/images/kalelinear.jpg b/docs/images/kalelinear.jpg index ffb3519..9f76423 100644 Binary files a/docs/images/kalelinear.jpg and b/docs/images/kalelinear.jpg differ diff --git a/docs/source/api_embed.rst b/docs/source/api_embed.rst index a3dd0dd..db8f038 100644 --- a/docs/source/api_embed.rst +++ b/docs/source/api_embed.rst @@ -35,3 +35,15 @@ the same transformer classes as :mod:`kalelinear.transformer`. :undoc-members: :show-inheritance: :no-index: + +.. autoclass:: CIFE + :members: + :undoc-members: + :show-inheritance: + :no-index: + +.. autoclass:: AJIVE + :members: + :undoc-members: + :show-inheritance: + :no-index: diff --git a/docs/source/api_transformers.rst b/docs/source/api_transformers.rst index b7db20e..f73863b 100644 --- a/docs/source/api_transformers.rst +++ b/docs/source/api_transformers.rst @@ -27,3 +27,13 @@ Transformers :members: :undoc-members: :show-inheritance: + +.. autoclass:: CIFE + :members: + :undoc-members: + :show-inheritance: + +.. autoclass:: AJIVE + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/source/index.rst b/docs/source/index.rst index 343bc78..063967f 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -1,4 +1,4 @@ -kalelinear Documentation +KaleLinear Documentation ========================= Getting Started @@ -17,15 +17,13 @@ kalelinear API .. toctree:: :maxdepth: 2 - api_embed api_transformers - api_predict api_estimators - api_utilities -kalelinear APIs above are ordered following the machine learning pipeline, -i.e., feature embedding transformers, predictive estimators, and reusable -utilities, rather than alphabetically. +The API pages above are ordered following the machine learning pipeline, +i.e., feature embedding transformers and predictive estimators, rather than +alphabetically. `kalelinear.embed` and `kalelinear.predict` are PyKale-compatible +aliases of `kalelinear.transformer` and `kalelinear.estimator`, respectively. Project ------- diff --git a/docs/source/installation.rst b/docs/source/installation.rst index aad43fd..f16e79d 100644 --- a/docs/source/installation.rst +++ b/docs/source/installation.rst @@ -13,7 +13,7 @@ Install a local checkout for development: pip install -e ".[dev]" -Kale-Linear requires Python 3.10 or later. Core dependencies include NumPy, +KaleLinear requires Python 3.10 or later. Core dependencies include NumPy, SciPy, scikit-learn, pandas, tensorly, cvxopt, and osqp. To build the documentation locally, install the documentation requirements and diff --git a/docs/source/introduction.rst b/docs/source/introduction.rst index 1611d71..a21aa52 100644 --- a/docs/source/introduction.rst +++ b/docs/source/introduction.rst @@ -1,7 +1,7 @@ Introduction ============ -Kale-Linear is a Python library for non-deep, knowledge-aware machine learning +KaleLinear is a Python library for non-deep, knowledge-aware machine learning from multiple sources, domains, or views. It provides NumPy-based implementations of transfer learning, domain adaptation, manifold regularization, and group-aware linear learning methods with a scikit-learn @@ -14,8 +14,8 @@ covariates, side information, or unlabeled target samples. Main Features ------------- -* Transformer models for learning feature embeddings: MPCA, TCA, JDA, BDA, and - MIDA. +* Transformer models for learning feature embeddings: 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/kalelinear/__init__.py b/kalelinear/__init__.py index 73855ce..e54af09 100644 --- a/kalelinear/__init__.py +++ b/kalelinear/__init__.py @@ -1,7 +1,7 @@ """ kalelinear. -Learning harmonized or individualized models from multi-source/multi-view data in linear or reproducing kernel Hilbert spaces (RKHS). +Non-deep machine learning that learns transferable, shared, or group-specific models from data across multiple sources, groups, blocks, or views. """ from importlib import import_module diff --git a/kalelinear/embed.py b/kalelinear/embed.py index a9f59b3..3cb98b5 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 BDA, JDA, MIDA, MPCA, TCA +from kalelinear.transformer import AJIVE, BDA, CIFE, JDA, MIDA, MPCA, TCA -__all__ = ["TCA", "JDA", "BDA", "MIDA", "MPCA"] +__all__ = ["TCA", "JDA", "BDA", "MIDA", "MPCA", "CIFE", "AJIVE"] diff --git a/kalelinear/transformer/__init__.py b/kalelinear/transformer/__init__.py index 7a68902..b0c44bb 100644 --- a/kalelinear/transformer/__init__.py +++ b/kalelinear/transformer/__init__.py @@ -1,3 +1,5 @@ +from kalelinear.transformer._ajive import AJIVE +from kalelinear.transformer._cife import CIFE from kalelinear.transformer._jda import BDA, JDA from kalelinear.transformer._mida import MIDA from kalelinear.transformer._mpca import MPCA @@ -9,4 +11,6 @@ "BDA", "MIDA", "MPCA", + "CIFE", + "AJIVE", ] diff --git a/kalelinear/transformer/_ajive.py b/kalelinear/transformer/_ajive.py new file mode 100644 index 0000000..98db9c4 --- /dev/null +++ b/kalelinear/transformer/_ajive.py @@ -0,0 +1,253 @@ +# ============================================================================= +# @author: Shuo Zhou, Lalu Muhammad Riza Rizky, The University of Sheffield +# @contact: shuo.zhou@sheffield.ac.uk +# ============================================================================= +"""Angle-based Joint and Individual Variation Explained (AJIVE) implementation. + +References +---------- +Feng, Q., Jiang, M., Hannig, J. and Marron, J.S., 2018. Angle-based joint and +individual variation explained. Journal of Multivariate Analysis, 166, +pp.241-265. + +The algorithm follows the authors' reference MATLAB implementation +(MeileiJiang/AJIVE_Project) with the same perturbation-bound and rank +selection steps. +""" + +from numbers import Integral, Real + +import numpy as np +from sklearn.utils._param_validation import Interval + +from kalelinear.transformer._multiblock import _check_per_block_ranks, BaseCommonIndividualTransformer + +_FERROR = 1e-10 + + +def _jive_rand_null_norm(data, basis, n_sim, random_state): + """Energy of data on random directions orthogonal to ``basis``.""" + n_ambient_dims = basis.shape[0] + if basis.shape[1] >= n_ambient_dims: + return np.zeros(n_sim) + null_norms = np.empty(n_sim) + for i in range(n_sim): + current = basis.copy() + directions = [] + for _ in range(basis.shape[1]): + direction = random_state.randn(n_ambient_dims) + direction = direction - current @ (current.T @ direction) + norm = np.linalg.norm(direction) + while norm == 0: + direction = random_state.randn(n_ambient_dims) + direction = direction - current @ (current.T @ direction) + norm = np.linalg.norm(direction) + direction /= norm + directions.append(direction) + current = np.column_stack((current, direction)) + directions = np.column_stack(directions) + null_norms[i] = np.linalg.norm(data @ directions) + return null_norms + + +def _wedin_angle_bound(block, n_sim, U, S, V, random_state): + """Resampled Wedin perturbation-angle bound for one data block.""" + row_bound = _jive_rand_null_norm(block, V, n_sim, random_state) + column_bound = _jive_rand_null_norm(block.T, U, n_sim, random_state) + delta = S[-1] + ratio = np.maximum(row_bound, column_bound) / delta + ratio = np.clip(ratio, 0.0, 1.0) + return np.rad2deg(np.arcsin(ratio)) + + +def _random_direction_ssv(n, ranks, n_sim, random_state): + """Largest squared singular values of random stacked subspaces.""" + stacked = np.zeros((int(np.sum(ranks)), n)) + values = np.empty(n_sim) + for i in range(n_sim): + row = 0 + for rank in ranks: + Q, _ = np.linalg.qr(random_state.randn(n, rank)) + stacked[row : row + rank] = Q.T + row += rank + values[i] = np.linalg.norm(stacked, 2) ** 2 + return values + + +class AJIVE(BaseCommonIndividualTransformer): + """Angle-based Joint and Individual Variation Explained (AJIVE). + + AJIVE decomposes multiblock data into a joint (common) subspace shared by + all blocks and block-specific individual subspaces. It estimates the joint + rank from the Wedin perturbation bound of the stacked per-block row spaces + and reconstructs the joint and individual structures with low-rank SVDs. + + ``X`` can be a single matrix with a ``groups`` array of block ids, or a + list of block matrices that all share the same feature space. + :meth:`transform` returns the common feature scores of the samples while + :meth:`transform_individual` returns the block-specific individual scores. + + Parameters + ---------- + n_common_components : int or None, default=None + Number of joint components. When None, the joint rank is selected from + the Wedin perturbation bound. + n_individual_components : int, array-like or None, default=None + Number of individual components per block. When None, components above + the per-block singular-value threshold are kept. + initial_ranks : int, array-like or None, default=None + Initial signal rank of each block. When None, ranks are estimated from + the fraction of variance explained by ``variance_threshold``. + variance_threshold : float, default=0.95 + Cumulative explained-variance fraction used to estimate ``initial_ranks`` + when it is None. + n_resamples : int, default=1000 + Number of re-samples for the Wedin perturbation bound. + percentile : float, default=5 + Percentile of the Wedin bounds used for the joint rank selection. + random_state : int, RandomState or None, default=None + Random seed for the perturbation-bound re-sampling. + + Attributes + ---------- + common_components_ : ndarray of shape (n_features, n_common_components_) + Orthonormal joint feature-space basis. + individual_components_ : list of ndarray of shape (n_features, rank) + Per-block orthonormal individual feature-space bases. + common_scores_ : list of ndarray of shape (n_samples_in_block, n_common_components_) + Common scores of the training samples in each block. + individual_scores_ : list of ndarray of shape (n_samples_in_block, rank) + Individual scores of the training samples in each block. + individual_ranks_ : ndarray of shape (n_blocks,) + Number of individual components retained per block. + n_common_components_ : int + Number of joint components retained. + """ + + _parameter_constraints: dict = { + **BaseCommonIndividualTransformer._parameter_constraints, + "initial_ranks": ["array-like", Interval(Integral, 1, None, closed="left"), None], + "variance_threshold": [Interval(Real, 0, 1, closed="both")], + "n_resamples": [Interval(Integral, 1, None, closed="left")], + "percentile": [Interval(Real, 0, 100, closed="both")], + } + + def __init__( + self, + n_common_components=None, + n_individual_components=None, + initial_ranks=None, + variance_threshold=0.95, + n_resamples=1000, + percentile=5, + random_state=None, + ): + self.initial_ranks = initial_ranks + self.variance_threshold = variance_threshold + self.n_resamples = n_resamples + self.percentile = percentile + super().__init__( + n_common_components=n_common_components, + n_individual_components=n_individual_components, + random_state=random_state, + ) + + def _resolve_initial_ranks(self, blocks): + if self.initial_ranks is None: + ranks = [] + for block in blocks: + singular_values = np.linalg.svd(block, compute_uv=False) + if singular_values.size == 0: + ranks.append(0) + continue + explained = np.cumsum(singular_values**2) / np.sum(singular_values**2) + rank = int(np.searchsorted(explained, self.variance_threshold) + 1) + ranks.append(int(min(rank, singular_values.size))) + return np.asarray(ranks, dtype=int) + ranks = _check_per_block_ranks(self.initial_ranks, self.n_blocks_, "initial_ranks") + for n, (rank, block) in enumerate(zip(ranks, blocks)): + if rank < 1 or rank > min(block.shape): + raise ValueError( + f"`initial_ranks[{n}]` must be between 1 and min(samples, features) " + f"= {min(block.shape)}, got {rank}." + ) + return ranks + + def _fit_blocks(self, blocks): + D = self.n_features_in_ + ranks = self._resolve_initial_ranks(blocks) + if np.any(ranks < 1): + raise ValueError("`initial_ranks` must contain positive values for every block.") + + stacked = np.zeros((int(np.sum(ranks)), D)) + thresholds = np.empty(self.n_blocks_) + angle_bounds = [] + row = 0 + for n, block in enumerate(blocks): + rank = ranks[n] + U, s, Vt = np.linalg.svd(block, full_matrices=False) + if rank + 1 <= len(s): + thresholds[n] = 0.5 * (s[rank - 1] + s[rank]) + else: + thresholds[n] = 0.5 * s[-1] + U0, S0, V0 = U[:, :rank], s[:rank], Vt[:rank].T + stacked[row : row + rank] = V0.T + row += rank + angle_bounds.append(_wedin_angle_bound(block, self.n_resamples, U0, S0, V0, self.random_state_)) + angle_bounds = np.vstack(angle_bounds) + + _, s_stacked, Vt_stacked = np.linalg.svd(stacked, full_matrices=False) + s_stacked = s_stacked[: min(stacked.shape)] + wedin_ssv_bounds = np.maximum(np.sum(np.cos(np.deg2rad(angle_bounds)) ** 2, axis=0), 1.0) + wedin_ssv_bound = np.percentile(wedin_ssv_bounds, self.percentile) + random_ssvs = _random_direction_ssv(D, ranks, 100, self.random_state_) + random_ssv_bound = np.percentile(random_ssvs, 95) + # Take the more conservative (larger) of the two perturbation bounds, + # following the reference implementation: max(wedin, random). + joint_threshold = max(wedin_ssv_bound, random_ssv_bound) + joint_rank = int(np.sum(s_stacked**2 + _FERROR > joint_threshold)) + if self.n_common_components is not None: + joint_rank = min(int(self.n_common_components), len(s_stacked)) + + row_joint = Vt_stacked[:joint_rank] + drop_rows = set() + for n, block in enumerate(blocks): + projected = block @ row_joint.T + low_variance = np.flatnonzero(np.sqrt(np.sum(projected**2, axis=0)) <= thresholds[n] + _FERROR) + drop_rows.update(low_variance.tolist()) + if drop_rows: + keep_rows = [j for j in range(row_joint.shape[0]) if j not in drop_rows] + row_joint = row_joint[keep_rows] + joint_rank = row_joint.shape[0] + + common_components = row_joint.T + common_scores = [block @ common_components for block in blocks] + + ranks_spec = _check_per_block_ranks(self.n_individual_components, self.n_blocks_, "n_individual_components") + individual_components = [] + individual_scores = [] + individual_ranks = [] + for n, block in enumerate(blocks): + individual = block - block @ common_components @ common_components.T + s_individual = np.linalg.svd(individual, compute_uv=False) + if ranks_spec is None: + rank = int(np.sum(s_individual + _FERROR > thresholds[n])) + else: + rank = min(int(ranks_spec[n]), len(s_individual)) + rank = max(rank, 0) + if rank > 0: + U_i, S_i, Vt_i = np.linalg.svd(individual, full_matrices=False) + U_i, S_i, Vt_i = U_i[:, :rank], S_i[:rank], Vt_i[:rank] + individual_components.append(Vt_i.T) + individual_scores.append(U_i * S_i) + else: + individual_components.append(np.zeros((D, 0))) + individual_scores.append(np.zeros((block.shape[0], 0))) + individual_ranks.append(rank) + + self.common_components_ = common_components + self.common_scores_ = common_scores + self.individual_components_ = individual_components + self.individual_scores_ = individual_scores + self.individual_ranks_ = np.asarray(individual_ranks, dtype=int) + self.n_common_components_ = joint_rank diff --git a/kalelinear/transformer/_cife.py b/kalelinear/transformer/_cife.py new file mode 100644 index 0000000..b2236fd --- /dev/null +++ b/kalelinear/transformer/_cife.py @@ -0,0 +1,287 @@ +# ============================================================================= +# @author: Shuo Zhou, The University of Sheffield +# @contact: shuo.zhou@sheffield.ac.uk +# ============================================================================= +"""Common and Individual Feature Extraction (CIFE) implementation. + +References +---------- +Zhou, G., Cichocki, A., Zhang, Y. and Mandic, D., 2016. Group component +analysis for multiblock data: Common and individual feature extraction. +IEEE Transactions on Neural Networks and Learning Systems, 27(11), +pp.2426-2439. + +The common orthogonal basis extraction (COBE) steps follow the authors' +reference implementations (pycifa and the accompanying MATLAB code). +""" + +from numbers import Integral, Real + +import numpy as np +from sklearn.utils._param_validation import Interval + +from kalelinear.transformer._multiblock import _check_per_block_ranks, BaseCommonIndividualTransformer + + +def _column_space_basis(Y, pca_dim=None): + """Return an orthonormal basis of the column space of ``Y`` (D x J).""" + D, J = Y.shape + U, s, _ = np.linalg.svd(Y, full_matrices=False) + if s.size == 0 or s[0] == 0: + raise ValueError("Each block must have a non-zero column space.") + rank = int(np.sum(s > s[0] * max(D, J) * np.finfo(s.dtype).eps)) + if rank >= D: + if pca_dim is None: + raise ValueError( + "A block spans the whole feature space, so common and individual " + "subspaces cannot be separated. Reduce the dimensionality first " + "or set `pca_dim` to truncate the per-block column spaces." + ) + if 0 < pca_dim < 1: + rank = min(int(np.floor(D * pca_dim)), D - 1) + else: + rank = min(int(pca_dim), D - 1) + rank = max(rank, 1) + return U[:, :rank], rank + + +def _cobe_common_basis(blocks, c, max_iter, tol, epsilon, pca_dim, random_state): + """Extract a common orthogonal basis shared by all blocks. + + Parameters + ---------- + blocks : list of ndarray of shape (D, J_n) + Data blocks sharing the same feature dimension ``D``. + c : int or None + Number of common components. When None, the number is estimated from + the residual ``epsilon`` criterion. + max_iter : int + Maximum power iterations per common direction. + tol : float + Convergence tolerance for the power iterations. + epsilon : float + Residual threshold below which a direction counts as common. + pca_dim : int, float or None + Optional per-block dimensionality truncation for blocks spanning the + whole feature space. + random_state : RandomState + Random number generator for initializing the power iterations. + + Returns + ------- + common_basis : ndarray of shape (D, n_common) + Orthonormal common basis shared by all blocks. + """ + n_blocks = len(blocks) + D = blocks[0].shape[0] + bases = [] + ranks = [] + for Y in blocks: + basis, rank = _column_space_basis(Y, pca_dim=pca_dim) + bases.append(basis) + ranks.append(rank) + min_rank = min(ranks) + if min_rank == 0: + return np.zeros((D, 0)) + if c is not None and c <= 0: + return np.zeros((D, 0)) + + order = np.argsort(ranks) + projections = [np.zeros((bases[n].shape[1], min_rank)) for n in range(n_blocks)] + common_basis = np.zeros((D, min_rank)) + residuals = [] + + def _power_iteration(initial, column): + direction = initial / np.linalg.norm(initial) + for _ in range(max_iter): + previous = direction + update = np.zeros(D) + for n in range(n_blocks): + projections[n][:, column] = bases[n].T @ direction + update += bases[n] @ projections[n][:, column] + update_norm = np.linalg.norm(update) + if update_norm == 0: + break + direction = update / update_norm + if abs(previous @ direction) > 1 - tol: + break + return direction + + # Seek the first common direction. + initial = bases[order[0]] @ random_state.randn(bases[order[0]].shape[1]) + first = _power_iteration(initial, 0) + residual = 0.0 + for n in range(n_blocks): + projection = bases[n].T @ first + residual += 1 - projection @ projection + residual /= n_blocks + residuals.append(residual) + + if c is None and residual > epsilon: + return np.zeros((D, 0)) + + if c is not None: + c = min(c, min_rank) + common_basis = np.zeros((D, c)) + common_basis[:, 0] = first + residuals.extend([np.inf] * (c - 1)) + else: + common_basis[:, 0] = first + residuals.extend([np.inf] * (min_rank - 1)) + + # Seek the remaining common directions with deflation. + for j in range(1, min_rank): + if c is not None and j >= c: + break + for n in range(n_blocks): + basis = bases[n] + bases[n] = basis - np.outer(basis @ projections[n][:, j - 1], projections[n][:, j - 1]) + initial = bases[order[0]] @ random_state.randn(bases[order[0]].shape[1]) + direction = _power_iteration(initial, j) + residual = 0.0 + for n in range(n_blocks): + projection = bases[n].T @ direction + residual += 1 - projection @ projection + residual /= n_blocks + residuals[j] = residual + if c is None and residual > epsilon: + residuals[j] = np.inf + break + common_basis[:, j] = direction + + common_basis = common_basis[:, ~np.isinf(np.asarray(residuals))] + if common_basis.shape[1] > 0: + u, _, vt = np.linalg.svd(common_basis, full_matrices=False) + common_basis = u @ vt + return common_basis + + +class CIFE(BaseCommonIndividualTransformer): + """Common and Individual Feature Extraction (CIFE). + + CIFE decomposes multiblock data into a common feature subspace shared by + all blocks and block-specific individual subspaces. The common subspace is + extracted with the common orthogonal basis extraction (COBE) algorithm and + the individual subspaces are obtained from the residual of each block after + removing its common part. + + ``X`` can be a single matrix with a ``groups`` array of block ids, or a + list of block matrices that all share the same feature space. + :meth:`transform` returns the common feature scores of the samples while + :meth:`transform_individual` returns the block-specific individual scores. + + Parameters + ---------- + n_common_components : int or None, default=None + Number of common components to extract. When None, the number is + estimated automatically from the residual threshold ``epsilon``. + n_individual_components : int, array-like or None, default=None + Number of individual components per block. When None, all numerically + non-zero residual directions are kept. + max_iter : int, default=200 + Maximum power iterations for each common direction. + tol : float, default=1e-6 + Convergence tolerance for the power iterations. + epsilon : float, default=0.01 + Residual threshold used to decide whether a direction is common when + ``n_common_components`` is None. + pca_dim : int, float or None, default=None + Optional truncation of per-block column spaces, either as a relative + fraction in (0, 1) or an absolute number of components. Required when + a block spans the whole feature space. + random_state : int, RandomState or None, default=None + Random seed for initializing the power iterations. + + Attributes + ---------- + common_components_ : ndarray of shape (n_features, n_common_components_) + Orthonormal common feature-space basis. + individual_components_ : list of ndarray of shape (n_features, rank) + Per-block orthonormal individual feature-space bases. + common_scores_ : list of ndarray of shape (n_samples_in_block, n_common_components_) + Common scores of the training samples in each block. + individual_scores_ : list of ndarray of shape (n_samples_in_block, rank) + Individual scores of the training samples in each block. + individual_ranks_ : ndarray of shape (n_blocks,) + Number of individual components retained per block. + n_common_components_ : int + Number of common components retained. + """ + + _parameter_constraints: dict = { + **BaseCommonIndividualTransformer._parameter_constraints, + "max_iter": [Interval(Integral, 1, None, closed="left")], + "tol": [Interval(Real, 0, None, closed="left")], + "epsilon": [Interval(Real, 0, None, closed="left")], + "pca_dim": [ + Interval(Real, 0, 1, closed="right"), + Interval(Integral, 2, None, closed="left"), + None, + ], + } + + def __init__( + self, + n_common_components=None, + n_individual_components=None, + max_iter=200, + tol=1e-6, + epsilon=0.01, + pca_dim=None, + random_state=None, + ): + self.max_iter = max_iter + self.tol = tol + self.epsilon = epsilon + self.pca_dim = pca_dim + super().__init__( + n_common_components=n_common_components, + n_individual_components=n_individual_components, + random_state=random_state, + ) + + def _fit_blocks(self, blocks): + D = self.n_features_in_ + transposed = [block.T for block in blocks] + common_components = _cobe_common_basis( + transposed, + c=self.n_common_components, + max_iter=self.max_iter, + tol=self.tol, + epsilon=self.epsilon, + pca_dim=self.pca_dim, + random_state=self.random_state_, + ) + n_common = common_components.shape[1] + common_scores = [block @ common_components for block in blocks] + + ranks_spec = _check_per_block_ranks(self.n_individual_components, self.n_blocks_, "n_individual_components") + individual_components = [] + individual_scores = [] + individual_ranks = [] + for n, (X_n, Y_n) in enumerate(zip(blocks, transposed)): + residual = Y_n - common_components @ (common_components.T @ Y_n) + U_i, s_i, _ = np.linalg.svd(residual, full_matrices=False) + if ranks_spec is None: + if s_i.size == 0: + rank = 0 + else: + noise_tol = s_i[0] * max(residual.shape) * np.finfo(s_i.dtype).eps + rank = int(np.sum(s_i > noise_tol)) + else: + rank = min(int(ranks_spec[n]), U_i.shape[1]) + if rank > 0: + U_i = U_i[:, :rank] + individual_components.append(U_i) + individual_scores.append(X_n @ U_i) + else: + individual_components.append(np.zeros((D, 0))) + individual_scores.append(np.zeros((X_n.shape[0], 0))) + individual_ranks.append(rank) + + self.common_components_ = common_components + self.common_scores_ = common_scores + self.individual_components_ = individual_components + self.individual_scores_ = individual_scores + self.individual_ranks_ = np.asarray(individual_ranks, dtype=int) + self.n_common_components_ = n_common diff --git a/kalelinear/transformer/_multiblock.py b/kalelinear/transformer/_multiblock.py new file mode 100644 index 0000000..f0dacc0 --- /dev/null +++ b/kalelinear/transformer/_multiblock.py @@ -0,0 +1,201 @@ +# ============================================================================= +# @author: Shuo Zhou, Lalu Muhammad Riza Rizky, The University of Sheffield +# @contact: shuo.zhou@sheffield.ac.uk +# ============================================================================= +"""Shared base classes for multiblock common and individual feature transformers.""" + +from abc import abstractmethod +from numbers import Integral + +import numpy as np +from sklearn.base import BaseEstimator, ClassNamePrefixFeaturesOutMixin, TransformerMixin +from sklearn.utils._param_validation import Interval +from sklearn.utils.validation import check_is_fitted, check_random_state + + +def _check_multiblock_input(X, groups=None, min_blocks=2): + """Validate a multiblock input and return one matrix per block. + + Parameters + ---------- + X : array-like of shape (n_samples, n_features) or list of array-like + When ``X`` is a single matrix, ``groups`` must give the block id of + each row. When ``X`` is a list, each element is one block and all + blocks must share the same feature space (columns). + groups : array-like of shape (n_samples,), default=None + Block id for each sample when ``X`` is a single stacked matrix. + min_blocks : int, default=2 + Minimum number of blocks required. The common projection in + :meth:`transform` only needs one block, so callers that project new + samples can pass ``min_blocks=1``. + + Returns + ------- + blocks : list of ndarray of shape (n_samples_in_block, n_features) + groups : ndarray of shape (n_samples,) or None + """ + if isinstance(X, (list, tuple)): + if groups is not None: + raise ValueError("`groups` must be None when `X` is a list of blocks.") + if len(X) == 0: + raise ValueError("`X` must contain at least one block.") + blocks = [] + n_features = None + for block in X: + block = np.asarray(block, dtype=float) + if block.ndim != 2: + raise ValueError("Each block in `X` must be a 2D array.") + if block.shape[0] == 0: + raise ValueError("Each block in `X` must contain at least one sample.") + if n_features is None: + n_features = block.shape[1] + elif block.shape[1] != n_features: + raise ValueError("All blocks must share the same number of features.") + blocks.append(block) + if len(blocks) < min_blocks: + raise ValueError("At least two blocks are required for common and individual feature extraction.") + return blocks, None + + X = np.asarray(X, dtype=float) + if X.ndim != 2: + raise ValueError("`X` must be a 2D array or a list of 2D block arrays.") + if groups is None: + raise ValueError("`groups` must be provided when `X` is a single stacked matrix.") + groups = np.asarray(groups) + if groups.ndim != 1 or groups.shape[0] != X.shape[0]: + raise ValueError("`groups` must be a 1D array aligned with the rows of `X`.") + block_ids = np.unique(groups) + blocks = [X[groups == block_id] for block_id in block_ids] + if len(blocks) < min_blocks: + raise ValueError("At least two blocks are required for common and individual feature extraction.") + if any(block.shape[0] == 0 for block in blocks): + raise ValueError("Each block must contain at least one sample.") + return blocks, groups + + +def _check_per_block_ranks(n_components, n_blocks, name): + """Validate an integer or a per-block sequence of component counts.""" + if n_components is None: + return None + if isinstance(n_components, (Integral, np.integer)): + return np.full(n_blocks, int(n_components), dtype=int) + ranks = np.asarray(n_components) + if ranks.ndim != 1 or ranks.shape[0] != n_blocks: + raise ValueError(f"{name} must be an integer or a sequence with one value per block.") + if not np.issubdtype(ranks.dtype, np.number): + raise ValueError(f"{name} must contain numeric values.") + if np.any(np.isnan(ranks)): + raise ValueError(f"{name} must not contain NaN values.") + if np.any(ranks < 0): + raise ValueError(f"{name} must contain non-negative values.") + if not np.all(np.equal(ranks, np.floor(ranks))): + raise ValueError(f"{name} must contain integer values.") + return ranks.astype(int) + + +class BaseCommonIndividualTransformer(ClassNamePrefixFeaturesOutMixin, TransformerMixin, BaseEstimator): + """Base class for multiblock common and individual feature transformers. + + Subclasses learn a common feature subspace shared by all data blocks plus + block-specific individual subspaces, following the common and individual + feature extraction (CIFE) framework of Zhou and Cichocki (2016). + + The input convention follows the rest of ``kalelinear``: ``X`` is a sample + matrix whose rows are partitioned into blocks by ``groups``, or a list of + block matrices that all share the same feature space. + """ + + _parameter_constraints: dict = { + "n_common_components": [Interval(Integral, 0, None, closed="left"), None], + "n_individual_components": ["array-like", Interval(Integral, 0, None, closed="left"), None], + "random_state": ["random_state"], + } + + def __init__(self, n_common_components=None, n_individual_components=None, random_state=None): + self.n_common_components = n_common_components + self.n_individual_components = n_individual_components + self.random_state = random_state + + def fit(self, X, y=None, groups=None, **fit_params): + """Fit the transformer on multiblock data. + + Parameters + ---------- + X : array-like of shape (n_samples, n_features) or list of array-like + Stacked samples from all blocks, or a list of block matrices + sharing the same feature space. + y : array-like of shape (n_samples,), default=None + Ignored. Present for scikit-learn API consistency. + groups : array-like of shape (n_samples,), default=None + Block id for each sample, required when ``X`` is a single stacked + matrix. + + Returns + ------- + self : object + Fitted transformer. + """ + self._validate_params() + blocks, groups = _check_multiblock_input(X, groups) + self.n_features_in_ = blocks[0].shape[1] + self.n_blocks_ = len(blocks) + self.block_sizes_ = np.array([block.shape[0] for block in blocks]) + self.random_state_ = check_random_state(self.random_state) + self._fit_blocks(blocks) + self._n_features_out = self.n_common_components_ + return self + + @abstractmethod + def _fit_blocks(self, blocks): + """Run the algorithm on validated per-block matrices.""" + + def transform(self, X, groups=None): + """Project samples onto the learned common feature subspace. + + Parameters + ---------- + X : array-like of shape (n_samples, n_features) or list of array-like + New samples, either stacked or given as a list of blocks. A + single block is accepted for the common projection. + groups : array-like of shape (n_samples,), default=None + Ignored for the common projection. Present for API consistency. + + Returns + ------- + X_new : ndarray of shape (n_samples, n_common_components_) + Common feature scores shared by all blocks. + """ + check_is_fitted(self, "common_components_") + if isinstance(X, (list, tuple)): + blocks, _ = _check_multiblock_input(X, min_blocks=1) + X_stacked = np.vstack(blocks) + else: + X_stacked = np.asarray(X, dtype=float) + if X_stacked.ndim != 2: + raise ValueError("`X` must be a 2D array or a list of 2D block arrays.") + if X_stacked.shape[1] != self.n_features_in_: + raise ValueError(f"Expected {self.n_features_in_} features, got {X_stacked.shape[1]}.") + return X_stacked @ self.common_components_ + + def transform_individual(self, X, groups=None): + """Project samples onto the individual feature subspaces. + + Parameters + ---------- + X : array-like of shape (n_samples, n_features) or list of array-like + New samples, either stacked or given as a list of blocks. + groups : array-like of shape (n_samples,), default=None + Block id for each sample when ``X`` is a single stacked matrix. + + Returns + ------- + scores : list of ndarray + One array per block, of shape (n_samples_in_block, individual_ranks_[i]). + """ + check_is_fitted(self, "individual_components_") + blocks, _ = _check_multiblock_input(X, groups) + if len(blocks) != self.n_blocks_: + raise ValueError(f"Expected {self.n_blocks_} blocks, got {len(blocks)}.") + if any(block.shape[1] != self.n_features_in_ for block in blocks): + raise ValueError(f"Expected {self.n_features_in_} features in every block, got mismatched blocks.") + return [block @ components for block, components in zip(blocks, self.individual_components_)] diff --git a/tests/test_public_api.py b/tests/test_public_api.py index b431f3b..708c32a 100644 --- a/tests/test_public_api.py +++ b/tests/test_public_api.py @@ -10,6 +10,8 @@ def test_embed_module_exposes_transformers(): assert embed.BDA is transformer.BDA assert embed.MIDA is transformer.MIDA assert embed.MPCA is transformer.MPCA + assert embed.CIFE is transformer.CIFE + assert embed.AJIVE is transformer.AJIVE def test_predict_module_exposes_estimators(): diff --git a/tests/transformer/test_ajive.py b/tests/transformer/test_ajive.py new file mode 100644 index 0000000..79943fc --- /dev/null +++ b/tests/transformer/test_ajive.py @@ -0,0 +1,123 @@ +import numpy as np +import pytest +from numpy import testing +from sklearn.base import clone + +from kalelinear.transformer import AJIVE +from tests.utils.test_utils import make_common_individual_dataset + + +@pytest.fixture(scope="module") +def multiblock_data(): + return make_common_individual_dataset(random_state=0) + + +def _subspace_error(estimated, planted): + projection = planted @ planted.T + return np.linalg.norm(estimated - projection @ estimated) / np.sqrt(estimated.shape[1]) + + +def test_ajive_recovers_joint_and_individual_subspaces(multiblock_data): + X, groups, _, common_basis = multiblock_data + ajive = AJIVE(initial_ranks=[5, 6, 4], n_resamples=50, random_state=0) + ajive.fit(X, groups=groups) + + assert ajive.n_common_components_ == 2 + assert ajive.common_components_.shape == (X.shape[1], 2) + assert _subspace_error(ajive.common_components_, common_basis) < 1e-6 + testing.assert_array_equal(ajive.individual_ranks_, [3, 4, 2]) + + +def test_ajive_automatic_initial_ranks(multiblock_data): + X, groups, _, _ = multiblock_data + ajive = AJIVE(variance_threshold=1.0, n_resamples=50, random_state=0) + ajive.fit(X, groups=groups) + assert ajive.n_common_components_ == 2 + + +def test_ajive_percentile_uses_larger_perturbation_bound(multiblock_data): + X, groups, _, _ = multiblock_data + ranks = [ + AJIVE(percentile=percentile, n_resamples=50, random_state=0).fit(X, groups=groups).n_common_components_ + for percentile in (5, 50, 95) + ] + # A larger percentile gives a non-decreasing Wedin threshold, so the joint + # rank selected from max(wedin, random-direction) bound cannot increase. + assert ranks[0] >= ranks[1] >= ranks[2] + + +def test_ajive_manual_joint_rank(multiblock_data): + X, groups, _, _ = multiblock_data + ajive = AJIVE( + n_common_components=1, + initial_ranks=[5, 6, 4], + n_resamples=50, + random_state=0, + ) + ajive.fit(X, groups=groups) + assert ajive.n_common_components_ == 1 + + +def test_ajive_individual_rank_override(multiblock_data): + X, groups, _, _ = multiblock_data + ajive = AJIVE( + initial_ranks=[5, 6, 4], + n_individual_components=[2, 2, 2], + n_resamples=50, + random_state=0, + ) + ajive.fit(X, groups=groups) + testing.assert_array_equal(ajive.individual_ranks_, [2, 2, 2]) + + +def test_ajive_list_input_matches_stacked(multiblock_data): + X, groups, blocks, _ = multiblock_data + from_stacked = AJIVE(n_resamples=50, random_state=0).fit(X, groups=groups) + from_list = AJIVE(n_resamples=50, random_state=0).fit(blocks) + testing.assert_allclose(from_stacked.common_components_, from_list.common_components_) + + +def test_ajive_transform_consistency(multiblock_data): + X, groups, blocks, _ = multiblock_data + ajive = AJIVE(n_resamples=50, random_state=0) + ajive.fit(X, groups=groups) + + z = ajive.transform(X) + testing.assert_allclose(z, X @ ajive.common_components_) + testing.assert_allclose(ajive.transform(blocks), z) + single_block = ajive.transform([blocks[0]]) + testing.assert_allclose(single_block, blocks[0] @ ajive.common_components_) + testing.assert_allclose(ajive.transform((blocks[0],)), single_block) + + individual = ajive.transform_individual(X, groups=groups) + assert [scores.shape for scores in individual] == [(60, 3), (50, 4), (70, 2)] + individual_from_blocks = ajive.transform_individual(blocks) + for expected, actual in zip(individual, individual_from_blocks): + testing.assert_allclose(expected, actual) + + +def test_ajive_fit_transform_and_clone(multiblock_data): + X, groups, _, _ = multiblock_data + ajive = AJIVE(n_resamples=50, random_state=0) + z = ajive.fit_transform(X, groups=groups) + assert z.shape == (X.shape[0], 2) + + copied = clone(ajive) + copied.fit(X, groups=groups) + testing.assert_allclose(copied.common_components_, ajive.common_components_) + + +def test_ajive_feature_names(multiblock_data): + X, groups, _, _ = multiblock_data + ajive = AJIVE(n_resamples=50, random_state=0).fit(X, groups=groups) + testing.assert_array_equal(ajive.get_feature_names_out(), np.array(["ajive0", "ajive1"])) + + +def test_ajive_input_validation(): + X = np.ones((10, 4)) + with pytest.raises(ValueError, match="`groups` must be provided"): + AJIVE().fit(X) + with pytest.raises(ValueError, match="same number of features"): + AJIVE().fit([np.ones((5, 3)), np.ones((5, 4))]) + with pytest.raises(ValueError, match="initial_ranks"): + AJIVE(initial_ranks=[5, 6]).fit([np.ones((4, 3)), np.ones((4, 3))]) diff --git a/tests/transformer/test_cife.py b/tests/transformer/test_cife.py new file mode 100644 index 0000000..8419829 --- /dev/null +++ b/tests/transformer/test_cife.py @@ -0,0 +1,113 @@ +import numpy as np +import pytest +from numpy import testing +from sklearn.base import clone + +from kalelinear.transformer import CIFE +from tests.utils.test_utils import make_common_individual_dataset + + +@pytest.fixture(scope="module") +def multiblock_data(): + return make_common_individual_dataset(random_state=0) + + +def _subspace_error(estimated, planted): + projection = planted @ planted.T + return np.linalg.norm(estimated - projection @ estimated) / np.sqrt(estimated.shape[1]) + + +def test_cife_recovers_common_and_individual_subspaces(multiblock_data): + X, groups, _, common_basis = multiblock_data + cife = CIFE( + n_common_components=2, + n_individual_components=[3, 4, 2], + tol=1e-10, + random_state=0, + ) + cife.fit(X, groups=groups) + + assert cife.n_common_components_ == 2 + assert cife.common_components_.shape == (X.shape[1], 2) + assert _subspace_error(cife.common_components_, common_basis) < 1e-4 + testing.assert_array_equal(cife.individual_ranks_, [3, 4, 2]) + assert [components.shape[1] for components in cife.individual_components_] == [3, 4, 2] + + +def test_cife_automatic_common_rank(multiblock_data): + X, groups, _, _ = multiblock_data + cife = CIFE(tol=1e-10, random_state=0) + cife.fit(X, groups=groups) + assert cife.n_common_components_ == 2 + + +def test_cife_detects_no_common_structure(): + random_state = np.random.RandomState(1) + blocks = [random_state.randn(50, 4) @ random_state.randn(30, 4).T for _ in range(3)] + cife = CIFE(tol=1e-10, random_state=0) + cife.fit(blocks) + assert cife.n_common_components_ == 0 + assert cife.common_components_.shape == (30, 0) + + +def test_cife_list_input_matches_stacked(multiblock_data): + X, groups, blocks, _ = multiblock_data + from_stacked = CIFE(tol=1e-10, random_state=0).fit(X, groups=groups) + from_list = CIFE(tol=1e-10, random_state=0).fit(blocks) + testing.assert_allclose(from_stacked.common_components_, from_list.common_components_) + + +def test_cife_transform_consistency(multiblock_data): + X, groups, blocks, _ = multiblock_data + cife = CIFE(n_common_components=2, n_individual_components=[3, 4, 2], random_state=0) + cife.fit(X, groups=groups) + + z = cife.transform(X) + testing.assert_allclose(z, X @ cife.common_components_) + testing.assert_allclose(cife.transform(blocks), z) + single_block = cife.transform([blocks[0]]) + testing.assert_allclose(single_block, blocks[0] @ cife.common_components_) + testing.assert_allclose(cife.transform((blocks[0],)), single_block) + + individual = cife.transform_individual(X, groups=groups) + assert [scores.shape for scores in individual] == [(60, 3), (50, 4), (70, 2)] + individual_from_blocks = cife.transform_individual(blocks) + for expected, actual in zip(individual, individual_from_blocks): + testing.assert_allclose(expected, actual) + + +def test_cife_fit_transform_and_clone(multiblock_data): + X, groups, _, _ = multiblock_data + cife = CIFE(n_common_components=2, random_state=0) + z = cife.fit_transform(X, groups=groups) + assert z.shape == (X.shape[0], 2) + + copied = clone(cife) + copied.fit(X, groups=groups) + testing.assert_allclose(copied.common_components_, cife.common_components_) + + +def test_cife_feature_names(multiblock_data): + X, groups, _, _ = multiblock_data + cife = CIFE(n_common_components=2, random_state=0).fit(X, groups=groups) + testing.assert_array_equal(cife.get_feature_names_out(), np.array(["cife0", "cife1"])) + + +def test_cife_full_rank_block_requires_pca_dim(): + random_state = np.random.RandomState(2) + blocks = [random_state.randn(40, 10) for _ in range(3)] + with pytest.raises(ValueError, match="whole feature space"): + CIFE(random_state=0).fit(blocks) + cife = CIFE(pca_dim=0.5, random_state=0) + cife.fit(blocks) + assert cife.n_common_components_ == 0 + + +def test_cife_input_validation(): + X = np.ones((10, 4)) + with pytest.raises(ValueError, match="`groups` must be provided"): + CIFE().fit(X) + with pytest.raises(ValueError, match="same number of features"): + CIFE().fit([np.ones((5, 3)), np.ones((5, 4))]) + with pytest.raises(ValueError, match="(?i)at least two blocks"): + CIFE().fit([np.ones((5, 3))]) diff --git a/tests/utils/test_utils.py b/tests/utils/test_utils.py index c1776a7..614c373 100644 --- a/tests/utils/test_utils.py +++ b/tests/utils/test_utils.py @@ -108,3 +108,65 @@ def make_domain_shifted_dataset( domains = domains[idx] return X, y, domains + + +def make_common_individual_dataset( + n_blocks=3, + n_features=30, + n_common=2, + individual_ranks=(3, 4, 2), + n_samples=(60, 50, 70), + noise=0.0, + random_state=None, +): + """Create multiblock data with planted common and individual subspaces. + + Each block is generated as ``B_c A_c^T + B_i A_i^T (+ noise)`` where + ``A_c`` is a common feature-space basis shared by all blocks and ``A_i`` + is a block-specific basis orthogonal to ``A_c``. + """ + individual_ranks = np.asarray(individual_ranks) + n_samples = np.asarray(n_samples) + if individual_ranks.ndim != 1 or individual_ranks.shape[0] != n_blocks: + raise ValueError( + f"`individual_ranks` must be a sequence with one rank per block: " + f"expected {n_blocks} values, got {individual_ranks.size}." + ) + if n_samples.ndim != 1 or n_samples.shape[0] != n_blocks: + raise ValueError( + f"`n_samples` must be a sequence with one sample count per block: " + f"expected {n_blocks} values, got {n_samples.size}." + ) + random_state = check_random_state(random_state) + common_basis, _ = np.linalg.qr(random_state.randn(n_features, n_common)) + blocks = [] + group_lists = [] + for k in range(n_blocks): + individual_basis, _ = np.linalg.qr(random_state.randn(n_features, individual_ranks[k])) + individual_basis -= common_basis @ (common_basis.T @ individual_basis) + individual_basis, _ = np.linalg.qr(individual_basis) + block = random_state.randn(n_samples[k], n_common) @ common_basis.T + block += random_state.randn(n_samples[k], individual_ranks[k]) @ individual_basis.T + if noise: + block += noise * random_state.randn(n_samples[k], n_features) + blocks.append(block) + group_lists.append(np.full(n_samples[k], k)) + X = np.vstack(blocks) + groups = np.concatenate(group_lists) + return X, groups, blocks, common_basis + + +def test_make_common_individual_dataset_validates_block_parameters(): + with pytest.raises(ValueError, match="individual_ranks"): + make_common_individual_dataset(n_blocks=4) + with pytest.raises(ValueError, match="n_samples"): + make_common_individual_dataset(n_blocks=4, individual_ranks=(1, 2, 3, 4)) + + +def test_make_common_individual_dataset_accepts_custom_n_blocks(): + X, groups, blocks, _ = make_common_individual_dataset( + n_blocks=2, individual_ranks=(1, 2), n_samples=(20, 30), random_state=0 + ) + assert len(blocks) == 2 + assert [block.shape[0] for block in blocks] == [20, 30] + assert groups.shape == (50,)