Target-based Hyperspectral Demixing via Generalized Robust PCA

TL;DR

This paper introduces a target localization method in hyperspectral images using a generalized robust PCA model that separates the image into low-rank and sparse components based on known spectral signatures, with theoretical guarantees and experimental validation.

Contribution

It proposes a novel hyperspectral demixing technique leveraging a superposition model with spectral signatures, providing recovery guarantees and demonstrating effectiveness on real data.

Findings

Effective localization of targets based on spectral signatures.

Theoretical recovery guarantees for the proposed method.

Superior performance compared to related techniques in experiments.

Abstract

Localizing targets of interest in a given hyperspectral (HS) image has applications ranging from remote sensing to surveillance. This task of target detection leverages the fact that each material/object possesses its own characteristic spectral response, depending upon its composition. As $\textit{signatures}$ of different materials are often correlated, matched filtering based approaches may not be appropriate in this case. In this work, we present a technique to localize targets of interest based on their spectral signatures. We also present the corresponding recovery guarantees, leveraging our recent theoretical results. To this end, we model a HS image as a superposition of a low-rank component and a dictionary sparse component, wherein the dictionary consists of the $\textit{a priori}$ known characteristic spectral responses of the target we wish to localize. Finally, we analyze the performance of the proposed approach via experimental validation on real HS data for a classification task, and compare it with related techniques.

Tables

Table 1. TABLE I : Simulation results for the proposed approach ( 𝐗𝐩𝐑𝐀 𝐗𝐩𝐑𝐀 \mathbf{XpRA} ), robust-PCA based approach ( RPCA † ), matched filtering ( MF ) on original data 𝐘 𝐘 \mathbf{Y} , and matched filtering on transformed data 𝐑 † ​ 𝐘 superscript 𝐑 † 𝐘 \mathbf{R^{\dagger}Y} ( MF † ), across dictionary elements d 𝑑 d , and the regularization parameter for initial dictionary learning step ρ 𝜌 \rho . Threshold selects columns with column-norm greater than threshold such that AUC is maximized. Further, ` ​ ` ∗ " ` ` " ``*" denotes the case where ROC curve was “flipped” (i.e. classifier output was inverted to improve performance).

$d$	$\rho$	Method	Threshold	Performance at best operating point		AUC
				TPR	FPR
4	0.01	XpRA	0.30	0.979	0.023	0.989
		RPCA†	0.65	0.957	0.049	0.974
		MF∗	N/A	0.957	0.036	0.994
		MF${}_{}{}^{†}{}_{\ast }{}^{}$	N/A	0.914	0.104	0.946
	0.1	XpRA	0.8	0.989	0.017	0.997
		RPCA†	0.8	0.989	0.014	0.997
		MF	N/A	0.989	0.016	0.998
		MF†	N/A	0.989	0.010	0.998
	0.5	XpRA	0.6	0.968	0.031	0.991
		RPCA†	0.6	0.935	0.067	0.988
		MF	N/A	0.548	0.474	0.555
		MF${}_{}{}^{†}{}_{\ast }{}^{}$	N/A	0.849	0.119	0.939

Table 2. (a) Learned dictionary, d = 4 𝑑 4 d=4

$d$	$\rho$	Method	Threshold	Performance at best operating point		AUC
				TPR	FPR
4	0.01	XpRA	0.30	0.979	0.023	0.989
		RPCA†	0.65	0.957	0.049	0.974
		MF∗	N/A	0.957	0.036	0.994
		MF${}_{}{}^{†}{}_{\ast }{}^{}$	N/A	0.914	0.104	0.946
	0.1	XpRA	0.8	0.989	0.017	0.997
		RPCA†	0.8	0.989	0.014	0.997
		MF	N/A	0.989	0.016	0.998
		MF†	N/A	0.989	0.010	0.998
	0.5	XpRA	0.6	0.968	0.031	0.991
		RPCA†	0.6	0.935	0.067	0.988
		MF	N/A	0.548	0.474	0.555
		MF${}_{}{}^{†}{}_{\ast }{}^{}$	N/A	0.849	0.119	0.939

Table 3. (b) Learned dictionary, d = 10 𝑑 10 d=10

$d$	$\rho$	Method	Threshold	Performance at best operating point		AUC
				TPR	FPR
