Linear Independent Component Analysis via Optimal Transport

Summary§

This paper proposes a new approach to Linear Independent Component Analysis (ICA) by directly measuring non-Gaussianity using the squared 2-Wasserstein distance $W_2^2$ to a standard Gaussian. The authors prove that maximizing this distance over linear projections recovers an independent component. They introduce OT-ICA, a gradient-based optimization algorithm that avoids proxy contrast functions like kurtosis or log-likelihoods, which are limited to specific source distributions.

Core Methodology§

ICA assumes observed data $\mathbf{x} \in \mathbb{R}^d$ is a linear mixture of independent sources $\mathbf{s} \in \mathbb{R}^d$: $$ \mathbf{x} = \mathbf{A} \mathbf{s}, $$ with $\mathbf{A}$ invertible. The goal is to find a demixing matrix $\mathbf{W}$ such that $\mathbf{y} = \mathbf{W} \mathbf{x}$ recovers independent components. Classical ICA maximizes non-Gaussianity of $\mathbf{y}$, often via negentropy $J(\mathbf{y}) = H(\mathcal{N}) - H(\mathbf{y})$, where $H$ is entropy. Since negentropy is intractable, proxies like kurtosis or log-likelihoods are used.

The key insight: For a fixed projection $\mathbf{w}$, the squared Wasserstein distance between the distribution of $\mathbf{w}^\top \mathbf{x}$ (assumed standardized to zero mean and unit variance) and a standard Gaussian $\mathcal{N}(0,1)$ is: $$ W_2^2(\mathbb{P}_{\mathbf{w}^\top \mathbf{x}}, \mathcal{N}) = \int_0^1 (F^{-1}(t) - \Phi^{-1}(t))^2 dt, $$ where $F$ and $\Phi$ are CDFs. The authors prove that if the sources have unit variance and are independent, then $W_2^2$ is maximized when $\mathbf{w}$ recovers a single independent component (up to scaling and sign). This is shown using the fact that $W_2^2$ is strictly convex and decomposable across dimensions under independence.

OT-ICA Algorithm§

OT-ICA optimizes a contrast function based on $W_2^2$ between each output dimension and a standard normal. For a demixing matrix $\mathbf{W}$, the loss is: $$ \mathcal{L}(\mathbf{W}) = -\sum_{i=1}^d W_2^2(\mathbb{P}_{(\mathbf{W}\mathbf{x})_i}, \mathcal{N}(0,1)), $$ where the $y_i$ are standardized to have zero mean and unit variance. Minimizing $\mathcal{L}$ corresponds to maximizing non-Gaussianity.

To compute $W_2^2$ from samples $\{z_j = \mathbf{w}^\top \mathbf{x}_j\}_{j=1}^n$, the empirical CDF $\hat{F}$ is used. The squared 2-Wasserstein distance simplifies to the $L^2$ distance between sorted empirical quantiles and theoretical Gaussian quantiles: $$ \hat{W}_2^2 = \frac{1}{n} \sum_{j=1}^n (z_{(j)} - \Phi^{-1}(\frac{j-0.5}{n}))^2, $$ where $z_{(j)}$ are sorted values. This is differentiable w.r.t. $\mathbf{w}$ via the sorted indices, allowing backpropagation.

Implementation Details§

1. Preprocessing: Whitening via PCA to obtain zero-mean, unit-variance data. 2. Initialization: Random orthogonal matrix $\mathbf{W}$. 3. Training: Gradient descent on $\mathcal{L}$ with the empirical $\hat{W}_2^2$. 4. Regularization: Enforce orthogonality of $\mathbf{W}$ via projection steps or a penalty term.

The algorithm does not assume any parametric form for source distributions.

Code Snippet (PyTorch)§

import torch

def wasserstein_distance_projection(X, w, eps=1e-8):
    """
    X: (n, d) data matrix, w: (d,) projection vector
    Returns: squared 2-Wasserstein distance to standard normal
    """
    y = X @ w  # (n,)
    y = y - y.mean()  # center
    y = y / (y.std(unbiased=False) + eps)  # scale to unit variance
    y_sorted, _ = torch.sort(y)
    n = y.shape[0]
    # Gaussian quantiles: inverse CDF of N(0,1)
    q = torch.distributions.Normal(0, 1).icdf(
        torch.linspace(0.5/n, 1 - 0.5/n, n, device=X.device)
    )
    return ((y_sorted - q) ** 2).mean()

# Loss for one output dimension (each output is separate in ICA)
loss = -wasserstein_distance_projection(X, w)  # minimize negative

# For full demixing matrix W (d x d), sum over dimensions
# Use orthogonalization steps during optimization

Empirical Results§

• On simulated data with various source distributions (super-Gaussian, sub-Gaussian, multimodal), OT-ICA consistently recovers sources with lower Amari distance than FastICA (kurtosis) and Infomax (logistic).
• In EEG artifact removal, OT-ICA denoises signals comparably to state-of-the-art while being distribution-agnostic.
• For econometric price discovery, OT-ICA identifies latent factors without assuming Gaussianity.

Conclusion§

OT-ICA provides a principled, distribution-free method for ICA by leveraging optimal transport geometry. The squared Wasserstein distance serves as a direct measure of non-Gaussianity, enabling gradient-based optimization. The method is simple to implement and performs robustly across diverse settings.