otica 0.1.3

Optimal Transport Independent Component Analysis

📊 Contrast-Free ICA

otica is a Python package for linear independent component analysis (ICA) based on optimal transport. It recovers latent sources by maximizing their empirical squared 2-Wasserstein distances to the standard Gaussian, using a fixed non-Gaussianity criterion that requires no user-chosen contrast function or nonlinearity.

---

✨ Features

• Contrast-free source separation: Uses the squared 2-Wasserstein distance to the standard Gaussian as a fixed non-Gaussianity criterion.
• Exact empirical objective: Computes the one-dimensional Wasserstein criterion directly from ordered samples and Gaussian quantiles, without density estimation.
• Riemannian optimization: Optimizes the whitened ICA objective on the orthogonal group with a Picard-style limited-memory BFGS method and Armijo backtracking.
• Dimension reduction: Supports extraction of a specified number of components through principal-component whitening.
• Flexible initialization: Accepts FastICA, random, or user-provided initial unmixing matrices through w_init.
• scikit-learn integration: Native BaseEstimator integration with the standard transformer API, including fit, transform, fit_transform, and inverse_transform.

---

⚡ Method

For observations generated by the linear ICA model $X = S A^\top$, otica first centers and whitens the data. Let $Z \in \mathbb{R}^{d}$ be the whitened random vector. For an orthogonal unmixing matrix $W \in \mathbb{R}^{d \times d}$, the population objective is

$$ F(W) = \sum_{k = 1}^{d} \mathcal{W}_2\left( (Z W^\top)_k, \mathcal{N}(0, 1) \right)^2, \quad W W^\top = I_d. $$

Given $n$ whitened observations collected as the rows of $Z \in \mathbb{R}^{n \times d}$, let $Y = Z W^\top$. OTICA maximizes

$$ \widehat{F}n(W) = \sum{k = 1}^{d} \mathcal{W}2\left( \frac{1}{n} \sum{i = 1}^{n} \delta_{Y_{ik}}, \mathcal{N}(0, 1) \right)^2. $$

For each component $k$, the implementation sorts the entries of the $k$-th column as $Y_{(1)k} \leq \cdots \leq Y_{(n)k}$ and matches them with the corresponding standard-Gaussian rank statistics. Under the usual ICA assumptions, including mutually independent sources with at most one Gaussian component, the population objective identifies the sources up to permutation, sign, and scale.

---

🚀 Installation

pip install otica

🔧 Usage

Example

The following example generates three independent non-Gaussian signals, mixes them linearly, and recovers them with OTICA. Because ICA is identifiable only up to permutation and sign, recovery is evaluated using the best absolute correlation for each true source.

import matplotlib.pyplot as plt
import numpy as np
from otica import OTICA

rng = np.random.default_rng(42)
n_samples = 5_000
time = np.linspace(0.0, 8.0, n_samples)

# Generate independent, non-Gaussian latent sources.
sources = np.column_stack(
    [
        rng.laplace(size=n_samples),
        rng.uniform(-np.sqrt(3.0), np.sqrt(3.0), size=n_samples),
        rng.standard_t(df=5, size=n_samples) * np.sqrt(3.0 / 5.0),
    ]
)

# Mix the sources into the observed signals.
mixing = np.array(
    [
        [1.0, 0.5, -0.2],
        [0.2, 1.0, 0.4],
        [-0.4, 0.1, 1.0],
    ]
)
X = sources @ mixing.T

# Fit OTICA and recover the latent components.
model = OTICA(random_state=42)
estimated_sources = model.fit_transform(X)

correlations = np.corrcoef(sources.T, estimated_sources.T)[:3, 3:]
best_indices = np.abs(correlations).argmax(axis=1)
best_correlations = correlations[np.arange(3), best_indices]
recovered_sources = estimated_sources[:, best_indices] * np.sign(best_correlations)
print("Best absolute correlation per source:", np.abs(best_correlations))

fig, axes = plt.subplots(3, 2, sharex=True, figsize=(12, 6))
for component in range(3):
    axes[component, 0].plot(time[:500], sources[:500, component])
    axes[component, 1].plot(
        time[:500], recovered_sources[:500, component], color="tab:orange"
    )
    axes[component, 0].set_ylabel(f"Source {component + 1}")

axes[0, 0].set_title("True sources")
axes[0, 1].set_title("Recovered sources")
axes[-1, 0].set_xlabel("Time")
axes[-1, 1].set_xlabel("Time")
fig.tight_layout()
plt.show()

---

📖 Learn More

For the mathematical formulation, configuration details, and API reference, visit otica's documentation.