10	0.01	XpRA	0.6	0.935	0.060	0.972
		RPCA†	0.7	0.978	0.023	0.990
		MF∗	N/A	0.624	0.415	0.681
		MF${}_{\ast }{}^{}{}_{}{}^{†}$	N/A	0.569	0.421	0.619
	0.1	XpRA	0.5	0.968	0.029	0.993
		RPCA†	0.5	0.871	0.144	0.961
		MF∗	N/A	0.688	0.302	0.713
		MF†	N/A	0.527	0.469	0.523
	0.5	XpRA	1	0.978	0.031	0.996
		RPCA†	2.2	0.849	0.113	0.908
		MF	N/A	0.807	0.309	0.781
		MF${}_{\ast }{}^{}{}_{}{}^{†}$	N/A	0.527	0.465	0.539

Table 4. (c) Dictionary by sampling voxels, d = 15 𝑑 15 d=15

$d$	Method	Threshold	Performance at best operating point		AUC
			TPR	FPR
15	XpRA	0.3	0.989	0.021	0.998
	RPCA†	3	0.849	0.146	0.900
	MF	N/A	0.957	0.085	0.978
	MF†	N/A	0.796	0.217	0.857

Equations

$$\mathbf{Y}=\mathbf{X}+\mathbf{RA}.$$

$$\tilde{\mathbf{Y}}=\tilde{\mathbf{X}}+\mathbf{A},$$

$$\underset{\mathbf{X},\mathbf{A}}{\mathrm{minimize}}\Vert \mathbf{X}{\Vert }_{\ast }+\lambda \Vert \mathbf{A}{\Vert }_1\text{s.t.}\mathbf{Y}=\mathbf{X}+\mathbf{RA}.$$

$${\mathbf{F}}_L\Vert \mathbf{v}{\Vert }_2^2\le \Vert \mathbf{Rv}{\Vert }^2\le {\mathbf{F}}_U\Vert \mathbf{v}{\Vert }_2^2,$$

$$\mathbf{\Phi }:=\{\mathbf{U}{\mathbf{W}}_1^{\top }+{\mathbf{W}}_2{\mathbf{V}}^{\top },{\mathbf{W}}_1\in {\mathbb{R}}^{nm\times r},{\mathbf{W}}_2\in {\mathbb{R}}^{f\times r}\}.$$

$${\mathbf{\Omega }}_{\mathbf{R}}:=\{\mathbf{Z}=\mathbf{RH},\mathbf{H}\in \mathbf{\Omega }\}.$$

$${\mathcal{P}}_{\Phi }(\mathbf{X})$$

$$\mu :=\underset{\mathbf{Z}\in {\mathbf{\Omega }}_R\backslash \{0_{d\times nm}\}}{\max }\frac{\Vert {\mathcal{P}}_{\Phi }(\mathbf{Z}){\Vert }_F}{\Vert \mathbf{Z}{\Vert }_F},$$

$${\gamma }_{UR}:=\underset{i}{\max }\frac{\Vert {\mathbf{P}}_U\mathbf{R}{\mathbf{e}}_i{\Vert }^2}{\Vert {\mathbf{Re}}_i{\Vert }^2}\text{and}{\gamma }_V:=\underset{i}{\max }\Vert {\mathbf{P}}_V{\mathbf{e}}_i{\Vert }^2,$$

$$C:=\frac{c}{{\mathbf{F}}_L(1-\mu {)}^2-c},$$

