A data-driven approach to sampling matrix selection for compressive sensing

TL;DR

This paper introduces a data-driven method for selecting sampling orders in compressive sensing, optimizing variance capture with applications to various sensors and improved chemical detection at high compression rates.

Contribution

It presents a novel approach to basis ordering in compressive sensing that enhances data capture efficiency and sensor performance, applicable across different sampling bases.

Findings

Successful chemical detection at 90% compression rate

Improved peak SNR by over 30% in depth images

Outperforms standard sampling orderings in experiments

Abstract

Sampling is a fundamental aspect of any implementation of compressive sensing. Typically, the choice of sampling method is guided by the reconstruction basis. However, this approach can be problematic with respect to certain hardware constraints and is not responsive to domain-specific context. We propose a method for defining an order for a sampling basis that is optimal with respect to capturing variance in data, thus allowing for meaningful sensing at any desired level of compression. We focus on the Walsh-Hadamard sampling basis for its relevance to hardware constraints, but our approach applies to any sampling basis of interest. We illustrate the effectiveness of our method on the Physical Sciences Inc. Fabry-P\'{e}rot interferometer sensor multispectral dataset, the Johns Hopkins Applied Physics Lab FTIR-based longwave infrared sensor hyperspectral dataset, and a Colorado State University Swiss Ranger depth image dataset. The spectral datasets consist of simulant experiments, including releases of chemicals such as GAA and SF6. We combine our sampling and reconstruction with the adaptive coherence estimator (ACE) and bulk coherence for chemical detection and we incorporate an algorithmic threshold for ACE values to determine the presence or absence of a chemical. We compare results across sampling methods in this context. We have successful chemical detection at a compression rate of 90%. For all three datasets, we compare our sampling approach to standard orderings of sampling basis such as random, sequency, and an analog of sequency that we term `frequency.' In one instance, the peak signal to noise ratio was improved by over 30% across a test set of depth images.

Equations

$$y=S\tilde{x}\in {\mathbb{R}}^k$$

$${\text{argmin}}_{u\in {\mathbb{R}}^n}||u|{|}_{{\ell }_1}\text{such that }y=S{H}^{-1}u.$$

$$(1-{\delta }_p)\le \frac{||Az|{|}_{{\ell }_2}^2}{||z|{|}_{{\ell }_s}^2}\le (1+{\delta }_p).$$

$$\text{RMSE}:=\sqrt{\frac{1}{mn}\sum _{i=1}^m\sum _{j=1}^n{\left(X_{ij}-Y_{ij}\right)}^2}.$$

$$\text{PSNR}:=20{\log }_{10}\left(\frac{N}{\text{RMSE}}\right).$$

$$\frac{(s^T{\Gamma }^{-1}x{)}^2}{(s^T{\Gamma }^{-1}s)(x^T{\Gamma }^{-1}x)}.$$

$$1-\prod _{i=1}^M(1-c_i).$$

$$\left[\begin{array}{cccccccc}[r]1& 1& 1& 1& 1& 1& 1& 1\\ 1& -1& -1& 1& 1& -1& -1& 1\\ 1& -1& 1& -1& -1& 1& -1& 1\\ 1& 1& -1& -1& -1& -1& 1& 1\\ 1& 1& -1& -1& 1& 1& -1& -1\\ 1& -1& 1& -1& 1& -1& 1& -1\\ 1& -1& -1& 1& -1& 1& 1& -1\\ 1& 1& 1& 1& -1& -1& -1& -1\end{array}\right].$$

$$\left[\begin{array}{cccccccc}[r]1& 1& 1& 1& 1& 1& 1& 1\\ 1& 1& 1& 1& -1& -1& -1& -1\\ 1& 1& -1& -1& -1& -1& 1& 1\\ 1& 1& -1& -1& 1& 1& -1& -1\\ 1& -1& -1& 1& 1& -1& -1& 1\\ 1& -1& -1& 1& -1& 1& 1& -1\\ 1& -1& 1& -1& -1& 1& -1& 1\\ 1& -1& 1& -1& 1& -1& 1& -1\end{array}\right].$$

