Linear Independent Component Analysis via Optimal Transport
Ashutosh Jha3
Michel Besserve1,2 1
Max Planck Institute for Intelligent Systems, Tübingen 2 Institute of Artificial Intelligence, TU Braunschweig 3 Methods Center, University of Tübingen
arXiv:2607.14081v1 [cs.LG] 15 Jul 2026
Abstract Linear Independent Component Analysis (ICA) recovers jointly independent source signals from their linear mixtures. To achieve this, classical ICA algorithms attempt to maximize non-Gaussianity, measured by negentropy, which is linked to independence by information theory. Because exact negentropy optimization is intractable, they rely on proxy contrast functions, such as fourth-order cumulants, and parametric log-likelihoods. We propose instead to measure non-Gaussianity using the squared Wasserstein distance W22 to a standard Gaussian. We prove that the Wasserstein distance between a standard normal distribution and linear projections of the data is maximized when the projection recovers an independent component. Based on this observation, we propose the OT-ICA algorithm which finds this projection by gradient-based optimization. Empirical evaluation on simulated data shows that OT-ICA outperforms proxy-based methods for different distributions of the latent variables. Application to EEG artifact removal and econometric price discovery confirm OT-ICA can be used for applied ICA tasks without distributional assumptions.
1
Simon Buchholz1
INTRODUCTION
Independent Component Analysis (ICA) [Jutten and Herault, 1985, Comon, 1994] recovers mutually independent source signals Z from a linear mixture X = AS. Identifiability holds under the following conditions. Theorem 1 (Linear ICA [Comon, 1994]). Let X ∈ Rd be observed data and Z ∈ Rd latent sources. Assume (i) X = AZ with A invertible, (ii) z1 , . . . , zd are mutually independent, and (iii) at most one zi is Gaussian. Then A is identifiable up to column permutation and scaling.
Algorithms for learning A from data generally rely on the property that linear mixtures of independent non-Gaussian sources are more Gaussian than their constituents [Hyvärinen et al., 2001], indeed, in the limit of many sources, convergence to a Gaussian follows from the central limit theorem. A rigorous statement of this property was given in Cardoso [2022]. Theorem 2 (Independence and Non-Gaussianity [Cardoso, 2022]). Let Y ∈ Rd be whitened (Cov(Y ) = I). The mutual information I (Y ), which measures statistical dependence among components Y1 , . . . , Yd , satisfies: X I (Y ) = G(Y ) − G(Yi ) (1) i
where G(·) denotes non-Gaussianity defined as the minimal KL divergence to a Gaussian distribution. This result implies that minimizing I (Y ) is equivalent P to maximizing total marginal non-Gaussianity i G(Yi ). However, evaluating G(Yi ) requires the marginal density. Classical algorithms therefore replace it with proxy contrasts: FastICA [Hyvärinen et al., 2001, Amari et al., 1996] uses logcosh negentropy approximations, JADE [Cardoso and Souloumiac, 1993] fourth-order cumulants, InfoMax [Bell and Sejnowski, 1995] a parametric log-likelihood, and Picard [Ablin et al., 2018] an adaptive nonlinearity with a provably converging solver. These proxies exhibit shortcomings such as evaluating to zero on specific non-Gaussian distributions or inducing singular Hessian matrices that stall the fixed-point solver shown in Appendix E. In this paper we use Optimal Transport (OT) [Monge, 1781, Kantorovich, 1942] to define a contrast. This allows us to apply OT to linear ICA and extends the long list of applications of optimal transport techniques in machine learning such as generative modelling [Arjovsky et al., 2017], normality testing [del Barrio et al., 1999], and distribution comparison [Peyré and Cuturi, 2019]. Here we propose to quantify non-Gaussianity as the cost of transporting the (empirical) projection of the observations to a fixed standard Gaussian
reference. We justify this choice in Theorem 3 below which implies that (in the population setting) this contrast is maximized when an independent component is recovered. This then motivates our OT-ICA algorithm, which relies on gradient ascent on the optimal transport cost.
Corollary 1. Under the assumptions of Theorem 3 we have
2
b such that This corollary implies that (3) is maximized by b ⊤ ⊤ b X = b b AZ = Zi for some i, i.e., when the projecb tion recovers an independent component.OT-ICA optimizes the Wasserstein contrast W22 Pb⊤ X e , Γ through gradient ascent. This requires a fully differentiable implementation of the Wasserstein distance. This is straightforward in the one dimensional case, where an optimal coupling is given by matching the percentiles of the distributions, such that the Wasserstein distance is given by Z 1 W22 (µ, ν) = |Fµ−1 (t) − Fν−1 (t)|2 dt (6)
π∈Π(µ,ν)
where Π(µ, ν) denotes (for measures µ, ν on Rd ) the set of all measures π on Rd+d with first marginal µ and second marginal ν. OT-ICA for a single component. We consider the contrast given by W22 Pb⊤ X , Γ , where Γ ≡ N (0, 1) denotes the e standard Gaussian. Since W2 induces a proper metric on probability distributions [Villani, 2003], this contrast is nonzero on all non-Gaussian distributions. OT-ICA maximizes this contrast over the direction b b = arg max W 2 P ⊤ e , Γ . b 2 b X b∈Rd ,|b|=1
(3)
where Fµ denotes the CDF of µ. In particular, the optimal coupling for finitely many samples simply matches the ordered values from the two empirical distributions. OT-ICA for all components. After whitening the data, the projection directions associated to all independent components are mutually orthogonal. While components can be extracted iteratively by searching the shrinking orthogonal complement one by one, Theorem 3 extends to the following corollary 2 which justifies the joint estimation of the entire unmixing matrix B ∈ O(d). Corollary 2. Under the assumptions of Theorem 3, for any orthogonal matrix B ∈ O(d) with rows b⊤ i , the total
Baseline Validation
Amari Error: 0.015
Student-t ( =3)
Beta(0.5, 0.5)
0.6
b⊤X The following theorem establishes that the projection b from (3) indeed recovers an independent component.
W22 (α · Z, Γ) <
d X
αi2 W22 (Zi , Γ).
(4)
i=1
The proof of this result is in Appendix B. We find the following immediate corollary 1.
0.4 0.2 0.0
10 0
1.5 1.0 0.5 0.0
10
Amplitude
Uniform
0.6
0.3
Density
Theorem 3 (W22 ICA Contrast). Let Z = (Z1 , . . . , Zd )⊤ be centered independent sources with unit variance and at most one Zi Gaussian. Assume they all have a smooth and strictly positive density. Then the following bound holds for any normalized weight vector α with at least two nonzero components
Original OT-ICA
Density
We will use the square of the L2 -Wasserstein distance W2 [Villani, 2003], defined as a minimum-cost coupling between probability distributions (see Appendix A): Z 2 2 W2 (µ, ν) ≜ W2 (µ, ν) = inf ∥x−y∥2 dπ(x, y), (2)
0
Density
We consider the standard linear ICA setting X = AZ where A ∈ Rd×d is an invertible matrix. By preprocessing the observations X through centering and whitening, we can assume without loss of generality that E(X) = 0 and Cov[X] = Id×d (see Appendix A for details). Then we can constrain the unmixing matrix B to the Orthogonal Group O(d) (the set of d × d real matrices satisfying BB⊤ = I), and seek B ∈ O(d) such that BX = BAZ has independent components, or equivalently, B inverts A up to a scaled permutation matrix. This can be done by iteratively seeking optimal projection directions b ∈ Rd , corresponding to each row of B, maximizing a non-Gaussianity function G(b⊤ X) of the projected data.
(5)
i.e., the contrast is maximal only for independent components.
Density
THE OT-ICA FRAMEWORK
argmax|α|=1 W22 (α · Z, Γ) ⊂ {e1 , . . . , ed },
0.2 0.1 0.0
2
0
Amplitude
2
2
0
2
Amplitude
Laplace
0.4 0.2 0.0
5
0
Amplitude
5
Figure 1: OT-ICA convergence on four standard continuous source types. Estimated latent components are smoothed with kernel density estimation.
Table 1: IS recovery across 500 Monte Carlo runs without imposing Cholesky ordering.
marginal contrast satisfies: d X
W22 (b⊤ i Z, Γ) ≤
d X
W22 (Zj , Γ),
(7)
Market / Source
True IS
Estimated IS (Mean)
Std. Dev.
where equality holds if and only if B is a signed permutation matrix.
Market 1 Market 2 Market 3
0.1200 0.2400 0.6400
0.1212 0.2399 0.6389
0.0161 0.0245 0.0260
This corollary 2 (proof in Appendix B) theoretically justifies symmetric joint optimization [Hyvärinen et al., 2001]. Rather than iterative deflation, all d rows of B are optimized simultaneously by ascending the total contrast sum, restoring the O(d) constraint via symmetric decorrelation retraction at every step. The full details of the resulting algorithm can be found in Appendix C and pseudocode is given in Algorithm 1.
error 40–45% over FastICA across all three dimensions. The margin is larger on Full Hybrid mixtures (Laplace, Gaussian, Uniform, Student-t, Beta combined): OT-ICA attains E = 0.059, 0.261, 0.335 at d = 10, 20, 30 against FastICA’s 0.255, 0.451, 0.671, a factor of 2–4×.
i=1
j=1
Let us here briefly mention three important algorithmic components utilized in practice: (i) Analytical Gaussian Targets [Fournier and Guillin, 2015], replacing sampled quantile targets of the standard normal distribution with closed-form bin expectations to eliminate approximation noise; (ii) Riemannian Gradient Retraction [Absil et al., 2008], projecting gradients onto the tangent space of Orthogonal Group O(d) and retracting via symmetric decorrelation; and (iii) Gaussian Dithering [Schuchman, 1964], smoothing discrete CDFs via narrow Gaussian convolution.
3
SYNTHETIC EXPERIMENTS
We first evaluate our method on synthetic data1 . Separation quality is measured by the Amari Performance Index (E) [Amari et al., 1996] (defined in Appendix D), which measures the deviation of the global transfer matrix P = Best A from a generalized permutation matrix (E < 0.1: nearperfect; E < 0.3: good separation; E ≥ 0.5: failure). As a sanity check we first consider separation of four different sources with N = 10, 000 samples where OT-ICA achieves E = 0.016 against FastICA’s 0.019 on this mixture (Table 2, Appendix D). Figure 1 illustrates that the latent distributions are accurately recovered. We now study the empirical performance of the proposed algorithm systematically. We use FastICA [Hyvärinen et al., 2001], JADE [Cardoso and Souloumiac, 1993], InfoMax [Bell and Sejnowski, 1995], and Picard [Ablin et al., 2018] as baselines and consider five mixture regimes at d ∈ {10, 20, 30}, N = 10,000, 10 trials per condition. The results can be found in Figure 2. Continuous and heterogeneous regimes. OT-ICA achieves the lowest Amari error on all four continuous and mixedtype configurations at every dimension (Table 3, Appendix D). On Continuous Only sources, OT-ICA reduces 1
Code to reproduce the experiments is available under https://github.com/ashutoshjha3103/ot_in_ linear_ica
Discrete sources. On Discrete only mixtures, OT-ICA underperforms JADE at d = 10 (0.706 vs. 0.443) and d = 20 (1.331 vs. 1.267); note both values exceed E = 0.5, placing all methods in failure territory where differences reflect degree of failure rather than separation quality (a ddimensional random orthogonal matrix yields E ≈ (d − 1)/d, e.g. ≈ 0.97 at d = 30; Appendix D). Table 4 shows W22 retains a clear non-Gaussianity signal on count data, exceeding logcosh resolution by more than 10× on standard Poisson(λ = 3), so the failure is not a contrast deficiency. Step-function CDFs create gradient plateaus that stall the Riemannian solver (Appendix D.5). OT-ICA’s gradient-based Riemannian solver requires more time than FastICA’s fixed-point Newton iterations, with the gap widening as d increases; this cost is offset on heterogeneous mixtures where proxy methods fail to separate independent sources (Appendix F).
4
APPLICATIONS
Structural Identification in Price Discovery. In fragmented financial markets, the Information Share (IS) [Hasbrouck, 1995, Baillie et al., 2002] quantifies each venue’s contribution to price discovery from the structural decomposition ut = Bϵt , with ut the VECM reduced-form residuals and ϵt structural innovations [Engle and Granger, 1987, Johansen, 1991]. The covariance constraint Ω = BB⊤ leaves B unidentified up to an orthogonal rotation; Cholesky identification depends on an arbitrary market ordering. Factoring B = SC [Zema and Cordoni, 2025] reduces identification to finding C, which OT-ICA recovers from the non-Gaussianity of innovations without ordering assumptions (full model in Appendix G). We validate on a simulated three-market VECM [Zema and Cordoni, 2025]: √ √ √ Btrue = diag( 0.12, 0.24, 0.64) encodes true IS values [0.12, 0.24, 0.64] by construction, with Student-t innovations enabling ICA identification. IS estimates match true values to within 0.002 across 500 Monte Carlo runs (see Table 1). EEG Artifact Removal. The skull acts as a linear volume
d = 10
Amari Error ( )
1.5 1.0
Picard
d = 20
1.5 Amari Error ( )
JADE InfoMax
0.5 0.0
1.0 0.5 0.0
d = 30
1.5 Amari Error ( )
Method
OT-ICA FastICA
1.0 0.5 0.0
Continuous Only
Discrete Only
SuperGaussian
Zero Gaussian
Mixture Regime
Full Hybrid
Figure 2: Mean Amari error (↓ better, clipped at 1.5) across five mixture regimes and d ∈ {10, 20, 30} (N = 10,000, 10 trials, 95% CI). Dotted line marks E = 0.3, the good-separation threshold.
conductor [Hyvärinen et al., 2001]: ocular blink artifacts mix with neural signals before reaching scalp electrodes. Blink artifacts are impulsive and super-Gaussian; their high W22 to Gaussian separates them from neural components without specifying a density model (dataset, preprocessing, and configuration in Appendix H). Figure 3 shows OT-ICA concentrating the blink artifact into a single component with excess kurtosis substantially above the remaining four. Zeroing that component and inverting the total unmixing matrix produces the RMS reduction quantified in the panel title, without tuning a density model.
5
CONCLUSION
(a) Raw frontal EEG EEG 001 EEG 002 EEG 003 EEG 004 EEG 005
(b) OT-ICA components [Comp 1: =59.5; others mean =-0.0] Comp 1 Comp 2 Comp 3 Comp 4 Comp 5
(c) Cleaned EEG [92\% RMS reduction in ±250\,ms blink window] EEG 001 EEG 002
OT-ICA replaces proxy contrast functions with the W22 metric computed by quantile sorting; Theorem 3 guarantees the contrast is maximized at independent components, and empirical results confirm OT-ICA outperforms proxy-based algorithms on all continuous and heterogeneous distributions, with the margin widening on heterogeneous mixtures where proxy contrasts fail to extract independent components. Price discovery in simulated VECM and EEG artifact removal demonstrate successful application to ICA tasks without distributional assumptions. On purely discrete sources, all ICA methods fail as d increases (E > 0.5); for OT-ICA the W22 contrast retains signal throughout but
EEG 003 EEG 004 EEG 005
0
2
4
6
Time (s)
8
10
Figure 3: OT-ICA on five frontal EEG channels (MNE sample dataset). Panel (b): isolated blink component (orange) with excess kurtosis κ in the panel title; remaining components have mean κ near zero. Panel (c): signal reconstructed by zeroing the artifact component, with RMS reduction in the ±250 ms blink window in the panel title.
step-function CDFs stall the solver, and per-iteration cost grows with d (Appendix F). These limitations point to concrete improvements: Sinkhorn distances [Cuturi, 2013, Peyré and Cuturi, 2019] as a differentiable substitute for Gaussian dithering, and amortized optimal transport to reduce per-inference sorting cost. Beyond these, OT-ICA’s distribution-free contrast opens new directions in causal discovery, including extensions to LiNGAM [Shimizu et al., 2006], Causal Component Analysis [Liang et al., 2023], and identifiable nonlinear ICA [Buchholz et al., 2022]. References Pierre Ablin, Jean-François Cardoso, and Alexandre Gramfort. Faster independent component analysis by preconditioning with Hessian approximations. IEEE Transactions on Signal Processing, 66(15):4040–4049, 2018. P-A Absil, Robert Mahony, and Rodolphe Sepulchre. Optimization Algorithms on Matrix Manifolds. Princeton University Press, 2008. Shun-Ichi Amari, Andrzej Cichocki, and Howard H Yang. A new learning algorithm for blind signal separation. In Advances in neural information processing systems, pages 757–763, 1996. Martin Arjovsky, Soumith Chintala, and Léon Bottou. Wasserstein generative adversarial networks. In Doina Precup and Yee Whye Teh, editors, Proceedings of the 34th International Conference on Machine Learning, volume 70 of Proceedings of Machine Learning Research, pages 214–223. PMLR, 06–11 Aug 2017. Richard T Baillie, G Geoffrey Booth, Yiuman Tse, and Tatyana Zabotina. Price discovery and common factor models. Journal of Financial Markets, 5(3):309–321, 2002. Anthony J Bell and Terrence J Sejnowski. An informationmaximization approach to blind separation and blind deconvolution. Neural Computation, 7(6):1129–1159, 1995. Yann Brenier. Polar factorization and monotone rearrangement of vector-valued functions. Communications on pure and applied mathematics, 44(4):375–417, 1991. Simon Buchholz, Michel Besserve, and Bernhard Schölkopf. Function classes for identifiable nonlinear independent component analysis. In Advances in Neural Information Processing Systems, volume 35, 2022. Luis A Caffarelli. The regularity of mappings with a convex potential. Journal of the American Mathematical Society, 5(1):99–104, 1992.
Jean-François Cardoso. Independent component analysis in the light of information geometry. Entropy, 24(3):377, 2022. Jean-François Cardoso and Antoine Souloumiac. Blind beamforming for non-Gaussian signals. IEE Proceedings F – Radar and Signal Processing, 140(6):362–370, 1993. Pierre Comon. Independent component analysis, a new concept? Signal processing, 36(3):287–314, 1994. Marco Cuturi. Sinkhorn distances: Lightspeed computation of optimal transport. In Advances in Neural Information Processing Systems, volume 26, 2013. Eustasio del Barrio, Juan A Cuesta-Albertos, Carlos Matrán, and Jesús M Rodríguez-Rodríguez. Tests of goodness of fit based on the L2 -Wasserstein distance. Annals of Statistics, 27(4):1230–1239, 1999. Robert F Engle and Clive WJ Granger. Co-integration and error correction: Representation, estimation, and testing. Econometrica: Journal of the Econometric Society, 55(2): 251–276, 1987. Nicolas Fournier and Arnaud Guillin. On the rate of convergence in wasserstein distance of the empirical measure. Probability Theory and Related Fields, 162(3):707–738, 2015. Joel Hasbrouck. One security, many markets: Determining the contributions to price discovery. The Journal of Finance, 50(4):1175–1199, 1995. Aapo Hyvärinen, Juha Karhunen, and Erkki Oja. Independent Component Analysis. John Wiley & Sons, 2001. Søren Johansen. Estimation and hypothesis testing of cointegration vectors in gaussian vector autoregressive models. Econometrica: Journal of the Econometric Society, 59(6): 1551–1580, 1991. Christian Jutten and Jean Herault. Blind separation of sources, part I: An adaptive algorithm based on neuromimetic architecture. Signal processing, 24(1):1–10, 1985. Leonid V Kantorovich. On the translocation of masses. C. R. (Doklady) Acad. Sci. URSS (N.S.), 37:199–201, 1942. Harold W Kuhn. The hungarian method for the assignment problem. Naval Research Logistics Quarterly, 2(1–2): 83–97, 1955. Te-Won Lee, Mark Girolami, and Terrence J Sejnowski. Independent component analysis using an extended infomax algorithm for mixed subgaussian and supergaussian sources. Neural Computation, 11(2):417–441, 1999.
Wendong Liang, Armin Kekić, Julius von Kügelgen, Simon Buchholz, Michel Besserve, Luigi Gresele, and Bernhard Schölkopf. Causal component analysis. In Advances in Neural Information Processing Systems, volume 36, 2023. Gaspard Monge. Mémoire sur la théorie des déblais et des remblais. Histoire de l’Académie Royale des Sciences de Paris, 1781. Gabriel Peyré and Marco Cuturi. Computational optimal transport: With applications to data science. Foundations and Trends® in Machine Learning, 11(5-6):355– 607, 2019. Leonard Schuchman. Dither signals and their effect on quantization noise. IEEE Transactions on Communication Technology, 12(4):162–165, 1964. Shohei Shimizu, Patrik O Hoyer, Aapo Hyvärinen, and Antti Kerminen. A linear non-gaussian acyclic model for causal discovery. Journal of Machine Learning Research, 7(10), 2006. Cédric Villani. Topics in Optimal Transportation, volume 58. American Mathematical Society, 2003. Sebastiano Michele Zema and Francesco Cordoni. A unifying non-gaussian information share approach to price discovery. Available at SSRN 5234231, 2025.
Supplementary Material for: Linear Independent Component Analysis via Optimal Transport A
BACKGROUND AND MOTIVATION
A.1
CENTERING, WHITENING, AND THE ORTHOGONAL GROUP O(d) SEARCH SPACE
The ICA model X = AS has d2 free parameters in the mixing matrix A. Two lossless preprocessing steps reduce this to a d(d−1) -dimensional rotational search. 2 Centering. Subtracting the mean, X ← X − E[X], ensures the data has zero mean. If X = AS, centering X is equivalent to centering S, leaving the mixing geometry intact. e = QX Whitening. Given the eigendecomposition Cov(X) = VΛV⊤ , the whitening transform Q = Λ−1/2 V⊤ yields X e = I [Hyvärinen et al., 2001]. After whitening, the unmixing operation BX e must satisfy Cov(BX) e = BB⊤ = with Cov(X) I, so B lies in Orthogonal Group O(d). The original mixing matrix decomposes as A = Q−1 B, so ICA on whitened data is lossless [Cardoso, 2022]. A.2
THE CLT MOTIVATION FOR NON-GAUSSIANITY
The Central Limit Theorem (CLT) provides the following intuition: a normalized sum of independent non-Gaussian variables converges to a Gaussian as the number of terms grows [Hyvärinen et al., 2001]. In ICA, any linear mixture b · S with multiple nonzero weights is therefore more Gaussian than the individual sources. Recovering an independent source means finding the projection that is least Gaussian, motivating the maximization of non-Gaussianity as the ICA objective. A.3
INFORMATION GEOMETRY: CONNECTING INDEPENDENCE AND NON-GAUSSIANITY
Following Cardoso [2022], ICA can be understood geometrically via two Pythagorean identities. Let G(·) = DKL (· ∥ N (Cov ·)) denote non-Gaussianity and C(Y ) = DKL (N (Cov Y ) ∥ N (diag Cov Y )) denote correlation. Q Product manifold identity. For any product distribution Q = i qi : X DKL (PY ∥Q) = I (Y ) + DKL (PYi ∥qi ). (8) i
Gaussian manifold identity. For any Gaussian N (Σ): DKL (PY ∥N (Σ)) = G(Y ) + DKL (N (Cov Y )∥N (Σ)). Evaluating the KL divergence from PY to N (diag Cov Y ) via both paths yields Cardoso’s identity: X I (Y ) + G(Yi ) = G(Y ) + C(Y ).
(9)
(10)
i
P After whitening, Cov(Y ) = I forces C(Y ) = 0, giving I (Y ) = G(YP ) − i G(Yi ). Since G(Y ) is invariant under orthogonal rotation, minimizing I (Y ) (independence) equals maximizing i G(Yi ) (total marginal non-Gaussianity). This is Theorem 2. A.4
OPTIMAL TRANSPORT AS THE NON-GAUSSIANITY CONTRAST
Classical ICA algorithms approximate G(Yi ) with proxies. OT-ICA instead measures it directly as the cost of transporting the empirical marginal to a standard Gaussian. Definition 1 (Coupling and Push-Forward). A coupling of distributions µ and ν is a joint distribution π on Rd × Rd with marginals µ and ν; the set of all such couplings is Π(µ, ν). A measurable map T : Rd → Rd pushes forward µ to ν (written T# µ = ν) if ν(B) = µ(T −1 (B)) for every Borel set B.
The p-Wasserstein distance is the minimum transport cost over all couplings: 1/p Z p Wp (µ, ν) = inf ∥x − y∥ dπ(x, y) . π∈Π(µ,ν)
(11)
The structure of the optimal coupling π varies qualitatively depending on whether the marginals are discrete, continuous, or mixed (Peyré and Cuturi [2019]). In 1D, Brenier’s theorem guarantees the optimal coupling is induced by the monotone R1 quantile map T = Fν−1 ◦ Fµ , so W22 (µ, ν) = 0 |Fµ−1 (t) − Fν−1 (t)|2 dt. ⊤e Taking ν = Γ ≡ N (0, 1), the quantity W22 (Pb⊤ X e , Γ) measures how far the projection b X is from Gaussian. Unlike 2 negentropy proxies, W2 is a true metric, geometry-aware, and computable by sorting without density estimation. OT-ICA maximizes this over b in Orthogonal Group O(d); Theorem 3 shows the maximum is attained at pure source directions.
B
MATHEMATICAL PROOFS FOR OT-ICA BOUNDS
B.1
PROOF OF THEOREM 3: W22 ICA CONTRAST
Theorem 3 (W22 ICA Contrast). Let Z = (Z1 , . . . , Zd )⊤ be centered independent sources with unit variance and at most one Zi Gaussian. Assume they all have a smooth and strictly positive density. Then the following bound holds for any normalized weight vector α with at least two nonzero components d X
W22 (α · Z, Γ) <
αi2 W22 (Zi , Γ).
(12)
i=1
Strategy. We first establish a weak upper bound on W22 (α · Z, Γ) (used as an intermediate step), then show the bound is strict for any proper mixture under the non-Gaussianity conditions, proving Theorem 3. Step 1: Upper Bound (intermediate). Proof. We construct a specific joint probability distribution, or coupling, using a common source of randomness. Let N = (N1 , . . . , Nd )⊤ be a vector of d independent standard normal variables, N ∼ N (0, Id ), assumed to be statistically independent of the sources Z. For each independent latent source Zi , there exists an optimal transport map Ti that pushes forward the standard normal distribution to the source distribution, denoted as (Ti )# Γ = Zi . We construct a random variable X representing the mixture: X=
d X
αi Ti (Ni ).
(13)
i=1
To compare this mixture to a Gaussian distribution, we construct a target variable Y using the same underlying noise vector N: d X Y = αi Ni . (14) i=1
Given that the weight vector has unit norm ( i αi2 = 1), the resulting variable Y follows the standard normal distribution, Y ∼ Γ. The pair (X, Y ) creates a valid coupling. Because W22 is the infimum over all valid couplings, the optimal distance P
is less than or equal to the cost of our constructed coupling: !2 d X . W22 (α · Z, Γ) ≤ E |X − Y |2 = E αi Ti (Ni ) − Ni
(15)
i=1
Squaring the summation produces diagonal terms (i = j) and cross terms (i ̸= j): |X − Y |2 =
d X i=1
αi2 Ti (Ni ) − Ni
2
+
X i̸=j
αi αj Ti (Ni ) − Ni Tj (Nj ) − Nj .
(16)
Because the underlying noise variables Ni and Nj are independent, and the functions evaluate to a mean of zero, when we take the expectation, all cross terms vanish. We are left with the expectation of the diagonal terms: d h X 2 i E |X − Y |2 = αi2 E Ti (Ni ) − Ni .
(17)
i=1
Recognizing that E[(Ti (Ni ) − Ni )2 ] is the definition of the squared Wasserstein distance between the individual source Zi and Γ, we arrive at the bound: d X W22 (α · Z, Γ) ≤ αi2 W22 (Zi , Γ). (18) i=1
Step 2: Strict Inequality (proof of Theorem 3). In a one-dimensional setting, the W2 distance equals the expected squared difference of a specific coupling if and only if the random variables are comonotonic. By Brenier’s Theorem [Brenier, 1991], this requires the gradients to be parallel at every point in the probability space: ∇X(N) = λ(N)∇Y (N).
(19)
Assuming the source distributions possess smooth and strictly positive densities, Caffarelli’s regularity theory [Caffarelli, 1992] guarantees that the optimal transport maps Ti are continuously differentiable. 2
X We differentiate both sides to compute their Hessian matrices. On the LHS, the cross-derivatives ∂N∂i ∂N are zero, yielding a j ′′ diagonal matrix of αi Ti (Ni ). On the RHS, we differentiate the product λ(N)b, yielding a Rank-1 outer product matrix: ∂λ α1 T1′′ (N1 ) 0 ··· 0 ∂λ α1 ∂N1 · · · α1 ∂N ′′ d 0 α T (N ) · · · 0 2 2 2 .. .. . . (20) = .. .. .. .. . . . . . . . ∂λ ∂λ αd ∂N1 · · · αd ∂Nd 0 0 · · · αd Td′′ (Nd ) {z } | {z } |Outer Product b(∇λ) ⊤ (Rank-1) Hessian of X (Diagonal)
∂λ = 0. Because we assume a Examining any off-diagonal entry (i ̸= j), the LHS dictates it must be zero, establishing αi ∂N j mixture where multiple weights are non-zero (αi ̸= 0), it follows that the gradient of the scaling function is zero (∇λ = 0). Consequently, λ(N) is a global scalar constant, λ.
Equating the diagonal elements implies Ti′ (Ni ) = λ. Integrating this indicates the optimal transport maps are linear functions of the form: Ti (x) = λx + c. (21) Applying a purely linear transformation to Gaussian noise Ni produces another Gaussian distribution. The identifiability assumption of ICA dictates that the original latent sources Zi are non-Gaussian. Therefore, the linear map requirement contradicts the non-Gaussianity of the sources. Because the condition for equality cannot be met, the relationship is a strict inequality.
W22 (α · Z, Γ) <
d X
αi2 W22 (Zi , Γ).
(22)
i=1
B.2
PROOF OF COROLLARY 2: SYMMETRIC JOINT MAXIMIZATION
Proof. Let B ∈ O(d) be an orthogonal matrix with rows b⊤ i . Because B is orthogonal, its squared elements sum to 1 across Pd Pd 2 2 both rows ( j=1 Bij = 1) and columns ( i=1 Bij = 1). By the intermediate step of Theorem 3, the contrast for each row projection is bounded by: d X 2 W22 (b⊤ Z, Γ) ≤ Bij W22 (Zj , Γ). (23) i j=1
Summing this inequality over all d rows yields an upper bound on the total marginal contrast: d X
W22 (b⊤ i Z, Γ) ≤
i=1
d X d X
2 Bij W22 (Zj , Γ)
i=1 j=1
=
d d X X j=1
=
d X
! 2 Bij
W22 (Zj , Γ)
(24)
i=1
W22 (Zj , Γ).
j=1
The maximum possible sum equals the unweighted sum of the pure sources’ contrasts. Because Theorem 3 establishes that the inequality is strict whenever a row mixes multiple sources, the global equality holds if and only if every row bi contains exactly one non-zero entry (which must be ±1 due to the unit norm constraint). Thus, the sum is maximized if and only if B is a signed permutation matrix.
C
ALGORITHMIC ENHANCEMENTS AND METHODOLOGY
This appendix explains the additional ingredients used in the OT-ICA algorithm. C.1
ANALYTICAL GAUSSIAN TARGETS
Discrete approximation of the target Gaussian distribution introduces approximation noise. To eliminate this noise, we replace point-sampled targets with an analytical formulation, computing the expected value of a standard normal variable within each discrete quantile bin. Definition 2 (Analytical Gaussian Target). Let the uniform probability bin edges for N samples be defined as pi = Ni for i = 0, . . . , N , with corresponding Gaussian domain boundaries zi = Φ−1 (pi ), where Φ is the standard normal CDF. The analytical target value Ti for the i-th sorted sample is: Z zi Ti = N xϕ(x) dx = N ϕ(zi−1 ) − ϕ(zi ) . (25) zi−1
C.2
RIEMANNIAN GRADIENT AND RETRACTION
Because the input data is whitened, the unmixing matrix B must remain in Orthogonal Group O(d), satisfying BB⊤ = I. Computing the standard Euclidean gradient G = ∇B W22 points into the unconstrained ambient space. Definition 3 (Riemannian Gradient on Orthogonal Group O(d)). Given the Euclidean gradient G, the Riemannian gradient resulting from the projection onto the tangent space TB O(d) of Orthogonal Group O(d) is given by [Absil et al., 2008]: ∇O(d) B = G −
1 GB⊤ + BG⊤ B. 2
(26)
Definition 4 (Symmetric Decorrelation Retraction). To map the matrix back to the orthogonal surface after a tangent step (Bstep ), we compute the overlap covariance C = Bstep B⊤ step and apply the inverse square root: Bnew = C−1/2 Bstep . C.3
(27)
CONTINUOUS SMOOTHING VIA GAUSSIAN DITHERING
Applying continuous optimal transport directly to discrete mixtures introduces non-smooth step-functions. Adapted from 2 signal processing [Schuchman, 1964], we inject continuous, zero-mean Gaussian noise with variance σdither = 0.01 into the 1D projected components prior to sorting. This decorrelates deterministic quantization errors and convolves the discrete empirical distribution with a Gaussian kernel.
G = ∇ B W 22 ∇S B
B
TB S
S
B new
(Orthogonal Matrices Space)
Figure 4: Geometric sketch of optimization on Orthogonal Group O(d). The Euclidean gradient G is projected onto the tangent space to yield the Riemannian gradient ∇O(d) B. A retraction maps the estimate back to the orthogonal surface.
Cumulative Probability Fw (y)
1.0 0.8 2 ∼ N(0, σdither )
0.6 0.4 0.2 0.0
Discrete Empirical CDF (Staircase) Dithered CDF (Continuous)
0
1
2
Projected Value y = w > x
3
Figure 5: The non-differentiable staircase of a discrete empirical CDF (gray) is convolved with a narrow Gaussian kernel via dithering, yielding a continuous, differentiable curve (black).
C.4
TOTAL COMPLEXITY AND STATISTICAL EFFICIENCY
For a dataset of dimension d with N samples, K random restarts (batched into a tensor Bbatch ), and T optimization iterations, the per-iteration complexity is: O(K · (d · N log N + d3 ))
(28)
The N log N term represents the sorting of projections, while the d3 term represents the cost of the symmetric decorrelation used for retraction.
C.5
THE OT-ICA ALGORITHM
The mutual orthogonality of ICs in the whitened space admits an iterative extraction strategy in which each component is found in the orthogonal complement of previously extracted ones. We implement this as symmetric joint optimization: all rows of B are updated simultaneously under the O(d) constraint via symmetric decorrelation retraction. Symmetric optimization is preferred because it eliminates error accumulation across extraction steps and ensures mutual orthogonality of all components is maintained at every iteration rather than enforced only by projection. Algorithm 1 The Optimal Transport ICA (OT-ICA) Algorithm Require: Observed mixture X ∈ Rd×N , iterations K, learning rate η, batch size Nbatch , dither noise σdither Ensure: Estimated unmixing matrix Bfinal ∈ Rd×d 1: Phase 1: Preprocessing & Target Generation 2: Center the data: X ← X − E[X] e ← (XX⊤ )−1/2 X 3: Whiten the data: X 4: Compute analytical target T ∈ RNbatch via analytical Gaussian CDF integration 5: Phase 2: Initialization 6: Initialize B ∈ Rd×d randomly from N (0, 1) 7: Apply symmetric decorrelation: B ← (BB⊤ )−1/2 B 8: Phase 3: Stochastic Riemannian Optimization 9: for k = 1 to K do e e batch ∈ Rd×Nbatch from X 10: Sample stochastic mini-batch X e 11: Project data onto candidates: Y ← BXbatch 2 12: Apply Gaussian Dithering: Ydither ← Y + N (0, σdither ) 13: Sort each row of Ydither ascending to yield empirical quantiles Ysorted Pd (i) 1 2 14: Compute Wasserstein objective: L = Nbatch i=1 ||Ysorted − T||2 15: Backpropagate to compute Euclidean gradient: G = ∇B L 16: Project to tangent space of Orthogonal Group O(d): ∇O(d) B = G − 21 GB⊤ + BG⊤ B 17: Update matrix along tangent plane: Bstep = B + η∇O(d) B −1/2 18: Retract to Orthogonal Group O(d): B ← (Bstep B⊤ Bstep step ) 19: end for 20: Phase 4: Finalization 21: Bfinal ← B(XX⊤ )−1/2 22: return Bfinal
D
EXPERIMENTAL METHODOLOGY AND METRICS
D.1
THE GENERALIZED PERMUTATION MATRIX AND AMARI INDEX
To empirically evaluate the success of source separation, we quantify the distance between the estimated unmixing matrix Best and the true mixing matrix A. Because ICA is subject to scaling and permutation ambiguities, the target recovery is a generalized permutation matrix. Let P = Best A be the global transfer matrix. For a d × d matrix P with elements pij , the
Amari Performance Index [Amari et al., 1996] measures the structural divergence: ! d d d d X X 1 X X 1 |pij | |p | ij E(P) = − 1 + −1 . 2d i=1 j=1 maxk |pik | 2d j=1 i=1 maxk |pkj |
D.2
(29)
COMPUTATIONAL RESOURCE REGIMES
Because OT-ICA utilizes stochastic gradient descent while FastICA employs fixed-point Newton iterations, we standardize resources via two defined regimes to ensure fairness: • Low Compute Baseline: OT-ICA is allocated restarts min(4 × d, 150) and 150/200 deflation/symmetric phase iterations. FastICA is allocated 50 random restarts and 1,000 maximum iterations per component. • High Compute Regime: OT-ICA is allocated restarts min(15 × d, 600) and 300/500 iterations. FastICA is allocated 50 restarts and 5,000 maximum iterations per component.
D.3
METHODOLOGY VALIDATION: GLOBAL TRANSFER MATRICES
We validated OT-ICA on a 4D mixture of Laplace, Uniform, Student-t(3), and Beta(0.5, 0.5) sources with N = 10,000, mixed by the random matrix A below. The global transfer matrices P = B̂A for both algorithms are near-permutation, with Amari errors below 0.02. True Mixing Matrix A
1.058 0.042 −0.154 −0.006
1.520 −0.486 −1.572 0.503
OT-ICA Transfer Matrix P −0.012 −0.009 −1.000 −0.003 0.000 0.003 −0.001 1.000 −0.007 1.000 0.000 0.012 −1.000 −0.007 0.005 −0.010
−0.249 0.388 −0.304 0.053
1.012 −0.523 −1.234 0.928
FastICA Transfer Matrix P
−1.000 0.012 −0.002 −0.010
Amari Error: OT-ICA = 0.0155
−0.012 −1.000 −0.010 −0.010
0.003 0.001 −0.001 −1.000
−0.008 −0.004 −1.000 −0.001
FastICA = 0.0188
Table 2: Global transfer matrices for the 4D validation mixture. Off-diagonal entries are ≤ 0.012 in absolute value for both algorithms, confirming near-perfect source recovery.
D.4
FULL BENCHMARK RESULTS
Table 3 reports mean Amari errors for all five methods across all mixture regimes at d ∈ {10, 20, 30} (N = 10,000, 10 trials). Values at 1.500 indicate total separation failure (clip ceiling). Bold marks the lowest error per row. The threshold E < 0.3 separates meaningful separation from failure: for a d-dimensional random orthogonal matrix, the expected Amari error approaches d−1 d (e.g., ≈ 0.97 at d = 30), so differences among methods above E = 0.3 reflect degree of failure rather than quality of separation; fine-grained comparisons in that region are nonetheless preserved for completeness. D.5
DISCRETE DISTRIBUTIONS: CONTRAST VS. OPTIMIZATION
The Discrete Only failure is mechanistically distinct from the continuous failure modes above. Table 4 compares the non-Gaussianity resolution of W22 against the logcosh proxy negentropy (E[G(x)] − E[G(ν)])2 across standard and highly non-Gaussian parameterizations.
Regime
d
OT-ICA
FastICA
JADE
InfoMax
Picard
Continuous Only
10 20 30
0.050 0.106 0.160
0.092 0.183 0.271
0.146 0.321 0.489
0.100 0.188 0.277
0.093 0.184 0.271
Strictly Super-Gaussian
10 20 30
0.057 0.113 0.176
0.089 0.183 0.280
0.146 0.316 0.482
0.103 0.190 0.286
0.090 0.184 0.280
Zero Gaussian
10 20 30
0.042 0.103 0.207
0.118 0.327 0.510
0.152 0.366 0.686
0.105 0.315 0.549
0.121 0.403 0.628
Full Hybrid
10 20 30
0.059 0.261 0.335
0.255 0.451 0.671
0.216 0.571 0.938
0.297 0.558 0.765
0.277 0.555 0.917
Discrete Only
10 20 30
0.706 1.331 1.500
0.613 1.402 1.500
0.443 1.267 1.500
0.757 1.495 1.500
0.835 1.440 1.500
√ Pure Laplace (µ = 0, , b = 1/ 2)
10 20 30
0.077 0.165 0.280
0.080 0.170 0.261
0.136 0.285 0.439
0.083 0.177 0.272
0.080 0.170 0.262
√ √ Pure Uniform (− 3, , 3)
10 20 30
0.038 0.094 1.500
0.054 0.116 0.179
0.051 0.110 0.171
0.056 0.122 0.187
0.054 0.116 0.180
Pure Student-t (ν = 3)
10 20 30
0.068 0.142 0.219
0.065 0.135 0.207
0.134 0.284 0.436
0.064 0.136 0.208
0.065 0.136 0.208
Pure Chi-square χ2 (k = 2)
10 20 30
0.039 0.085 0.129
0.094 0.205 0.315
0.120 0.255 0.392
0.087 0.190 0.292
0.094 0.205 0.316
Pure Exponential (λ = 1)
10 20 30
0.039 0.085 0.129
0.094 0.205 0.315
0.120 0.255 0.392
0.087 0.190 0.292
0.094 0.205 0.316
Table 3: Mean Amari error (↓ better, clipped at 1.500) for all five methods and mixture regimes (N = 10,000, 10 trials). OT-ICA achieves the lowest error on all five mixed-source configurations. On pure single-distribution sources, OT-ICA leads on Chi-square and Exponential at all d; FastICA leads on Student-t and Laplace (d = 30); OT-ICA fails on Uniform at d = 30. On Discrete Only, JADE leads at d ≤ 20; all methods including OT-ICA saturate the clip ceiling at d = 30.
Distribution
W22
Logcosh Negentropy
Laplace (continuous baseline) Binomial (standard, n = 10) Poisson (standard, λ = 3.0) Binomial (non-Gaussian, n = 2) Poisson (non-Gaussian, λ = 0.5)
0.0398 0.0344 0.0513 0.2022 0.3384
0.001214 0.000031 ≈0 0.000267 0.000085
Table 4: W22 provides ≫10× greater non-Gaussianity resolution than logcosh on count-based data. Standard Poisson (λ = 3.0) evaluates to near-zero under logcosh yet retains a clear W22 signal. The failure to unmix discrete sources therefore cannot originate from an insufficient W22 contrast; the issue is in the optimization landscape.
Discrete Mixture Optimization Landscape: OT-ICA vs.FastICA Highly Non-Gaussian Poisson ( = 0.5): = 0.5
2
0.05
0.0001
0.00
0.0000
0
/2
Projection Angle (Radians) OT-ICA (W22 Distance)
W2 Distance (W22)
Proxy Negentropy
0.0002
0.10
1.00
0.03
0.75
0.02
0.50
0.01
0.25
0.00
0.0003
0.15
1.25
0.04
Standard Binomial (n = 10): n = 10
1e 5
0.03
4
0.02
3
Proxy Negentropy
Highly Non-Gaussian Binomial (n = 2): n = 2
0
1e 5
Proxy Negentropy
4
0.1
0.20
W2 Distance (W22)
6
0.2
0.0
W2 Distance (W22)
8
0.3
Standard Poisson ( = 3.0): = 3.0
0.05
Proxy Negentropy
W2 Distance (W22)
1e 5
2
0.01 0.00
1 0
/2
Projection Angle (Radians)
FastICA (Logcosh Negentropy)
Figure 6: W22 and logcosh optimization landscapes for two-source Poisson and Binomial mixtures across parameterizations. The W22 landscape retains a non-zero signal throughout (confirming the contrast is present), but the step-function geometry of discrete CDFs creates flat plateaus where the gradient evaluates to near-zero, stalling gradient-based solvers. The logcosh landscape degenerates to near-zero for standard count-based parameterizations, indicating a fundamental contrast failure.
E
PROXY CONTRAST BASED ICA ALGORITHMS
This appendix collects a brief overview of the baselines used in the paper. E.1
FASTICA: NEGENTROPY APPROXIMATION FAILURES
FastICA [Hyvärinen et al., 2001] approximates negentropy using a non-quadratic contrast function G(·), such as the logcosh function (E[G(x)] − E[G(ν)], where ν ∼ N (0, 1)). We demonstrate two limitations of this proxy approach. E.1.1
The Zero Negentropy Condition
If a latent independent source possesses a non-Gaussian distribution such that its expected value under the contrast function equals that of a Gaussian (E[G(si )] = E[G(ν)]), the approximated negentropy evaluates to zero. In this scenario, the objective function provides no gradient signal, and the algorithm fails to extract the source. As shown in Figure 7, this failure leads to increasing Amari error with increasing dimensions.
Amari Error vs. Dimension: The Zero-Negentropy Condition
Amari Error (Lower is Better)
4.0 3.5 3.0
Low Compute Regime
High Compute Regime
Method FastICA OT-ICA
2.5 2.0 1.5 1.0 0.5 0.0 5
Good Separation 10
15
Number of Dimensions
20
5
10
15
Number of Dimensions
20
Figure 7: Amari error comparison for an engineered zero-negentropy distribution across Low and High compute regimes. FastICA’s error scales with dimension as the proxy gradient vanishes, crossing the good separation threshold (E = 0.3). OT-ICA maintains good separation till higher dimensions depending on compute regime.
E.1.2
The Vanishing Curvature Condition
FastICA utilizes a Newton fixed-point iteration. Defining g(x) = G′ (x) and g ′ (x) = G′′ (x), the algorithm approximates the Hessian matrix. For a given component i, the diagonal entry scales according to the expectation Hii ≈ E si g(si ) − g ′ (si ) . The Newton update requires applying the inverse of this Hessian (H−1 ). If a non-Gaussian source distribution causes the expectation in the denominator to evaluate to zero (E[si g(si ) − g ′ (si )] → 0), the inverse becomes undefined [Hyvärinen et al., 2001, Chapter 8]. This analytical divide-by-zero error makes the fixed-point update step intractable. Figure 8 illustrates this failure mode, displaying the specific trimodal distribution that induces zero curvature and the resulting divergence in separation performance compared to our OT-ICA framework. E.2
JADE: SAMPLE COMPLEXITY OF FOURTH-ORDER CUMULANTS
JADE [Cardoso and Souloumiac, 1993] identifies sources by jointly diagonalizing a set of fourth-order cumulant matrices via Jacobi sweeps. Reliable estimation of the full cumulant tensor requires O(d4 ) samples [Cardoso and Souloumiac, 1993, Hyvärinen et al., 2001]. At d = 30, this implies approximately 810,000 samples; against the N = 10,000 available in our benchmark, the cumulant matrices are heavily undersampled. The resulting joint diagonalization converges on estimated
Zero Curvature: Engineered Trimodal Gaussian
Probability Density
2.5
Trimodal PDF (KDE) Side Peaks (b = ±1.44) Negative Curvature Zone Empirical Histogram (Data)
2.0 1.5 1.0 0.5 0.0 4
h(x) Amplitude
3 2 1 0 1
4
3
2
1
0
FastICA Curvature: h(x) 2 3 4
1
Signal Amplitude (x)
Amari Error vs. Dimension: The Vanishing Curvature Condition
Amari Error (Lower is Better)
4.0 3.5 3.0
Low Compute Regime
High Compute Regime
Method FastICA OT-ICA
2.5 2.0 1.5 1.0 0.5 0.0 5
Good Separation 10
15
Number of Dimensions
20
5
10
15
Number of Dimensions
20
Figure 8: The Vanishing Curvature failure mode. Top: A trimodal distribution engineered such that negative central curvature and positive outer curvature (side peaks at b = ±1.44) cancel (E[sg(s) − g ′ (s)] = 0) while maintaining unit variance. Bottom: Amari error comparison on mixtures containing this distribution. FastICA’s Newton solver diverges rapidly as dimension increases across both compute regimes. OT-ICA maintains good separation.
tensors reflecting sampling noise rather than true source structure, yielding E = 0.489 on continuous sources and E = 0.938 on hybrid sources at d = 30, both approaching the failure threshold despite JADE’s consistency on well-sampled data. E.3
INFOMAX: SCORE FUNCTION MISSPECIFICATION
InfoMax [Bell and Sejnowski, 1995] maximizes a log-likelihood under a fixed logistic nonlinearity, implicitly treating all sources as sub-Gaussian. On mixtures containing super-Gaussian components (Laplace, Student-t), the score function is misspecified: the gradient points away from the true unmixing direction. Lee et al. [1999] documented this failure mode explicitly, proposing extended InfoMax with online source-type switching. That approach partially addresses the misspecification but requires correct per-component classification, an assumption OT-ICA avoids entirely. In our benchmark, InfoMax produces E = 0.765 on hybrid sources and E = 0.549 on Zero Gaussian sources at d = 30, against OT-ICA’s 0.335 and 0.207 respectively, consistent with the misspecified-score explanation. E.4
PICARD: THE CONTRAST AS THE BINDING CONSTRAINT
Picard [Ablin et al., 2018] replaces InfoMax’s gradient ascent with an L-BFGS preconditioned solver on the same loglikelihood objective, achieving provably faster convergence per iteration. In our benchmark, FastICA and Picard produce Amari errors within 0.003 of each other across all continuous configurations at each dimension. This near-identical performance isolates the source of failure: solver speed is not the binding constraint. The shared tanh nonlinearity, and not the optimizer, determines the error floor. Replacing it with an exact, distribution-free contrast reduces error by 40–45% on continuous sources.
F
COMPUTATIONAL SCALING AND DIMENSIONALITY LIMITS
We benchmarked OT-ICA against FastICA using a linear mixture of continuous Laplace sources across increasing dimensions with a fixed sample size of N = 10, 000.
Amari Scores: FastICA vs OT-ICA Laplacian Independent Sources
Amari Error (Lower is Better)
3.0 2.5
Low Compute Regime
High Compute Regime
Method
FastICA OT-ICA
2.0 1.5 1.0 0.5 Good Separation 0.010
15
20
25
30
35
40
45
Number of Dimensions
50
55
60 10
15
20
25
30
35
40
45
Number of Dimensions
50
55
60
Figure 9: Amari error across dimensions for Laplacian sources at fixed N = 10,000 samples. Both OT-ICA and FastICA hit the same curse-of-dimensionality noise ceiling at d = 40, a finite-sample sparsity effect shared by all empirical ICA methods at this regime.
While OT-ICA maintains error rates comparable to FastICA (E < 0.3) in lower dimensions, both algorithms exhibit a loss in unmixing quality at d = 40. This limit is a manifestation of the curse of dimensionality. The expected error of the empirical Wasserstein distance scales as O(N −1/d ). As d grows, a fixed finite sample fails to densely populate the state space, creating massive empty regions that swallow the gradient signal. The execution time for OT-ICA scales more steeply than FastICA. FastICA costs O(dN ) per iteration, making it computationally superior for exclusively homogeneous, continuous mixtures that do not trigger its proxy blind spots. OT-ICA incorporates parallel sorting across multiple restarts (O(K · dN log N )), trading execution speed for reliability on complex topologies.
Time Complexity: FastICA vs OT-ICA Laplacian Independent Sources Low Compute Regime
Total Time (Seconds)
1000
High Compute Regime
Method
FastICA OT-ICA
800 600 400 200 0 10
15
20
25
30
35
40
45
Number of Dimensions
50
55
60 10
15
20
25
30
35
40
45
50
Number of Dimensions
55
60
Figure 10: Total execution time (in seconds) for FastICA and OT-ICA across dimensions.
G
PRICE DISCOVERY APPLICATION: DATA GENERATING PROCESS AND IS DEFINITION
G.1
DATA GENERATING PROCESS
We simulate a three-market VECM following the non-Gaussian information-share framework of Zema and Cordoni [2025]. The three markets share a single common efficient price, with long-run impact vector ψ = (1, 1, 1)⊤ (each market moves one-for-one with the common trend). The true structural mixing matrix is diagonal: √ √ √ Btrue = diag 0.12, 0.24, 0.64 , (30) which, as shown below, yields the target Information Shares [0.12, 0.24, 0.64] by construction. Structural innovations ϵt are drawn from Student-t distributions with degrees of freedom ν1 = 5, ν2 = 6, ν3 = 7, normalized to unit variance. The VECM reduced-form residuals are ut = Btrue ϵt . Heavy tails are a realistic feature of high-frequency financial innovations and are what enables ICA identification: the non-Gaussianity of ut carries the rotation information that a Gaussian model would lose. G.2
INFORMATION SHARE DEFINITION
Under Hasbrouck [1995], the Information Share of venue i measures its contribution to price discovery: ISi =
(ψ ⊤ B)2i . ψ ⊤ BB⊤ ψ
2 With ψ = (1, 1, 1)⊤ and diagonal Btrue , the numerator equals Bii and the denominator equals 2 0.64 = 1, so ISi = Bii = {0.12, 0.24, 0.64}.
G.3
(31) 2 j Bjj = 0.12 + 0.24 +
P
ESTIMATION PROCEDURE
The reduced-form covariance is Ω = BB⊤ . Cholesky or eigendecomposition of Ω yields the whitening transform S = Ω1/2 ; the orthogonal rotation C satisfying B = SC is then identified from the non-Gaussianity of S−1 ut by OT-ICA. In this simulation Ω is treated as known; in practice it is estimated from VECM residuals by OLS. Two economic constraints are imposed post-estimation: the Hungarian algorithm [Kuhn, 1955] resolves the permutation ambiguity by assigning each structural shock to the market for which it explains maximum variance; sign normalization enforces a positive own-market impact. IS estimates are then computed from B̂ = SĈ using the Hasbrouck formula above.
H
EEG APPLICATION DETAILS
H.1
DATASET AND PREPROCESSING
We use the MNE sample dataset (sample_audvis_raw.fif), a publicly available MEG/EEG recording with auditory and visual stimuli. Five frontal EEG channels (EEG 001–005) are extracted and cropped to a 10 s segment (10–20 s from onset). A zero-phase FIR bandpass filter (1–40 Hz, firwin design) is applied before ICA to remove DC drift and high-frequency noise. Signals are then z-score normalized per channel: x̃i = (xi − µi )/σi . H.2
OT-ICA CONFIGURATION
Deflation phase: 50 random restarts, 200 iterations per restart, dither σ = 0.01. Symmetric refinement phase: 400 gradient steps, learning rate η = 0.05, mini-batch size 512, dither σ = 0.01, symmetric decorrelation retraction. No distributional assumptions are imposed; the algorithm operates on the 5-channel whitened signal. H.3
ARTIFACT IDENTIFICATION AND EVALUATION
Independent components are ranked by excess kurtosis κ = E[(z − µ)4 ]/σ 4 − 3. The component with the highest κ is designated the ocular artifact; blink artifacts are impulsive (super-Gaussian) and are expected to dominate on frontal channels. Reconstruction quality is evaluated as the RMS amplitude reduction in a ±250 ms window around the maximum-amplitude blink event, computed as: ∥xclean ∥RMS RMS reduction = 1 − . (32) ∥xraw ∥RMS The cleaned signal is obtained by zeroing the artifact component in the ICA source space and applying the inverse total unmixing matrix (Weeg Wwhite )−1 .