$${\gamma }_{UR}\le \{\begin{array}{l}\frac{(1-\mu {)}^2-2s{\gamma }_V}{2s(1+{\gamma }_V)},\text{ for }s\le \text{min}(d,s_{\text{max}})\\ \frac{(1-\mu {)}^2-2s{\gamma }_V}{2(d+s{\gamma }_V)},\text{ for }d<s\le s_{\text{max}}\end{array}.$$

Full text

I Introduction

Hyperspectral (HS) imaging is an imaging modality which senses the intensities of the reflected electromagnetic waves corresponding to different wavelengths of the electromagnetic spectra, often invisible to the human eye. As the spectral response associated with an object/material is dependent on its composition, HS imaging lends itself very useful in identifying target objects/materials via their characteristic spectra or signature responses, often referred to as "endmembers" in the literature. Typical applications of HS imaging range from monitoring agricultural use of land, catchment areas of rivers and water bodies, food processing and surveillance, to detecting various minerals, chemicals, and even presence of life sustaining compounds on distant planets. However, often, these spectral signatures are highly correlated, making it difficult to detect region of interest based on these endmembers. Analysis of these characteristic responses to localize a target material/object serves as the primary motivation of this work.

I-A Model

We begin by formalizing the problem. Each HS image $\mathbf{I}\in\mathbb{R}^{n\times m\times f}$, is a stack of $f$ 2-D images each of size $n\times m$, as shown in Fig. 1(a). Here, $f$ is determined by the number of frequencies or frequency bands across which we measure the reflectances. Overall, each volumetric element or voxel , of a HS image is a vector of length $f$, and represents the response of the material in a pixel of the 2-D grid to $f$ frequencies.

As a particular scene is composed of only a small number of objects/materials, the corresponding characteristic responses of each voxel are low-rank [1]. For example, while imaging an agricultural area, we would expect to record responses from materials like biomass, farm vehicles, roads, houses and water bodies and so on. Further, the spectra of complex materials can be assumed to be a linear mixture of the constituent materials [1, 2], i.e. HS image voxels can be viewed as being generated by a linear mixture model [3]. Since, a particular scene does not contain a large number of distinct materials, we can decompose the responses into a low-rank part and a component which is sparse in a known dictionary – a dictionary sparse part.

Formally, let $\mathbf{Y}\in\mathbb{R}^{f\times nm}$ be formed by unfolding the HS image $\mathbf{I}$, such that, each column of $\mathbf{Y}$ corresponds to a voxel. Now, let $\mathbf{Y}$ arise as a result of a superposition of a low-rank component $\mathbf{X}\in\mathbb{R}^{f\times nm}$ with rank $r$, and a dictionary sparse component expressed here as $\mathbf{RA}$, i.e.,

[TABLE]

Here, $\mathbf{R}\in\mathbb{R}^{f\times d}$ is an a priori known dictionary composed of normalized characteristic responses of the material/object, or of the constituents of the material we wish to localize, and $\mathbf{A}\in\mathbb{R}^{d\times nm}$ is the sparse coefficient matrix (also referred to as "abundances" in the literature) with $s$ non-zero elements globally, i.e. $\|\mathbf{A}\|_{0}=s$. Note that, $\mathbf{R}$ can also be constructed by learning a dictionary based on the known spectral signatures of a target[4, 5, 6].

I-B Prior Art

The model shown in (1) is a close counterpart of a number of well-known problems. To start, in the absence of the dictionary sparse part, (1) reduces to the popular problem of principal component analysis (PCA) [7]. This problem also shares its structure with variants of PCA such as robust-PCA[8, 9] where $\mathbf{R}$ is absent, outlier pursuit [10] where $\mathbf{R}=\mathbf{I}$ and $\mathbf{A}$ is column sparse, and others [11, 12, 13, 14, 15, 16, 17, 18, 19]. On the other hand, the problem can be identified as that of sparse recovery [20, 21, 22], in the absence of the low-rank part $\mathbf{X}$. Following which, sparse recovery methods for analysis of HS images have been explored in [23, 24, 25, 26]. In addition, in a recent work [27], the authors further impose a low-rank constraint on the coefficient matrix $\mathbf{A}$. Further, applications of compressive sampling have been explored in [28], while [3] analyzes the case where HS images are noisy and incomplete.

The model described in (1) was introduced in [29] as a means to detect traffic anomalies in a network, wherein, the authors focus on a case where the dictionary $\mathbf{R}$ is overcomplete, i.e., fat, and the rows of $\mathbf{R}$ are orthogonal, e.g., $\mathbf{RR^{\top}}=\mathbf{I}$. Further, the coefficient matrix $\mathbf{A}$ is assumed to possess at most $k$ nonzero elements per row and column and $s$ nonzero elements overall.

In a recent work [30], we analyze the extension of [29] to a case where the dictionary has more rows than columns, i.e., is thin, while removing the orthogonality constraint for both the thin and the fat dictionary cases, when $s$ is small. The thin dictionary case is amenable to a number of applications, as often, we are looking for only a small number of signatures. For example, in our motivating application of HS imaging, we are interested in detecting the presence of a small number of spectral signatures corresponding to the target of interest. To this end, we focus our attention on the thin case, although a similar analysis applies for the fat case [30].

Note that, when the known dictionary $\mathbf{R}$ is thin, we can alternatively multiply (1) on the left by the pseudo-inverse of $\mathbf{R}$, say $\mathbf{R}^{\dagger}$, in which case, the model shown in (1) reduces to that of robust PCA, i.e.,

[TABLE]

where $\mathbf{\tilde{Y}}=\mathbf{R}^{\dagger}\mathbf{Y}$ and $\mathbf{\tilde{X}}=\mathbf{R}^{\dagger}\mathbf{X}$. Therefore, in this case, we can adopt a two-step procedure: (1) recover the sparse matrix $\mathbf{A}$ by robust PCA and, (2) estimate the low-rank part using the estimate of $\mathbf{RA}$. 111We denote this method of estimating $\mathbf{A}$ as RPCA† in our experiments. Models considered in (1) and (2) are two alternate methods of recovering the components for thin dictionaries. These two methods inherently solve different optimization problems. The proposed method is a one-shot procedure with guarantees to solve the problem, while the conditions under which the pseudo-inverse based procedure succeeds will depend on the interaction between $\mathbf{R}^{\dagger}$ and the low-rank part. Further, our experiments indicate that solving for model in (1) is more robust as compared to RPCA† for the classification problem at hand, across different choice of dictionaries.

I-C Our Contribution

In this work, we present a one-step method for target detection in a hyperspectral image with the underlying model of (1). This is done by building upon the results of [30] for the underlying application. In particular, we show the conditions under which target detection is possible with thin dictionaries. In general, the method is also applicable to fat dictionaries and for other applications, which follow the underlying model. The rest of the paper is organized as follows. We formulate the problem, introduce the notation and present the theoretical result for the case at hand, i.e., the thin case, in section II. Next, we present the description of the dataset along with the experimental evaluations in section III. Finally, we conclude this discussion in section IV.

II Problem Formulation

Our aim is to recover the low-rank component $\mathbf{X}$, and the sparse coefficient matrix $\mathbf{A}$, given the dictionary $\mathbf{R}$, and samples $\mathbf{Y}$ generated according to the model described in (1). We utilize the structure of the components, $\mathbf{X}$ and $\mathbf{A}$, to arrive at the following convex problem for $\lambda\geq 0$,

[TABLE]

where, $\|.\|_{*}$ denotes the nuclear norm, and $\|.\|_{1}$ refers to the $l_{1}$- norm, which serve as convex relaxations of rank and sparsity (i.e. $l_{0}$-norm), respectively. For the application at hand, we present the theoretical result for the case when the number of dictionary elements $d\leq f$. Analysis of a more general case can be found in [30]. For the this case, we assume the dictionary $\mathbf{R}$ to be a frame, i.e. for any vector $\mathbf{v}\in\mathbb{R}^{d}$, we have

[TABLE]

where $\mathbf{F}_{L}$ and $\mathbf{F}_{U}$ are the lower and upper frame bounds, respectively, with $0<\mathbf{F}_{L}\leq\mathbf{F}_{U}$.

We begin by defining a few relevant subspaces, similar to those used in [29].

Let the pair $\{\mathbf{X_{0},A_{0}}\}$ be the solution to the problem shown in (3). We define $\mathbf{\Phi}$ as the space of matrices spanning either the row or the column space of the low-rank component $\mathbf{X}_{0}$. Specifically, let $\mathbf{U\Sigma V^{\top}}$ denote the singular value decomposition of $\mathbf{X}_{0}$, then the space $\mathbf{\Phi}$ is defined as

[TABLE]

Next, let $\mathbf{\Omega}$ be the space spanned by the matrices which have the same support (locations of non-zero elements) as $\mathbf{A}_{0}$, and let $\mathbf{\Omega_{R}}$ be defined as

[TABLE]

Next, let $\mathcal{P}_{\Phi}(.)$, $\mathcal{P}_{\Omega}(.)$ and $\mathcal{P}_{\Omega_{R}}(.)$ be the orthogonal projection operator(s) onto the space of matrices defined above. In addition, we will use $\mathbf{P}_{U}$ and $\mathbf{P}_{V}$ to denote the projection matrices for the column and row spaces of $\mathbf{X}_{0}$, respectively, i.e., implying the following for any matrix $\mathbf{X}$,

[TABLE]

Indeed, there are situations under which we cannot hope to recover the matrices $\mathbf{X}$ and $\mathbf{A}$. To characterize these scenarios, we employ various notions of incoherence, out of which, the first is the incoherence between the low-rank part, $\mathbf{X}_{0}$, and the dictionary sparse part, $\mathbf{RA_{0}}$,

[TABLE]

where $\mu\in[0,1]$ is small when these components are incoherent (which is good for recovery). The next two measures of incoherence can be interpreted as a way to avoid the cases where for $\mathbf{X}_{0}=\mathbf{U\Sigma V^{\top}}$, (a) $\mathbf{U}$ resembles the dictionary $\mathbf{R}$, and (b) $\mathbf{V}$ as the sparse coefficient matrix $\mathbf{A}_{0}$. To this end, we define respectively, the following to measure these properties,

[TABLE]

where $\mathbf{\gamma}_{V}{\small\in[r/nm,1]}$. Also, we define $\xi:=\|\mathbf{R^{\top}}\mathbf{UV^{\top}}\|_{\infty}$. Finally, for a matrix $\mathbf{M}$, we use $\|\mathbf{M}\|:=\sigma_{\text{max}}(\mathbf{M})$ for the spectral norm, where $\sigma_{\text{max}}(\mathbf{M})$ denotes the maximum singular value of the matrix $\mathbf{M}$, and $\|\mathbf{M}\|_{\infty}:=\underset{\{i,~{}j\}}{\text{max}}|M_{ij}|$.

II-A Theoretical Result

We now present our main theoretical result for the thin case, i.e., when $d\leq f$, which was developed in [30]. We begin by introducing the following definitions and assumptions.

Definition D.1**.**

$\lambda_{\text{min}}=\tfrac{1+C}{1-C}~{}\xi$, where $C$ is defined as

[TABLE]

$c:=\tfrac{\mathbf{F}_{U}}{2}[(1+2\gamma_{U\!R})(\text{min}(s,d)+s\gamma_{V})+2s\gamma_{V}]-\tfrac{\mathbf{F}_{L}}{2}[\text{min}(s,d)+s\gamma_{V}]$.

Definition D.2**.**

\lambda_{\text{max}}:=\tfrac{1}{\sqrt{s}}\big{(}\sqrt{\mathbf{F}_{L}}{(}~{}1-\mu~{}{)}-\sqrt{r\mathbf{F}_{U}}\mu\big{)}.

Assumption A.1**.**

$\lambda_{\text{max}}\geq\lambda_{\text{min}}$

Assumption A.2**.**

Let $s_{\text{max}}:=\tfrac{(1-\mu)^{2}}{2}\tfrac{nm}{r}$, then

[TABLE]

The following theorem establishes the conditions under which we can recover the pair $\{\mathbf{X_{0},A_{0}}\}$ successfully.

Theorem 1**.**

Consider a superposition of a low-rank matrix $\mathbf{X}_{0}\in\mathbb{R}^{f\times nm}$ of rank $r$, and a dictionary sparse component $\mathbf{RA_{0}}$, wherein the dictionary $\mathbf{R}\in\mathbb{R}^{f\times d}$, with $d\leq f$, obeys the frame condition with frame bounds $[\mathbf{F}_{L},\mathbf{F}_{U}]$, the sparse coefficient matrix $\mathbf{A}_{0}$ has at most $s$ non-zeros, i.e., $\|\mathbf{A}_{0}\|_{0}=s$, where $s\leq s_{\text{max}}$, and $\mathbf{Y}=\mathbf{X}_{0}+\mathbf{RA_{0}}$, with parameters $\gamma_{U\!R}$, $\xi$, $\mathbf{\gamma}_{V}\in[r/nm,1]$ and $\mu\in[0,1]$. Then, solving the formulation shown in (3) will recover matrices $\mathbf{X}_{0}$ and $\mathbf{A}_{0}$ if assumptions A.1 and A.2 hold for any $\lambda\in[\lambda_{\text{min}},\lambda_{\text{max}}]$, as defined in D.1 and D.2, respectively.

In other words, Theorem 1 establishes the sufficient conditions for the existence of $\lambda$s to guarantee recovery of $\{\mathbf{X_{0},A_{0}}\}$. We see that these conditions are closely related to the various incoherence measures $\gamma_{U\!R}$, $\mathbf{\gamma}_{V}$ and $\mu$, between $\mathbf{X}$, $\mathbf{R}$ and $\mathbf{A}$ . Note that, we have here an upper-bound on the global sparsity, i.e., $s\leq s_{\text{max}}$, however, our simulation results in [30] on synthetic data reveal that the recovery is indeed possible even for higher $s$; See [30] for the details of the analysis. To this end, we propose a randomized analysis of the problem in (3) to improve these results as a future work.

III Experimental Evaluation

In this section, we evaluate the performance of the proposed technique on real HS data. We begin by describing the dataset used for the simulations, following which we describe the experimental set-up and conclude by presenting the results.

III-A Data

The Airborne Visible/Infrared Imaging Spectrometer (AVIRIS) [31] sensor is one of the popular choices for collecting HS images for various remote sensing applications. In this present exposition, we consider the “Indian Pines” dataset [32], which was collected over the Indian Pines test site in North-western Indiana in the June of 1992. This dataset consists of spectral reflectances across $224$ bands in wavelength of ranges $400-2500$ nm from a scene which is composed mostly of agricultural land, along with two major dual lane highways, a rail line and some built structures. This dataset is further processed by removing the bands corresponding to those of water absorption, and the resulting HS data-cube with dimensions $\{145\times 145\times 200\}$ is as visualized in Fig. 1. This modified dataset is available as “corrected Indian Pines” dataset [32], with the ground-truth containing $16$ classes; Henceforth, referred to as the “Dataset".

We will analyze the performance of the proposed technique for the identification of the stone-steel towers (class $16$ in the Dataset), shown as ground-truth in Fig. 2(b), which constitutes about $93$ voxels in the Dataset. We consider two cases for choosing the known dictionary $\mathbf{R}$ – (1) when a (thin) dictionary is learned[4, 5, 6] based on the voxels (class-$16$), and (2) when the dictionary is formed by randomly sampling voxels (from class-$16$). This is to emulate the two ways in which we can arrive at the characteristics of a target – (1) the case where the exact signatures are not available, and/or there is noise, and (2) where we have access to the exact signatures we are looking for. Next, we form the data matrix $\mathbf{Y}\in\mathbb{R}^{f\times nm}$ by stacking each voxel of the image side-by-side, which results in a $\{200\times 145^{2}\}$ data matrix $\mathbf{Y}$.

Normalization: For normalization, we divide the data matrix $\mathbf{Y}$ and the dictionary $\mathbf{R}$ by $\|\mathbf{Y}\|_{\infty}$ to preserve the inter-voxel scaling. Next, we normalize the columns of $\mathbf{R}$, such that they are unit-norm. In practice, we can store the un-normalized dictionary $\mathbf{R}$, consisting of actual signatures of the target material, and can normalize it, as described above, after the HS image has been acquired.

III-B Experimental Setup

We evaluate and compare the performance of the proposed method (XpRA) with RPCA† using transformed data $\mathbf{\tilde{Y}=R^{\dagger}Y}$ (as described in section I-B), matched filtering (MF) based approach on data $\mathbf{Y}$ ( denoted by MF), and a MF based approach on $\mathbf{\tilde{Y}}$, denoted by MF†, via receiver operating characteristic (ROC), shown in Table I(a)-(d).

For the proposed technique, we employ the accelerated proximal gradient (APG) algorithm presented in [29] to solve the optimization problem shown in (3). Similarly, for RPCA† we employ the APG with transformed data matrix $\mathbf{\tilde{Y}}$, while setting $\mathbf{R=I}$; See (2). With reference to selection of tuning parameters for the APG solver for XpRA (RPCA†, respectively), we choose $v=0.95$, $\nu=\|\mathbf{Y}\|$ ($\nu=\|\mathbf{\tilde{Y}}\|$), $\bar{\nu}=10^{-4}$, and scan through $100$ $\lambda$s in the range $\lambda\in(0,{|\mathbf{R^{\top}Y}\|_{\infty}}/{\|\mathbf{Y}\|}]$ ($\lambda\in(0,{\|\mathbf{\tilde{Y}}{\|_{\infty}}/{\|\mathbf{\tilde{Y}}\|}}]$), to generate the ROCs. The output of XpRA and RPCA† is a data-cube corresponding to the low-rank part and the dictionary sparse part (which contains the region of interest). We threshold the resulting estimates of the sparse part $\mathbf{A}$, based on their column norms. We choose the threshold such that the area under the curve (AUC) is maximized for both cases (XpRA and RPCA†).

In case of MF, we form the inner-product of the column-normalized data matrix $\mathbf{Y}$, say $\mathbf{Y}_{n}$, with the dictionary $\mathbf{R}$, i.e., $\mathbf{R^{\top}Y}_{n}$, and select the maximum absolute inner-product per column. For MF†, matched filtering corresponds to finding maximum absolute entry for each column of the column-normalized $\mathbf{\tilde{Y}}$. Next, in both cases we scan through $1000$ threshold values between (0,{\color[rgb]{0,0,0}1}] to generate the ROC results shown in Table I.

Table I(a)-(b) show the ROC characteristics for the aforementioned methods for a dictionary learned using all voxels of class-$16$. We evaluate the performance of the aforementioned techniques for two dictionary sizes, $d=4$ and $d=10$, for three values of the regularization parameter $\rho$, used for the initial dictionary learning step, i.e., $\rho=0.01,~{}0.1$ and $0.5$. (The parameter $\rho$ controls the sparsity during the initial dictionary learning step.)

Table I(c) shows the ROC results for the case where we randomly select $15$ voxel from class-$16$ to form the known dictionary. Further, we show the overall performance in terms of the mean and standard deviation of ROC metrics (TPR, FPR and AUC), for each algorithm in Table I(d), i.e., across different dictionaries considered in Table I(a)-(c).

III-C Analysis

As described above, Table I shows the ROCs for the classification performance of the proposed method as compared to RPCA†, MF and MF†. We note that when we use $\mathbf{R}$ which has been learned based on samples of class-$16$, matched filtering based approaches under perform (with an exception for $d=4$ and $\rho=0.1$). RPCA† performs well for $d=10$ and $\rho=0.01$. For all other choices the proposed method performs the best. Even for cases in which it is not the best, it still is comparable to the best performer. Following this, we conclude that the performance of XpRA is most reliable across different the choices of the dictionaries. This is also evident by Table I(d) which shows the mean and and standard deviation of ROC metrics (TPR, FPR and AUC). We observe that the mean and standard deviation metrics for the proposed method are superior to other methods across the different dictionary choices.

There are other interesting recovery results which warrant our attention. For example, Fig. 2 shows the low-rank component and the dictionary sparse component recovered by the proposed method for two different values of $\lambda$, for the case where we form the dictionary by randomly sampling the voxels (Table I(c)). Interestingly, we recover the rail tracks/roads running diagonally on the top-right corner, along with some low-density housing; See Fig 2 (f). This is because the signatures we seek (stone-steel towers) are similar to those of the materials used in these structures. This further corroborates the applicability of the proposed approach in detecting the presence of a particular spectral signature in a HS image. However, this also highlights potential drawback of this technique. As XpRA is based on identifying materials with similar composition, it may not be effective in distinguishing between very closely related classes, say two agricultural crops.

IV Conclusions

We present a generalized robust PCA-based technique to localize a target in a HS image, based on the a priori known spectral signature of the material we wish to localize. We model the data as being composed of a low-rank component and a dictionary sparse component. Here, the dictionary contains the a priori known spectral signatures of the target. In this work, we adapt the theoretical results of [30] and present the conditions under which such a decomposition recovers the two components for the HS demixing task. Further, we evaluate and compare the performance of the proposed method via experimental evaluations on real HS image data for a classification task.

V Acknowledgement

The authors graciously acknowledge support from the DARPA YFA, Grant N66001-14-1-4047.