$$M=\{s_{{\rho }^{-1}(i)}{\}}_{i=1}^{\lceil (1-\alpha )m\rceil }.$$

$$\left[1,1,-1,-1,1,1,-1,-1,1,1,-1,-1,1,1,-1,-1\right]$$

Full text

1 INTRODUCTION

Many promising kinds of imaging currently have limited application because the sensors required are either too expensive, too slow, or too large for a given application. One of the attractive aspects of the compressive sensing (CS) framework is that it allows for alternative sampling routines that can sometimes overcome these limitations. The single pixel camera proposed by Duarte et. al. [1] is a good example of this; the use of CS in MRI is another [2].

Much of the theory behind CS attempts to understand which sampling techniques maximize accuracy in a reconstructed signal. Foundational work in the field provides reconstruction guarantees based on whether the sampling matrix $S$ (which controls how the signal is sampled) satisfies the restricted isometry property [3]. When sampling is actually done via a product of matrices $SD$ where $D$ is fixed (as often happens in practice when the signal is not sparse in its natural basis), then $S$ can be chosen via its incoherence [4] with $D$. Incoherence itself can (with high probability) be obtained by drawing entries of $S$ randomly from some known distribution. This scheme has been used for many specific applications of CS [5, 2, 6]. It has also been shown to be flexible enough to give reconstruction guarantees when the specific structure is expected between samples [7].

In this paper we take a different approach. Rather than designing a sampling scheme which is guaranteed to be able to reconstruct any signal, we instead choose $S$ based on its ability to capture the maximum amount of variance on a set of training signals (or in other words, a set of signals representative of the signals we expect to later sample and reconstruct). We show that when $S$ is constructed in such a way it performs similarly or better than other popular sampling schemes. Based on these results we propose this method for CS applications where a set of representative example signals is known before imaging begins.

To simplify the analysis (and also to make the results applicable to the single-pixel-camera framework) we focus on sampling matrices $S$ whose entries consist of [math]’s and $1$’s. Specifically, we study sampling matrices generated from subsets of rows of the shifted Walsh-Hadamard matrix. In order to show the flexibility of our method we present results for reconstruction of both LIDAR and hyperspectral data (the Physical Sciences Inc. Fabry-Pérot interferometer sensor multispectral dataset, the Johns Hopkins Applied Physics Lab FTIR-based longwave infrared sensor hyperspectral dataset, and the Colorado State University Swiss Ranger depth image dataset) [8, 9]. One particular advantage is that since rows of the matrix $S$ are ranked by the amount of variance they capture in the data, $S$ does not need to be recalculated when decreasing the sampling level. Instead we can simply omit the appropriate rows of $S$.

This paper is organized as follows. In Secs. 2.1-2.2 we review the basics of CS and state the optimization problem that we are solving in order to reconstruct signals. In 2.3 we review the ideas that underlie the standard methods for generating a sampling matrix. We state and describe the metrics that we used to analyze signal reconstruction accuracy in Sec. 2.4. Since two of the datasets used to measure the performance of our method consist of hyperspectral data, in Sec. 2.5 we describe the algorithms that we used to detect signals in this context. In Sec. 3, we describe our maximal-variance based method of generating a sampling matrix. In this section we also give performance results for our method on both hyperspectral and LIDAR datasets. Finally, we conclude and suggest directions for future work in Sec. 4.

2 BACKGROUND

2.1 Compressive Sensing

Compressive sensing (CS) is a collection of methods that, in certain cases, allows for highly accurate reconstruction of signals even when they have been sampled well under the Nyquist-rate [6, 10]. Framed in terms of linear algebra (as much of this paper is), CS can be equivalently formulated as being a set of methods for solving the ill-posed problem of finding $\tilde{x}$ from

