GP-unmix: Unsupervised Source Separation

• GP-unmix is a class of unsupervised algorithms that fuse gradient projection methods with Gaussian process regularization to achieve efficient source separation.
• These methods are applied in nonlinear hyperspectral unmixing and time-domain audio separation, enforcing physical constraints like nonnegativity and simplex adherence.
• By incorporating scalable optimization techniques such as variational sparse inference and block-term tensor decomposition, GP-unmix delivers state-of-the-art performance and significant speedups.

GP-unmix denotes a class of unsupervised algorithms for source separation—often referred to as unmixing—whose central computational paradigm involves gradient projection methods or Gaussian-process (GP) regularization. These methods are principally found in two application domains: nonlinear hyperspectral unmixing for remote sensing, and time-domain audio source separation for acoustics and signal processing. Additionally, GP-unmix refers to recent accelerated algorithms for block-term tensor decomposition (specifically, LL1 decomposition) in structured hyperspectral unmixing. Across these domains, GP-unmix frameworks balance expressivity of source models (nonlinearities and data priors) with computational scalability, incorporating both statistical inference and constraint projection. This overview systematically describes three representative GP-unmix methodologies as developed in (Altmann et al., 2012, Alvarado et al., 2018), and (Ding et al., 2022).

1. Nonlinear Mixing and Gaussian Process Models in Hyperspectral Unmixing

In hyperspectral image (HSI) analysis, GP-unmix methods address the scenario in which the observed spectrum $\mathbf{x}_n\in\mathbb{R}^L$ at each pixel arises from an unknown, potentially nonlinear function $f$ of abundance vectors:

$$\mathbf{x}_n = f(\mathbf{a}_n) + \boldsymbol{\epsilon}_n,\quad \boldsymbol{\epsilon}_n \sim \mathcal{N}(\mathbf{0}, \sigma^2 I_L).$$

The abundance vector $\mathbf{a}_n$ is required to obey the simplex (nonnegativity and sum-to-one) constraint: $a_{n,r} \geq 0,\quad \sum_{r=1}^R a_{n,r} = 1.$ The nonlinear mapping $f$ is modeled as a multivariate GP, parameterized via a bilinear feature expansion $\psi(\mathbf{a})$ of dimension $D = R(R+1)/2$, resulting in the kernel

$$k(\mathbf{a}_i, \mathbf{a}_j) = \psi(\mathbf{a}_i)^\top \psi(\mathbf{a}_j).$$

Joint Bayesian inference is performed over the latent abundances, the GP function $f$, and the noise hyperparameters. The full (marginal) likelihood of the data is

$$p(X|A, \sigma^2) = \prod_{\ell=1}^L \mathcal{N}(y_{\ell} | 0, K + \sigma^2 I_N),$$

where each $y_\ell$ is an $N$-vector of $\ell$-th band values and $K$ is the Gram matrix over all abundances.

Physical constraints are enforced via a locally linear embedding (LLE) prior on the latent arrangement, yielding: $$p(X) \propto \exp\left[ -\frac{\gamma}{2} \sum_{i=1}^N \left\| x_i - \sum_{j \in \nu_i} \lambda_{ij} x_j \right\|^2 \right] \prod_{n=1}^N \mathbb{1}_\Delta(x_n).$$

MAP estimation over this highly structured posterior employs a scaled conjugate-gradient (SCG) method, exploiting the Woodbury matrix identity to ensure computational complexity scales as $O(ND^2 + D^3)$ with typically $D \ll N, L$. The final abundance solutions are projected into the positive simplex via a minimal-volume simplex fit, and endmember spectra are recovered using GP regression conditioned on the optimized latent variables.

2. GP-unmix in Time-Domain Audio Source Separation

