{'fraction_non_alphanumeric': 0.030461847389558233, 'fraction_numerical': 0.012931726907630522, 'mean_word_length': 4.996508127633955, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Differentiable Programming for Hyperspectral Unmixing using a Physics-based Dispersion Model

Hyperspectral unmixing is an important remote sensing task with applications including material identification and analysis. Characteristic spectral features make many pure materials identifiable from their visible-to-infrared spectra, but quantifying their presence within a mixture is a challenging task due to nonlinearities and factors of variation. In this paper, spectral variation is considered from a physics-based approach and incorporated into an end-to-end spectral unmixing algorithm via differentiable programming. The dispersion model is introduced to simulate realistic spectral variation, and an efficient method to fit the parameters is presented. Then, this dispersion model is utilized as a generative model within an analysis-by-synthesis spectral unmixing algorithm. Further, a technique for inverse rendering using a convolutional neural network to predict parameters of the generative model is introduced to enhance performance and speed when training data is available. Results achieve state-of-the-art on both infrared and visible-to-near-infrared (VNIR) datasets, and show promise for the synergy between physics-based models and deep learning in hyperspectral unmixing in the future.


Introduction
Hyperspectral imaging is a method of imaging where light radiance is densely sampled at multiple wavelengths. Increasing spectral resolution beyond a traditional camera's red, green, and blue spectral bands typically requires more expensive detectors, optics, and/or lowered spatial resolution. However, hyperspectral imaging has demonstrated its utility in computer vision, biomedical imaging, and remote sensing [42,8]. In particular, spectral information is critically important for understanding material reflectance and emission properties, important for recognizing materials.
Spectral unmixing is a specific task within hyperspectral imaging with application to many land classification problems related to ecology, hydrology, and arXiv:2007.05996v1 [cs.CV] 12 Jul 2020 mineralogy [29,32]. It is particularly useful for analyzing aerial images from aircraft or spacecraft to map the abundance of materials in a region of interest. While pure materials have characteristic spectral features, mixtures require algorithms to identify and quantify material presence.
A common model for this problem is linear mixing, which assumes electromagnetic waves produced from pure materials combine linearly and are scaled by the material abundance. Mathematically this is expressed as b = Ax + η where b is the observed spectra, A is a matrix whose columns are the pure material spectra, η is the measurement noise, and x is the abundance of each pure material. The model assumes that the pure material spectra, referred to as endmember spectra, is known before-hand. Nonlinear effects are known to occur when photons interact with multiple materials within a scene for which we refer readers to the review by Heylen et al. for techniques to account for the non-linear mixing [29].
A key challenge that affects both linear and nonlinear mixing models is that pure materials have an inherent variability in their spectral signatures, and thus cannot be represented by a single characteristic spectrum. Spectral variability of endmembers is caused by subtle absorption band differences due to factors such as different grain sizes [43,52,47,48] or differing ratios of molecular bonds [7,56] as shown in Figure 1. Since variability causes significant errors in unmixing algorithms, it is an active area of research [66,18,71]. Recently, differentiable programming has become a popular research area due to its potential to bridge gaps between physics-based and machine learningbased techniques for computer vision and graphics [59,25,1]. Our key insight is to leverage differentiable programming by modelling the variation of spectra with a physics-based dispersion model, and incorporating this differentiable model into an end-to-end spectral unmixing algorithm. Such an approach has the capacity to unmix scenes with a large amount of variability, while constraining the predictions to be physically plausible. These physically plausible variations of endmember spectra also provide additional science data as the variation of absorption bands can reveal properties about the composition and history of the material. To our knowledge, we are the first to use a generative physics model to account for spectral variability in an unmixing algorithm.
Contributions: Our specific contributions in this paper are the following: -We introduce a physics-based dispersion model (first presented in [53,34,61]) to generate and render spectral variation for various pure materials. We provide an efficient optimization method via gradient descent to find dispersion model parameters for this spectral variation. -We incorporate this dispersion model into an end-to-end spectral unmixing algorithm utilizing differentiable programming to perform analysis-bysynthesis optimization. Analysis-by-synthesis is solved via alternating minimization optimization and requires no training data. -We further design an inverse rendering algorithm consisting of a convolutional neural network to jointly estimate dispersion model parameters and mineral abundances for spectral unmixing. This method requires training data, but is computationally efficient at test time and outperforms analysisby-synthesis and other state-of-the-art methods.
We provide extensive analysis of our proposed methods with respect to noise and convergence criteria. To validate our contributions, we test on both synthetic and real datasets using hyperspectral observations in the visible and near infrared (VNIR), and mid to far infrared (IR). The datasets also span three different environments from laboratory, aircraft, and satellite based spectrometers. Our methods achieve state-of-the-art across all datasets, and we compare against several baselines from literature. Our code is openly available and accessible here: https://github.com/johnjaniczek/InfraRender. We hope this work inspires more fusion between physics models and machine learning for hyperspectral imaging and computer vision more generally in the future.

Related Work
Optimization-based Approaches. Standard optimization techniques for linear unmixing include projection, non-negative least squares, weighted least squares, and interior point methods [47,28,51,13,27]. Further, sparsity-based optimization can improve abundance prediction [69,9]. However, most optimization have not leveraged physics priors as we do in our model. Spectral Variability. Spectral variability has been a topic of recent interest [66,6]. One approach is to augment A with multiple variations or spectra for each endmember. To do this, multiple endmember spectral mixture analysis (MESMA) [50] and multiple-endmember linear spectral unmixing (MEL-SUM) [17] both require labeled data of the spectral variation for each endmember. In contrast, unsupervised techniques learn endmember sets from unlabelled hyperspectral images, including semi-automated techniques [3], k-means clustering [3], and the sparsity promoting iterated constrained endmember algorithm (SPICE) [64,65] which simultaneously finds endmember sets while unmixing for material abundances. These techniques are limited by the amount of sets in the endmember library, and computational complexity increases with more additions. Our method by contrast finds an efficient parameter set to physically model the spectral variation.
Another category of endmember variability techniques models the endmember spectral variation as samples from a multivariate distribution P(e|θ) where e is the endmember spectra, and θ are the distribution parameters. Common statistical distributions proposed include the normal compositional model [54], Gaussian mixture models [71], and the beta compositional model [18]. These distribution models have large capacity to model spectral variations, however sometimes they can render endmember spectra that are not physically realistic.
Deep Learning for Hyperspectral Classification and Unmixing. Deep learning has recently improved many hyperspectral imaging tasks [68,37]. In particular, networks process hyperspectral pixel vectors using both deep belief networks [40] and CNNs [30,10,62]. For spatial hyperspectral data, CNNs [12], joint spectral-spatial feature extraction [70], and 3D CNNs [39] are used. All these methods require large hyperspectral datasets that are annotated correctly. One of our methods uses analysis-by-synthesis and differentiable programming to avoid low training data issues, but our technique can be made complementary to deep learning architectures as we show in our inverse rendering CNN.
Differentiable Programming and Rendering. Differentiable programming refers to the paradigm of writing algorithms which can be fully differentiated end-to-end using automatic differentiation for any parameter [59,60,4]. This has been applied for audio [20] and 3D geometry processing [49]. In graphics, differentiable rendering has improved ray tracing [38,44,41,67], solved analysisby-synthesis problems in volumetric scattering [26,25], estimated reflectance and lighting [1], and performed 3D reconstruction [58]. In our paper, we write a forward imaging model utilizing the physics of dispersion in spectral variation to allow our pipeline to be differentiable end-to-end.

Method
Our approach to hyperspectral unmixing features two main components: (1) use of a physically-accurate dispersion model for pure endmember spectra, and (2) a differentiable programming pipeline to perform spectral unmixing. This approach has synergistic benefits of leveraging prior domain knowledge while learning from data. Our first algorithm solves spectral unmixing in a self-supervised fashion using analysis-by-synthesis optimization with the dispersion model as the synthesis step. Further, we show how inverse rendering via a convolutional neural network (CNN) can learn parameters of this model to help speed up our end-to-end pipeline and improves performance when training data is available.

Dispersion Model
We first describe the dispersion model for generating endmember spectra. Endmember and/or endmember spectra is what we call the spectral curve for emissivity as a function of wavenumber ω. Each pure material has a characteristic endmember spectrum, although spectra can vary, which is the problem we are trying to solve/disambiguate. Let ε measured (ω) be endmember spectra we have measured, typically in a lab or in the field, whose emissivity is sampled at different wavenumbers: ε measured (ω 1 ), · · · ε measured (ω N ) T . Our goal is to propose a model ε model (Λ; ω) with parameters Λ such that the following loss is minimized: That is, we fit the model emissivity of an endmember spectrum to the measured spectrum. In practice, we need to add regularization and constraints to this endmember loss for better fitting which we describe after the derivation of the dispersion model.
Derivation of the Dispersion model: Our model of endmember spectra is derived from an atomistic oscillator driven by electromagnetic waves impinging on the molecular structure of the pure material [53,34]. In Figure 2, we show a conceptual diagram of this model, and how it generates emissivity curves as a function of wavelength. For the full derivation of the model from first principles, we refer the reader to Appendix A in the supplementary material. Instead, we outline the model below based on the equations derived from that analysis.  Figure: The insight of the dispersion model is that optical properties can be related to molecular structure through first principles via an atomistic oscillator model. We use this generative model for the formation of spectral variation in our spectral unmixing algorithm Let Λ = [ρ, ω o , γ, r ] be a matrix of parameters, where ρ, ω 0 , γ, r ∈ R K and K is a model hyperparameter corresponding to the number of distinct massspring equations used to model the emissivity. ρ is the band strength, ω o is the resonant frequency, γ is the frictional force (dampening coefficient), and r is relative dielectric permeability. Please see the supplemental material for the physical significance of these parameters to the atomistic oscillator model, and their control over the shape of spectral absorption bands. Note: usually r is a constant vector which does not vary with K. Thus Λ ∈ R K×4 . The refractive index terms n, k are given as follows [53,34]: where the expressions for θ, b, φ are given as follows: We note that subscript k denotes the k-th coordinate of the corresponding vector. Also there is another useful relation (derived in Appendix A) that n 2 − k 2 = θ, nk = φ. We then define the complex refractive index asn(Λ; ω) = n(Λ; ω) − i · k(Λ; ω), where i = √ −1 is the imaginary number. Hence, we can calculate the emissivity as follows: When considering minerals, we introduce M ∈ N, the number of optical axes of symmetry in crystal structures, (eg. 2 axes of symmetry in quartz [53,61]), to define the full model: where we use a different parameter matrix Λ m and weight α m for each optical axis of symmetry. The dispersion model has been primarily used to analyze optical properties of materials to determine n and k, which then can be subsequently applied to optical models like radiative transfer [53,61]. After n and k are found, spectra such as reflectance, emissivity, and transmissivity can be generated. In particular, we notice that fine-grained control of the dispersion model parameters can realistically render spectral variation that occurs in hyperspectral data. Our contribution is to leverage these properties in a differentiable programming pipeline for spectral unmixing.
Endmember Fitting. Using the dispersion model presented above, we want to robustly estimate the model parameters to fit the spectra ε measured captured in a lab or in the field. To fit the model, we wish to perform gradient descent to efficiently find these parameters. Using chain rule on the loss function, we see that ∂L ∂Λij = ∂L ∂ε model ∂ε model ∂Λij , where (ij) corresponds to that element of the parameter matrix, and all expressions are scalars once the coordinate is specified. While the partial derivatives can be calculated explicitly via symbolic toolboxes, the resulting expressions are too long to be presented here. For simplicity and ease of use, we use the autograd function [4,45] in PyTorch [46] to automatically compute derivatives for our model as we are performing backpropagation.
One main challenge in performing endmember fitting is that the dispersion model is not an injective function, and hence is typically not identifiable, that is more than one Λ can result in the same fit. This can be solved, in-part, through regularization to enforce sparsity, especially since a preference for fewer dispersion parameters has been suggested in the literature [53,61]. In our implementation, we initialize our model with K = 50 rows of the parameter matrix. Since the parameter ρ controls the strength of the absorption band, small values of ρ do not contribute much energy to the spectra (unnecessary absorption bands), and can be pruned. After performing sparse regression by penalizing the L 1 norm of ρ, K is typically around 10-15 in our experiments.
Thus, our modified sparse regression problem may be written as arg min where Λ min and Λ max restrict the variation of the dispersion parameters to a plausible range. In addition, endmembers (particularly minerals) can have multiple optical axes of symmetry described by separate spectra, which has been noted in the literature [53,61]. Without prior knowledge of the number of axes for every material we encounter, we run this optimization for a single and double axes, and pick the one with the lowest error. See Section 4 for results on endmember fitting and Figure 6 for examples of modelled vs. measured spectra. Despite the fact that this regression problem is non-convex, we solve it using gradient descent with a random initialization; this is known to converge to a local minimum with probability 1 [35]. A global minimum is not necessary at this stage, since we use endmember fitting to provide a good initialization point for the subsequent alternating minimization procedure introduced in the next subsection.

Differentiable Programming for End-to-End Spectral Unmixing
Analysis-by-Synthesis Optimization. In Figure 3, we show our full endto-end spectral unmixing pipeline. Here, ε model (Λ; ω), which is initially fit to ε measured , is then aggregated into the columns of A. Then, the observed spectra b is linearly unmixed by solving a regularized least-squares optimization: 2 +λ x p subject to sum-to-one and non-negativity constraints x 1 = 1, x ≥ 0. Given these constraints, one cannot impose sparsity with the usual L 1 norm. Instead, we use the L p norm to induce sparsity for the predicted abundances; this has been proposed before for spectral unmixing [9].
The key to our pipeline is that everything is fully differentiable, and thus we can actually minimize the following equation: arg min Fig. 3: Analysis-by-Synthesis: The analysis-by-synthesis algorithm uses differentiable programming to find optimal dispersion parameters and abundances. The initial dispersion parameters and the target spectra are fed as inputs, and the algorithm alternates between optimizing the abundances in the least squares sense and updating the dispersion parameters with respect to the gradient with respect to both the parameters of the dispersion model Λ and the abundances x. This gives us our recipe for hyperspectral unmixing: first, perform endmember fitting to initialize A(Λ), then, solve Equation 7 in an alternating fashion for x and Λ. One could also solve this equation jointly for both unknowns, however, we found that the alternating optimization was faster and converged to better results. The optimization problem established in equation (7) is an alternating minimization problem and is unfortunately not convex [31]. One popular approach to tackle nonconvex problems is to find a good initialization point [19,5], and then execute a form of gradient descent. Inspired by this, we first initialize A(Λ) by performing endmember fitting using Equation 6 as described in the previous subsection. Our experiments indicate that this provides a useful initialization for our subsequent step. We then perform alternating minimization on Equation 7 for x and Λ. Note that each iterate of the resulting alternating minimization involves the solution of a subproblem which has a convergence rate which depends on the condition number of the matrix A(Λ). For more details on this, we refer the reader to Appendix C where we discuss on the properties of A(Λ) across multiple runs.
In the ideal scenario, this initial matrix A(Λ) would consist of the endmember spectra that fully characterizes the mixed spectra b. However, since spectra for the same material can significantly vary [43,52,47,48,7,56] (see Figure 1), the initialization can be slightly off and we follow up with (7) to obtain a better fit. Note that this optimization problem is solving for the maximum likelihood estimator under a Gaussian noise model. Our optimization technique is performing analysis-by-synthesis, as given a single observation b, the dispersion model synthesizes endmember variation until a good fit is achieved.
Inverse Rendering of Dispersion Model Parameters. The previous analysis-by-synthesis optimization does not require training data (labeled abun-dances in spectral mixtures) in order to perform spectral unmixing. However, there is room for even more improvement by using labeled data to help improve the parameter fitting of the model in the synthesis step. We train a CNN to predict the parameters for a generative model, known as inverse rendering in other domains [63]. In Figure 4, we show this inverse rendering conceptually, and how it can be fed into our differentiable programming pipeline for end-to-end spectral unmixing.
Our CNN architecture consists of convolutional layers followed by a series of fully-connected layers. We refer the reader to the supplemental material for the exact network structure and implementation details. Using a CNN for inverse rendering is significantly faster at test time as compared to the analysis-bysynthesis optimization. However, it does have a drawback of requiring training data which is unavailable for certain tasks/datasets. Fig. 4: Inverse Rendering: A CNN is trained to "inversely render" pixels of the hyperspectral image, by predicting both the dispersion parameters that control the spectral variability, and the abundances that control the mixing model. During training, the reconstruction error is back-propagated through the differentiable dispersion model to boost the performance of the network at making physically realistic predictions

Experimental Results
Datasets. We utilize three separate datasets to validate our spectral unmixing algorithms. In Figure 5, we visually represent these datasets and their exemplar data. For specific implementation details and dataset pre-processing, please see Appendix D in the supplemental material.
Feely et al. Dataset. We utilize 90 samples from the Feely et al. dataset [21] of thermal emission spectra in the infrared for various minerals measured in the lab. Ground truth was determined via optical petrography [21], and a labeled endmember library is provided. The limited amount of data is challenging for machine learning methods, so we utilize the dispersion model to generate 50,000 additional synthetic spectra for dataset augmentation.
Gulfport dataset. The Gulfport dataset from Gader et al. [22] contains hyperspectral aerial images in the VNIR along with ground truth classification labels segmenting pixels into land types (e.g. grass, road, building). Although the dataset is for spectral classification, it can also be used to benchmark unmixing algorithms by creating synthetic mixtures of pure pixels from the Gulfport dataset with random abundances as done by [18,71]. We perform both spectral classification (Gulfport) and unmixing (Gulfport synthetic) tasks in our testing. Both datasets are split into a train and test set (although some methods do not require training data), and the training data is augmented with 50,000 synthetically generated mixtures from the dispersion model.
One main difficulty of this dataset is the endmembers identified correspond to coarse materials such as grass and road as opposed to pure materials. Such endmembers can significantly vary across multiple pixels, but this spectral variation is not physically described by the dispersion model. To solve this problem, we utilize K-means clustering to learn examplar endmembers for each category (e.g. grass, road, etc). Then the resulting centroid endmember can be fit to the dispersion model to allow further variation such as absorption band shifts in the spectra. We found that K = 5 worked the best for the Gulfport dataset.
TES Martian Dataset. The Thermal Emission Spectrometer (TES) [14] uses Fourier Transform Infrared Spectroscopy to measure the Martian surface. We utilize pre-processing from Bandfield et al. [2], and the endmember library used by Rogers et al. to analyze Mars TES data [51]. There is no ground truth for this dataset, as the true abundance of minerals on the Martian surface is unknown, so other metrics such as reconstruction error of the spectra are considered. Fig. 5: Datasets: This figure shows representative data and instrumentation for the three datasets considered in this paper. Data includes laboratory, aircraft, and satellite measurements, and ground truth ranges from detailed abundance analysis under a microscope (Feely [21]) to pure pixel labels of land type for spectral classification (Gulfport [22]) to no ground truth for the Martian data (TES [14]) Baselines. We compare against several state-of-the-art baselines in the literature. The basic linear unmixing algorithm is Fully Constrained Least Squares (FCLS) [28] which solves least squares with sum-to-one and non-negativity constraints on the abundances. We also implement two state-of-the-art statistical methods for modelling endmember variability as distributions: the Normal Compositional Model (NCM) [54] and the Beta Compositional Model (BCM) [18]. NCM and BCM use a Gaussian and Beta distribution respectively, perform expectation-maximization for unmixing, and require a small amount of training data to determine model parameters.
We also compare against two state-of-the-art deep learning networks by Zhang et al. [69]. The first network utilizes a 1D CNN (CNN-1D) architecture, while the second network utilizes a 3D CNN (but with 1D convolutional kernels) (CNN-3D). CNN-3D is only applicable to datasets with spatial information, and not testable on the Feely and Gulfport synthetic data. We further created a modified CNN architecture (CNN-1D Modified) to maximize the performance on our datasets by changing the loss function to MSE, removing max-pooling layers, and adding an additional fully connected layer before the output. In the supplemental material, we provide information about the parameters, network architectures, and training procedures we used for these baselines as well as details for our own methods. We also have all of our code available here: https://github.com/johnjaniczek/InfraRender. (Right) Cluster centroids found for pixels labelled as grass in the Gulfport dataset, and the model fit to these centroids. Note the high fidelity of fit via the dispersion model for both these cases Endmember Fitting Results. To bootstrap both the analysis-by-synthesis and inverse rendering algorithms, good initial conditions for the dispersion parameters need to be input to the model. Determining dispersion parameters typically required detailed molecular structure analysis or exhaustive parameter searching methods [53,61,36]. One main advantage of our method is that we utilize gradient descent to efficiently find parameter sets for different materials. In Appendix B of the supplemental material, we share some of these parameter sets and our insights using the dispersion model for the scientific community.
In Figure 6, we show qualitative results of our endmember fitting by minimizing the loss in Equation 6 using gradient descent. The reconstructed spectra achieves a low MSE with the measured spectra with an average MSE of 0.016 for the TES library, 0.0019 for the Feely library, and an MSE of 2.6e-5 on the Gulfport cluster centroids. Note that there is noise in the measurements, and so MSE is not an absolute metric of the fit to the true unknown spectra.
Spectral Unmixing Results. In Table 1, we show results on the Feely, Gulfport, and the Gulfport synthetic mixture datasets. For Feely, the analysisby-synthesis method achieved a MSE of 0.052, with the next closest method (NCM) achieving 0.119. Due to the Feely dataset only containing 90 test samples, the machine learning methods were trained on synthetic data which explains their lower performance as data mismatch. Thus, the low error of analysis-bysynthesis shows the utility of the dispersion model for modelling endmember variability, particularly in cases with low training data.
For the Gulfport dataset, the task was to predict the material present since the labeled data is for single coarse materials (e.g. road, grass, etc) at 100% abundance per pixel. Here, the deep learning methods of CNNs and Inverse Rendering have the highest performance. This is expected as there exists a large amount of training data to learn from. Note that Inverse Rendering performs the best at 0.272 MSE, demonstrating that the addition of a generative dispersion model to the output of the CNN improves performance over purely learned approaches. Also note that our analysis-by-synthesis method still has relatively high performance (0.45 MSE) without using any training data at all.
For the Gulfport synthetic mixture dataset, Inverse Rendering achieves the lowest MSE of 0.059, leveraging both physics-based modeling for spectral mixing as well as learns from available training data. The BCM and the analysis-bysynthesis methods both outperform the CNN methods, even though they do not have access to the training data. In fact, BCM even slightly outperforms the analysis-by-synthesis method, which could be because the sources of variation in this data are well-described by statistical distributions. Speed of Methods. The additional capacity of adding statistical and physical models usually has a cost of speed in implementation. Averaged over 90 mixtures, the convergence for a single operation was FCLS -10ms, BCM -1.23s, NCM -18ms, CNN -33ms, Inverse Rendering -39ms, and analysis-by-synthesis -10.2s. Future work could potentially increase the speed of analysis-by-synthesis with parallel processing.
Noise analysis. Prior to spectral unmixing, emissivity is separated from radiance by dividing out the black-body radiation curve at the estimated temperature [47,16]. In general, a Gaussian noise profile in the radiance space with variance σ 2 radiance results in wavenumber dependent noise source in the emissivity space with the profile σ 2 (ω) = σ 2 radiance · 1/B(ω, T ) where B is the black-body function given by Planck's law. In our noise experiments we use a black-body radiation curve for a 330K target, which is the approximate temperature the Feely dataset samples were held at. In Figure 7 left, we see that the emissivity noise is higher where the radiance signal is lower.
We simulated varying the noise power to determine the methods' robustness tested on 30 samples from the Feely dataset. In Figure 7, you can see that analysis-by-synthesis still has the best performance in the presence of noise, and is relatively flat as noise increases compared to other methods. We note that statistical methods, while having higher average error, seem to be robust to increased noise as they can handle random perturbations of each spectral band statistically. CNN and Inverse Rendering methods perform the worst for high noise, as these methods were trained on data without noise. Fig. 7: The left plot shows the radiance profile of a spectra perturbed by Gaussian noise and the resulting emissivity profile after separating out the blackbody radiance. The right figure shows the robustness of the algorithms to increasing amounts of noise TES Data. The Mars TES data was unmixed using our analysis-by-synthesis method to demonstrate it's utility on tasks where zero training data is available. The method produces mineral maps which correctly finds abundances of the mineral hematite at Meridiani Planum in Figure 8. This is an important Martian mineral which provides evidence for liquid water having existed at some point on Mars, and has been verified by NASA's Opportunity Rover [33]. Note how FCLS predicts many sites for hematite, while our method narrows down potential sites on the Martian surface, which is useful for planetary scientists. By allowing for spectral variation through our physics-based approach, our method has lower RMS reconstruction error than previous analysis of TES data. FCLS, which was previously used on TES because of the zero training-data problem, has an average RMS reconstruction error of 0.0043 while analysis-by-synthesis has an average of 0.0038. This is an exciting result as our methods could provide a new suite of hyperspectral analysis tools for scientists studying the Martian surface.

Discussion
This paper incorporated generative physics models into spectral unmixing algorithms via differentiable programming. We adopt a physics-based dispersion model to simulate spectral variability, and show how this model can realistically fit several real measured spectra via gradient descent. We further show how to jointly optimize for the dispersion parameters and material abundances with an analysis-by-synthesis optimization. A second algorithm is introduced for tasks where additional data is available by training a CNN to "Inversely Render" a hyperspectral image with the differentiable dispersion model in the loop.
We validate these contributions extensively on three datasets ranging from mid to far IR and VNIR, and compared against state-of-the-art optimization, statistical and deep learning benchmarks. From these experiments we observe that analysis-by-synthesis has the best performance when training data is not available, and that Inverse Rendering has the best performance when training data is available. We also see that analysis-by-synthesis is noise resilient, and reconstructs Mars spectra with lower error than previous techniques.
There are still limitations for the methods proposed. First, analysis-by-synthesis has a large computational cost compared to other methods, although this could be mitigated through parallelization of the algorithm. Secondly, the spectral unmixing community is limited by the training data available. This is difficult to overcome, because it is expensive to produce quality datasets, and it is not easy for experts to label remote sensing datasets from prior knowledge alone. Future work could investigate generating realistic synthetic data suitable for training machine learning based algorithms for better performance. We hope using generative physics-based models inspires others to produce realistic synthetic data as well as differentiable programming methods which require low training data.
Acknowledgements. This work was supported by NSF grant IIS-1909192 as well as GPU resources from ASU Research Computing. We thank Dr. Alina Zare, Christopher Haberle, and Dr. Deanna Rogers for their helpful discussions, and Kim Murray (formerly Kim Feely) for providing the laboratory measurements and analysis contributing to this paper. In this section, we derive the dispersion model from first principles, modeling the generation of spectral radiance as a sum of mass-spring oscillations driven by an electromagnetic wave. This induces changes in the index of refraction, which governs the reflectance of the material with respect to light wavelength/frequency. This is based on earlier work by Garbuny and by Spitzer et al. [23,53]. We first start with the equation for a mass-spring oscillator driven by an external force: For a charged particle, F = qE, where q is charge and E = E 0 e iωt for a propagating electromagnetic wave. Thus we can substitute these in to get: which has the solution: where ω 2 0 = G/m and γ = R/m. At the same time, we can also relate x to E via the band strength: where α is the polarizability. Using the identity = 1 + 4πα, we can derive the following band strength equation: Combining Eq. 3 and Eq. 5, we get Solving for : Relating to the refractive index, µ =n 2 where µ = 1 andn = n − ik, we get This yields the refractive index equations: Using the Lorentz-Lorenz formula, we can get where ω 2 1 = ω 2 0 − 4πq 2 3m where ω 1 < ω 0 . So plugging in ω 1 for ω 0 yields: This entire derivation was for a single oscillator, but in practice, there are multiple oscillators that interact. We write this as a linear superposition given as follows: where f i is the strength of each individual oscillator. Using these equations, we have two equations for two unknowns (n and k), which we showed in Section 3 of the main paper is the basis of calculating reflectance and emission.

Appendix B: Spectral Variation
As stated in the main paper, the motivation for incorporating the dispersion model into a differentiable program for spectral unmixing is to allow for physically plausible spectral variation of pure materials. Because it is known that absorption bands can shift in frequency and strength between different endmember samples, the goal was to find a model that described these changes in a physically plausible way. That is we wanted a generative model for the formation of spectra with parameters that have "dials" to tune the frequency, strength, and shape of absorption bands. From the literature on analysis of the formation of spectra from an atomistic perspective [53,61] we find that the Lorentz-Lorenz dispersion model is the correct approach to take. However, unlike previous works we go further than using the model to derive optical properties of materials, we also incorporate the dispersion model into an end-to-end spectral unmixing pipeline that allows the parameters to be fine-tuned via differentiable programming to account for spectral variability.
In Appendix A, the dispersion model is derived from first principles and each absorption band is described by the parameters ρ, ω o , γ, r . ρ is the band strength and as it increases the absorption band becomes deeper. ω o is the resonant frequency and as it increases the absorption band shifts in wavenumber (and also slightly shifts the shape). γ is the frictional force (dampening coefficient) and controls the shape/width of the absorption bands. r is relative dielectric permeability and as it increases the entire emissivity curve is shifted downwards. Also note that absorption bands which are close to each other interact in highly non-linear ways. Fig. 1: A single absorption band is initialized with r = 2.356, ω 0 = 1161, γ = 0.1, ρ = 0.67. Then the parameters are perturbed such that ω 0 is increased by 100, γ is increased to 120%, and ρ is increased to 120%. The plots show the effect of changing each parameter individually to show it's control over the shape and width of the absorption band.
The importance of initializing the alternating optimization with good initial dispersion parameters was emphasized in the main paper, as the problem is non-convex and good initialization is essential. It also makes intuitive sense to initialize with parameters fit to an endmember sample to give physical significance to the generated spectra. As shown in the results of the main paper, we achieve good fits with low MSE on endmember libraries used to analyze the Martian surface as well as a semi-urban university scene. The endmember libraries used to fit the minerals to analyze the Mars TES data are of high quality from the Arizona State University Thermal Emission Spectral Library [15]. The resulting parameters from a few of the important materials from this endmember library are provided in the following tables.

Appendix C: Alternating Optimization Convergence
Our ultimate goal is to solve the spectral unmixing problem which can be formulated as min A,x b − Ax 2 2 , where the minimization occurs over both the matrix A and the unmixing vector x. This is a standard case of alternate minimization which is known to be nonconvex [31]. In practice, alternating minimzation are particularly hard to tackle due to the presence of suboptimal local minimas.  Recent progress on tackling nonconvex problems involves either characterizing the optimization landscape [24,57,55] or providing initialization to descent algorithms [5,11] to assure convergence to the global optimum. It is known that gradient descent applied to alternate minimization problem faces the issue of getting stuck at local minimas [31] and hence initialization plays an important role in solving Equation 7 in the main paper. With that in mind, in this paper we provide a mechanism to provide good initialization to gradient descent algorithm with the hope of tackling the alternating minimization problem effectively.
Initialization using dictionary A( model ): We investigate the properties of matrix A as relates to the convergence of the alternating optimization. We denote the measured the emissivity spectrum of various materials in the lab as measured , and the physics-based dispersion model as model . We then use these emissivity spectrum to construct a dictionary A( ) which servers as the initialization for A in the alternating minimization approach in Equation 7 of the main paper. The intuition behind this revolves around the ability of matrix A as a dictionary of known emissivity spectra and we expect that the unknown spectra unknown would be described as a linear combination of columns from matrix A. Consider the following subproblem of the alternating minimization: without the regularization terms. In order to ensure the uniqueness of the solution x * , we need to ensure that the matrix A is full rank. The rate of convergence for the above minimization is inversely dependent on the condition number of the matrix A(Λ). While it is difficult to analyze this matrix analytically, we perform an experimental characterization of the rank, condition number, and eigenvalues of the matrix for several different runs of the optimization algorithm with random initializations. From the plots, we can note that the matrix A(Λ) has full rank with condition number of around 6000. The minimum and maximum eigenvalues are not showing drastic difference which goes well with our motive to incorporate small changes using alternate minimization to fit the spectral differences due to geographic differences. The high condition number is the reason for the relatively slow performance for running the alternating minimization framework, with our method taking tens of seconds to converge.

Appendix D: Implementation Details
Fully Constrained Least Squares We compare against Fully Constrained Least Squares (FCLS) [28] which is a popular classical unmixing algorithm. FCLS solves for the the (aerial) abundances, x:x = arg min x ||b − Ax|| 2 2 subject to ||x|| 1 = 1, x ≥ 0. The constraints, referred to as the sum-to-one constraint and the non-negativity constraint, are enforced since abundances are interpreted as percentages.   [54] is one of the most popular methods for modelling endmember variability via statistical methods. The method requires a small amount of training data (roughly 50 samples per endmember) to learn the mean and variance of reflec-tivity (or emissivity) of each spectral band, and modelling the variation as a Gaussian distribution. During unmixing, Expectation Maximization is used to simultaneously learn the abundances and the endmember variation, subject to the abundance sum-to-one and non-negativity constraints.
The NCM is run using the Matlab code provided by Du et al. [18]. Training data of about 50 samples of endmember variation were provided to the NCM for each dataset. There are no hyperparameters needed for this method.
Beta Compositional Model. The Beta Compositional Model (BCM) [18] is a more recent method for modelling spectral variability via a statistical method. Similar to the NCM, a small amount of training data is used to learn the beta parameters of each spectral band, and an Expectation Maximization algorithm is used during unmixing. The beta parameters allow each spectral band to be modelled as a more complex distribution than the NCM and has been shown to increase performance.
The BCM is run using the Matlab code provided by Du et al. [18]. Training data of about 50 samples of endmember variation were provided to the BCM for each dataset. For datasets without sufficient endmember samples, we generated synthetic endmember variation with the dispersion model. We search for the optimal hyperparameters through repeated experiments and report the best results. The optimal BCM across all datasets was run with K = 3, σ V = 100, and σ M = 0.001.
CNN for Spectral Classification and Unmixing. We baseline against the architecture recently proposed by Zhang et al. [69] for hyperspectral unmixing using both a 1D and 3D CNN. The main difference between the 1D and 3D CNN is that the 3D CNN is operates on 3×3 bundles of pixels while the 1D CNN is provided a single pixel. However, in both architectures the convolutional kernal is 1D and operates along the spectral dimension. Both architectures have four convolutional layers with Rectified Linear Unit (ReLU) and max-pooling nonlinear operations. These layers are followed by 2 fully connected layers, where the last layer is the output abundance predictions. The ReLU of the fully connected layers ensure non-negativity, and the sum-to-one constraint is ensured by normalizing the output layer. The network is trained to minimize −x log(x), wherê x and x are the predicted and ground truth abundances respectively. While the 1D architecture performs well on the dataset it was designed for in [69], we found that modifications were necessary to maximize the performance on our datasets. Namely, we found that training the network with respect to the mean squared error (MSE) loss function, removing the max-pooling layers, and adding an additional fully connected layer before the output improved performance of the 1D CNN. The 3D CNN is only applicable to datasets with spatial information, thus is ignored for the Feely and Synthetic datasets.
CNN Baselines [69]. The architecture from Zhang et al. was used to baseline against our method (CNN-1D and CNN-3D), as well as introducing our own modified baseline (CNN-1D modified). All CNNs were trained for 100 epochs using the Adam optimizer with learning = 0.0005, betas = (0.9, 0.999) and weight decay = 0.
CNN-1D: A 1x1 hyperspectral pixel is input into the network, four convolutional layers with alternating 1x5 and 1x4 kernels and a depth of 3, 6, 12, and 24 kernels in each layer respectively. All convolutional layers have ReLU activations a 1x2 maxpooling layer. The convolutional layers are followed by a fully connected layer with 150 hidden units, and an output fully connected layer with a size that depends on the number of abundances. Normalization is used to enforce the abundance sum-to-one constraint and the ReLU activation enforces the non-negativity constraint. The CNN is trained to minimize the log loss between the predicted and ground truth abundances. The network converges in about 100 epochs with a learning rate of 1e-3. CNN-3D: CNN-3D has an almost identical architecture, although it accepts a 3x3 set of pixels at the input. Although a spatial dimension exists at the input, the convolutions only occur in the spectral dimension. four convolutional layers with alternating 1x5 and 1x4 kernels and a depth of 16, 32, 64, and 128 kernels in each layer respectively. The convolutional layers are followed by a fully connected layer with 150 hidden units, and an output fully connected layer with a size that depends on the number of abundances. Normalization is used to enforce the abundance sum-to-one constraint and the ReLU activation enforces the non-negativity constraint. The CNN is trained to minimize the log loss between the predicted and ground truth abundances. The network converges in about 100 epochs with a learning rate of 1e-3. CNN-1D Modified: Finally, a modified version of CNN-1D is baselined against to try to find the optimal architecture for performance on our datasets. The first 2 max-pooling layers are removed, an additional hidden fully connected layer with 150 units is added before the output, and a softmax operation is applied to the output to enforce the abundance sum-to-one constraint. Also, the network is trained to minimize the mean squared error between the predicted and ground truth abundances. The network converges in about 100 epochs with a learning rate of 1e-3.
Analysis-by-Synthesis Optimization. For analysis-by-synthesis, the sparse regularization was set with p = 0.95 and λ p = 0.0001. Dispersion parameters were constrained within a tolerance of their initial conditions with ρ tol = 0.05, γ tol = 0.005, tol = 0.001, and ω tol = 0.0001. On the Gulfport datasets γ tol and tol were increased to 0.05 to compensate for increased variation. Analysis-by-Synthesis alternates between finding optimal abundances (solving a regularized least squares problem), and updating the dispersion parameters for 100 iterations using the Adam optimizer with learning rate = 0.01, betas = (0.9, 0.999), and weight decay = 0.
Inverse Rendering CNN. The inverse rendering CNN uses the same CNN architecture as CNN-1D modified. The input to the CNN is the spectrum (or batch of spectra). That is there are four convolutional layers with alternating 1x5 and 1x4 kernels and a depth of 3, 6, 12, and 24 kernels in each layer respectively. The convolutional layers also have ReLU activations and the last 2 layers have a 1x2 maxpooling layer. The convolutional layers are followed by fully connected layers with 150 hidden units. The final fully connected layer has enough units for the amount of dispersion parameters and abundances depending on the size of the endmember library and number of dispersion parameters per endmember. Then, the dispersion parameters are used to render endmember spectra and the mixture is reconstructed under the linear mixing model with the predicted abundances as inputs. The network only needs the input spectra and the abundances as inputs for training, as the reconstruction error of the spectra is used to back-propagated through the differentiable dispersion model to teach the network to predict good dispersion parameters. Real data (when available) and synthetic data (around 50,000 samples) are used to train the network, which converges after about 100 epochs. An Adam optimizer is used with learning rate set to 1e-3, betas set to (0.9, 0.999), and weight decay set to 0.

Appendix E: Additional Results
Due to size limitations on ArXiv, we cannot display all the mineral maps of Mars TES data here. We encourage the reader to go to Dr. Jayasuriya's webpage https://web.asu.edu/sites/default/files/imaging-lyceum/files/eccv2020_ hyperspectral_paperwithsupplement-min.pdf to see the full supplemental material.