[TABLE]

when $S$ is a $k\times n$ matrix, $\tilde{x}\in\mathbb{R}^{n}$, and $k<n$. CS methods used to solve such problems utilize additional assumptions about the signal $\tilde{x}$. These additional assumptions allow one to choose the correct signal $\tilde{x}$ out of the infinite number of solutions $x^{\prime}$ to $Sx^{\prime}=y$. For example, it is often the case that $\tilde{x}$ will be sparse (or compressible) in some particular basis. Formally this means that there is some $n\times n$ basis transformation matrix $D$ such that $\tilde{x}=D\tilde{u}$ and $\tilde{u}$ is sparse.

The matrix $S$ is known as the sampling matrix and $D$ is known as the sparsity-promoting matrix. As its name suggests, the sampling matrix defines the way in which samples are collected from a signal. In many CS applications, one has some latitude to choose $S$ up to context-specific constraints, and hence understanding which matrices $S$ will capture the maximum amount of information about a signal from a minimum number of samples is a question of central importance to the CS community. As a consequence, rich theory has been developed to explain the relative performance of different sampling matrices (see Sec. 2.3 for a brief overview). In this paper we will present a new data-driven approach to the choice of $S$.

2.2 Reconstruction of signal via optimization

Given a sampled signal $y=SD\tilde{u}$, CS theory shows that, provided that $S$ has been well-chosen and $\tilde{u}$ has certain basic characteristics (such as sparsity), a close approximation $x$ to $\tilde{x}$ can be reconstructed by solving an optimization problem. In this paper, the optimization problem that was solved in order to approximate $\tilde{x}$ from the set of measurements $y=S\tilde{x}\in\mathbb{R}^{k}$ with $k\ll n$ is

[TABLE]

Here $H,H^{-1}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ is the Haar wavelet transform and its inverse, respectively. Informally, this optimization problem requires that we find a solution $u$ which satisfies $SH^{-1}u=y$ such that $||u||_{\ell_{1}}$ is minimized. Minimizing the term $||u||_{\ell_{1}}$ typically leads to a solution which is sparse in the wavelet domain. At an intuitive level, this makes sense as many types of signals (particularly various kinds of images) have been shown to be naturally sparse in wavelet bases. The attentive reader will notice that the $\ell_{1}$-norm is used rather than the $\ell_{0}$ pseudo-norm (which measures sparsity directly). This substitution is common practice in CS since the $\ell_{1}$-norm is convex whereas the $\ell_{0}$-pseudo norm is not and results in problems which are impractical to solve.

Driven in part by the interest in CS, there are now many methods of solving $\ell_{1}$-type optimization problems. We choose to use the split Bregman method [11] to solve (1) because it provides fast, lightweight reconstruction with reliable convergence properties.

2.3 Restricted Isometry Property and Incoherence

In this section we review one of the most popular approaches to producing sampling matrices that guarantee good reconstructions of signals. These ideas hinge on a characterization of matrices called the restricted isometry property (RIP) [12]. Our introduction of this topic follows the outline given by Baraniuk [6].

In the unrealistic situation where we know that $\tilde{u}$ is $p$-sparse and we know the exact coefficients in which $\tilde{u}$ is sparse, then we only need to consider the corresponding $p$ columns of $A=SD$ (or for optimization problem (1), the corresponding $p$ columns of $SH^{-1}$). In particular, for $A$ to have a well-conditioned inverse on all $p$-sparse vectors $z$ following the sparsity pattern of $\tilde{u}$ we need to ensure that for all such $z$ the following inequality is satisfied:

[TABLE]