For single-channel audio source separation, GP-unmix operates wholly in the time domain, avoiding the phase reconstruction ambiguities of time–frequency approaches: $$y(t)=f(t)+\epsilon(t),\quad f(t) = \sum_{j=1}^J s_j(t),\quad \epsilon(t)\sim\mathcal{N}(0,\nu^2).$$ Each latent source $s_j(t)$ is modeled as an independent GP with zero-mean and a spectral mixture (SM) kernel reflecting the temporal spectrum of the source: $$k_{\mathrm{sm}}(\tau) = \sum_{q=1}^Q w_q\, \exp\!\bigl(-2\pi^2\,v_q\,\tau^2\bigr) \cos\!\bigl(2\pi\,\mu_q\,\tau\bigr).$$ The sum $f(t)$ is thus a GP with covariance $\sum_{j=1}^J k_j(t,t')$.

To reduce cubic complexity in the sequence length $N$, the model employs a variational sparse GP framework, introducing $M \ll N$ inducing variables. The evidence lower bound (ELBO) is optimized per time frame (e.g., 125 ms segments), leading to per-frame complexity $\mathcal{O}(\hat{n}\,M^2)$ and total $\mathcal{O}(N M^2)$. Kernel parameters can be pre-trained on isolated recordings for improved spectral priors.

Frame-wise GP posteriors are recomposed to yield separated, phase-correct time-domain source estimates. Quantitative evaluation demonstrates improved source-to-distortion ratio (SDR, up to 24.1 dB) and source-to-interference ratio (SIR, up to 31.4 dB) relative to baseline methods (NMF variants, tensor factorization).

3. Block-Term Tensor Decomposition via Gradient Projection

In the linear mixing regime for HSI, GP-unmix also refers to an efficient algorithm for LL1 block-term tensor decomposition, formalized as: $$\underline{\mathbf{Y}} = \sum_{r=1}^R \mathbf{S}_r \circ \mathbf{c}_r,$$ where $\mathbf{S}_r\in \mathbb{R}^{I\times J}$ (abundance map, rank $\le L$), and $\mathbf{c}_r \in \mathbb{R}^K$ (endmember spectrum). Aggregating all abundance maps as $\mathbf{S} \in \mathbb{R}^{R \times IJ}$, the objective is

$$\min_{\mathbf{C}, \mathbf{S}} \;\; \frac{1}{2} \|\mathbf{Y} - \mathbf{C}\mathbf{S}\|_F^2 + \sum_{r=1}^R \theta_r \varphi(\mathbf{S}_r)$$

subject to nonnegativity, per-pixel simplex, and per-map rank constraints: $$\begin{align*} \mathbf{C} &\ge 0, \;\; \mathbf{S} \ge 0, \;\; \mathbf{1}^\top \mathbf{S} = \mathbf{1}^\top,\ \text{rank}(\mathbf{S}_r) &\le L, \;\; r=1, \ldots, R. \ \end{align*}$$ The iterative GP-unmix (also termed GradPAPA‐LL1) alternates between projected gradient steps on $\mathbf{C}$ and $\mathbf{S}$, followed by explicit projections:

• $\mathbf{C}$: Clamp at zero following a gradient step.
• $\mathbf{S}$: Alternate projections onto the positive simplex (per-pixel abundance at each location), and per-map low-rank constraints (via SVD truncation) or nuclear norm surrogates.

A spatial regularizer (e.g., smoothed $\ell_q$-TV) $\varphi(\mathbf{S}_r)$ enforces abundance map smoothness, and all feasibility constraints are maintained by analytic projection, ensuring $100\%$ feasible iterates. The method achieves per-iteration complexity $\mathcal{O}(IJKR + m[IJR\log R + IJL])$ and sublinear $O(1/t)$ convergence rate.

4. Optimization, Algorithmic Outline, and Scalability

Across domains, GP-unmix methods share a reliance on scalable optimization. For nonlinear hyperspectral unmixing (Altmann et al., 2012), SCG with analytic gradients is employed for the joint GP/abundance posterior, leveraging low-rank structure and PCA subspace priors on the endmember spectra. For time-domain audio separation (Alvarado et al., 2018), blockwise (per-frame) variational inference breaks a potentially intractable global Bayes problem into parallelizable, GPU-friendly GP subproblems.

In LL1 tensor unmixing (Ding et al., 2022), the two-block gradient projection scheme bypasses the high per-iteration cost and slow convergence endemic to three-factor ALS–MU LL1 algorithms, leading to an order-of-magnitude wall-time speedup (e.g., from $1-2$ hours to $5-6$ minutes on a $500\times 307\times 166$ scene, $R=5$). Fast closed-form projection solvers and alternating projections enable the exact imposition of hard physical constraints.

5. Experimental Evaluation and Quantitative Outcomes

Empirical performance summaries demonstrate the effectiveness of GP-unmix variants:

Application	Domain	Metrics/Results
Nonlinear HSI unmixing	Imaging	GP-unmix outperforms PCA/linear/nonlinear baselines by ARE and SAM, notably w/o pure pixels (Altmann et al., 2012)
Audio source separation	Acoustics	SDR up to 24.1 dB, SIR up to 31.4 dB, 98% faster than full GP; state-of-the-art against NMF/tensor baselines (Alvarado et al., 2018)
Block-term tensor unmixing	Imaging	10–100× speedup, endmember MSE $10^{-3}$–$10^{-2}$, 100% feasibility vs. $\sim$10% in three-factor ALS–MU (Ding et al., 2022)

Key observations include the capability of GP-unmix to recover endmembers when pure pixels are absent, robust abundance estimation under moderate noise, and the enforceability of exact physical constraints (nonnegativity, simplex, low-rank) across all solutions.

6. Physical Priors, Limitations, and Extension Prospects

Hard projection onto the simplex and nonnegativity constraints ensures that GP-unmix solutions remain physically interpretable. Physics-motivated priors—local isometry (LLE), spatial smoothness ($\ell_q$-TV), spectral kernel priors from isolated source recordings—enhance generalization and identifiability. However, GP-unmix methods may depend on pre-trained kernels (audio), heuristic hyperparameter tuning (number/mode of inducing points or spectral mixture components), and frame-based processing that may not fully capture long-range dependencies.

Possible extensions cited include end-to-end kernel and waveform learning, adaptive time-varying or nonstationary kernels, and the integration of deep kernel or neural feature maps to further enhance the modeling power.

7. Summary and Significance

GP-unmix methods realize a fusion between expressive probabilistic modeling (Gaussian processes, block-term tensor factorization) and rigorous projection-based optimization, delivering accurate, physically plausible source separation in both hyperspectral and audio data. Their combination of flexible nonlinear modeling with computational efficiency and constraint enforcement underlies their success in unmixing tasks where classic linear or unconstrained methods falter. Empirical evidence shows superior accuracy, feasibility, and speed relative to standard alternatives, attesting to their suitability as state-of-the-art approaches in the respective domains.