Scalable Approximate Algorithms for Optimal Transport Linear Models

Optimal Transport (OT) has become a fundamental tool for measuring distribution differences, with applications spanning image processing, machine learning, and spectral analysis. Cuturi (2013) introduced entropy regularization, enabling the Sinkhorn algorithm to efficiently compute OT distances at scale. Subsequent research expanded OT's scope to non-probability measures, unbalanced OT, and transfer learning. Despite these advances, applying OT to large-scale linear regression remains computationally challenging due to high complexity of classical solvers like LP and QP, especially with high-dimensional data. Recent efforts focus on scalable algorithms—Sinkhorn iterations, gradient-based methods, and proximal algorithms—but often lack flexibility to incorporate diverse penalties and loss functions. This paper addresses this gap by proposing a unified, efficient framework that extends Sinkhorn's scalability to generalized linear models with convex penalties, supporting large datasets in spectral, mass spectrometry, and transfer learning applications.

The core challenge is developing an efficient, scalable algorithm for solving entropy-regularized OT problems within large-scale non-negative linear regression models. Traditional solvers like LP and QP become computationally infeasible as data size grows, often exceeding memory and time limits. Existing scalable OT algorithms are limited to probability distributions or simple cost structures, restricting their applicability to complex regression tasks involving diverse penalties (L1, L2, Elastic Net) and loss functions (squared error, total variation, Poisson). Additionally, balancing computational efficiency with theoretical convergence guarantees remains unresolved. Addressing these issues is crucial for deploying OT-based models in real-world high-dimensional problems such as spectral unmixing, mass spectrometry, and transfer learning, where data sizes can reach hundreds of thousands of samples.

This work introduces a novel scalable algorithm that combines Sinkhorn-like scaling with Majorization-Minimization (MM) techniques to support a broad class of convex penalties and loss functions in OT linear models. Key innovations include: 1) extending the Sinkhorn algorithm to handle non-negative linear regression with convex penalties, 2) deriving proximal operators for squared error and Poisson likelihood losses within the OT framework, 3) supporting both balanced and unbalanced OT formulations, and 4) ensuring linear computational complexity with convergence guarantees. The framework allows efficient parallelization and straightforward implementation, enabling large-scale problems to be solved with high accuracy and speed, surpassing traditional LP/QP solvers in both performance and scalability.

• �� Formulate the non-negative linear regression as an OT problem with source marginal as a non-negative linear combination of dictionary basis X and weights w.
• �� Use Sinkhorn-like iterative updates for the transport plan Q, alternating between updating the marginal distributions and the weights.
• �� Incorporate convex penalties (L2, L1, Elastic Net) via proximal operators, designed through Majorization-Minimization (MM) schemes, ensuring each iteration reduces the objective.
• �� Derive specific proximal operators for data fitting functions like squared error, total variation, and negative Poisson likelihood, enabling flexible loss modeling.
• �� Extend the scaling algorithm to handle unbalanced OT by relaxing marginal constraints with convex penalty functions, maintaining convergence.
• �� Implement the algorithm with parallelizable multiplicative updates, ensuring scalability to large datasets.
• �� Theoretical analysis guarantees convergence under convexity assumptions, with empirical validation on spectral and mass spectrometry data.

The experiments involve synthetic spectral datasets with sample sizes from 10^2 to 10^5, comparing the proposed algorithm against LP/QP solvers in terms of runtime and accuracy. Regularization parameters (ϵ, α, λ) are tuned via cross-validation. The datasets include Gaussian basis functions and skewed target distributions, simulating real spectral data. Performance metrics include RMSE, convergence speed, and robustness under noise. Ablation studies analyze the impact of penalty choices and regularization strength. Additional tests on spectral unmixing and mass spectrometry data demonstrate the method's practical utility, with results showing significant speedups and error reductions, validating scalability and robustness.

The proposed algorithm achieves over tenfold faster convergence than classical solvers for problems with 10^5 samples, with RMSE below 0.05. It effectively supports L2, L1, and Elastic Net penalties, maintaining stability across different regularization strengths. Incorporating negative Poisson likelihood improves spectral fit accuracy, reducing errors by about 30%. The method scales linearly with problem size, handling high-dimensional spectral data efficiently. Results confirm robustness against noise and bias, outperforming existing OT models in both speed and precision, especially in unbalanced scenarios where traditional methods struggle.

The framework is directly applicable to spectral unmixing, mass spectrometry, and high-dimensional transfer learning, where large datasets and diverse penalties are common. It enables real-time spectral analysis, robust signal separation, and domain adaptation in complex environments. The approach can be integrated into existing spectral analysis pipelines, providing scalable, accurate, and interpretable models. Long-term, it supports the development of intelligent systems for environmental monitoring, medical diagnostics, and industrial quality control, where rapid, large-scale data processing is essential.

The method's performance depends on careful tuning of the entropy regularization parameter, which can be challenging in practice. Its robustness diminishes under extreme imbalance or high noise levels, requiring further enhancements. Currently, the framework primarily supports convex penalties; extending to non-convex or non-smooth penalties remains an open challenge. Additionally, while scalable, the method's efficiency may degrade with very high-dimensional feature spaces or highly complex loss functions, necessitating further algorithmic refinements.