In general of course, while we may have a guess for the natural sparsity of our signal $\tilde{u}$, we rarely know which specific $p$ coefficients of $\tilde{u}$ will be nonzero (and in most real-world situations, $u$ will not be sparse in the sense of having many terms be exactly zero, but instead will be compressible in the sense of many terms being very close to zero). To generalize to such a situation, $A$ is said to have the $p$-restricted isometry property provided that for any $p$-sparse vector $z$, (2) holds. It has been shown that recovery of signals that are either sparse or compressible is guaranteed when $A$ satisfies a specific RIP depending on the assumed sparsity of the signal.

However, even verifying that a matrix satisfies an RIP property is a combinatorially difficult problem. Furthermore, in our approach we can only change $S$, not $H^{-1}$ which is already assumed to be fixed. In order to try to obtain a product, $SH^{-1}$, which guarantees good reconstruction properties, another (related) option is to try to maximize the incoherence between $S$ and $H^{-1}$. Roughly, this means constructing $S$ such that the columns of $S$ cannot sparsely represent the columns of $H^{-1}$.

The standard way to either find a matrix $S$ satisfying an RIP property or $S$ such that $S$ and $H^{-1}$ are incoherent is to utilize randomness. It has been shown that if $S$ is drawn randomly from a number of standard distributions, then with high probability $S$ and $H^{-1}$ will be incoherent and $SH^{-1}$ will satisfy the RIP property [13, 14, 15]. Of course, the methods above produce $S$ with good properties in general. In some cases however, characteristics of the signal we wish to sample are known sufficiently well that it makes sense to try to adapt $S$ to best capture information from this specific class of signals. Additionally, hardware constraints sometimes make the methods above inaccessible. Presenting a procedure that addresses this setting is the purpose of this paper.

2.4 Image Quality Metrics

In order to compare reconstructed depth images against original uncompressed data in Sec. 3.3.1, we rely on two common image quality metrics: root-mean-square error (RMSE) and peak signal-to-noise ratio (PSNR). We will introduce them here in the context of image comparison. Let $X$ be an $m\times n$ image and let $Y$ be an approximation of $X.$ RMSE is defined as its name suggests: specifically, it is the square root of the average of the set of squared differences in pixel values. That is,

[TABLE]

Note that a consequence of this definition is that large differences in pixel values are exaggerated in their contributions to RMSE and small values ($<1$) have diminished contributions. The definition of PSNR can be written in terms of RMSE, though interpretations of the two metrics differ. If $N$ is the maximum possible pixel value that can occur in $X,$ then we define PSNR for the approximation $Y$ to be

[TABLE]

There are many image quality metrics [16]; we choose to focus on two of the most widely-used for the purpose of this work.

2.5 Chemical Detection: ACE, MPACE, Thresholds

The results we present in Sec. 3.3.2 for reconstructed hyperspectral data involve videos of chemical releases. Thus our measure of quality is with respect to quantitative comparisons of chemical detection on reconstructed vs. raw data. We therefore present background on several relevant chemical detection techniques here.

A standard approach to chemical detection is via the adaptive coherence/cosine estimator (ACE) algorithm [17, 18]. Given a signature of a target chemical $s$ and a signature in a pixel $x$ in a hyperspectral cube, the ACE algorithm calculates the square of the cosine of the angle between $s$ and $x$ relative to background data. Specifically, if $\Gamma$ is the maximum likelihood estimator for the covariance matrix of the background, then the ACE statistic is

[TABLE]

The signal can be strengthened in this setting by additionally applying a technique called bulk coherence or multipulse coherence (MPACE) [19, 20, 21]. This statistic is intended to enhance the signal in neighborhoods in which several pixels have high ACE values. Let $c_{i}$ be the value of the ACE statistic in pixel $i$ and consider a neighborhood of $M$ pixels. The bulk coherence statistic is then defined as

[TABLE]

If several of the ACE values in the neighborhood are high (close to 1), then the product of the $(1-c_{i})$ terms will be close to zero, resulting in a bulk coherence value close to 1. We notice significant improvement of signal strength as a result of the incorporation of bulk coherence. We further add noise reduction by including what we refer to as persistence. Specifically, we add a filter that associates a value of zero to a pixel if the bulk coherence value in that pixel doesn’t remain above a chosen threshold for five consecutive time frames.

