I Problem statement
In many signal processing applications, the signals of interest do not span the entire observation space and a relevant and frequently used assumption is that they evolve in a low-dimensional subspace . Subspace modeling is accurate when the signals consist of a linear combination of modes in a
-dimensional space, and constitute a good approximation for example when the signal covariance matrix is close to rank-deficient. As a consequence, subspace estimation plays a central role in recovering these signals with maximum accuracy. An ubiquitous solution to this problem is to resort to the singular value decomposition (SVD) of the data matrix. The SVD emerges naturally as the maximum likelihood estimator in the classical model, where stands for the observation matrix, is the (deterministic) matrix, with , whose columns span the -dimensional subspace of interest, is the (deterministic) waveform matrix and is the additive noise. The
principal left singular vectors ofprovide very accurate estimates of a basis for the range space of , and have been used successfully, e.g., in estimating the frequencies of damped exponentials or the directions of arrival of multiple plane waves, see [2, 3]
among others. However, the SVD can incur some performance loss in two main cases, namely when the signal-to-noise ratio (SNR) is very low and thereof the probability of a subspace swap or subspace leakage is high[4, 5, 6, 7]. A second case occurs when the number of samples is lower than the subspace dimension : indeed, is at most of rank and information is lacking about how to complement in order to estimate .
Under such circumstances, a Bayesian approach might be helpful as it enables one to assist estimation by providing some statistical information about . We investigate such an approach herein and assign to the unknown matrix an appropriate prior distribution, taking into account the specific structure of . The paper is organized as follows. In section II, we propose an approach based on minimizing a natural distance on the Grassmann manifold, which yields a new estimator of . The theory is illustrated in section III where the new estimator is derived for some specific examples. In section IV its performance is assessed through numerical simulations, and compared with conventional approaches. Section V studies an application to the analysis of interactions between pure materials contained in hyperspectral images.
Ii Minimum mean square distance estimation
In this section, we introduce an alternative to the conventional minimum mean square error (MMSE) estimator, in the case where a subspace is to be estimated. Let us consider that we wish to estimate the range space of
from the joint distributionwhere stands for the available data matrix. Usually, one is not interested in per se but rather in its range space , and thus we are operating in the Grassmann manifold , i.e., the set of -dimensional subspaces in . It is thus natural to wonder whether the MMSE estimator (which is the chief systematic approach in Bayesian estimation ) is suitable in . The MMSE estimator of a vector minimizes the average Euclidean distance between and , i.e., . Despite the fact that this distance is natural in an Euclidean space, it may not be the more natural metric in . In fact, the natural distance between two subspaces and is given by  where are the principal angles between these subspaces, which can be obtained by SVD of where and denote orthonormal bases for these subspaces . The SVD of is defined as , where and are two unitary matrices. Therefore, it seems more adequate, rather than minimizing as the MMSE estimator does, to minimize the natural distance between the subspaces spanned by and . Although this is the most intuitively appealing method, it faces the drawback that the cosines of the angles and not the angles themselves emerge naturally from the SVD. Therefore, we consider minimizing the sum of the squared sine of the angles between and , since for small , . As argued in [8, 10], this cost function is natural in the Grassmann manifold since it corresponds to the Frobenius norm of the difference between the projection matrices on the two subspaces, viz . It should be mentioned that our approach follows along the same principles as in  where a Bayesian framework is proposed for subspace estimation, and where the author considers minimizing . Hence the theory presented in this section is similar to that of , with some exceptions. Indeed the parameterization of the problem in  differs from ours and the application of the theory is also very different, see the next section.
Given that , we define the minimum mean-square distance (MMSD) estimator of as
it follows that
Therefore, the MMSD estimate of the subspace spanned by is given by the largest eigenvectors of the matrix , which we denote as
In other words, MMSD estimation amounts to find the best rank- approximation to the posterior mean of the projection matrix on . For notational convenience, let us denote . Except for a few cases where this matrix can be derived in closed-form (an example is given in the next section), there usually does not exist any analytical expression for . In such situation, an efficient way to approximate the matrix is to use a Markov chain Monte Carlo simulation method whose goal is to generate random matrices drawn from the posterior distribution , and to approximate the integral in (4) by a finite sum. This aspect will be further elaborated in the next section. Let
denote the eigenvalue decomposition ofwith and . Then the average distance between and is given by
The latter expression constitutes a lower bound on and is referred to as the Hilbert-Schmidt bound in [11, 12]. As indicated in these references, and similarly to , this lower bound may be difficult to obtain analytically.
The MMSD approach can be extended to the mixed case where, in addition to , a parameter vector which can take arbitrary values in needs to be estimated jointly with . Under such circumstances, one can estimate and as
Doing so, the MMSD estimator of is still be given by (4) while the MMSD and MMSE estimators of coincide.
The MMSD approach differs from an MMSE approach which would entail calculating the posterior mean of , viz . Note that the latter may not be meaningful, in particular when the posterior distribution depends on only through , see next section for an example. In such a case, post-multiplication of by any unitary matrix yields the same value of . Therefore averaging over does not make sense while computing (4) is relevant. On the other hand, if depends on directly, then computing the posterior mean of can be investigated: an example where this situation occurs will be presented in the next section. As a final comment, observe that is not necessarily unitary but its range space can be used to estimate .
We open a parenthesis here regarding the framework of this paper. Although it is not directly related to this paper (we do not address optimization problems here) it is interesting to note the recent growing interest in optimization problems on special manifolds, especially on the Stiefel manifold (the set of matrices such that ) and the Grassmann manifold, see the excellent tutorial paper by Edelman et al.  as well as [13, 14], and [15, 16, 17] for signal processing applications. These references show the interest of taking into account the underlying geometry of the problem, as we attempt to do herein.
Iii Illustration examples
In this section we illustrate the previous theory on some examples, including the conventional linear Gaussian model (conditioned on ) and a model involving the eigenvalue decomposition of the data covariance matrix. As a first step, we address the issue of selecting prior distributions for and then move on to the derivation of the MMSD estimator.
Iii-a Prior distributions
A crucial step in any Bayesian estimation scheme consists of selecting the prior distribution for the variables to be estimated. We focus here on distributions on the Stiefel or Grassmann manifold, depending whether we consider the matrix itself or its range space. There exist only a few distributions on the Stiefel or Grassmann manifolds, the most widely accepted being the Bingham or von Mises Fisher (vMF) distributions [18, 19], which are given respectively by
where stands for the exponential of the trace of the matrix between braces, is an symmetric matrix, is an arbitrary matrix, and , are hypergeometric functions of matrix arguments, see e.g.  for their definitions. We will denote these distributions as and , respectively. Observe that the Bingham distribution depends on only, and can thus be viewed as a distribution on the Grassmann manifold [18, 19] while the vMF distribution depends on and is a distribution on the Stiefel manifold. In our case, in order to introduce some knowledge about , we assume that it is “close” to a given subspace spanned by the columns of an orthonormal matrix , and hence we consider two possible prior distributions for , namely
where means “proportional to”. The distribution in (9) is proportional to the sum of the squared cosine angles between and while is proportional to the sum of the cosine angles between and . Note that is a concentration parameter: the larger the more concentrated around are the subspaces . The difference between the two distributions is the following. In the Bingham distribution only and are close (at least for large values of ) since is invariant to post-multiplication of by any unitary matrix . Hence is not necessarily close to . In contrast, under the vMF prior distribution, and are close. For illustration purposes, Figure 1 displays the average fraction of energy of in defined as
As can be observed from these figures, both distributions allow the distance between and to be set in a rather flexible way. Their AFE is shown to be identical for small values of the concentration parameter but, when increases, the AFE of the vMF distribution increases faster.
Iii-B Linear model
In order to illustrate how the previous theory can be used in practice, we first consider a simple example, namely a linear Gaussian model (conditioned on ), i.e., we assume that the data follows the model where the columns of are independent and identically distributed (i.i.d.) Gaussian vectors with zero-mean and (known) covariance matrix . We assume that no knowledge about is available and hence its prior distribution is set to . Therefore, conditioned on we have
When follows the Bingham prior distribution, the posterior distribution of , conditioned on is given by
which is recognized as a Bingham distribution with parameter matrix , i.e., . For such a Bingham distribution, it turns out that the eigenvectors of coincide with those of , with the same ordering of their eigenvalues, see Appendix A for a proof. Therefore the MMSD estimator is obtained in closed-form as
Therefore, the MMSD estimator has a very simple form in this case. It consists of the principal subspace of a (weighted) combination of the a priori projection matrix and the information brought by the data through . Observe that, in this particular case of a Bingham posterior, the MMSD estimator coincides with the maximum a posteriori (MAP) estimator.
Let us now consider the case where the prior distribution of is vMF, and contrast it with the previous example. Using (III-B) along with along with (10), it follows that the posterior distribution now writes
which is referred to as the Bingham-von-Mises-Fisher (BMF) distribution with parameter matrices , and respectively111The matrix is said to have a distribution -where is an symmetric matrix, is a diagonal matrix and is an matrix- if .. Although this distribution is known , to our knowledge, there does not exist any analytic expression for the integral in (4) when has the BMF distribution (15). Therefore, the MMSD estimator cannot be computed in closed-form. In order to remedy this problem, a Markov chain Monte Carlo simulation method can be advocated [20, 21] to generate a large number of matrices drawn from (15), and to approximate (4) as
In (16), is the number of burn-in samples and is the number of samples used to approximate the estimator. An efficient Gibbs sampling scheme to generate random unitary matrices drawn from a distribution with arbitrary full-rank matrix was proposed in . It amounts to sampling successively each column of by generating a random unit norm vector drawn from a (vector) BMF distribution. In our case, whose rank is and hence is rank-deficient whenever , a case of most interest to us. Note also that to generate matrices drawn from the Bingham distribution in (9), we need to consider which has rank . Therefore, the scheme of  needs to be adapted in order to generate random matrices drawn from (15). In Appendix B, we review the method of  and show how it can be modified to handle the case of a rank-deficient matrix .
Interestingly enough, the above estimator in (16) is the so-called induced arithmetic mean (IAM)  of the set of unitary matrices , . It differs from the Karcher mean of the set , , which truly minimizes the sum of the distances to all . However, the Karcher mean may not exist and requires iterative schemes to be computed  while the IAM is straightforward to compute.
In the particular case where has a Bingham prior distribution, the MAP estimator of and its MMSD estimator are equal. This is no longer true when has a vMF prior distribution, and hence a BMF posterior distribution. The mode of the latter is not known in closed-form either. However, it can be approximated by selecting, among the matrices generated by the Gibbs sampler, the matrix which results in the largest value of the posterior distribution.
Iii-C Covariance matrix model
We now consider a more complicated case where , conditioned on and
, is Gaussian distributed with zero-mean and covariance matrix
where is an orthonormal basis for the signal subspace, is the diagonal matrix of the eigenvalues and
stands for the white noise power which is assumed to be known here. As it will be more convenient and more intuitively appealing, we re-parametrize the covariance matrix as follows. The inverse ofcan be written as
where , with and
The idea is to parametrize the problem in terms of and rather than and . The interest of this transformation is twofold. First, it enables one to express all eigenvalues with respect to the white noise level. Indeed, one has where is an orthonormal basis for and hence the s are representative of the scaling between the “signal” eigenvalues and the noise eigenvalues. In fact, they carry information about the signal-to-noise ratio since and represents the SNR of the -th signal component. Second, this new parametrization will facilitate derivation of the conditional distributions required for the Gibbs sampler.
Since conditioned on and is Gaussian, it follows that
From , it ensues that and hence
Let us now consider the prior distributions for and . We will consider either a Bingham or vMF distribution for . As for , we assume that, i.e.,
The value of [respectively ] can be set to [respectively ] if a non-informative prior is desired. Otherwise, if some information is available about the SNR, and can be chosen so as to reflect this knowledge since : [resp. ] rules the lowest [resp. highest] value of the SNR, say [resp. ].
With the Bingham assumption for , the joint posterior distribution of and is
where is the incomplete Gamma function. Unfortunately, the above distribution does not belong to any known family and it is thus problematic to generate samples drawn from it. Instead, in order to sample according to (III-C), we propose to use a Gibbs sampler drawing samples according to and for . From (III-C), the conditional distribution of is
which is recognized as a (modified) Bingham distribution222
Let us now turn to the conditional distribution of . From (III-C) one has
which is the product of independent gamma distributions with parametersand , truncated in the interval . We denote this distribution as . Random variables with such a distribution can be efficiently generated using the accept-reject scheme of .
The above conditional distributions can now be used in a Gibbs sampler, as described in Table I. When has a vMF prior distribution, it is straightforward to show that , conditioned on and , follows a BMF distribution while the posterior distribution of is still given by (III-C). Therefore line 2 of the Gibbs sampler in Table I just needs to be modified in order to handle this case.
In this section we illustrate the performance of the approach developed above through Monte Carlo simulations. In all simulations , and . The matrix is generated from a Gaussian distribution with zero-mean and covariance matrix and the signal-to-noise ratio is defined as . The matrix is generated from the Bingham distribution (9) or the vMF distribution (10) and, for the sake of simplicity, . The number of burn-in iterations in the Gibbs sampler is set to and . The MMSD estimator (4) is compared with the MAP estimator, the MMSE estimator, the usual SVD-based estimator and the estimator that discards the available data and use only the a priori knowledge. The latter is referred to as “Ubar” in the figures. The estimators are evaluated in terms of the fraction of energy of in , i.e., .
Iv-a Linear model
We begin with the linear model. Figures 4 to 7 investigate the influence of and onto the performance of the estimators. Figures 4 and 5 concern the Bingham prior while the vMF prior has been used to obtain Figures 6 and 7. From inspection of these figures, the following conclusions can be drawn:
the MMSD estimator performs better than the estimator , even at low SNR. The improvement is all the more pronounced that is large. Therefore, the MMSD estimator makes a sound use of the data to improve accuracy compared to using the prior knowledge only.
the MMSD estimator performs better than the SVD, especially at low SNR. Moreover, and this is a distinctive feature of this Bayesian approach, it enables one to estimate the subspace even when the number of snapshots is less than the size of the subspace .
for a Bingham prior, the MMSE performs very poorly since the posterior distribution of conditioned on depends on only. Hence, averaging the matrix itself does not make sense, see our remark 1. In contrast, when has a vMF prior, the posterior depends on both and : in this case, the MMSE performs well and is close to the MMSD. Note however that the vMF prior is more restrictive than the Bingham prior.
the MMSD estimator also outperforms the MAP estimator.
As a conclusion, the MMSD estimator performs better than most other estimators in the large majority of cases.
Iv-B Covariance matrix model
We now conduct simulations with the covariance matrix model. The simulation parameters are essentially the same as in the previous section, except for the SNR. More precisely, the random variables are drawn from the uniform distribution in (22) where and are selected such that dB and dB. The results are shown in Fig. 8 for the Bingham prior and Fig. 9 for the vMF prior. They corroborate the previous observations made on the linear model, viz that the MMSD estimator offers the best performance over all methods.
V Application to hyperspectral imagery
In this section, we show how the proposed subspace estimation procedure can be efficiently used for an application to multi-band image analysis. For several decades, hyperspectral imagery has received considerable attention because of its great interest for various purposes: agriculture monitoring, mineral mapping, military concerns, etc. One of the crucial issue when analyzing such image is the spectral unmixing which aims to decompose an observed pixel into a collection of reference signatures, (called endmembers) and to retrieve the respective proportions of these signatures (or abundances) in this pixel . To describe the physical process that links the endmembers and their abundances to the measurements, the most widely admitted mixing model is linear
where is the pixel spectrum measured in spectral bands, () are the endmember spectra and () are their corresponding abundances. Due to obvious physical considerations, the abundances obey two kinds of constraints. Since they represent proportions, they must satisfy the following positivity and additivity constraints
Let now consider pixels
of an hyperspectral image induced by the linear mixing model (LMM) in (28) with the abundance constraints (29). It is clear that the dataset formed by these pixels lies in a lower-dimensional subspace . More precisely, in this subspace , the dataset belongs to a simplex whose vertices are the endmembers to be recovered. Most of the unmixing strategies developed in the hyperspectral imagery literature are based on this underlying geometrical formulation of the LMM. Indeed, the estimation of the endmembers is generally conducted in the lower-dimensional space
, previously identified by a standard dimension reduction technique such as the principal component analysis (PCA). However, it is well known that the model linearity is a simplifying assumption and does not hold anymore in several contexts, circumventing the standard unmixing algorithms. Specifically, non-linearities are known to occur for scenes including mixtures of minerals or vegetation. As a consequence, evaluating the suitability of the LMM assumption for a given hyperspectral image is a capital question that can be conveniently addressed by the approach introduced above.
V-a Synthetic data
First, we investigate the estimation of the subspace when the image pixels are non-linear functions of the abundances. For this purpose, a synthetic hyperspectral image is generated following a recently introduced non-linear model referred to as generalized bilinear model (GBM). As indicated in , the GBM is notably well adapted to describe non-linearities due to multipath effects. It assumes that the observed pixel spectrum can be written
where stands for the Hadamard (termwise) product and the abundances () satisfy the constraints in (29). In (30), the parameters (which belong to ) characterize the importance of non-linear interactions between the endmembers and in the -th pixel. In particular, when (), the GBM reduces to the standard LMM (28). Moreover, when (), the GBM leads to the non-linear model introduced by Fan et al. in . In this simulation, the synthetic image has been generated using the GBM with endmember signatures extracted from a spectral library. The corresponding abundances have been uniformly drawn in the set defined by the constraints (29). We have assumed that there is no interaction between endmembers and , and between endmembers and resulting in , . Moreover, the interactions between endmembers and are defined by the map of coefficients displayed in Fig. 10 (top, left panel) where a black (resp. white) pixel represents the lowest (resp. highest) degree of non-linearity. As can be seen in this figure, of the pixels (located in the bottom and upper right squares of the image) are mixed according to the LMM resulting in . The remaining image pixels (located in the upper left square of the image) are mixed according to the GBM with nonlinearity coefficients radially increasing from to ( in the image center and in the upper left corner of the image). Note that this image contains a majority of pixels that are mixed linearly and belong to a common subspace of . Conversely, the non-linearly mixed pixels do not belong to this subspace333Assuming there is a majority of image pixels that are mixed linearly is a reasonable assumption for most hyperspectral images.. We propose here to estimate the local subspace where a given image pixel and its nearest spectral neighbors live ( denotes the set of the ()-nearest neighbors of ).
Assuming as a first approximation that all the image pixels are linearly mixed, all these pixels are approximately contained in a common -dimensional subspace that can be determined by performing a PCA of (see  for more details). The corresponding principal vectors spanning are gathered in a matrix . This matrix is used as a priori knowledge regarding the -dimensional subspace containing . However, this crude estimation can be refined by the Bayesian estimation strategy developed in the previous sections. More precisely, for each pixel , we compute the MMSD estimator of the matrix , whose columns are supposed to span the subspace containing and its -nearest neighbors . The Bayesian estimator is computed from its closed-form expression (14), i.e., using the Bingham prior where has been introduced above. Then, for each pixel, we evaluate the distance between the two projection matrices and onto the subspaces and , respectively. As stated in Section II, the natural distance between these two projection matrices is given by . The resulting distance maps are depicted in Fig. 10 (bottom panels) for non-zero values of (as it can be noticed in (14
), this hyperparameterbalances the quantity of a priori knowledge included in the estimation with respect to the information brought by the data). For comparison purpose, the subspace has been also estimated by a crude SVD of (top right panel). In this case, simply reduces to the associated principal singular vectors and can be considered as the MMSD estimator of obtained for .
These figures show that, for the of the pixels generated using the LMM (bottom and right parts of the image), the subspace estimated by an SVD of the whole dataset
is very close to the hyperplaneslocally estimated from through the proposed approach (for any value of ). Regarding the remaining pixels resulting from the GBM (top left part of the image), the following comments can be made. When a crude SVD of is conducted, i.e., when no prior knowledge is taken into account to compute the MMSD (, top right panel), the distance between the locally estimated subspace and the a priori assumed hyperplane does not reflect the non-linearities contained in the image. Conversely, when this crude SVD is regularized by incorporating prior knowledge with and (bottom left and right panels, respectively), leading to the MMSD estimator, the larger the degree of non-linearity, the larger the distance between and . To summarize, evaluating the distance between the MMSD estimator and the a priori given matrix allows the degree of non-linearity to be quantified. This interesting property is exploited on a real hyperspectral image in the following section.
V-B Real data
The real hyperspectral image considered in this section has been acquired in 1997 over Moffett Field, CA, by the NASA spectro-imager AVIRIS. This image, depicted with composite true colors in Fig. 11 (top, left panel), has been minutely studied in  assuming a linear mixing model. The scene consists of a large part of a lake (black pixels, top) and a coastal area (bottom) composed of soil (brown pixels) and vegetation (green pixels), leading to endmembers whose spectra and abundance maps can be found in . A simple estimation of a lower-dimensional space where the pixels live can be conducted through a direct SVD of the whole dataset, providing the a priori matrix . As in the previous section, this crude estimation can be refined by computing locally the MMSD estimators spanning the subspaces (bottom panels). These estimators have been also computed with , corresponding to an SVD of (top, right figure). The distances between and have been reported in the maps of Fig. 11. Again, for (top, right panel), a simple local SVD is unable to locate possible non-linearities in the scene. However, for two444Additional results obtained with other values of are available online at http://dobigeon.perso.enseeiht.fr/app_MMSD.html. non-zero values and (bottom left and right panels, respectively), the distances between the a priori recovered subspace and the MMSD-based subspace clearly indicate that some non-linear effects occur in specific parts of the image, especially in the lake shore. Note that the non-linearities identified by the proposed algorithm are very similar to the ones highlighted in  where the unmixing procedure was conducted by using the GBM defined in (30). This shows the accuracy of the proposed MMSD estimator to localize the non-linearities occurring in the scene, which is interesting for the analysis of hyperspectral images.
This paper considered the problem of estimating a subspace using some available a priori information. Towards this end, a Bayesian framework was advocated, where the subspace is assumed to be drawn from an appropriate prior distribution. However, since we operate in a Grassmann manifold, the conventional MMSE approach is questionable as it amounts to minimizing a distance which is not the most meaningful on the Grassmann manifold. Consequently, we revisited the MMSE approach and proposed, as an alternative, to minimize a natural distance on the Grassmann manifold. A general framework was formulated resulting in a novel estimator which entails computing the principal eigenvectors of the posterior mean of . The theory was exemplified on a few simple examples, where the MMSD estimator can either be obtained in closed-form or requires resorting to an MCMC simulation method. The new approach enables one to combine efficiently the prior knowledge and the data information, resulting in a method that performs well at low SNR or with very small sample support. A successful application to the analysis of non-linearities contained in hyperspectral images was also presented.
Appendix A The eigenvalue decomposition of
The purpose of this appendix is to prove the following proposition which can be invoked to obtain the MMSD estimator whenever the posterior distribution is a Bingham distribution.
Let be an orthogonal matrix -
be an orthogonal matrix -- drawn from a Bingham distribution with parameter matrix
with . Let denote the eigenvalue decomposition of where the eigenvalues are ordered in descending order. Let us define . Then the eigenvalue decomposition of writes
with and where .
For notational convenience, let us work with the projection matrix whose distribution on the Grassmann manifold is 
We have then that
Moreover is diagonal since, for any orthogonal diagonal matrix ,
where, to obtain the third line, we made use of the fact that . It follows that the eigenvectors of and coincide, and that the eigenvalues of are , for . Moreover, it is known that and, from (32), one has
Differentiating the latter equation with respect to and denoting , one obtains