Paper deep dive
Robust blind unmixing: A geometric approach to overcoming basis variation
Dumitru Mirauta, Vladimir V. Gusev, Michael W. Gaultois, Matthew J. Rosseinsky, Yannis Goulermas
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 93%
Last extracted: 10/6/2026, 10:44:09 AM
Summary
The paper proposes a novel geometric approach to robust blind unmixing that addresses basis variation in signal separation problems. Unlike traditional parametric methods, this approach utilizes optimal transport distances, specifically Wasserstein distances, to model variations in basis patterns within a metric space. The method minimizes in-class variance of candidate solutions and is demonstrated on Gaussian mixture models, simulated powder X-ray diffraction, and hyperspectral imaging datasets.
Entities (7)
Relation Signals (6)
Basis Variation → addressedby → Robust blind unmixing
confidence 95% · One of the core challenges of this task is the variation of the basis... We propose a novel geometric approach to unmixing problems
Robust blind unmixing → uses → Optimal Transport
confidence 95% · We take advantage of the optimal transport distances to capture commonly occurring basis variations
Robust blind unmixing → utilizes → Wasserstein Distance
confidence 95% · We take advantage of the optimal transport distances... Our approach to unmixing is novel in... Allowing for non-Euclidean perturbation
Robust blind unmixing → appliedto → Powder X-ray Diffraction
confidence 90% · The effectiveness of our approach is demonstrated on a range of unmixing tasks using ... simulated powder X-ray diffraction
Robust blind unmixing → appliedto → Hyperspectral Imaging
confidence 90% · The effectiveness of our approach is demonstrated on a range of unmixing tasks using ... laboratory hyperspectral imaging datasets.
Robust blind unmixing → comparedwith → Non-negative Matrix Factorisation
confidence 85% · In Section 5 we compare the effectiveness of our method vs Non-negative Matrix Factorisation (NMF)
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Signal separation problems are common in science. A prominent example of this occurs during the use of diffraction or spectroscopy to identify the individual components of a mixture by measuring it. In the simplest case, the measured signal is a linear combination of basis patterns corresponding to the constituent parts. The unmixing problem is to infer all or some of these basis patterns and abundances of components from measurements of distinct mixtures. One of the core challenges of this task is the variation of the basis from mixture to mixture due to noise and the exact physics of the measurement process. This is usually addressed with tailored model-based and parametric methods that are then limited in use to specific application domains by the nature of the assumptions made. We propose a novel geometric approach to unmixing problems which views the generation of data during measurement through a metric space lens, thereby shifting the focus from parametrised models to a general relationship between basis transformations and the corresponding geometry. We take advantage of the optimal transport distances to capture commonly occurring basis variations, and use minimisation of in-class variance of candidate solutions to drive the optimisation. We pay special attention to the one-dimensional case due to its practical importance and availability of efficient distance and transport map routines. The effectiveness of our approach is demonstrated on a range of unmixing tasks using random Gaussian mixture models, simulated powder X-ray diffraction, and laboratory hyperspectral imaging datasets.
Tags
Links
- Source: https://arxiv.org/abs/2610.04091v1
- Canonical: https://arxiv.org/abs/2610.04091v1
Trouble viewing inline? Open PDF directly →
Full Text
123,118 characters extracted from source content.
Expand or collapse full text
Robust blind unmixing: A geometric approach to overcoming basis variation Dumitru Mirauta Affiliation: School of Computer Science and Informatics, University of Liverpool, Ashton Street, Liverpool, L69 3BX, Merseyside, United Kingdom Vladimir V. Gusev †thanks: Corresponding author: vladimir.gusev@liverpool.ac.uk. Affiliation: School of Computer Science and Informatics, University of Liverpool, Ashton Street, Liverpool, L69 3BX, Merseyside, United Kingdom Michael W. Gaultois Affiliation: Materials Innovation Factory, University of Liverpool, Oxford Street, Liverpool, L7 3NY, Merseyside, United Kingdom Matthew J. Rosseinsky Affiliation: Materials Innovation Factory, University of Liverpool, Oxford Street, Liverpool, L7 3NY, Merseyside, United Kingdom Yannis Goulermas †thanks: This author has passed away. Affiliation: School of Computer Science and Informatics, University of Liverpool, Ashton Street, Liverpool, L69 3BX, Merseyside, United Kingdom Abstract Signal separation problems are common in science. A prominent example of this occurs during the use of diffraction or spectroscopy to identify the individual components of a mixture by measuring it. In the simplest case, the measured signal is a linear combination of basis patterns corresponding to the constituent parts. The unmixing problem is to infer all or some of these basis patterns and abundances of components from measurements of distinct mixtures. One of the core challenges of this task is the variation of the basis from mixture to mixture due to noise and the exact physics of the measurement process. This is usually addressed with tailored model-based and parametric methods that are then limited in use to specific application domains by the nature of the assumptions made. We propose a novel geometric approach to unmixing problems which views the generation of data during measurement through a metric space lens, thereby shifting the focus from parametrised models to a general relationship between basis transformations and the corresponding geometry. We take advantage of the optimal transport distances to capture commonly occurring basis variations, and use minimisation of in-class variance of candidate solutions to drive the optimisation. We pay special attention to the one-dimensional case due to its practical importance and availability of efficient distance and transport map routines. The effectiveness of our approach is demonstrated on a range of unmixing tasks using random Gaussian mixture models, simulated powder X-ray diffraction, and laboratory hyperspectral imaging datasets. Keywords: Unmixing, Basis variation, Optimal transport, Wasserstein distance, Earth movers distance, Error-in-variables 1 Introduction The objects and processes of investigation in many scientific disciplines require advanced characterisation, often with multiple instrumental techniques to examine different aspects. In many cases, measurements generate a nonnegative response across some domain (e.g., across a range of angles, frequencies, energies, etc.). Examples of such responses include Powder X-ray Diffraction (PXRD) patterns, Hyperspectral Imaging (HI) pixels, and Raman spectra. For techniques producing real valued but periodic traces (e.g. infra-red spectroscopy, ion cyclotron resonance mass spectrometry, or even of audio waves) frequency distributions may be considered. These response profiles can be useful as a pseudoidentifier of the underlying object, where a similar profile is likely to refer to a similar—though not necessarily the same—object. A response may be characterised in itself, by a “pattern” of highs and lows, large and small features in different locations. As such, in this work, these response profiles will generally be referred to as patterns. For identification purposes, the magnitude of features relative to one another is more instructive than the absolute magnitude of any individual feature. As such, we will treat patterns as normalised distributions. Crucially, the investigation of what may be considered a single underlying object, may produce a range of patterns. For instance, the PXRD pattern of a “pure phase” is composed of a series of peaks; these peaks can individually change position or shape, due to small variation in crystal structure induced as external parameters—such as temperature or pressure—change. Astronomical HI pixels may exhibit Doppler shift as relative speed to an object changes, or Earth HI may exhibit transformations due to atmospheric scattering, subsurface scattering, and illumination. In audio, the same structure (e.g. chord or overtone series of a voice) can change pitch and will have an associated (log) frequency shift. As such, multiple distinct points in pattern space can represent the same basis pattern, we call this basis variation. Further, it is not always feasible to probe an object of interest alone. For instance, a single PXRD sample may contain responses from many “pure phases”. In such cases, the available response profiles will be a combination of responses from multiple objects. In order for the pattern(s) of interest to be isolated, it/they must first be separated. The set of all distinct patterns that can play a role in mixing will be referred to as the basis. In linear mixing, there is no basis variation, and mixtures are produced through linear combination of the basis. Though linear combination is a common feature in many types of mixed data, across disciplines, so is basis variation. Here, it is variations of the basis patterns that are combined. With this, a set of patterns corresponding to the same set of objects, combined in the same proportions, may produce two different mixtures. As such, in this work a mixing model is defined by providing both a basis variation mechanism and a combination mechanism. This flexible umbrella model covers data encountered in materials discovery Ermon et al. (2015); Bai et al. (2018a); CHERNI et al. (2020), Remote sensing Borsoi et al. (2021); Heylen et al. (2014); Ma et al. (2014); Aharoni et al. (2022), Biology Hosseini-Asl et al. (2015); Yan et al. (2019); Glatz et al. (2011), Astrophysics Boulais, A. et al. (2021); Rameau, J. et al. (2021); Nakhostin et al. (2016); Liu et al. (2018), Food testing Holland et al. (1998), Music transcription Smaragdis (2004); Flamary et al. (2016) and other forms of audio separation (e.g. cocktail party problem). Image segmentation in computer vision Ho et al. (2017), multiple curve clustering/alignment Gaffney and Smyth (2004) and distribution template estimation Boissard et al. (2015) are also related. While mixing produces a set of combined patterns from a basis, blind unmixing aims to recover a basis solely from a given set of combined patterns. While there is of course information loss in forming a single mixture, we may be able to patch together a good estimate of a basis from many mixed examples. Here domain-wise (cross-position) variation works in our favor, as different parts of a basis may be visible or obstructed in different mixtures depending on how the bases have transformed. Blind unmixing can of course be seen as an Inverse Problem (IP) Tarantola (2005); Vito et al. (2005), where the basis is a model parameter which essentially fixes a mixing model, which otherwise depends on per mixture variables. Otherwise, nonblind unmixing seeks to perform separation of the constituents in each mixture, according to a given basis. With sufficiently complex basis variation, there may be many ways of generating a given mixed example from a given basis, making even a nonblind “inversion” rather difficult. Consider the importance of tackling the aforementioned basis variation in materials discovery; as sufficiently flexible unmixing methods are not available, much expert time can be directed towards manual analysis, instead of being able to focus on interpretation of higher level data and further experiment planning. The impact from the reduction of such a bottleneck may be felt far, as new materials enable new technologies and techniques for many other fields, including towards the improvement of the very hardware that could be helping accelerate the materials discovery process. A mixing method that is not too specific to any single kind of characterisation technique is particularly beneficial here, as many kinds of characterisation techniques are in use. We classify previous solutions to this problem as either (parametric) model fitting based, or geometric. Fitting based methods of course require the specification of a model to be fitted, however, it may be difficult to capture the multitude of transformation mechanisms required across the application domains, in a single, reasonably constrained model. A geometric method instead views mixing and its inverse problem geometrically. For instance, the linear unmixing case translates to a problem of recognising vertices in a point cloud that is expected to span a simplex. As many families of parametric mixing model can produce geometries with the same quality, they could all be covered by one geometric approach targeting this shared quality. The more general form of the algorithm we suggest operates in an arbitrarily chosen metric space (,d)(P,d_P). If the spectra of interest exhibit transformations that are Lipschitz in this space, then this space has a key property. Namely, that up to some transformation parameter limit (which depends on the distinctiveness of the true basis), we show that the basis variations form clean, separable clusters (see Appendix A). This expected regularity can be exploited when searching for a way to pull apart our mixtures such that the resultant components correspond to plausible basis variants. The set of patterns (or representations thereof) P is a choice as well, for instance 1D peak patterns might be expressed through a space of vectors, of Gaussian Mixture Model (GMM) parameters, of abstract measures, etc. This in turn may affect the kind of distance that can be defined between patterns. After this our focus shifts to the choice of Wasserstein space, and generating a practical algorithm that can execute in time scales feasible for interactive laboratory exploration. Wasserstein space can represent distribution-like data that is continuous or discrete, with weights on uniform or non-uniform grids, when grids are shared between discrete distributions or when not. Wasserstein distances can compare any combination of the above, and deal with transformations that would not be continuous if we were required to continue respecting a shared uniform grid for instance, such as when shifting data provided on a uniform grid. Wasserstein distances are Lipschitz to shifting, broadening and more Engquist et al. (2016). Overall, our approach to unmixing is novel in: 1. Allowing for non-Euclidean perturbation rather than explicitly fitting a fixed basis variation model. 2. Focusing on how a given mixture can be broken apart rather than on how a similar mixture might be put together. 3. Using the geometry of mixing to inform plausible decompositions, rather than a regressive cost. 4. Allowing for sparse and dynamic descriptors through signatures, rather than fixed grids, which can more accurately represent Wasserstein centroids and enable efficient separation. This paper is structured as follows. In Section 2.1 we present our unmbrella model for mixing problems, establishing what is shared among them, and what can arbitrarily change. Section 2.2 gives an overview of Wasserstein distances and their place within Optimal Transport (OT) distances. In Section 2.3 we discuss existing fitting based and geometric approaches to unmixing. Section 3 contains our main contributions, our general approach, followed by the Wasserstein specialisation, and a practical implementation for 1D. A fast implementation hinges on efficient calculation of both Wasserstein distances, and crucially, transport plans, as such Section 4 discusses efficient computaton of these. In Section 5 we compare the effectiveness of our method vs Non-negative Matrix Factorisation (NMF) and Overlapping NMF (oNMF), when dealing with data with basis variation in three contexts; GMM, PXRD and Hyperspectral data. Finally in Section 6 we will conclude and discuss how the approach could be further extended. 2 Preliminaries Lower indices will correspond to vector indices while upper indices may sometimes enumerate non-vector quantities, such as discrete measures. By X⊆X we will implicitly mean measurable subsets. Glossary Acronyms AFD Agile Factor Decomposition ED Euclidean Distance EMD Earth Movers Distance GC Ground Cost GMM Gaussian Mixture Model GS Ground Space HI Hyperspectral Imaging IP Inverse Problem NMF Non-negative Matrix Factorisation NWCA North-west Corner Algorithm oNMF Overlapping NMF OT Optimal Transport PCA Principal Component Analysis PXRD Powder X-ray Diffraction RMSE Root Mean Square Error 2.1 Formal definition of mixing and unmixing An unmixing problem seeks to invert the outputs of a mixing process for the corresponding inputs, as such, let us first define mixing. Let P be the space of allowed (pure or mixed) patterns, such as =ℝNP=R^N as in many cases, or the space of distributions over the real line =P(ℝ)P=P(R) as we will later consider. The basis B=(b1,…,bK)∈K B=(b^1,…,b^K) ^K is a tuple of such patterns, representing prototypes to be drawn from in the mixing process. A combination function κ:ℝK→κ:R^K will describe the mechanism which prototypes may be combined through, such that each mixture is a function of its mixing proportions (here we take the basis to be fixed). Figure 1: A mixing model operating in the space of pattern P, incorporating independent variation of basis patterns prior to combination. Left: Basis patterns bib^i. Middle: The varied basis copies b~ij b^ij that correspond to the mixture cjc^j. Right: The resulting mixture cjc^j. Transitions from left to middle and from middle to right are provided through the variation ϕ and combination κ functions respectively (illustrated in red and blue respectively). The simplest and most common example is linear superposition κlin(,B)=∑i=1Kaibiκ^lin(a; B)= _i=1^Ka_ib^i. In the case of P being a space of distributions, mixing fractions will need to belong to a simplex Sd−1=∈ℝ+d:‖1=1S^d-1=\x ^d_+:\|x\|_1=1\, for valid combinations, and so κ:SK−1→κ:S^K-1 . Other examples include bilinear mixing κ(,B)=∑m,namnbm⊙bnκ(A; B)= _m,na_mnb^m b^n in some hyperspectral unmixing methods (where =ℝNP=R^N and ⊙ is the Hadamard product) Heylen et al. (2014) or weighted Wasserstein centroids as in Schmitz et al. (2018). While in the linear case a mixture might be constructed as cj=κlin(j,B~j)c^j=κ^lin(a^j; B^j), in the presence of independent variation of each basis pattern, we could interpret each mixture to be constructed according to its own, subtly different basis cj=κlin(j,B~j)c^j=κ^lin(a^j; B^j), with B~j=(b~1j,…,b~Kj) B^j=( b^1j,…, b^Kj). The basis variation model, if parameterised, can be captured by a function ϕ:Λ→ϕ:Λ in the variation parameters ∈Λ λ∈ Λ. For example, letting P be continuous functions over the real line, we could have shifting and broadening b~=ϕ(,b)=b(λ2⋅(x−λ1)) b=ϕ( λ;b)=b(λ_2·(x-λ_1)). For brevity, let this also act elementwise over a tuple of patterns, such that we may write B~j=ϕ(j,B) B^j=ϕ( Λ^j; B) where j Λ^j collects the individual parameters of each basis for mixture j. In full, we say that a mixture/combined pattern cjc^j arises according to the mixing model ψ(j,j,B)=κ(j,ϕ(j,B))ψ(a^j, Λ^j; B)=κ(a^j;ϕ( Λ^j; B)), generating the dataset C=cj:j∈[M] C=\c^j:j∈[M]\. The blind unmixing problem seeks to reverse this process, generating the basis B B from the dataset C C alone. In the framework of IPs Tarantola (2005); Vito et al. (2005), B B could be seen to parameterise a family of mixing models ψB:SK−1×ΛK→ψ_ B:S^K-1× Λ^K , in this context recovering the correct B B is equivalent to recovering the correct mixing model. 2.2 Optimal transport primer We give a brief overview of a few concepts in OT, as they form a basis for the intuition, analysis and the practical implementation of our approach. It will be useful to highlight the overlap between the Monge and Kantorovich models, discrete and continuous, as well as the Earth Movers Distance (EMD) and Wasserstein costs or their minimising maps. When equivalent, it can be instructive to switch between these different perspectives. For a deeper overview see Peyré and Cuturi (2019); Santambrogio (2015); Villani (2009). Figure 2: A distance between two distributions μ and ν can often be written as a sum of distances over parts of μ and ν that it takes to correspond. Here, DTW is dynamic time warping. Take μ to be an initial pattern (blue) and ν a transformed (orange) counterpart, with corresponding points (red markers) linked by dotted lines. In such a sum, change in one part is measured independently of change in another part, as if a point along μ had moved along its dotted line to become the corresponding point in ν. The problem of OT is intuitively expressed in terms of the earth movers analogy. In this we have a mound of dirt, the shape of which is described by one distribution, and a hole described by another. We seek a set of movements from the former to the latter which would expend the least overall energy11 1 The energy/work done in this analogy is proportional to both the “mass” of some piece of the mound and to the distance it needs to slide over.. While the ‘transport’ in this analogy is literal, transport in the general problem should be understood more as association or matching22 2 And can in fact be thought of as a generalised matching problem.. More formally, the OT problem aims to find minimum cost associations between two distributions (defined by measures) μ∈P()μ∈ P(X) and ν∈P()ν∈ P(Y) over Ground Spaces (GSs) ,X,Y. Given a GS cost c:×→ℝ+c:X×Y _+ of transferring mass from x to y, an OT cost can be formulated as the average GS cost when following a particular transport map. Similar costs appear in information theory as measures of distortion Pollard (1982); Lloyd (1982). There are generally two types of maps that are considered, leading to two ways to write this cost. Figure 3: While a Monge map seeks a one to one correspondence, with a Kantorovich plan, mass can be broken up and sent to multiple locations. This splitting of mass can be necessary in the discrete setting, as a larger “stick” (Dirac) may need to be broken into parts and contribute towards multiple others for a mass preserving matching. Definition 1 (Pushforward). The pushforward of μ∈ℳ()μ (X) through T:→T:X , results in a measure from ℳ()M(Y) denoted T#μT_\#μ. This is defined as (T#μ)(Y)=μ(T−1(Y))(T_\#μ)(Y)=μ(T^-1(Y)). That is to say, though T redistributes some mass such that it now is spread over Y rather than over X, the amounts of mass before and after this redistribution should agree. Definition 2 (Monge cost). Provided a Monge map T:→T:X that satisfies mass conservation T#μ=νT_\#μ=ν, the average Ground Cost (GC) cost may be written Mc(T)=∫c(x,T(x))μ(x).M^c(T)= _Xc(x,T(x))dμ(x). (1) Definition 3 (Space of couplings). Space of distributions π with μ,νμ,ν as marginals (μ,ν)=π∈P(×):π=μ,π=ν,U(μ,ν)=\π∈ P(X×Y):π_X=μ,π_Y=ν\, (2) where we denote the marginal π(X)=∫X×dπ(x,y)π_X(X)= _X×Ydπ(x,y) for X⊂X and similarly π(Y)=∫×Ydπ(x,y)π_Y(Y)= _X× Ydπ(x,y) for Y⊂Y . Definition 4 (Kantorovich cost). Provided a Kantorovich plan π∈P(×)π∈ P(X×Y) that satisfies mass conservation π∈(μ,ν)π (μ,ν), the average GC may be written Kc(π)=∫×c(x,y)π(x,y).K^c(π)= _X×Yc(x,y)dπ(x,y). (3) When =Y=X and c(x,y)c(x,y) is a distance in X, the problem can represent literal transport (movement across this common GS), as in an EMD. Here both the Monge and Kantorovich costs can be interpreted either as the average displacement during a rearrangement of μ into ν or equivalently the energy required by this rearrangement. In this case, an OT distance on P()P(X) can be defined by taking the minimum cost mass preserving maps dMc(μ,ν)=minT:T#μ=νMc(T)ordKc(μ,ν)=minπ∈(μ,ν)Kc(π).d^c_M(μ,ν)= _T:T_\#μ=νM^c(T) d^c_K(μ,ν)= _π (μ,ν)K^c(π). (4) The Monge variant can be written in PDE form by considering the change of variables formula for distributions, where T acts as a reparameterisation. The Kantorovich variant is a linear functional (in π), where the dual problem—max∫λ,γλ(x)μ(x)+∫γ(y)ν(y) _λ,γ λ(x)dμ(x)+ γ(y)dν(y) such that c(x,y)≥λ(x)+γ(y)c(x,y)≥ λ(x)+γ(y)—can be also be used. If an optimal Monge map exists, it is equivalent to its Kantorovich counterpart in the sense that π∗=(id,T∗)#μπ^*=(id,T^*)_\#μ, (see Fig. 4) where Kc(π∗)=∫c(x,y)dπ∗(x,y) K_c(π^*)= c(x,y)dπ^*(x,y) =∫c(x,y)d((id,T∗)#μ)(x,y) = c(x,y)d((id,T^*)_\#μ)(x,y) (5) =∫c(id(x),T∗(x))dμ(x)=Mc(T∗). = c(id(x),T^*(x))dμ(x)=M_c(T^*). Or loosely, the mass of π∗π^* around the position (x,T∗(x))(x,T^*(x)) is equal to that of μ around x, and 00 elsewhere, i.e. all of the mass of π∗π^* is on the graph of the function T∗T^*. Figure 4: Overlap in Monge maps T and Kantorovich plans π. The density correponding to bivariate distributions π are represented here through heatmaps, it is only noteworthy here where the density is zero (white) and where it is not. Nonzero density at (x,y)(x,y) signifies transport from x to y, these (x,y)(x,y) are in the “support” of the Kantorovich plan, supp(π∗)supp(π^*). While a Kantorovich plan allows for mass splitting if necessary, it is generally more economical to avoid splitting if possible. Therefore if possible supp(π∗)⊆graph(T∗)supp(π^*) (T^*). Here, the pattern ν was generated by shifting and changing the breadth of the peaks in μ. Valid distances can also be constructed in the format Wpc(μ,ν)=dMcp(μ,ν)porWpc(μ,ν)=dKcp(μ,ν)p,W^c_p(μ,ν)= [p]d^c^p_M(μ,ν) W^c_p(μ,ν)= [p]d^c^p_K(μ,ν), (6) these are called Wasserstein distances. Let the omission of the superscript indicate the usual choice c(x,y)=‖x−y‖pc(x,y)=\|x-y\|_p. 2.2.1 Discrete optimal transport Let the measure μ∈P()μ∈ P(X) have a finite support, so that μ=∑i=1Nμwμiδxμiμ= _i=1^N_μw^i_μδ_x^i_μ, where we may wish to refer solely to its positions Xμ=xμi∈:i∈[Nμ] X_μ=\x^i_μ :i∈[N_μ]\ or its corresponding weights Wμ=wμi:i∈[Nμ] W_μ=\w_μ^i:i∈[N_μ]\. The measure μ will be specified via its signature, a tuple of positions and weights (Xμ,Wμ)( X_μ, W_μ). The subscripts will be necessary later to distinguish the positions and weights of many different measures. Sometimes it may be more convenient to refer instead to the weight vector μ∈ℝ+Nμw_μ ^N_μ_+, which we will usually take to be normalised, i.e. μ∈SNμ−1w_μ∈ S^N_μ-1. Let33 3 Here we chose to treat =X=Y for brevity, but this trivially applies to measures over separate GSs. ν∈P()ν∈ P(X) also be finite with signature (Xν,Wν)( X_ν, W_ν). The transport plan between these can only be nonzero on Xμ×Xν X_μ× X_ν, taking the form π=∑i=1,j=1Nμ,Nνpijδxμi,xνjπ= _i=1,j=1^N_μ,N_νp_ijδ_x^i_μ,x^j_ν, hence Kc(π)=∑i=1,j=1Nμ,Nνpij∫×c(x,x′)dδxμi,xνj(x,x′)=∑i=1,j=1Nμ,Nνc(xμi,xνj)pijK^c(π)= _i=1,j=1^N_μ,N_νp_ij _X×Xc(x,x )dδ_x^i_μ,x^j_ν(x,x )= _i=1,j=1^N_μ,N_νc(x^i_μ,x^j_ν)p_ij (7) and so we arrive at the common form of EMD dKc(μ,ν)=min∈ℝ+Nμ×Nν⟨,⟩s.t.T=μand=ν gatheredd_K^c(μ,ν)= _P ^N_μ× N_ν_+ ,P \\ s.t. ^T1=w_μ 1=w_ν gathered (8) where cij=c(xμi,xνj)c_ij=c(x^i_μ,x^j_ν). Remark 1. Other applications of OT include: a distance between signed measures Engquist et al. (2016); a distance between combinations of elements in any set, as combinations may be treated as distribution, such as in Hargreaves et al. (2020); bootstrapping from intra as opposed to inter GS costs, with Gromov-Wasserstein distances; and more Peyré and Cuturi (2019); Flamary (2019). Use of the Wasserstein distance has been noted to produce a more stable optimisation in GANs Arjovsky et al. (2017). OT has been used in complex IPs, such as the inverse design of glass geometries that produces specific images with its caustics Schwartzburg et al. (2014). 2.3 Existing foundational techniques There is a strong association between unmixing and factor based models, and indeed many solutions in which the factor representation is central. The factor representation of mixing is possible in linear mixing, where the basis is combined linearly, i.e. solely via κlinκ^lin, with no basis variation and taking patterns to be expressed as vectors ⊆ℝNP ^N. This representation usually corresponds to sampling on a shared uniform grid of positions. Here, multiple patterns (and similarly fractions) can be packed into the columns of a matrix. Noting that j=κlin(j,B)=jc^j=κ^lin(a^j; B)=Ba^j, the factor representation follows immediately =(1…M)=C=(Ba^1…Ba^M)=BA, where =(1…k)B=(b^1…b^k) and =(1…M)A=(a^1…a^M). Since linear mixing can be expressed through matrix multiplication (in ℝNR^N), linear unmixing, can here be seen as a matrix factorisation problem. Matrix factorisation techniques such as Principal Component Analysis (PCA), Independent component analysis, and NMF can be therefore be seen to perform the required inversion, in the cases of orthogonal basis, independent basis, and nonnegative factor conditions being required respectively Mu et al. (2012); Agrawal et al. (2022); Lee and Seung (2001). Accordingly, NMF often forms a base for unmixing where nonnegative data is involved. Of course, while the model being fit may be linear, the fitting itself generally is a nonconvex optimisation. The model utilised by oNMF Cichocki et al. (2009), diverges from the ordinary two factor model in NMF, by allowing bases to be varied via discrete convolution44 4 As defined here, it is actually a cross correlation, though the kernel λ can be learned reversed. ϕ(,)=∗S=∑s=1Sλs↓s−1ϕ( λ;b)= λ*_Sb= _s=1^Sλ_sb s-1. By ↓su s we mean the shift-down operation where (↓s)i=ui−s(u s)_i=u_i-s if i−s>1i-s>1 else 00. Correspondingly, ↓sU s represents downward shifting each column vector. This allows for addition of many overlapping shifted copies of each basis in each mixture, hence the name55 5 Sometimes this may also be referred to by the name of convolutional NMF (cNMF), which is a closely related algorithm Smaragdis (2004), though relates to a different composition, of 2D objects.. Similarly to NMF, a solution is sought by minimising a matrix dissimilarity between the mixtures C and their reconstruction according to the corresponding mixing model. By observing that the reconstructed mixture can be expressed as ~j c^j =∑i=1Kaijij∗Si=∑i=1Kaij∑s=1Sλsij(i)↓s−1=∑i=1K∑s=1Sλ¯ijs(i)↓s−1 = _i=1^Ka^j_i λ^ij*_Sb^i= _i=1^Ka^j_i _s=1^Sλ^ij_s(b^i) s-1= _i=1^K _s=1^S λ_ij^s(b^i) s-1 (9) =∑s=1S()↓s−1¯∗js = _s=1^S(B) s-1 Λ^s_*j with matrices (¯s)ij=aijλsij( Λ^s)_ij=a^j_iλ^ij_s, the usual matrix product sum form of the reconstructed mixtures ~=∑s=1S↓s−1¯s C= _s=1^SB s-1 Λ^s is obtained. This sum of products can be written as a single product of concatenated factors, hence this model can be seen through the lens of the linear mixing ~=¯¯ C= B A, where ¯=(↓0…↓S−1) B=(B 0…B S-1) and similarly ¯=(¯0…¯S−1) A=( Λ^0… Λ^S-1). Here the basis may be considered over-complete, holding multiple alternative representations of the same basis pattern. In mixing, only a single representation is desired for each basis pattern and as such sparse fractions are desired. An instance of this is Agile Factor Decomposition (AFD) Bai et al. (2018b), where varied bases are generated by a shift. Though additional constraints would be required so that different elements of the dictionary may correspond to different representations of the same basis, this is not dissimilar to the sparse dictionary learning problem. Accounts of methods more directly applying a sparse dictionary learning approach can be found in Borsoi et al. (2021). The IPs of NMF and oNMF are usually tackled through reconstruction error minimisation, through costs such as ∑jD(ψ(j,j,B),cj) _jD(ψ(a^j, Λ^j; B),c^j) where D is a distance, dissimilarity or divergence, for example, ∑j‖j−j‖22 _j\|Ba^j-c^j\|_2^2. Recall that while all three metric axioms are satisfied by a distance, dissimilarities and divergences may violate the triangle inequality and divergences may also violate symmetry. A reconstructed data matrix ~ C can also be defined, where the jthj^th column is given by the reconstructed mixture ψ(j,j,B)ψ(a^j, Λ^j; B) and with the cost D(~,)D( C,C) being minimised instead, for some matrix divergence, for example ‖~−‖22\| C-C\|_2^2, or column-wise EMD such as in EMD-NMF66 6 Note, while EMD-NMF performs a more robust fit of the linear mixing model, it still only fits a linear/two factor model. Sandler and Lindenbaum (2011). The Wasserstein barycenter fitting approach Schmitz et al. (2018), also forms such a cost. A cost of course needs to act as a notion of proximity to a desirable solution, however, there are two reasons why a cost of the format described above may not do this as well as possible. A good fit of the mixture may not require a good basis. If the bases may be arbitrarily contorted, as is the case in using an OT model for basis transformation, a perfect fit may be achieved with any choice of basis, and so does not serve to guide us to a plausible basis. Figure 5: Distances after offsetting an original density f (a-c, solid blue), producing shifted counterparts (a-c, dashed) for various shifts s, the markers in panels (d) and (e) are colored according to corresponding values of s. Panels (d) and (e) track a pointwise distance (Euclidean) and an OT distance (Wasserstein) respectively, between the original and a shifted copy, as a function of shift. Note, since a pointwise distance only compares features that align vertically (see Fig. 2), we obtain local minima whenever non-corresponding peaks happen to align (such as in the yellow case). Note also the plateau in (d) when no features align (green case), where the gradient is near 00, and provides little direction in optimisation (even after perturbation). Avoidable local minima. As well as local minima corresponding to alternative solutions, a typical reconstruction cost can also be stationary at misaligned copies of the same solution. This is due to the use of pointwise distances and the many locally optimal comparisons that they can produce. Suppose A and B B are fixed. Even if there exists a unique globally optimal ∗j Λ^j_* with ψ(j,∗j,B)=cjψ(a^j, Λ^j_*; B)=c^j, there can still be other j Λ^j that locally minimise Fj(j)=D(ψ(j,j,B),cj)F^j( Λ^j)=D(ψ(a^j, Λ^j; B),c^j). For instance, consider registering a single uniformly shifted basis in a pointwise distance (see Fig. 5 - bottom left). Any combination of local minima of FjF^j is a local minimum of ∑j=1MFj(j) _j=1^MF^j( Λ^j), hence if each FjF^j has L local minima, the sum has LML^M local minima. If we have 100100 mixed examples, and the alignment of each is akin to Fig. 5, then each feasible decomposition has roughly an additional 3.23×10843.23×10^84 misregistered counterparts, entirely due to the choice of distance used in construction of a fitting cost. 2.3.1 Geometric approaches Not only can the image of a mixing model ψ, and hence the point cloud of generated mixtures, form a specific structure within pattern space P, but inputs such as the basis may be reflected through this structure as geometric features. Geometric approaches recognise this and look for such geometric features. We give these two classifications: those which consider the spatial structure arrising from a particular form of the mixing model ψ, which strongly informs which geometric features should be sought, and those which seek some analogue to features present in the linear case, without necessarily considering whether the presence of such analogues will be preserved under the change in context. In the linear case, taking all convex combinations of bib^i simply forms a convex hull77 7 This is a simplex given a number of bases not more than one above the dimensionality of the patterns. of these patterns, this contains all possible mixtures in the normalised case, with the bases at the vertices of the hull (Fig. 6 - left). Supposing a dataset provides a dense enough sampling of the possible mixtures, then to find the basis one needs only identify vertices. Simplex structure forms the foundations for simplex methods for linear mixing Ma et al. (2014). A modified structure may be obtained if one first allows bases to be perturbed prior to combination (Fig. 6 - middle) Stein (2003); Du et al. (2014); Zhou et al. (2018). While this modification considers variation of basis patterns, the distributions considered only assign significant likelihood to basis variation that does not drastically increase Euclidean distance (which can be large for non-pointwise variation). Some approaches seek a simplex embedded within a manifold, for instance through search of vertices via counting graph geodesics Heylen et al. (2014). This offers an analogue for simplex fitting in the non-linear case, however, as basis variation is not treated independently in this model, it would not be able to inspect how each basis has changed individually, in the context of a given mixture. It is also not clear what kinds of mixing might produce such structure. Manifold analogues of PCA have also been explored Fletcher et al. (2004), however a manifold analogue of NMF might be more appropriate for non-negative data. 3 Latent basis variant clustering Figure 6: Geometric views of mixing under three kinds of basis variation (see Fig. 1) and linear combination, with the action of the variation ϕ and combination κ functions denoted with red and blue arrows respectively. Dark grey represents possible variations b~ b of basis pattern b, lighter grey - possible combinations, and lightest - the whole space of patterns; the simplex Sd−1=∈ℝ+d:‖1=1 S^d-1=\x ^d_+:\|x\|_1=1\ or an arbitrary pattern space P. Consider the space of mixing inputs (fractions and basis variations) that combine to produce a given arbitrary pattern p∈p exactly ℐ(p)=(,B~)∈SK−1×K:κ(,B~)=p.I(p)=\(a, B)∈ S^K-1×P^K:κ(a; B)=p\. (10) A decomposition that perfectly reconstructs all of our given mixtures can be obtained by simply selecting one element from each ℐ(cj)I(c^j). Such a selection of b~ij b^ij may have no particular spatial structure (or connection to a shared basis), however under mild conditions (see Appendix A), we expect clustering of basis variants. There may be many solutions in ∏jℐ(cj) _jI(c^j) which appear to be clustered, a simple way to find one would be to optimise a clustering objective such as minimising in class variation. This has the desirable side effect of favouring bases requiring less modification on average to explain the data, which we might say provide a more natural explanation. To measure variation, we will employ the Fréchet variance. The weighted Fréchet variance of a cluster of patterns P~⊂ P is defined as follows FV(P~,)=min∑j=1|P~|p∈wjd(p,p~j)2.FV( P,w)= _p _j=1^| P|w_jd_P(p, p^j)^2. (11) Much like ordinary variance, it measures average squared deviation from the mean (defined by the minimiser), with the modification that deviation is given by d_P rather than the Euclidean distance. Our problem is therefore min(j,b~∗j)∈ℐ(cj)∀j∑i=1KFV(b~i∗,ai∗)=argmin(j,B~j)∈ℐ(cj)∀jbi∈∀i∑i=1K∑j=1Maijd(bi,b~ij)2 _(a^j, b^*j) (c^j)\>∀ j _i=1^KFV( b^i*,a_i^*)= _ subarrayc(a^j, B^j) (c^j)\>∀ j\\ b^i \>∀ i subarray _i=1^K _j=1^Ma_i^jd_P(b^i, b^ij)^2 (12) where we weight according to abundances, as we wish components that play a lesser role in a given mixture, to have less overall influence. This is of course much too large a problem in practice, however, it can be decoupled into many smaller subproblems through alternating optimisation. Averaging - given b~∗∗ b^**, update cluster centroids ∀i∈[K]bi=argminb∈∑j=1Maijd(b,b~ij)2.∀ i∈[K] b^i= _b _j=1^Ma_i^jd_P(b, b^ij)^2. (13) Separation - given b∗b^*, if possible, select feasible mixing inputs that are less spread out from the cluster centers ∀j∈[M]j,B~j=argmin(,B~)∈ℐ(cj)∑i=1Kaid(bi,b~i)2.∀ j∈[M] ^j, B^j= _(a, B) (c^j) _i=1^Ka_id_P(b^i, b^i)^2. (14) Clearly Eq. 12 is non-increasing under these alternating minimisations, as each updates a subsets of the overall sum in a non-increasing fashion. Note, either the cost is strictly decreasing, in which case we have convergence as the cost is bound below by 0, or neither averaging or separation have updated any terms at some iteration, in which case no further updates can be made, also converging. For data such as that produced by PXRD; perturbation mechanisms such as shifting and broadening are of particular importance. While pointwise distances do not behave well with respect to these (see Fig. 5), Wasserstein distances are convex (and therefore locally Lipschitz) to these and more Engquist et al. (2016), and so we can expect clustering of basis variations (see Appendix A). Accordingly, from here we consider the set of all distributions =P()P=P(X) on a geodesic GS (X, d_X). The distance on P will be the Wasserstein distance d=W2d_P=W_2^d_X. 3.1 Averaging Let (S,dS)(S,d_S) be some metric space. For convenience, let the weighted centroid of two points μ,ν∈Sμ,ν∈ S be expressed by BdS(μ,ν,λ)=argminβ∈S[(1−λ)dS(β,μ)2+λdS(β,ν)2],B_d_S(μ,ν,λ)= _β∈ S [(1-λ)d_S(β,μ)^2+λd_S(β,ν)^2 ], (15) this can in fact be obtained by walking a fraction of the way along the geodesic connecting the two points gμν:[0,1]→Sg_μ^ν:[0,1]→ S, or BdS(μ,ν,λ)=gμν(λ)B_d_S(μ,ν,λ)=g_μ^ν(λ) (see Appendix C). In the Euclidean Distance (ED) case, where dS(μ,ν)=‖μ−ν‖2d_S(μ,ν)=\|μ-ν\|_2, the centroids lie on the line (1−λ)μ+(1-λ)μ+λν. Note, this connection between centroids and geodesics comes in handy both in the ground space and in the overarching Wasserstein space. Definition 5 (The general McCann interpolation). Let Tλ(x)=Bd(x,T∗(x),λ)T^λ(x)=B_d_X(x,T^*(x),λ) where T∗T^* is the optimal transport map from μ to ν, then the general McCann interpolant is the pushforward T#λμT^λ_\#μ (for λ∈[0,1]λ∈[0,1]). A Monge map can tell us that mass at x should go to T(x)T(x), but does not specify how this should happen. In the McCann, or displacement interpolation, we continuously (and at constant speed) move mass along the GS geodesic connecting x and T(x)T(x), obtaining a measure somewhere between μ and ν at each instant (see Fig. 7). Since the general McCann interpolant follows the Wasserstein space geodesic (see Appendix D), just as the centroid does (see Appendix C), one can use the McCann interpolant to calculate centroids in Wasserstein space and vice versa BW2d(μ,ν,λ)=T#λμ.B_W_2^d_X(μ,ν,λ)=T^λ_\#μ. (16) When a Monge map does not exist, an approximation—in the sense W2(T¯#μ,ν)≈0W_2( T_\#μ,ν)≈ 0—can be generated via the barycentric projection, T¯(x)=∫ydπ∗(x,y) T(x)= _Xydπ^*(x,y) in the ED GC case. Thanks to the displacement interpolation, there is no need to solve a coupled transport problem between in π1∈(μ,β(λ))π^1 (μ,β(λ)) and π2∈(β(λ),ν)π^2 (β(λ),ν) for each weight λ, knowing just π∈(μ,ν)π (μ,ν) gives us β(λ)β(λ) for all λ∈[0,1]λ∈[0,1]. Algorithm 1 Pairwise barycenter function B¯W2d B_W_2^d_X(μ, ν, λ) ⊳ Approximation due to projection better for lower λ. ←P← deriveOTPdderiveOTP_d_X(μ, ν) ⊳ Choose deriveOTP as appropriate, we consider the North-west Corner Algorithm (NWCA) later for i from 11 to N do ν′i←Bd(Xν,i∗)x^i_ν ← B_d_X( X_ν,p_i*) ⊳ Barycentric projection, ∑j=1Mpijνj/wμi _j=1^Mp_ijx^j_ν/w^i_μ for ED βi←Bd(μi,ν′i,λ)x^i_β← B_d_X(x^i_μ,x^i_ν ,λ) ⊳ Displacement interpolation, λμi+(1−λ)ν′iλx^i_μ+(1-λ)x^i_ν for ED end for return (Xβ,Wμ)( X_β, W_μ) end function When averaging more than two patterns, we approximate barycenters hierarchically, as in Boissard et al. (2015), see Algorithm 2. Figure 7: The displacement interpolation (OT barycenter construction) from three perspectives (one per column). The barycenter captures aspects of multiple patterns by striving to be ‘in-between’ them, attempting to be similar to both. With an appropriate distance, this can combine information in a meaningful way. The barycenter β will be more μ-like or ν-like according to the weighting coefficient λ. Ordinarily a transport map associates some mass at position x in μ to mass at position y in ν (see Fig. 3). With this association, we can consider continuously moving this mass from its location within μ to its location within ν, along the shortest path between x and y. With λ=0λ=0, all mass elements are at their source position and constitute μ. With λ=1λ=1, all mass elements are at their target position and constitute ν. Mass elements in intermediate positions, constitute an intermediate pattern β. In latter columns, a segment of one colour in μ corresponds to that of the same color in β and in ν. The blue curve in the Kantorovich column illustrates the Monge map from μ to β. Algorithm 2 Hierarchical Wasserstein centroid approximation function HWCA(W=(w1,…,wNM) W=(w^1,…,w^N_M), M=(μ1,…,μNM) M=(μ^1,…,μ^N_M)) ⊳ For the sake of brevity, let NM=|M|N_M=| M| be a power of two. ⊳ Take W W, M M to be sorted such that weights are decreasing. while |M|>1| M|>1 do M′←(),W′←() M ←(), W ←() for k from 11 to |M|2 | M|2 do ⊳ Can be run in parallel wk′←w2k−1+w2kw_k ← w_2k-1+w_2k μk′←μ_k ← B¯W2d B_W_2^d_X(μ2k−1,μ2kμ^2k-1,μ^2k, w2k−1/wk′w_2k-1/w_k ) end for M←M′,W←W′ M← M , W← W ⊳ Note, these can be edited in place end while return μ1μ^1 end function 3.2 Separation Recall the notation of Section 2.2.1. Discretising and writing in terms of abundance scaled versions of the transport maps ¯i=aii P^i=a_iP^i between bib^i and b~i b^i, the discrete W2W_2 separation subproblem (see Eq. 14) becomes min∑i=1K¯1,…,¯K∈ℝ+N×N⟨i,¯i⟩s.t.(¯i)T=‖¯i‖1bi∀i∈[K]and∑i=1K¯i=c gathered _ P^1,…, P^K ^N× N_+ _i=1^K ^i, P^i \\ s.t. ( P^i )^T1=\| P^i\|_1w_b^i ∀ i∈[K] _i=1^K P^i1=w_c gathered (17) where the entries of iD^i are given by dkli=d(bik,cl)2d^i_kl=d_X (x^k_b^i,x^l_c )^2 and all measures have been taken to have the same size support for brevity88 8 Uneven measures can be trivially padded with zero weight Diracs.. The combination constraint ∑i=1Kaib~i=c⟹∑i=1Kai∑k=1Nwb~ikδb~ik=∑k=1Nwckδck _i=1^Ka_i b^i=c _i=1^Ka_i _k=1^Nw^k_ b^iδ_x^k_ b^i= _k=1^Nw^k_cδ_x^k_c (18) gives that the varied basis patterns cannot have weight where the mixture does not or supp(b~i)⊂supp(c)supp( b^i) (c), hence we can set b~ik=ckx^k_ b^i=x^k_c. And so, during the separation subproblem, minimisation with respect to weights alone is sufficient. The weight vector b~iw_ b^i is obtained from the marginal ¯i/∥¯i∥1 P^i1/ P^i _1. (a) Coupled multi-source transport problem, with individual transport plans ijP^ij between bases bib^i and variations b~ij b^ij, in each mixed datapoint cjc^j. (b) Single transport plan ~j P^j between mixture cjc^j and concatenated sources c~j=∑j=1Kbj c^j= _j=1^Kb^j. a Figure 8: The segmentation problem as it emerges from the variation cost in the Wasserstein context (a), minimising with respect to multiple coupled plans, and an equivalent, single transport plan/distance evaluation simplification (b). Note, in the former case we inherently discover the abundances a when optimising for unnormalised bjb^j, but in the latter case we must evaluate the distance/plan for each choice of a. However, since we can construct exact plans for the ordinary 1D optimal transport problem easily, now we only have to optimise for a. Concatenating ~ D =(12…K) =(D^1D^2…D^K) (19) ~ P =(¯1¯2…¯K) =( P^1 P^2… P^K) c~T w_ c^T =(a1b1T…akbkT), =(a_1w_b^1^T… a_kw_b^k^T), we can now express the separation as min∈SK−1⟨~,~⟩s.t.~T=c~and~=c gathered _a∈ S^K-1 D, P \\ s.t. P^T1=w_ c P1=w_c gathered (20) where c~()=κlin(,B) c(a)=κ^lin(a; B), as supposed to c=κlin(,ΦB)c=κ^lin(a;Φ B), so we can in fact write min∈SK−1W2d(c~(),c)2 _a∈ S^K-1W_2^d_X( c(a),c)^2 (21) which is a convex optimisation problem that benefits from being able to efficiently evaluate W22W_2^2 and its optimal transport map. For the initial retrieval of abundances, this optimisation is indeed equivalent to W2dW_2^d_X regression in the linear model, such as in Sandler and Lindenbaum (2011), however, we do not stop there (see Algorithm 3), as we are not after a linear representation, we are after the transformed basis patterns inferred from the unconstrained marginals. If the initial distance computation cannot expose the concatenated transport plan, we must use another method to recover it afterwards. We then separate this into individual plans (see Fig. 8). Note the advantage in non-overlapping supports here, where some positions can be associated with a particular basis pattern and therefore transport to/from such positions in the combined map can be identified as transport to/from the corresponding pattern. Algorithm 3 OT separation (unmixing with known basis) function separate(B B, c) ←argmin¯W2d(c~(¯),c)2a← _ aW_2^d_X ( c( a),c )^2 (¯1¯2…¯K)←( P^1 P^2… P^K)← deriveOTPdderiveOTP_d_X(c~() c(a), c) for i∈1,…,Ki∈\1,…,K\ do b~iT←¯i/∥¯i∥1w_ b^i^T← P^i1/ P^i _1 end for end function 3.3 Proposed unmixing algorithm Given a method for efficiently computing transport plans (see Section 4.1 for 1D) and provided that centroids can be computed (or well approximated) in the given GS (in a feasible timescale), such as via a hierarchical approach, we arrive at the general unmixing algorithm Algorithm 4. This algorithm is highly parallelisable99 9 The averaging step is broken in to a series of pairwise reduction passes, these can be run in parallel., as due to the degree of decoupling in the problem, there are many independent subproblems and for large datasets, GPUs can be leveraged to solve a larger subset of these concurrently. We provide a concrete implementation for the (P(ℝ),W2)(P(R),W_2) case. Algorithm 4 Warping robust unmixing for discrete measures function W2unmix(C, NiterN^iter, K) Generate random initial B~ B (random weights over each x∈supp(c)x (c)) for k∈1,…,Niterk∈\1,…,N^iter\ do for i∈1,…,Ki∈\1,…,K\ do ⊳ Averaging step bi←b^i← HWCA(ia_i, (b~i1,…,b~iM)( b^i1,…, b^iM)) end for for c∈Cc∈ C do ⊳ Can be performed independently per mixture separate(B B, c) end for end for return B~,B B, B end function 3.4 Variations of the unmixing problem Depending on what information is originally known, the problem of interest may in fact be more constrained. For instance our separation step provides a solution to the problem of separating components when the basis is known. For partially blind unmixing, the known patterns are included in the initial basis instead of fully random data and can be kept fixed by skipping the corresponding basis update. However, if we wish to detect novel patterns, as one often may in materials discovery, these will not have comparably weighted clusters in our geometry, as these patterns are effectively outliers and are drowned out by the more common bases. As such, we cannot simply use the partially blind approach in this case. A dedicated anomaly detection approach is necessary. Problem Amount of basis known Primary goal Separation (non-blind) All For each mixture, determine parts likely arising from each basis Partially blind unmixing Some Find remaining basis patterns Fully blind unmixing None Find basis B + Anomaly detection - Identify datapoints containing trace uncommon patterns Table 1: Table of variations of the unmixing problem. The anomaly detection relaxation can be optionally made from any of the first three variants, it is treated separately from partially blind unmixing since a sought additional component may be too statistically insignificant to estimate (though may be indirectly detected). The presence of sufficiently novel patterns can still be detected through our usual formulation, via an uncharacteristically large segmentation cost for a particular mixture. This can be seen as attempting to fit a mixture of K+1K+1 distinct basis patterns with only a K subset, will involve additional deformation of these K basis patterns. By constraining the deformation of these basis patterns, to less than what would minimally be required for fully fitting with a K subset, we can attempt to estimate the unknown through the residual, after a best partial fit. Again, residuals need only be enabled for those mixtures where anomalies are reasonably suspected. This changes the separation problem to the following linear program min∑i=1Kπ¯i∈M+(ℝ×ℝ)η¯∈M+(ℝ)⟨d,π¯i⟩+β‖η¯‖1s.t.(projx)#π¯i=‖π¯i‖1bi∀i∈[K],(projy)#∑i=1Kπ¯i+η¯=cand⟨d,π¯i⟩≤Δ2¯ gathered _ subarrayc π^i∈ M_+(R×R)\\ η∈ M_+(R) subarray _i=1^K d, π^i +β\| η\|_1\\ s.t. (proj_x)_\# π^i=\| π^i\|_1b^i ∀ i∈[K], (proj_y)_\# _i=1^K π^i+ η=c\\ d, π^i ≤ Δ^2 gathered (22) where the ‖η¯‖1=η¯(ℝ)=ε\| η\|_1= η(R)=ε penalty discourages unnecessarily attributing mass to η, and Δ2¯ Δ^2 is a specified maximum average squared displacement, set to a value less than the minimum of the unconstrained segmentation cost. In the discrete linear program, we can alternatively (or additionally) enforce a maximum displacement constraint, this can be achieved efficiently by disregarding variables p¯kli p^i_kl that correspond to dkli>Δ2d^i_kl>Δ^2, reducing the amount of variables that need to be solved for. 4 Fast Wasserstein Here we will cover two ways in which to exploit structure in the OT problem, for a fast exact and an approximate distance calculation. 4.1 The north west corner algorithm Definition 6 (Monotone increasing). A transport plan π is monotone increasing if ⟨x′−x,y′−y⟩≥0∀(x,y),(x′,y′)∈supp(π). x -x,y -y ≥ 0 ∀(x,y),(x ,y ) (π). (23) In the 1D case, the monotonicity of transport plans fully constrains them, such that no minimisation is required. Further, for GCs that admit monotone increasing solutions, such as the ED, this unique feasible map can be constructed through the NWCA (in (N)O(N) time). Though this fact is known in OT literature Bai et al. (2023), it does not appear to be always exploited downstream. Lemma 1. Though a Monge map is not sufficient to convey optimal transport between a pair of general measures (such that supp(π∗)⊆graph(T∗)supp(π^*) (T^*)), a parameterised trace (x,y)=(Tx(z),Ty(z))(x,y)=(T_x(z),T_y(z))—with monotone Tx,TyT_x,T_y—is in 1D with the ED GC. Similarly, in the discrete case, a sequence (xμi,xνj)(x^i_μ,x^j_ν) with monotone xμix^i_μ and xνjx^j_ν is sufficient. Proof. In the discrete case we have the following. Take the points in the support of the optimal (Kantorovich) map in any initial order, forming the sequence zi=(xi,yi)∈supp(π∗)z^i=(x^i,y^i) (π^*). Suppose there are no ziz^i with the same x value, then we can order these such that xi\x^i\ are strictly increasing. Since the optimal plan is monotone increasing π∗π^*, (xi+1−xi)(yi+1−yi)≥0(x^i+1-x^i)(y^i+1-y^i)≥ 0 must hold, where xi+1−xi>0x^i+1-x^i>0 by ordering and so yi+1−yi≥0y^i+1-y^i≥ 0, hence yi\y^i\ are non-decreasing. More generally, we may order such that xi\x^i\ are non-decreasing, containing both strictly increasing and constant subsequences. Over strictly increasing subsequences of x, we have non-decreasing y as before. By uniqueness of ziz^i (points in the support), y cannot also be equal for constant x subsequences, so over these, ziz^i can be ordered such that yiy^i are strictly increasing. Finally, note that for z in a preceding strictly increasing sequence, and z′z the next subsequence, x′−x>0x -x>0, hence all y in the preceding subsequence are lower than those in the next. Hence the support of any optimal transport map (for discrete, euclidean gd, 1D) can be expressed through two non-decreasing subsequences over the GS. Similarly, in the limit we require monotone Tx,TyT_x,T_y in the continuous case. ∎ Corollary 1. If xμix^i_μ and xνjx^j_ν are ordered, then it follows from the monotonicity of the support trace that pijpi′j′>0withi≤i′⇒j≤j′.p_ijp_i j >0 i≤ i j≤ j . (24) Definition 7 (Slacks). sμi=wμi−∑j=1Mpijandsνj=wνj−∑i=1Npij, gathereds_μ^i=w^i_μ- _j=1^Mp_ij s_ν^j=w^j_ν- _i=1^Np_ij, gathered (25) where sμi=0s_μ^i=0 signals the ithi^th row constraint being met. Definition 8 (North-west corner rule (NWCR)). Set the north-west corner of the transport plan to p11=min(wμ1,wν1)p_11= (w^1_μ,w^1_ν), so that at least one of sμ1=0s_μ^1=0 or sν1=0s_ν^1=0 is satisfied. Definition 9 (North-west corner algorithm (NWCA)). Having satisfied at least one row or column constraint, we move on to determine the remainder of the transport plan, which is structured (μ,ν)=(p11(μ1,ν1))or(μ,ν)=(p11(μ1,ν1))P(μ,ν)= pmatrixp_11 0\\ P(μ^1,ν^1) pmatrix (μ,ν)= pmatrix matrixp_11\\ 0 matrix&P(μ^1,ν^1) pmatrix (26) where μ1=μ−p11δxμ1μ^1=μ-p_11δ_x^1_μ and ν1=ν−p11δxν1ν^1=ν-p_11δ_x^1_ν. To determine (μ1,ν1)P(μ^1,ν^1), and therefore (μ,ν)P(μ,ν), we iterate this rule until μK()=νK()=0μ^K(X)=ν^K(X)=0. This yields a feasible transport plan (as all slacks are set to 00) in at most N+M−1N+M-1 steps1010 10 Too see why, consider that each step satisfies at least one additional constraint (than the previous). There are a total of N+MN+M constraints, and the final step must satisfy both the last row and column constraints. Alternatively, the longest monotonic trace through the matrix (μ,ν)P(μ,ν) visits N+M−1N+M-1 entries. (see 3.4.2 in Peyré and Cuturi (2019)). Algorithm 5 North-west corner algorithm function NWCA(μw_μ, νw_ν) i1←1,j1←1,sμ←wμ1,sν←wν1,←N×Mi_1← 1,j_1← 1,s_μ←w^1_μ,s_ν←w^1_ν,P 0_N× M, KC=0KC=0 for k=1k=1 to k=N+M−1k=N+M-1 do pikjk←min(sμ,sν)p_i_kj_k← (s_μ,s_ν) sμ←sμ−pikjks_μ← s_μ-p_i_kj_k sν←sν−pikjks_ν← s_ν-p_i_kj_k KC=KC+c(xμik,xνjk)pikjkKC=KC+c(x^i_k_μ,x^j_k_ν)p_i_kj_k if sμ=0s_μ=0 then ik+1←ik+1i_k+1← i_k+1, jk+1←jkj_k+1← j_k if ik+1=Ni_k+1=N then break end if sμ←wμik+1s_μ←w^i_k+1_μ end if if sν=0s_ν=0 then ik+1←iki_k+1← i_k, jk+1←jk+1j_k+1← j_k+1 if jk+1=Mj_k+1=M then break end if sν←wνjk+1s_ν←w^j_k+1_ν end if end for end function (a) Optimal Kantorovich plan. (b) Edges corresponding to pij∗>0p^*_ij>0 do not cross. Figure 9: The North-west corner algorithm (Algorithm 5) and the monotonic trace produced, corresponding to monotonic parameterisations i(k)i(k) and j(k)j(k), or in the continuous case x(t)x(t) and y(t)y(t). As the name suggests, we start in the north-west corner of the transport plan, and at each step k, compare how much mass is left at each marginal index i and j, transfer the most mutually available, and increment the saturated index/indices (i, j, or both). Note, the horizontal stepping of the Kantorovich plan in the middle panels of Fig. 7, how a Monge function cannot do this, and how a function pair (x(t),y(t))(x(t),y(t)) can (see Lemma 1), as one can linger at the same position for some time, e.g. i(k)=3i(k)=3 for k∈[3,6]k∈[3,6]. 4.1.1 Optimality for 1D sorted measures For μ∈P(ℝ)μ∈ P(R), its signature may be sorted so that xμ1<xμ2…<xμN−1<xμNx^1_μ<x^2_μ…<x^N-1_μ<x^N_μ (27) noting that if xμi=xμi′x^i_μ=x^i _μ, we can simplify wμiδxμi+wμjδxμi′=(wμi+wμi′)δxμiw^i_μδ_x^i_μ+w^j_μδ_x^i _μ=(w^i_μ+w^i _μ)δ_x^i_μ. Lemma 2. Taking 1D μ and ν to be sorted without loss of generality, the NWCR is a necessary condition for a monotone plan. Proof. Supposing this is not true, if neither the mass constraint at xμ1x^1_μ or xν1x^1_ν is satisfied by p11p_11 alone then sμ1>0⇒∃pi1>0andsν1>0⇒∃p1j′>0s_μ^1>0\>\> \>\>∃ p_i1>0 s_ν^1>0\>\> \>\>∃ p_1j >0 (28) with i,j′>1i,j >1, but by monotonicity (recall Corollary 1), since p1j′pi1>0p_1j p_i1>0 (where 1<i1<i) we must have j′≤1j ≤ 1, which contradicts 1<j′1<j . ∎ Corollary 2. If (μ,ν)P(μ,ν) is to be monotone, (μk,νk)P(μ^k,ν^k) must be monotone for all k, hence the application of the NWCA is necessary for (μ,ν)P(μ,ν) to be monotone. Theorem 1. In 1D, application of the NWCA on sorted signatures produces an optimal map, provided the GC satisfies c(x,y)+c(x′,y′)≤c(x,y′)+c(x′,y)c(x,y)+c(x ,y )≤ c(x,y )+c(x ,y) for all ⟨x′−x,y′−y⟩≥0 x -x,y -y ≥ 0. Proof. As by Proposition 2.5 of Rieger (2006), there is only one monotone increasing transport plan with marginals μ,νμ,ν in 1D and as we require a monotonic increasing plan for optimality, we have the unique optimal solution. ∎ 4.2 Optimal distance approximation quantisers Though the computation of a complex distance calculation might scale badly in the size of the original patterns compared, we may be able to economically approximate this via a much smaller but nearby quantised counterparts. In other words, there may exist much smaller but similar patterns that are much faster to compare. For instance the ‘signatures’ used in Rubner et al. (2000), generated via k-means clustering from some original distribution (which in practice will also be discrete but larger). It is well known that k-means clustering in a distribution is equivalent to looking for its W2W_2 optimal k point quantiser Lloyd (1982); Pollard (1982); Canas and Rosasco (2012) (see Appendix B for a quick justification). For the OT case this can be thought of as approximating the cost of delivering some goods between two sets of consumers, by the cost of the transport between major warehouses/hubs, where so long as the hubs are appropriately placed, this will represent the bulk of the overall transport (scale separation Shirdhonkar and Jacobs (2008); Bronstein et al. (2021)). Figure 10: Quantisation of various continuous measures μ with density ρμρ_μ (one per column), into a discrete measure μ¯ μ of various numbers of Diracs NδN_δ (one per row). The colored region surrounding a dashed stick denotes the continuous area that is condensed into that stick, with taller sticks (larger weight) having absorbed larger areas. With lower quantisation cost (higher NδN_δ), we of course have better proxies for estimating the Wasserstein (p=2) distance between different μs. Figure 11: Approximation of d(μ,ν)d(μ,ν) through d(μ¯,ν¯)d( μ, ν) where μ¯,ν¯ μ, ν are within balls of radius τ (solid black boundary), where d(μ¯,ν¯)→d(μ,ν)d( μ, ν)→ d(μ,ν) as τ→0τ→ 0 by the triangle inequality. For a geodesic distance, the distance d(μ,ν)d(μ,ν) is the length of the geodesic connecting μ and ν (solid blue line, while dashed connects μ¯ μ and ν¯ ν). Let the closest quantised pattern of at most size N be μ¯N=argminm∈d(μ,m)size(m)≤N, gathered μ_N= _m d_P(μ,m)\\ size(m)≤ N, gathered (29) where for a space of measures we can take size(m)=supp(m)size(m)=supp(m). So long as limN→∞d(μ¯N,μ)=0 _N→∞d_P( μ_N,μ)=0, then ∃N′∃ N for which |d(μ,ν)−d(μ¯N′,ν¯N′)|≤d(μ¯N′,μ)+d(ν¯N′,ν)<2τ,|d_P(μ,ν)-d_P( μ_N , ν_N )|≤ d_P( μ_N ,μ)+d_P( ν_N ,ν)<2τ, (30) from the triangle inequality. Additionally, with monotonically improving proxies d(μ¯N+1,μ)<d(μ¯N,μ)d_P( μ_N+1,μ)<d_P( μ_N,μ), the distance approximations also improve monotonically. Appendix Section B.1 shows that these conditions do indeed hold for W2W_2. Figure 12: Average error (evaluated over 1000 random GMM pattern pairs) vs variables required for two quantisation choices. GMMs are generated as in Section 5.1, their density is initially sampled on a uniform grid of 1024 points. Signatures require twice the amount of variables for each grid size as we are specifying positions as well as weights (plot is adjusted accordingly). The signature quantisation is achieved through hierarchical agglomerative centroid-linkage clustering. When using this approximation along with our unmixing framework, a parallel may be noted with multilevel clustering Ho et al. (2017), as we are quantising/clustering both between and within our patterns. Though of course in our case the varied basis distribution must also be learned and our discrete measures have weighted atoms. In practice, hierarchical agglomerative centroid-linkage provides good proxies, while allowing a quantiser to be built iteratively. Intuitively, similarity to k-means can be understood from noting that repeatedly merging Diracs that cause the smallest change in W2W_2 is likely to merge measures from the same Voronoi partition/cluster, when all points from a cluster are merged, we obtain a hierarchical approximation of the W2W_2 average of that cluster (see Appendix B). The merging continues until a criterion is met, such as target size or a target distance error (resulting in signatures of various sizes per input measure). Fig. 12 suggests that one can represent a distribution more accurately even with half as many “sticks” (Diracs) when compared to a uniform grid downsampling, such that there is better value in controlling positions as well as weights (compared to just a greater number of weights). In the context of the anomaly detection linear program, marginals of half the size results in transport plans of a quarter the size, which may result in solve time speed up factors of 16−6416-64 based on (n2)−(n3)O(n^2)-O(n^3) scaling for solving a naive linear program (while using more accurate marginals). 5 Experimental evaluation To demonstrate the flexibility of our method, we performed three experiments, each representing a different mixing context. In Section 5.1, we explore blind unmixing performance across 12001200 randomly generated large datasets, featuring various levels of basis variation. The independent variation of each peak represents a more extreme version of what might otherwise occur in PXRD patterns. In Section 5.2, we explore non-blind unmixing within a dataset corresponding to a more realistic, physically informed, simulation of PXRD patterns. This features a large basis, with some basis patterns being hard to distinguish after variation. In Section 5.3, we explore blind unmixing across three smaller laboratory hyperspectral imaging datasets corresponding to different physical mixing scenarios, each with different resulting mixing models. We compare against two other general1111 11 General in that they do not make any assumptions corresponding to a particular discipline or kind of spectra. unmixing methods. These are NMF and oNMF (with Kullbach-Leibler and Frobenius divergence based costs respectively). During blind unmixing, each method is tested for a constant Niter=60N_iter=60. In our visualisations, our derived signatures will be represented via kernel density estimation1212 12 If our signatures are taken to approximate continuous objects, this process approximates them by spreading out mass that was previously concentrated during discretisation., using a Gaussian kernel of specified width. As it is standard to provide Root Mean Square Error (RMSE) between ground truth and inferred quantities to measure performance, we will do so as well, wherever possible, however, as we are working with signatures, in order to measure these we must first approximate densities and evaluate them on a common grid. In such a case B refers to this further approximation, while B B is a tuple of measures as usual. RMSE(,′)=1NM∑i=1N∑j=1M(bij−bij′)2RMSE(B,B )= 1NM _i=1^N _j=1^M(b_ij-b_ij )^2 (31) To facilitate a more direct comparison, we will define a dissimilarity based a shift invariant modification of W2W_2, which we will refer to as “centered Wasserstein” W¯2(μ,ν)=W2(T#−μbμ,T#−νbν)TΔx(x)=x+Δxwithμb=∫xdμ(x). gathered W_2(μ,ν)=W_2(T_\#^-μ_bμ,T_\#^-ν_bν)\\ T^Δx(x)=x+Δx μ_b= xdμ(x). gathered (32) When ground truth is available we will use the following metrics to gauge performance, ground truth quantities will be decorated with a “t” subscript. 1. Inter Basis Centered Wasserstein (IBCW) IBCW(B,B′)=1K∑i=1KW¯2(bi,b′i)IBCW( B, B )= 1K _i=1^K W_2(b^i,b i) (33) We can similarly define inter component and basis-component variants, though these seemed to follow similar trends in our experiments. 2. Average absolute difference in percentage abundances (AADPA). 100KM∑i=1K∑j=1M|aij−(at)ij| 100KM _i=1^K _j=1^M|a_i^j-(a_t)_i^j| (34) Though the RMSE between abundances is perfectly appropriate, this is offered as a more interpretable alternative. As a standard oNMF implementation can only shift basis patterns to the right, the learned basis patterns will be offset, we center these in our plots and for RMSE calculations. Also note that with discrete convolutions we can only broaden and not thin patterns, so the basis patterns learned correspond to thinnest, rather than average variations. 5.1 Random GMM datasets Suppose our patterns are 1D GMM’s. We express our random measures with respect to the Gaussian density g(x,μ,σ)g(x;μ,σ) for mean μ and standard deviation σ ρbi(x)=∑k=1NGIkig(x,μki,σki),ρ_b^i(x)= _k=1^N_GI^i_kg(x;μ^i_k,σ^i_k), (35) where a continuous measure ν can be related to its density via ν(S)=∫Sρν(x)xν(S)= _Sρ_ν(x)dx for S⊂ℝS . This distribution is itself randomly drawn by randomly selecting Iki∼|(1,0.25)|,μki∼U(0.2,0.8)andσki∼0.02⋅2(0,0.5).I^i_k |N(1,0.25)|, μ^i_k U(0.2,0.8) σ^i_k 0.02· 2^N(0,0.5). (36) Peak intensity is denoted by I (where we must normalise ∑k=1NGIki=1 _k=1^N_GI_k^i=1), location by μ and breadth by σ. Variation is then induced by subsequently perturbing the parameters of each peak k away from those of the basis patterns, independently, per each warped basis j, per each mixture i ρb~ij(x)=∑k=0NGIkig(x,μkij,σkij)ρ_ b^ij(x)= _k=0^N_GI^i_kg(x;μ^ij_k,σ^ij_k) (37) where μkij=μki+Δμkij,σkij=σki⋅γkij.μ^ij_k=μ^i_k+Δμ^ij_k, σ^ij_k=σ^i_k· γ^ij_k. (38) with Δμkij∼(0,σΔμ)andγkij∼2(0,σp).Δμ^ij_k (0,σ_Δμ) γ^ij_k 2^N(0,σ_p). (39) The exponential is used as we would like a broadening such as σ=σ×2σ=σ× 2 to be roughly as likely as thinning by the same factor σ=σ/2σ=σ/2. Mixed examples can then be created by combining cj=∑i=1Khijb~ijc^j= _i=1^Kh^j_i b^ij where jh^j is drawn from a uniform distribution on SK−1S^K-1. We perform our assessment across three different basis variation strengths summarised in Table 2. Preset σΔμσ_Δμ σpσ_p 95% of ΔμkijΔμ^ij_k are in 95% of σkijσ^ij_k are in Low 0.015 0.15 [−0.03,+0.03][-0.03,+0.03] [0.81σki,1.23σki]~[0.81σ^i_k,1.23σ^i_k] Medium 0.025 0.3 [−0.05,+0.05][-0.05,+0.05] [0.66σki,1.52σki]~[0.66σ^i_k,1.52σ^i_k] High 0.035 0.5 [−0.07,+0.07][-0.07,+0.07] [0.50σki,2.00σki]~[0.50σ^i_k,2.00σ^i_k] Table 2: Basis variation strength presets used, where two variations of the same peak can be found as far as 14% of the domain width apart. (a) The basis attained with each decomposition method tested, one method per row. Each column corresponds to one basis pattern, which also has an associated color (e.g. b1b^1 in blue), while the ground truth basis patterns are overlayed in red. (b) Mixed datapoints cjc^j (black line), one per column, as a sum of transformed bases b~ij b^ij (area with coloring corresponding to bib^i), according to ground truth, and the tested methods, one per row. Figure 13: Visual comparison of decompositions of one randomly generated GMM dataset (high variation preset), according to each tested method. Signatures are visualised by smoothing with a Gaussian kernel with σKDE=0.004 _KDE=0.004. Note that the mixture in the first column of (b), has a component of the first basis (blue) in which some peaks are found further to the left and others are found further to the right, when compared to the left column of (a). Also, note that the mixture in the third column of (b), has a component with a different set of non-uniform shifts. Therefore, for a decomposition method to be able to explain both of these mixtures with respect to the same basis, it must have the ability to perturb each peak independently. Metric Variability Approach AADPA(,t)(A,A_t) RMSE(,t)(A,A_t) RMSE(,t)(B,B_t) IBCW(B,Bt)( B, B_t) ×100× 100 ×1000× 1000 ×100× 100 NMF 5.783±2.4695.783± 2.469 7.909±3.4217.909± 3.421 1.742±0.5691.742± 0.569 1.294±0.6021.294± 0.602 Low oNMF 8.300±2.9218.300± 2.921 10.882±3.72210.882± 3.722 4.611±1.4544.611± 1.454 3.559±1.5073.559± 1.507 Our approach 4.274±1.9734.274± 1.973 5.544±2.7575.544± 2.757 1.668±0.4851.668± 0.485 1.756±0.8941.756± 0.894 NMF 10.247±3.80010.247± 3.800 14.142±5.12014.142± 5.120 3.130±0.9263.130± 0.926 2.262±1.0562.262± 1.056 Medium oNMF 11.955±3.13411.955± 3.134 15.539±3.84615.539± 3.846 7.138±1.7107.138± 1.710 5.220±1.7715.220± 1.771 Our approach 6.019±2.0646.019± 2.064 8.154±2.9328.154± 2.932 2.574±0.5872.574± 0.587 1.837±0.8531.837± 0.853 NMF 14.167±4.02314.167± 4.023 19.276±5.18019.276± 5.180 4.271±1.0104.271± 1.010 3.299±1.2433.299± 1.243 High oNMF 14.413±2.92814.413± 2.928 18.743±3.55018.743± 3.550 8.723±1.6368.723± 1.636 5.769±1.5895.769± 1.589 Our approach 8.371±2.4938.371± 2.493 11.509±3.45111.509± 3.451 3.625±0.7323.625± 0.732 2.154±0.8002.154± 0.800 Table 3: Mean decomposition metrics over 400 random GMM datasets per variability preset. These were generated with K=3K=3 and M=800M=800, where the density of each cjc^j is initially sampled uniformly with N=256N=256 then discretised with Nδ=64N_δ=64. The best of eight decompositions with different initialisation (as measured by final cost), is taken for each dataset/method. The ground truth basis (Bt B_t) is also similarly sampled (producing tB_t), as is the kde density of the inferred basis signature for basis RMSEs. The oNMF is performed with a maximum shift of 25 gridpoints and the inferred basis patterns are centered by shifting for 13 gridpoints before measuring RMSE(,t)(B,B_t). We can see from Table 3 that sometimes the oNMF model can perform worse than NMF in extracting the information we are after, this is despite the fact that our oNMF solutions are initialised with NMF solutions. Even if the overall fit can improve, the accuracy of inferred constituents can simultaneously suffer (see Fig. 13). In fact, the more flexible a composition model, the more trivially a profile can be fit, regardless of basis. For instance, if one were to use OT as a basis variation model ϕ(T,b)=T#bϕ(T;b)=T_\#b, one could arbitrarily warp any given basis to perfectly fit any mixture, not giving us any indication as to weather one basis is more appropriate over another. This supports our original reasoning for defining a model driven by a non-reconstruction cost. 5.2 PXRD unmixing As an indicator of suitability to PXRD patterns, we can generate a benchmarking dataset with the Synthetic instance generator Bras et al. (2014), which provides a more physically informed model of such data (compared to GMMs). Similarly to the previous, this too features non-uniform shifting and broadening. The default preset of this generator (the AlLiFe system) features a relatively large basis (6 patterns, corresponding to pure phases). As we found that none of the methods we compare here could consistently blindly discover all 6 of the basis patterns, we assess separation/abundance estimation in the non-blind (see Table 1) regime here, the performance in this task is summarised in Table 4 and illustrated in Fig. 14. (a) The 6 ground truth basis patterns. (b) Mixed datapoints cjc^j (black line), one per column, as a sum of transformed bases b~ij b^ij (area with coloring corresponding to bib^i), according to ground truth, and the tested methods, one per row. Figure 14: Visual comparison of basis separation (non-blind unmixing) according to the tested methods (one per row), of a Synthetic instance generator dataset with K=6K=6, m=M=1275m=M=1275, Nδ=64N_δ=64. To visualise signatures, σKDE=0.188 _KDE=0.188 has been used, x domain (2θ2θ) has been truncated to [16,45][16,45] in these plots. Metric Approach AADPA(,t)(A,A_t) RMSE(,t)(A,A_t) NMF 10.467 0.203 oNMF 6.825 0.118 Our approach 3.115 0.055 Table 4: Abundance derivation metrics for the compared approaches (in the separation-only/fixed-basis regime), over the synthetic instance generator dataset. Table 4 shows significantly lower error in abundance prediction for our approach. Fig. 14 again demonstrates the merit of a non-reconstruction based cost, in that, though the sums of the components come close to matching the profile of the mixture for NMF and oNMF in the middle column of Fig. 14(b), they are not able to correctly relate features of the mixture with features from the basis patterns, and therefore wrongly partition the area under a mixture, and derive less accurate varied bases and abundances. In the right column of Fig. 14(b), where the purple pattern cannot be adequately varied by NMF and oNMF, they over assign mass to the similar, pink pattern instead. 5.3 Hyperspectral imaging This dataset Zhao et al. (2019) contains 3 scenes containing hyperspectral images of colored powders prepared in a laboratory in known amounts to showcase some different mixing scenarios in HI, these are as follows. Scene 1 (Linear), 3 distinct powders are arranged in a grid of squares smaller than that of the hyperspectral instruments pixels, the subpixel pattern is repeated periodically. A light source points directly at the powders, light generally travels from source, to a single powder (without significant scattering), to a hyperspectral pixel, though rays passing by each of the subpixel squares can make it to the same pixel, adding to the measurement independently. Scene 2 (Intimate), 4 distinct powders are mixed uniformly, the same ray of light can scatter off of multiple differently colored grains before reaching the sensor, as such this contains second-order reflections (and possibly higher) involving every combination of basis patterns. Scene 3 (Limited bilinear), 3 distinct powders are mixed as in the first case, though now the light source points at an adjacent plastic board, such that most rays will go from source, to board, to a subpixel square, to a pixel. This contains second-order reflections of the board with each of the first three basis patterns. (a) The basis attained with NMF vs our method. Each column corresponds to one basis pattern, which also has an associated color (e.g. b1b^1 in blue), while the ground truth basis patterns are overlayed in red. (b) Mixed datapoints cjc^j (black line), one per column, as a sum of transformed bases b~ij b^ij (area with coloring corresponding to bib^i), according to NMF in the top row and our approach in the bottom row. Figure 15: Visual comparison of decompositions of scene 1. To visualise signatures, σKDE=7.53 _KDE=7.53 has been used. Metric Scene Approach AADPA(,t)(A,A_t) RMSE(,t)(A,A_t) RMSE(,t)(B,B_t) IBCW(B,Bt)( B, B_t) ×10× 10 ×1000× 1000 1 NMF 4.5614.561 0.5490.549 2.0722.072 45.26945.269 Our approach 4.4374.437 0.5950.595 1.9471.947 41.68341.683 2 NMF 11.69711.697 1.4861.486 1.8301.830 35.21235.212 Our approach 9.7879.787 1.3841.384 1.6481.648 40.19840.198 3 NMF 12.78812.788 1.6021.602 2.4582.458 49.15749.157 Our approach 3.0653.065 0.4590.459 1.9851.985 40.40340.403 Table 5: Performance summary of NMF vs our method (in the fully blind regime), on a HI dataset. The best of 16 decompositions with different initial starts (as measured by final cost), is taken for each scene/method. oNMF has been omitted in this instance due to instability. Each scene contains a number of images (10, 14 and 8 respectively), where each image corresponds to different mixing ratios. While each pixel in each image could contain unique basis variation, due to the previously mentioned redundancy, we treat the Euclidean average1313 13 The Euclidean average is used simply to avoid resampling for NMF, which requires the inputs to be on the same grid. pixel of each mixed image as our dataset. Each pixel is a vector with 256 entries representing intensity at each point in the linearly spaced spectral band from 400400 to 10001000nm, we generate signatures with Nδ=128N_δ=128 for our method. For each scene we perform fully blind unmixing. Though a reconstruction based method might have otherwise used any additional degrees of freedom in order to try to provide a closer fit, injecting unnecessary variation (at the expense of ground truth recovery), scene 1 shows (see Table 5) that our method does not suffer from this issue (as we minimise basis variation). Thus we provide similar decompositions in purely linear cases, as the appropriately constrained model, while having flexibility should more complex cases show up. In cases featuring bilinear contributions κjbinB=∑p,qapqjbp⊙bqκ_A^j^bin B= _p,qa_pq^jb^p b^q, we can still make progress with a linear combination function, noting that our method has the flexibility to implicitly construct basis variations b~ij b^ij that are themselves some combination of the first-order and the second-order reflections of bib^i (with each of the other basis patterns). The more second-order contributions in cjc^j, the more b~ij b^ij deviates from bib^i. Our approach somewhat improves on a linear decomposition on scene 2, but improves greatly on a linear decomposition of scene 3 where the presence of second-order reflections is more limited, but still seemingly significant. 6 Conclusion Many scientific disciplines work with “mixed data,” where each measurement results from a combination of underlying basis patterns that can themselves vary between samples (see Fig. 1). Prevailing unmixing methods largely derive from a specific choice of variation function, meaning that datasets that feature different forms of variation likely need to be processed with different specific methods. If no such method is available, significant adaptation is then required. These methods typically optimize for reconstruction accuracy, that is, how closely the model’s output matches the measured signal. Reconstruction accuracy alone does not provide sufficient guidance, as multiple combinations of inputs might achieve a perfect fit for some mixtures. Not all inputs are equally plausible, yet reconstruction does not distinguish between them. Lastly, pointwise distances are often used in constructing a cost, which can vary erratically with respect to non-pointwise basis variation (see Fig. 5). Our work introduces a different paradigm for unmixing based on geometric approaches (see Section 2.3.1). One can imagine a set of seed inputs which transform into a latent point cloud after passing through the variation function ϕ, and into the final measurements after passing through the combination function κ. Geometric approaches typically draw conclusions from geometric features of the point cloud of measurements, which can have a particular structure as a result of the mixing model. We take this further by considering all latent point clouds that exactly produce the given measurements, and weigh these by considering the kinds of point cloud that the variation function ϕ can feasibly produce. We show that for any ϕ that is Lipchitz with respect to some distance d_P, one can expect clustering in the metric space (d,)(d_P,P), where a smaller bound on maximum variation parameters (e.g. shifting/broadening) produces tighter clusters (see Appendix A). As such, we propose a (in-class) Fréchet variance based cost to search for a latent point cloud that correspond to the least basis variation possible. To handle common non-pointwise basis variations like shifting and broadening, a corresponding non-pointwise distance metric is required. The Wasserstein distance is a natural choice when signals are treated as distributions. Our practical implementation leverages several key features of OT to ensure efficiency: quantization of distributions into sparse, representative “signatures”; a fast algorithm (NWCA) to compute exact 1D Wasserstein distances and transport maps in linear time; and displacement interpolation for the efficient calculation of Wasserstein centroids (averages). These innovations make our method computationally feasible for exploratory data analysis. As demonstrated through extensive experiments on random GMM, simulated PXRD, and laboratory HI data, our approach consistently outperforms other general-purpose methods in a range of unmixing scenarios. This work provides a strong foundation for a wide range of applications and opens several avenues for future research. These include clustering with a full Fisher discriminant by including between-class variance terms, exploring application to higher-dimensional data (this requires fast nD transport map calculation), bidirectional partial matching during anomaly detection, and, finally, exploring other ground distances, including learned distances, or even non-OT distances all together. Declarations These results were obtained as part of a PhD funded by Johnson Matthey. The authors acknowledge support from the Leverhulme Trust via the Leverhulme Research Centre for Functional Materials Design (RC-2015-036), UK Engineering and Physical Sciences Research Council (EPSRC) through grant number EP/V026887 as well as through AI for Chemistry: AIchemy Hub (EPSRC grant EP/Y028775/1 and EP/Y028759/1). V. V. Gusev is supported by the Leverhulme Research Fellowship. Appendix A Clustering and separability of metric space perturbations Theorem 2 gives us that, so long as ϕ(λ,b)ϕ(λ;b) is Lipchitz with respect to d_P, it is possible to guarantee well defined basis variant clusters for any basis B∈K B ^K (according to the criteria of Lemma 3), up to some (basis dependent) bound on the parameter λ. Lemma 3 (Voronoi separable). With a maximum radius among clusters maxijd(bi,b~ij)<R _ijd_P(b^i, b^ij)<R and adequate basis pattern spacing mini≠jd(bi,bj)>2R _i≠ jd_P(b^i,b^j)>2R (or maximum similarity), the corresponding basis of some varied pattern can be easily identified, from the Voronoi diagram of the bases b~jk∈Xi⇒j=i b^jk∈ X^i j=i for Voronoi cell XiX^i (see Eq. 45), and in fact d(bi,b~jk)<R⇒j=i,d_P(b^i, b^jk)<R j=i, where each cluster is enclosed in a ball which is enclosed in its own Voronoi cell, bi∗⊂BR(bi)⊂Xib^i*⊂ B_R(b^i)⊂ X^i. Proof. We require that no variations spill into a neighbouring Voronoi cell, i.e. no two clusters touch/share a common point min∀i≠k∀j,ld(b~ij,b~kl)>0. _ subarrayc∀ i≠ k\\ ∀ j,l subarrayd_P( b^ij, b^kl)>0. (40) From the triangle inequality d(bi,bk)≤d(bi,b~ij)+d(b~ij,b~kl)+d(b~kl,bk)⇒d(b~ij,b~kl)≥d(bi,bk)−2R, gatheredd_P(b^i,b^k)≤ d_P(b^i, b^ij)+d_P( b^ij, b^kl)+d_P( b^kl,b^k)\\ d_P( b^ij, b^kl)≥ d_P(b^i,b^k)-2R,\\ gathered (41) holds for all i,j,k,li,j,k,l, hence holds also for the choice min∀i≠k∀j,ld(b~ij,b~kl)≥mini≠jd(bi,bj)−2R>0.∎ _ subarrayc∀ i≠ k\\ ∀ j,l subarrayd_P( b^ij, b^kl)≥ _i≠ jd_P(b^i,b^j)-2R>0. Corollary 3. Clearly if d(ϕ(λ,bi),bi)≤R∀bi∈Bd_P(ϕ(λ;b^i),b^i)≤ R\;\;∀ b^i∈ B then maxijd(bi,b~ij)<R _ijd_P(b^i, b^ij)<R. Theorem 2. Let ϕ be Lipschitz continuous s.t. ∃Lb∃ L_b with d(ϕ(λ,b),ϕ(λ′,b))≤LbdΛ(λ,λ′)d_P(ϕ(λ;b),ϕ(λ ;b))≤ L_bd_Λ(λ,λ ), and let ϕ(λ0,b)=bϕ(λ_0;b)=b at some λ0λ_0. Voronoi separability can be guaranteed through a parameter limit dΛ(λ,λ0)≤dΛmaxd_Λ(λ,λ_0)≤ d_Λ , so long as dΛmax<mini≠jd(bi,bj)/2maxiLbid_Λ < _i≠ jd_P(b^i,b^j)/2 _iL_b^i. Proof. We have d(ϕ(λ,bi),bi)≤LbidΛ(λ,λ0)≤maxiLbidΛ(λ,λ0)≤maxiLbidΛmaxd_P(ϕ(λ;b^i),b^i)≤ L_b^id_Λ(λ,λ_0)≤ _iL_b^id_Λ(λ,λ_0)≤ _iL_b^id_Λ , giving R=maxiLbidΛmaxR= _iL_b^id_Λ . ∎ Remark 2. With a compact pattern space, mini≠jd(bi,bj) _i≠ jd_P(b^i,b^j) cannot be greater than some dmaxd_P , so it is necessary that d(ϕ(λ,b),b)≤R<dmax/2d_P(ϕ(λ;b),b)≤ R<d_P /2. Figure 16: Left: Simplex S2 S^2, where weights on a static grid such as =(0,0.5,1)x=(0,0.5,1) can be represented. The repeated application of the shift-down operation (initially on 1e^1) is illustrated. Right: The Dirac measures δx∈P(ℝ)δ_x∈ P(R) representing distributions with weights concentrated at a single point x, which we can continuously interpolate between. For instance, with Variations generated by discrete shifting where ϕ(s,)=↓s∈Hull(i:i∈[K]∪)ϕ(s;p)=p s (\e^i:i∈[K]\∪\0\) (accounting for mass loss from shifting past the end of a vector) with s∈0,1,2,…s∈\0,1,2,…\ and ∈SN−1p∈ S^N-1, we can have d(↓1,)=dmax≮dmax/2d_P(p 1,p)=d_P d_P /2 (see Fig. 16) when d_P is the Euclidean distance, so basis variation as a result of shifting, is not Voronoi separable over =Hull(i:i∈[K]∪)P=Hull(\e^i:i∈[K]\∪\0\) with the Euclidean distance. Appendix B Wasserstein quantisation and weighted k-means Minimising μ¯N=argminm∈P()W2d(μ,m)2|supp(m)|≤N, gathered μ_N= _m∈ P(X)W_2^d_X(μ,m)^2\\ |supp(m)|≤ N, gathered (42) is generally a semidiscrete OT problem1414 14 Where μ could be discrete as well, but only m is constrained to be. for the discrete measure m=∑i=1Nwmiδxmim= _i=1^Nw^i_mδ_x^i_m, where transport maps between μ and m are necessarily of the form T(x)=xmiT(x)=x^i_m for x∈Xix∈ X^i where XiX^i are disjoint so that they partition X. The weights are determined from mass conservation m(xmi)=wmi=μ(T−1(xmi))=μ(Xi)m(\x^i_m\)=w^i_m=μ(T^-1(\x^i_m\))=μ(X^i), hence Eq. 42 can be stated as minxmi∈Xi⊂∑i=1N∫Xid(x,xmi)2dμ(x). _ subarraycx^i_m \\ X^i subarray _i=1^N _X^id_X (x,x^i_m )^2dμ(x). (43) This can be solved in an alternating fashion, where given regions Xi∀i\X^i\>∀ i\, the corresponding optimal xmix^i_m are centroids1515 15 Where xmi=∫Xixμ(x)/∫Xidμ(x)x^i_m= _X^ixdμ(x)/ _X^idμ(x) when d_X is the Euclidean distance of each region xmi=argminx′∈∫Xid(x,x′)2dμ(x)x^i_m= _x _X^id_X (x,x )^2dμ(x) (44) and given supp(m)=xmi∀isupp(m)=\x^i_m\>∀ i\, the optimal regions are then given by the Voronoi tessellation Xi=x∈:d(x,xmi)<d(x,xmj)∀j≠i.X^i=\x :d_X(x,x^i_m)<d_X(x,x^j_m)\>∀ j≠ i\. (45) To see this, suppose XiX^i violated this condition, then for some j, X^=x∈Xi:d(x,xmi)>d(x,xmj) X=\x∈ X^i:d_X(x,x^i_m)>d_X(x,x^j_m)\ is non-empty, then one can reassign each x∈X^x∈ X according to X¯i=Xi∖X X^i=X^i X and X¯j=Xj∪X X^j=X^j∪ X and by doing so improve the cost ∫X¯id(x,xmi)2dμ(x)+∫X^d(x,xmi)2dμ(x)+∫Xjd(x,xmj)2dμ(x) _ X^id_X (x,x^i_m )^2dμ(x)+ _ Xd_X (x,x^i_m )^2dμ(x)+ _X^jd_X (x,x^j_m )^2dμ(x) (46) > > ∫X¯id(x,xmi)2dμ(x)+∫X^d(x,xmj)2dμ(x)+∫Xjd(x,xmj)2dμ(x), _ X^id_X (x,x^i_m )^2dμ(x)+ _ Xd_X (x,x^j_m )^2dμ(x)+ _X^jd_X (x,x^j_m )^2dμ(x), hence violating Eq. 45 is suboptimal. Regardless of the strict optimality of this step, the sequence of updates produced by alternating these steps are non-increasing and the cost is bound below by 00, hence this approach converges to a local minimiser of Eq. 42. This alternating approach is a generalised, weighted version of the k-means algorithm which can be seen as a generalisation of Lloyds algorithm Lloyd (1982), for which =ℝX=R, Xi=x:xmidi−1<x<xmidiX^i=\x:x_mid^i-1<x<x_mid^i\ and xmidi=(xmi+xmi+1)/2x_mid^i=(x^i_m+x^i+1_m)/2 with xmid0=−∞x_mid^0=-∞ and xmidN=∞x_mid^N=∞. B.1 Wasserstein quantiser monotonicity and convergence Theorem 3 (Monotonicity). Let μ∈P()μ∈ P(X), its optimal W2W_2 quantisers of size N (μ¯N μ_N) have monotonically decreasing error W22(μ,μ¯N)>W22(μ,μ¯N+1).W_2^2(μ, μ_N)>W_2^2(μ, μ_N+1). (47) This follows from being able to construct a suboptimal quantiser of size N+1N+1 that is better than the optimal one of size N. Proof. Let μ∈P()μ∈ P(X) and its optimal quantiser of size N be μ¯N μ_N. Pick x′∈XN∖xNx ∈ X^N \x^N\ and let A,BA,B be the corresponding Voronoi sets of xn,x′x^n,x partitioning XNX^N, we may construct μ~N+1=μ¯N−μ(XN)δxN+μ(A)δxN+μ(B)δx′, μ_N+1= μ_N-μ(X^N)δ_x^N+μ(A)δ_x^N+μ(B)δ_x , (48) where we have let xμ¯Ni=xix^i_ μ_N=x^i for brevity. We have that ∫XNd(x,xN)2dμ(x) _X^Nd_X (x,x^N )^2dμ(x) =∫Ad(x,xN)2dμ(x)+∫Bd(x,xN)2dμ(x) = _Ad_X (x,x^N )^2dμ(x)+ _Bd_X (x,x^N )^2dμ(x) (49) >∫Ad(x,xN)2dμ(x)+∫Bd(x,x′)2dμ(x), > _Ad_X (x,x^N )^2dμ(x)+ _Bd_X (x,x )^2dμ(x), since by definition B=x∈XN:d(x,x′)<d(x,xN)B=\x∈ X^N:d(x,x )<d(x,x^N)\. Hence W22(μ,μ¯N)−W22(μ,μ~N+1)>0W_2^2(μ, μ_N)-W_2^2(μ, μ_N+1)>0 and so we have monotonicity W22(μ,μ¯N)>W22(μ,μ~N+1)≥W22(μ,μ¯N+1).∎W_2^2(μ, μ_N)>W_2^2(μ, μ_N+1)≥ W_2^2(μ, μ_N+1). Theorem 4 (Convergence). Let μ∈P()μ∈ P(X), its optimal W2W_2 quantisers of size N (μ¯N μ_N) weakly converge to it. W22(μ,μ¯N)→0asN→∞W_2^2(μ, μ_N)→ 0 N→∞ (50) Proof. Consider a sequence of suboptimal quantisers μ~N μ_N that converge to μ, then optimal quantisers must also converge. Construct such a sequence by noting that at each N, one can pick the largest diameter Voronoi region, where diam(X)=maxx,y∈Xd(x,y)diam(X)= _x,y∈ Xd_X (x,y ), over which one has a 1 Dirac distribution, and by monotonicity we can generate a 2 Dirac quantisation over this region with a better cost. Clearly the diameter of all regions shrinks to 00 in the limit. W22(μ,μ¯N)<W22(μ,μ~N) W_2^2(μ, μ_N)<W_2^2(μ, μ_N) =∑i=1N∫Xμ~Nid(x,xμ~Ni)2dμ(x) = _i=1^N _X^i_ μ_Nd_X (x,x^i_ μ_N )^2dμ(x) (51) ≤∑i=1Nμ(Xμ~Ni)diam(Xμ~Ni)2 ≤ _i=1^Nμ(X^i_ μ_N)diam(X^i_ μ_N)^2 ≤[∑i=1Nμ(Xμ~Ni)]maxidiam(Xμ~Ni)2→0 ≤ [ _i=1^Nμ(X^i_ μ_N) ] _idiam(X^i_ μ_N)^2→ 0 where ∑i=1Nμ(Xμ~Ni)=1 _i=1^Nμ(X^i_ μ_N)=1. ∎ Appendix C Weighted centroids parameterise the geodesic Lemma 4. Provided that ∀p1,p2∈∀ p^1,p^2 and ∀d1b,d2b≥0∀ d_1b,d_2b≥ 0 with d1b+d2b≥d(p1,p2)=d12d_1b+d_2b≥ d(p^1,p^2)=d_12 ∃pb∈∃ p^b s.t. d(p1,pb)=d1bd(p^1,p^b)=d_1b and d(p2,pb)=d2bd(p^2,p^b)=d_2b, then the centroids pb(λ)=Bd(p1,p2,λ)p^b(λ)=B_d(p^1,p^2,λ) lie on the geodesic between p1p^1 and p2p^2. Further, the curve given by pb(λ)p^b(λ) is constant speed and is therefore the geodesic between p1p^1 and p2p^2. Proof. B(p1,p2,λ)B(p^1,p^2,λ) is attained when the barycentric cost takes the value mind1b,d2b∈ℝ+(1−λ)d1b2+λd2b2s.t.d1b+d2b≥d12 gathered _d_1b,d_2b _+(1-λ)d_1b^2+λd_2b^2\\ s.t. d_1b+d_2b≥ d_12 gathered (52) with Lagrangian ℒ(d1b,d2b,γ)=λd1b2+(1−λ)d2b2+γ(d12−d1b−d2b).L(d_1b,d_2b,γ)=λd_1b^2+(1-λ)d_2b^2+γ(d_12-d_1b-d_2b). (53) It is necessary for optimality that we have stationarity (Karush-Kuhn-Tucker conditions), ∂d1bℒ=∂d2bℒ=0 _d_1bL= _d_2bL=0, hence γ=2(1−λ)d1b=2λd2bγ=2(1-λ)d_1b=2λd_2b. Outside the trivial cases p1=p2p^1=p^2 or λ∈0,1λ∈\0,1\, it follows that γ>0γ>0. Since γ≠0γ≠ 0, the bound is tight, i.e. d1b+d2b=d12d_1b+d_2b=d_12 indicating that ∀λ∈(0,1)∀ λ∈(0,1), pb(λ)p^b(λ) sits on the geodesic connecting p1,p2p^1,p^2. Further, d2b=d1b(1−λ)/λd_2b=d_1b(1-λ)/λ, so that d12=[1+(1−λ)/λ]d1b=d1b/λd_12= [1+(1-λ)/λ ]d_1b=d_1b/λ. Since we know pb(λ)p^b(λ) at worst reparameterises the geodesic, d(p1,pb(λ))=λd(p1,p2)d(p^1,p^b(λ))=λd(p^1,p^2) is enough to know it is constant speed. ∎ Appendix D General McCann interpolation The following is related to the Monge OT minimisation minTt∫C[(Tt(x))0≤t≤1]dμ(x):T0=IdandT#1μ=ν, _T^t \ C[(T^t(x))_0≤ t≤ 1]dμ(x): T^0=Id T^1_\#μ=ν \, (54) that is to say, rather than evaluating the final displacement of x to y, evaluate the whole trajectory taken between these. This is one dynamic formulation of the OT problem. At time t, the mass at x is instead at Tt(x)T^t(x) and only at t=1t=1 does it arrive at y. These two problems are equivalent (with T∗1=T∗T^*1=T^*) provided that the static cost agrees with the cost of the best trajectory, i.e. that c(x,y)=min(Tt(x))0≤t≤1C[(Tt(x))0≤t≤1]:T0(x)=xandT1(x)=yc(x,y)= _(T^t(x))_0≤ t≤ 1 \C[(T^t(x))_0≤ t≤ 1]:T^0(x)=x T^1(x)=y \ (55) is met Villani (2003). Definition 10 (Length of a trajectory). For the trajectory τ:[0,1]→τ:[0,1] , the length is defined as Ld[τ]=sup∑i=1Nd(τ(ti−1),τ(ti)),L_d_X[τ]= _i=1^Nd_X(τ(t_i-1),τ(t_i)), (56) with the supremum taken over all N and sequences 0=t0≤…≤tN=10=t_0≤…≤ t_N=1. Definition 11 (Length space). A metric space (,d)(X,d_X) with d(x,y)=infτL[τ]d_X(x,y)= _τL[τ] over all τ:[0,1]→τ:[0,1] with τ(0)=xτ(0)=x and τ(1)=yτ(1)=y. Definition 12 (Geodesic). A geodesic from x∈x to y∈y , is a trajectory τ:[0,1]→τ:[0,1] with τ(0)=xτ(0)=x, τ(1)=yτ(1)=y which has (1.) minimal length and (2.) constant speed. (1.) can be ensured by L[(τ)t1≤t≤t2]=d(τ(t1),τ(t2))L[(τ)_t^1≤ t≤ t^2]=d_X(τ(t^1),τ(t^2)) for all 0≤t1<t2≤10≤ t^1<t^2≤ 1, and (2.) by L[(τ)t1≤t≤t2]=(t2−t1)λL[(τ)_t^1≤ t≤ t^2]=(t^2-t^1)λ for some λ, again for all 0≤t1<t2≤10≤ t^1<t^2≤ 1. By Definition 10, (1.) is given by (2.) when λ=d(x,y)λ=d_X(x,y). Definition 13 (Geodesic space). A metric space (,d)(X,d_X) in which ∀x,y∈∀ x,y , there exists a geodesic connecting x and y. Geodesic spaces are also length spaces. Suppose c=dpc=d_X^p where the ground space (,d)(X,d_X) is a geodesic space, if we let C[τ]=Ld[τ]pC[τ]=L_d_X[τ]^p and without loss of correspondence to the Monge problem, let Eq. 54 also require that all considered paths are constant speed, then clearly the optimal trajectories (T∗t(x))0≤t≤1(T^*t(x))_0≤ t≤ 1 are ground space geodesics between x and T∗(x)T^*(x). By Lemma 4, this optimal trajectory can then be expressed in terms of ground space barycenters T∗t(x)=Bd(x,T∗(x),t)T^*t(x)=B_d_X(x,T^*(x),t), where T∗T^* is the optimal map in the static problem. Lemma 5. Let the OT problem between μ and ν admit a Monge solution and let T∗tT^*t be the (constant speed) dynamic OT solution. Let T12T^12 be the (static) optimal transport map between T#∗t1μT^*t^1_\#μ and T#∗t2μT^*t^2_\#μ. Provided that T12(T∗t1(x))=T∗t2(x)T^12(T^*t^1(x))=T^*t^2(x) for all 0≤t1<t2≤10≤ t^1<t^2≤ 1, then (T#∗tμ)0≤t≤1(T^*t_\#μ)_0≤ t≤ 1 is a Wasserstein space geodesic between μ and ν. Proof. Note that the shortest path between T∗t1(x)T^*t^1(x) to T∗t2(x)T^*t^2(x) is given by the geodesic segments (T∗t(x))t1≤t≤t2(T^*t(x))_t^1≤ t≤ t^2, then Wpd(T#∗t1μ,T#∗t2μ)p W_p^d_X(T^*t^1_\#μ,T^*t^2_\#μ)^p =∫d(x,T12(x))pd(T#∗t1μ)(x) = d_X(x,T^12(x))^pd(T^*t^1_\#μ)(x) =∫d(T∗t1(x),T∗t2(x))pμ(x) = d_X(T^*t^1(x),T^*t^2(x))^pdμ(x) =∫Ld[(T∗t(x))t1≤t≤t2]pμ(x) = L_d_X[(T^*t(x))_t^1≤ t≤ t^2]^pdμ(x) =∫[(t2−t1)d(x,T∗(x))]pμ(x)=(t2−t1)pWpd(μ,ν)p.∎ = [(t^2-t^1)d_X(x,T^*(x))]^pdμ(x)=(t^2-t^1)^pW_p^d_X(μ,ν)^p. Where Lemma 5 applies, such as in 1D, where feasible monotonic plans are unique, Lemma 4 can be employed again, allowing one to calculate Wasserstein barycenters via B(μ,ν,λ)=T#∗λμB(μ,ν,λ)=T^*λ_\#μ, this is the displacement interpolation. This is usually stated for an Euclidean ground space where T∗λ=(1−λ)Id+λT∗T^*λ=(1-λ)Id+λT^*. References Agrawal et al. (2022) J. Agrawal, M. Gupta, and H. Garg Blind source separation in perspective of ica algorithms: a review. In 2022 International Conference on Computational Intelligence and Sustainable Engineering Solutions (CISES), Vol. , p. 78–85. External Links: Document Cited by: §2.3. Aharoni et al. (2022) R. Aharoni, A. Zuck, D. Peri, and S. Kendler Reflectance spectra analysis algorithms for the characterization of deposits and condensed traces on surfaces. In Colorimetry, Prof. A. K. Samanta (Ed.), External Links: Document Cited by: §1. Arjovsky et al. (2017) M. Arjovsky, S. Chintala, and L. Bottou Wasserstein generative adversarial networks. In Proceedings of the 34th International Conference on Machine Learning - Volume 70, ICML’17, Sydney, NSW, Australia, p. 214–223. Cited by: Remark 1. Bai et al. (2018a) J. Bai, S. Ament, G. Perez, J. Gregoire, and C. Gomes An efficient relaxed projection method for constrained non-negative matrix factorization with application to the phase-mapping problem in materials science. In Integration of Constraint Programming, Artificial Intelligence, and Operations Research, W. van Hoeve (Ed.), Cham, p. 52–62. External Links: ISBN 978-3-319-93031-2 Cited by: §1. Bai et al. (2018b) J. Bai, Y. Xue, J. Bjorck, R. L. Bras, B. Rappazzo, R. Bernstein, S. K. Suram, R. B. van Dover, J. M. Gregoire, and C. P. Gomes Phase‐mapper: accelerating materials discovery with ai. AI Mag. 39 (1), p. 15–26. External Links: ISSN 0738-4602, Link, Document Cited by: §2.3. Bai et al. (2023) Y. Bai, B. Schmitzer, M. Thorpe, and S. Kolouri Sliced optimal partial transport. In 2023 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), Vol. , p. 13681–13690. External Links: Document Cited by: §4.1. Boissard et al. (2015) E. Boissard, T. L. Gouic, and J. Loubes Distribution’s template estimate with Wasserstein metrics. Bernoulli 21 (2), p. 740 – 759. External Links: Document, Link Cited by: §1, §3.1. Borsoi et al. (2021) R. A. Borsoi, T. Imbiriba, J. C. M. Bermudez, C. Richard, J. Chanussot, L. Drumetz, J. Tourneret, A. Zare, and C. Jutten Spectral variability in hyperspectral data unmixing: a comprehensive review. IEEE Geoscience and Remote Sensing Magazine 9 (4), p. 223–270. External Links: Document Cited by: §1, §2.3. Boulais, A. et al. (2021) Boulais, A., Berné, O., Faury, G., and Deville, Y. Unmixing methods based on nonnegativity and weakly mixed pixels for astronomical hyperspectral datasets. A&A 647, p. A105. External Links: Document, Link Cited by: §1. Bras et al. (2014) R. Bras, R. Bernstein, J.M. Gregoire, S. Suram, C. Gomes, B. Selman, and R. van Dover A computational challenge problem in materials discovery: synthetic problem generator and real-world datasets. Proceedings of the National Conference on Artificial Intelligence 1, p. 438–443. Cited by: §5.2. Bronstein et al. (2021) M. M. Bronstein, J. Bruna, T. Cohen, and P. Veličković Geometric deep learning: grids, groups, graphs, geodesics, and gauges. arXiv. External Links: Document Cited by: §4.2. Canas and Rosasco (2012) G. Canas and L. Rosasco Learning probability measures with respect to optimal transport metrics. In Advances in Neural Information Processing Systems, F. Pereira, C.J. Burges, L. Bottou, and K.Q. Weinberger (Eds.), Vol. 25, Lake Tahoe, Nevada, USA, p. . Cited by: §4.2. CHERNI et al. (2020) A. CHERNI, S. ANTHOINE, and C. CHAUX β-NMF and sparsity promoting regularizations for complex mixture unmixing - application to 2d hsqc nmr. In ICASSP 2020 - 2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Vol. , p. 5795–5799. External Links: Document Cited by: §1. Cichocki et al. (2009) A. Cichocki, R. Zdunek, A. Phan, and S. Amari Nonnegative matrix and tensor factorizations: applications to exploratory multi-way data analysis and blind source separation. External Links: ISBN 978-0-470-74666-0, Document Cited by: §2.3. Du et al. (2014) X. Du, A. Zare, P. Gader, and D. Dranishnikov Spatial and spectral unmixing using the beta compositional model. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing 7 (6), p. 1994–2003. External Links: Document Cited by: §2.3.1. Engquist et al. (2016) B. Engquist, B. Hamfeldt, and Y. Yang Optimal transport for seismic full waveform inversion. Communications in Mathematical Sciences 14, p. . External Links: Document Cited by: §1, §3, Remark 1. Ermon et al. (2015) S. Ermon, R. Le Bras, S. Suram, J. Gregoire, C. Gomes, B. Selman, and R. van Dover Pattern decomposition with complex combinatorial constraints: application to materials discovery. Proceedings of the AAAI Conference on Artificial Intelligence 29 (1). External Links: Link, Document Cited by: §1. Flamary (2019) R. Flamary Optimal transport for machine learning. Ph.D. Thesis, Université Cote d’Azur. Cited by: Remark 1. Flamary et al. (2016) R. Flamary, C. Févotte, N. Courty, and V. Emiya Optimal spectral transportation with application to music transcription. In Proceedings of the 30th International Conference on Neural Information Processing Systems, NIPS’16, Red Hook, NY, USA, p. 703–711. External Links: ISBN 9781510838819 Cited by: §1. Fletcher et al. (2004) P.T. Fletcher, C. Lu, S.M. Pizer, and S. Joshi Principal geodesic analysis for the study of nonlinear statistics of shape. IEEE Transactions on Medical Imaging 23 (8), p. 995–1005. External Links: Document Cited by: §2.3.1. Gaffney and Smyth (2004) S. Gaffney and P. Smyth Joint probabilistic curve clustering and alignment. In Proceedings of the 17th International Conference on Neural Information Processing Systems, NIPS’04, Cambridge, MA, USA, p. 473–480. Cited by: §1. Glatz et al. (2011) J. Glatz, N. C. Deliolanis, A. Buehler, D. Razansky, and V. Ntziachristos Blind source unmixing in multi-spectral optoacoustic tomography. Opt. Express 19 (4), p. 3175–3184. External Links: Link, Document Cited by: §1. Hargreaves et al. (2020) C. J. Hargreaves, M. S. Dyer, M. W. Gaultois, V. A. Kurlin, and M. J. Rosseinsky The earth mover’s distance as a metric for the space of inorganic compositions. Chemistry of Materials 32 (24), p. 10610–10620. External Links: Document Cited by: Remark 1. Heylen et al. (2014) R. Heylen, M. Parente, and P. Gader A review of nonlinear hyperspectral unmixing methods. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing 7 (6), p. 1844–1868. Cited by: §1, §2.1, §2.3.1. Ho et al. (2017) N. Ho, X. L. Nguyen, M. Yurochkin, H. H. Bui, V. Huynh, and D. Phung Multilevel clustering via wasserstein means. In Proceedings of the 34th International Conference on Machine Learning - Volume 70, ICML’17, Sydney, NSW, Australia, p. 1501–1509. Cited by: §1, §4.2. Holland et al. (1998) J. K. Holland, E. K. Kemsley, and R. H. Wilson Use of fourier transform infrared spectroscopy and partial least squares regression for the detection of adulteration of strawberry purées. Journal of the Science of Food and Agriculture 76 (2), p. 263–269. External Links: Document Cited by: §1. Hosseini-Asl et al. (2015) E. Hosseini-Asl, J. M. Zurada, and A. El-Baz Automatic segmentation of pathological lung using incremental nonnegative matrix factorization. In 2015 IEEE International Conference on Image Processing (ICIP), Vol. , p. 3111–3115. Cited by: §1. Lee and Seung (2001) D. Lee and H. Seung Algorithms for non-negative matrix factorization. Adv. Neural Inform. Process. Syst. 13, p. . Cited by: §2.3. Liu et al. (2018) J. Liu, B. Luo, S. Douté, and J. Chanussot Exploration of planetary hyperspectral images with unsupervised spectral unmixing: a case study of planet mars. Remote Sensing 10 (5). External Links: Link, ISSN 2072-4292, Document Cited by: §1. Lloyd (1982) S. Lloyd Least squares quantization in pcm. IEEE Transactions on Information Theory 28 (2), p. 129–137. External Links: Document Cited by: Appendix B, §2.2, §4.2. Ma et al. (2014) W. Ma, J. M. Bioucas-Dias, T. Chan, N. Gillis, P. Gader, A. J. Plaza, A. Ambikapathi, and C. Chi A signal processing perspective on hyperspectral unmixing: insights from remote sensing. IEEE Signal Processing Magazine 31 (1), p. 67–81. Cited by: §1, §2.3.1. Mu et al. (2012) T. Mu, J. Y. Goulermas, J. Tsujii, and S. Ananiadou Proximity-based frameworks for generating embeddings from multi-output data. IEEE Transactions on Pattern Analysis and Machine Intelligence 34 (11), p. 2216–2232. External Links: Document Cited by: §2.3. Nakhostin et al. (2016) S. Nakhostin, H. Clenet, T. Corpetti, and N. Courty Joint anomaly detection and spectral unmixing for planetary hyperspectral images. IEEE Transactions on Geoscience and Remote Sensing 54 (12), p. 6879–6894. External Links: Document Cited by: §1. Peyré and Cuturi (2019) G. Peyré and M. Cuturi Computational optimal transport. Foundations and Trends in Machine Learning 11 (5-6), p. 355–602. External Links: Link Cited by: §2.2, Definition 9, Remark 1. Pollard (1982) D. Pollard Quantization and the method of k-means. IEEE Transactions on Information Theory 28 (2), p. 199–205. External Links: Document Cited by: §2.2, §4.2. Rameau, J. et al. (2021) Rameau, J., Chanussot, J., Carlotti, A., Bonnefoy, M., and Delorme, P. Spectral unmixing for exoplanet direct detection in hyperspectral data. A&A 649, p. A143. External Links: Document, Link Cited by: §1. Rieger (2006) M. O. Rieger Monotonicity of transport plans. Note: cvgmt preprint External Links: Link Cited by: §4.1.1. Rubner et al. (2000) Y. Rubner, C. Tomasi, and L. Guibas The earth mover’s distance as a metric for image retrieval. International Journal of Computer Vision 40, p. 99–121. External Links: Document Cited by: §4.2. Sandler and Lindenbaum (2011) R. Sandler and M. Lindenbaum Nonnegative matrix factorization with earth mover’s distance metric for image analysis. IEEE Transactions on Pattern Analysis and Machine Intelligence 33 (8), p. 1590–1602. Cited by: §2.3, §3.2. Santambrogio (2015) F. Santambrogio Optimal transport for applied mathematicians. Springer International Publishing, Cham. External Links: Document Cited by: §2.2. Schmitz et al. (2018) M. A. Schmitz, M. Heitz, N. Bonneel, F. Ngolè, D. Coeurjolly, M. Cuturi, G. Peyré, and J. Starck Wasserstein dictionary learning: optimal transport-based unsupervised nonlinear dictionary learning. SIAM Journal on Imaging Sciences 11 (1), p. 643–678. External Links: Document Cited by: §2.1, §2.3. Schwartzburg et al. (2014) Y. Schwartzburg, R. Testuz, A. Tagliasacchi, and M. Pauly High-contrast computational caustic design. ACM Trans. Graph. 33 (4). External Links: ISSN 0730-0301, Link, Document Cited by: Remark 1. Shirdhonkar and Jacobs (2008) S. Shirdhonkar and D. W. Jacobs Approximate earth mover’s distance in linear time. In 2008 IEEE Conference on Computer Vision and Pattern Recognition, Vol. , p. 1–8. External Links: Document Cited by: §4.2. Smaragdis (2004) P. Smaragdis Non-negative matrix factor deconvolution; extraction of multiple sound sources from monophonic inputs. In Independent Component Analysis and Blind Signal Separation, C. G. Puntonet and A. Prieto (Eds.), Berlin, Heidelberg, p. 494–499. External Links: ISBN 978-3-540-30110-3 Cited by: §1, footnote 5. Stein (2003) D. Stein Application of the normal compositional model to the analysis of hyperspectral imagery. In IEEE Workshop on Advances in Techniques for Analysis of Remotely Sensed Data, 2003, Vol. , p. 44–51. External Links: Document Cited by: §2.3.1. Tarantola (2005) A. Tarantola Inverse problem theory and methods for model parameter estimation. Society for Industrial and Applied Mathematics, Philadelphia. External Links: Document Cited by: §1, §2.1. Villani (2003) C. Villani Topics in optimal transportation. In Graduate Studies in Mathematics, External Links: Link Cited by: Appendix D. Villani (2009) C. Villani Optimal transport. Springer Berlin Heidelberg, Berlin, Heidelberg. External Links: Document Cited by: §2.2. Vito et al. (2005) E. D. Vito, L. Rosasco, A. Caponnetto, U. D. Giovannini, and F. Odone Learning from examples as an inverse problem. Journal of Machine Learning Research 6 (30), p. 883–904. External Links: Link Cited by: §1, §2.1. Yan et al. (2019) J. Yan, C. Deng, L. Luo, X. Wang, X. Yao, L. Shen, and H. Huang Identifying imaging markers for predicting cognitive assessments using wasserstein distances based matrix regression. Frontiers in Neuroscience 13, p. 668. External Links: Link, Document, ISSN 1662-453X Cited by: §1. Zhao et al. (2019) M. Zhao, J. Chen, and Z. He A laboratory-created dataset with ground truth for hyperspectral unmixing evaluation. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing 12 (7), p. 2170–2183. External Links: Document Cited by: §5.3. Zhou et al. (2018) Y. Zhou, A. Rangarajan, and P. D. Gader A gaussian mixture model representation of endmember variability in hyperspectral unmixing. IEEE Transactions on Image Processing 27 (5), p. 2242–2256. External Links: Document Cited by: §2.3.1.