To make comparison of results objective, we algorithmically defined a threshold for determining whether a chemical is present or absent in a pixel based on the ACE statistic [22]. We construct this threshold to be responsive to the device in use, the chemical of interest, and the reconstruction approach (e.g. $\ell_{1}$ optimization). Experimentally, this threshold appears to faithfully detect the presence or absence of chemicals in hyperspectral cubes and makes objective comparison possible.

3 DATA-DRIVEN SAMPLING

As described in Sec. 2.3, the classical approach to sampling within the CS framework is to use a random sampling set. In this section we will present an alternative data-driven method. One of several potential benefits of our method is that it is flexible enough to accommodate hardware constraints that might make the use of traditional random sampling difficult. For this paper, we consider a setting in which hardware constraints specifically require either the use of entries of $\pm 1$ only or the use of entries of [math] and $1$ only (the latter is relevant to the single pixel camera architecture). A natural choice in this context is to use the sampling basis provided by a Walsh-Hadamard transform or its ‘shifted’ version. From the matrix for the transform, we select a subset of rows in order to perform CS at the desired level. Thus when we refer to sampling, we often refer to an ordering - the ordering of the rows determines the choice of subset since we simply select the first $k$ rows for sampling. The shifted version of a Walsh-Hadamard transform simply replaces $-1$ entries in the corresponding matrix with [math] entries. Walsh-Hadamard transforms were defined as early as 1867 by Sylvester; they can be defined in various ways and the corresponding matrices have been given different names and orderings [23, 24, 25, 26]. We consider the standard ordering (‘Hadamard order’) of the rows of a Walsh-Hadamard matrix, the sequency ordering, an analog of sequency that we term frequency, random ordering, and the maximal-variance order as defined in terms of a relevant dataset. We also consider the corresponding orderings of the shifted Walsh-Hadamard matrices.

We briefly review some properties of Walsh-Hadamard matrices along with some important orderings of the rows of the Walsh-Hadamard matrix transform. Walsh-Hadamard matrices are square $(\pm 1)$-matrices and satisfy $W_{m}W_{m}^{T}=mI_{m},$ where $W_{m}$ is an $m\times m$ Walsh-Hadamard matrix and $I_{m}$ is the $m\times m$ identity. We provide examples to illustrate various row orderings here. Let $m=8.$ Then the standard order produces

[TABLE]

The sequency order, in contrast, is defined so that the rows are in increasing order of zero-crossings between consecutive entries. Thus, the sequency order gives

[TABLE]

Finally, inspired by the success of the use of randomness in standard CS settings, we will also compare against random orderings of the rows of Walsh-Hadamard matrices.

3.1 Maximal-Variance Sampling

In this section, we present the main definition of this paper: a data-driven method for selecting a subset of a basis for optimal sampling.

Definition: Maximal-Variance Sampling Choose a parameter $\alpha\in[0,1]$ that is the desired percent of compression. Let $X\in\mathbb{R}^{m\times n}$ be a set of training data of $n$ samples of points in $\mathbb{R}^{m}.$ Given a set of sampling vectors $S=\{s_{i}\}_{i=1}^{p},$ with $p\geq\lceil(1-\alpha)m\rceil$, where $\lceil\cdot\rceil$ is the ceiling function, we compute the variance captured by each vector $s_{i}.$ Let $\sigma_{i}^{2}$ be the variance of the set $\{s_{i}\cdot x_{j}\}_{j=1}^{n};$ equivalently $\sigma_{i}^{2}$ is the variance of the vector $s_{i}^{T}X.$ Define a permutation $\rho$ of $(1,2,\dots,p)$ so that $(\sigma_{\rho^{-1}(1)}^{2},\sigma_{\rho^{-1}(2)}^{2},\dots,\sigma_{\rho^{-1}(p)}^{2})$ is decreasing. Then we define the maximal-variance sampling subset $M$ of $S$ to be

[TABLE]

The maximal-variance sampling algorithm proposed in the definition above is independent of the size of the sampling set. In the case where $S$ is defined by the Walsh-Hadamard matrix, the set of sampling vectors forms a basis for $\mathbb{R}^{m},$ so $p=m.$ In general, the set $S$ need not have this property. The only requirement is that $p\geq\lceil(1-\alpha)m\rceil.$

3.2 Frequency Ordering

While we propose maximal-variance sampling and demonstrate in Secs. 3.3.1 and 3.3.2 that it outperforms other orderings of rows from the Walsh-Hadamard matrix in terms of accuracy of reconstruction, it is sometimes preferable to have an ordering that does not require training data. Thus we also propose the following Walsh-Hadamard specific ordering, which is data independent and experimentally yields similarly strong results (see Secs. 3.3.1 and 3.3.2). To sample scenes at a resolution of $m_{1}\times m_{2}$, where $m=m_{1}m_{2}$ is a power of $2$, we define an order on the rows of the Walsh-Hadamard matrix $W_{m}$ that we call the frequency order. The frequency order is so called because it is a generalization of the sequency order that incorporates the 2-dimensional nature of the sampling. In the frequency order, the rows of the Walsh-Hadamard matrix are ordered in increasing order based on the number of connected components when each row is reshaped (via column-major format) to form an $m_{1}\times m_{2}$ matrix. For example, when $m=16,$ the frequency-ordered Walsh-Hadamard matrix has as its first row the vector of all ones; the second row is

[TABLE]

since it reshapes to a matrix with two connected components. Similarly, the remaining rows are ordered to produce increasing numbers of connected components in matrix form; the last row results in 16 components and is $\left[1,0,1,0,0,1,0,1,1,0,1,0,0,1,0,1\right].$ Unlike the sequency order, the frequency order is not unique; for example, in this case two different rows have two connected components. When there is a tie for the number of connected components, we order the tied rows arbitrarily.

As we demonstrate in Sec. 3.3, the frequency order produces good performance consistently on various datasets. It is worth noting, however, that the frequency order is specific to the Walsh-Hadamard (or shifted Walsh-Hadamard) matrix because reshaped vectors naturally have meaningful connected components. If one wished to design a frequency ordering for a sampling matrix other than Walsh-Hadamard, it would likely be necessary to do thresholding or binning in order to produce a similar ordering. Even so, it remains to be determined whether such analogs of the frequency ordering are effective for sampling and reconstruction. In contrast, the maximal-variance ordering generalizes in a natural way for other sampling matrices.

3.3 Performance on Data

We demonstrate the effectiveness of maximal-variance sampling on several datasets of interest. In each case, we simulate CS in software and compare accuracy of reconstructed data against the raw data. We begin in Sec. 3.3.1 with results on a set of LIDAR depth images.

3.3.1 Depth imaging

In this section, we present results that compare maximal-variance sampling against other sampling orders of Walsh-Hadamard matrices on a dataset of over 4000 depth images collected with a Swiss Ranger SR4000 Time of Flight camera at Colorado State University. Each depth image has a spatial resolution of $64\times 64.$

In Fig. 1, we show the number of connected components in the sampling basis vectors corresponding to various orders: standard, maximal-variance for the full Swiss Ranger data, and frequency. Note that the maximal-variance order defined for the Swiss Ranger depth data seems in general to preference sampling vectors with lower numbers of connected components.

In the remainder of this section, we split the Swiss Ranger depth dataset into a training and testing set. The training set was taken to be the first 80% of the images and the testing set was taken to be the remaining 20%. Note that we intentionally split the data in this way as opposed to taking a random sample of 80%. The data is collected sequentially as the camera moves around a scene, so data points that appear within a small time range are generally similar to each other. Thus a random sample would actually provide an unfair advantage in training.

We compare reconstructions after simulated CS against an original test set depth image in Fig.  2. Note that the maximal-variance and frequency orders perform best. The corresponding PSNR values are: standard order (49.4), random order (49.6), sequency order (58.0), frequency order (62.7), and maximal-variance order (62.8).

In Figs. 3 and 4, we provide quantitative comparisons of the performance of sampling and reconstruction with the various sampling orders on the Swiss Ranger depth image test set. Fig. 3 provides a histogram of PSNR values for each of the five sampling orders on the test set. Note that frequency and maximal-variance orders have higher occurrences of relatively large PSNR values. Fig. 4 provides a comparison of PSNR and RMSE on each test image for each of the different sampling orders. Again we see the best performance arising from frequency and maximal-variance sampling. The mean PSNR value across the test set for the maximal-variance order is 67.7, whereas that for standard and random orders is 51.3 and 50.1, representing a 32% and 35% increase, respectively. The mean PSNR value for the frequency order is similar to that for the maximal-variance order, demonstrating that both the maximal-variance order and the frequency order produce strong results on this particular dataset.

3.3.2 Hyperspectral imaging for chemical detection

In this section, we present results for multi/hyper-spectral datasets that contain chemical explosions. The two datasets we consider are the Physical Sciences Inc. Fabry-Pérot interferometer sensor multispectral dataset and the Johns Hopkins Applied Physics Lab FTIR-based longwave infrared sensor hyperspectral dataset. We present analysis of detection of a chemical release in each dataset: a GAA-release in Fabry-Pérot and a release of SF6 in Johns Hopkins. We select a $64\times 64$ spatial subset of the cubes and include all associated spectral bands. In both datasets, the impact of sampling is apparent in the results. Chemical detection is the main motivation for this particular CS application; thus our performance comparison emphasizes detection over other potential metrics like PSNR.

In Fig. 5, we provide a plot of the number of pixels over the associated threshold (defined based on ACE values in background cubes - see Sec. 2.5). The three subfigures show the different results for different variants of the ACE chemical detection algorithm. Moving from left to right, the subfigures show the number of pixels over the threshold for ACE, ACE with bulk coherence, and ACE with bulk coherence and persistence. The rightmost subfigure (which includes the most post-processing) is the easiest to visually examine given the large amount of noise in this particular video.

It is apparent from Fig. 5 that the standard, random, and sequency orders produce artificially high numbers of pixels with ACE values over the threshold. The reason for this is made clear from an examination of an example reconstruction. In Fig. 6, we apply ACE to reconstructions of Fabry-Pérot GAA cube 100 after sampling with various orders. Note that the standard, random, and sequency orders have artifacts and make it impossible to localize the signal. The signal appears strongest in the reconstruction based on the frequency order, but the localization is strongest with the maximal-variance order.

We next consider a release of SF6 in the Johns Hopkins hyperspectral dataset (specifically the SF6 27 Romeo dataset). In this analysis, we include both the maximal-variance order defined on Johns Hopkins SF6 data and that defined on Fabry-Pérot GAA data to emphasize the robustness of the orders to different settings. Even though the chemical and device are both different, the GAA maximal-variance order performs well in the SF6 context.

We chose the SF6 dataset because of the dynamic behavior of the chemical within the $64\times 64$ spatial field of view we selected. First, the chemical is released, then it effectively disappears from the scene, and then it reappears (presumably due to a quick change in wind direction). This made for an interesting performance test - specifically, we were curious about which orderings would lead to successful detection of the chemical in the time frames when the chemical reappears, necessarily in a weaker, dissipated form.

In Fig. 7, we show the number of pixels over the threshold after applying ACE, ACE with bulk coherence, and ACE with bulk coherence and persistence to the reconstructed cubes that were sampled with various orderings of the shifted Walsh-Hadamard basis. As in the GAA example, it is the case that some orderings detect the chemical in a high number of pixels and consequently appear to be producing good performance. However, from example reconstructions (e.g. those in Figs. 8 and 9), it is apparent that the standard, random, and sequency orders produce reconstructions that fail to localize the signal.

From Figs. 8 and 9, we observe that the orderings that produce reconstructions that both detect the presence of SF6 and do a reasonable job of localizing that presence are frequency, maximal-variance (from GAA data), and maximal-variance (from SF6 data). Combining this with the information displayed in Fig. 7, we observe that during the time when the chemical reappears in the scene, the SF6 maximal-variance order does best, then the GAA maximal-variance order, then frequency. (Recall that while standard, random, and sequency have many pixels above the threshold, the reconstructions themselves are not informative.)

3.3.3 Example: maximal-variance order sampling vectors

In this section, we provide example sampling basis vectors in the order defined by maximal-variance for the Swiss Ranger depth data. Note that there is a clear preference in the basis ordering for large-scale features.

We note here that the maximal-variance ordering for the shifted Walsh-Hadamard matrix is more effective if the data is mean-subtracted prior to defining the ordering. To be precise, we subtract the mean of each data point from itself individually; this is somewhat different than what is often meant by mean-subtraction for datasets. In Fig. 10, we show the cumulative percent of variance explained by the ordered shifted Walsh-Hadamard basis for the two cases in which the data has and has not been mean-subtracted. Note that we see a significant improvement in how quickly variance is captured in the case in which data is mean-subtracted. It is also worth pointing out here that it is possible to incorporate this form of mean-subtraction into CS by using the sample corresponding to a vector of all $1$’s.

Mean-subtraction has little effect on the maximal-variance ordering for the Walsh-Hadamard matrix. Additionally, in our explorations, we have seen quite similar orderings produced for the Walsh-Hadamard matrix and the shifted Walsh-Hadamard matrix on mean-subtracted data.

In Fig. 11, we show examples of (re-shaped) sampling vectors that appear early in the order and late in the order to provide some context for the nature of the maximal-variance order. Recall that the vectors are ordered by variance of the corresponding set of dot products of samples. The associated percent of the total variance across all sampling vectors is shown above each reshaped sampling vector. Note that the sampling vectors that show up early appear to capture relatively large-scale features and those that show up late appear to generally capture high frequency information.

4 CONCLUSION

.

In this paper, we presented a method for selecting CS sampling vectors from a sampling basis in a manner that is responsive to training data. Specifically, we defined the maximal-variance sampling order and we demonstrated the benefits of this method of sampling in the context of both depth images and multi/hyper-spectral images. We further proposed and compared results for a ‘frequency’ order that is appropriate for sampling bases such as Walsh-Hadamard bases and requires no training. These sampling methods are particularly valuable in settings with hardware constraints and settings in which context-specific sampling is desired.

We suggest the following questions as directions for further study:

• •

It would be interesting to study methods in which the choice of sampling matrix described in this paper is continually adapted on the fly as new data is acquired. Will this adaptive form of sampling produce better results in applications in which we receive a continuous stream of data?

• •

In this paper we built $S$ starting from a Walsh-Hadamard matrix whose entries had been shifted to only be [math]’s and $1$’s. It would be interesting to explore how this method performs when $S$ is initialized with another matrix.

• •

The method proposed in this paper was only tested with respect to a single optimization problem and a single sparsity-promoting basis. In the future we would like to expand our study to include explorations of how our method performs when other optimization problems are used (for example, minimizing total variation or using custom dictionaries as sparsity-promoting bases).

Acknowledgements.

The authors would like to thank Louis Scharf for insightful discussions related to this work, especially with regard to content involving ACE and MPACE. This research was partially supported by Department of Defense Army STTR Compressive Sensing Flash IR 3D Imager contract W911NF-16-C-0107 and Department of Energy STTR Compressive Spectral Video in the LWIR contract W911SR-17-C-0012.