Sparse Estimation of the Precision Matrix and Plug-In Principle in Linear Discriminant Analysis for Hyperspectral Image Classification

ABSTRACT

In this paper, a new method for supervised classification of hyperspectral images is proposed for the case in which the size of the training sample is small. It consists of replacing in the Mahalanobis distance the maximum likelihood estimator of the precision matrix by a sparse estimator. The method is compared with two other existing versions of LDA sparse, both in real and simulated images.

Keywords:
hight dimensionality; linear discriminant analysis; hyperspectral image

1 INTRODUCTION

Hyperspectral sensors collect information from the earth’s surface in hundreds of bands of the electromagnetic spectrum whit relatively narrow bandwidths. Each pixel of a hyperspectral image can be regarded as a high-dimensional vector, $\mathit{x}=\left(x_1,\dots ,x_p\right)\text{'}$, called feature vector, whose entries x _j , 1 ≤ j ≤ p corresponds to the spectral reflectance collected in the band j of the electromagnetic spectrum.

Image classification is an essential step wich consists in classification rules that assign every data point x into one of K classes based on a features vector. In supervised image classification, a classifier is trained on a training data set of size n. Since the collection of such training data is either expensive or time-consuming, there is usually only a small amount of labeled data available.

The huge number of features (p) of hyperspectral data and the limited availability of the training sample data (n) greatly worsen the traditional classification methods performance such as linear discriminant analysis (LDA) as ²¹ emphasize. This problem is typical in high-dimensional data analysis.

For example, if p > n the sample covariance matrix S is singular and LDA cannot be applied. If n is not much greater than p, S can be calculated but due to the large number of estimated parameters (p(p + 1)/2), it will have a significant amount of sample error, and its inverse will be a highly biased estimator of Σ^-1(¹). On the other hand, because the bandwidth decreases as p grows, many bands are highly correlated and hyperspectral data can be contain a significant amount of redundant spectral information. Recently, two algorithms have been proposed called sparse discriminant analysis (SDA) (⁵) and penalized LDA (²⁰) which simultaneously solve the problem of singularity of the covariance matrix (when p > n) and allow variable selection by imposing a penalty of type ℓ₁ to the optimal scoring approach from which one can derive the LDA rule and the discriminant vectors in Fisher’s discriminant problem respectively. In the context of high dimensionality, the estimation of the covariance matrix and its inverse, the precision matrix, is already a research area of important results. Most of the proposals rely on sparsity assumptions and requires some kind of regularization strategy. In this framework, it is assumed that many of the matrix entries to be estimate are zeros. (¹⁰; ⁸; ¹⁵; ⁹). Using estimation of sparse precision matrices and plug-in principle in the Mahalanobis distance allows to extend LDA to high dimensional problems ⁷. This strategy, will be called plug-in LDA.

On the other hand, we are interested in non-paametric tree-based classification methods introduced by Breiman, Friedman, Olshen and Stone in mid 1980’s. Therefore, we can expect this approach outperforms LDA when the decision boundary is highly non-linear. However these methods can lead to a phenomenon known as overfitting the data, which essentially means that on the training data set they follow the errors too closely having a very poor perfomance on new observations. For this reason trees can be very sensitive to small changes in the training sample dramatically affecting the final estimated tree ¹¹.

The aim of this paper is to compare the performance of SDA and penalized - LDA versus plugin-LDA using Kashlak proposal based on precision matrix estimation ¹³. Since random forest (RF) is a very popular tree-based classification methods it will be also included in this paper for comparative purposes ²; ⁶; ¹⁴.

This paper is organized as follows. We review LDA, SDA and penalized LDA in Section 2 and in Section 3 we present our proposal, the plug-in-LDA. A simulation study and applications to real hyperspectral images are presented in Section 5. Section 6 contains a discussion.

2 REVIEW OF LDA AND PENALIZED LDA

There are three distinct arguments that result in the LDA classifier: the multivariate Gaussian model, Fisher’s discriminant problem, and the optimal scoring problem.

Consider K populations (classes) where a p-dimensional random vector x is defined and assume that $\mathit{x}~N\left({\mathit{\mu }}_k,{\Sigma }_w\right)$ where µ _k is the mean vector for class k and Σ_w is a pooled within-class covariance matrix common to all K classes. The Bayesian approach classifies a new observation x* to the class k that maximizes Mahalanobis’s distance

${\delta }_k\left(\mathit{x}*\right)=\left(\mathit{x}*-{\hat{\mathbf{\mu }}}_k\right)\text{'}{\hat{\Sigma }}_w^{-1}\left(\mathit{x}*-{\hat{\mathit{\mu }}}_k\right)-2\log {\pi }_k$ (2.1)

where ${\hat{\mathbf{\mu }}}_k$ and ${\hat{\Sigma }}_w$ denote maximum likelihood estimates for µ _k and Σ_w , obtained from a training sample $\left({\mathit{x}}_1^{\text{'}},\dots {\mathit{x}}_n^{\text{'}}\right)$

LDA, considered as arising from Fisher’s discriminant problem, consists in finding K vectors, called discriminant vectors, ${\hat{\mathit{\beta }}}_1,\dots {\hat{\mathit{\beta }}}_K$ such that sequentially solve

$$\begin{gathered} {\text{maximize}}_{{\mathit{\beta }}_k\in {\mathrm{ℝ}}^p}\left\{{\mathit{\beta }}_k^{\text{'}}{\hat{\Sigma }}_b{\mathit{\beta }}_k\right\} \\ \text{subject to} \\ {\mathit{\beta }}_k^{\text{'}}{\hat{\Sigma }}_w{\mathit{\beta }}_k=1\text{and}{\mathit{\beta }}_j^{\text{'}}{\hat{\Sigma }}_w{\mathit{\beta }}_k=0,\forall j<k, \end{gathered}$$ (2.2)

where ${\hat{\Sigma }}_b$ is the standard estimate for the between-class covariance matrix. Since ${\hat{\Sigma }}_b$ has rank at most K − 1 then there exits at most K − 1 non trivial solutions to (2.2), ${\hat{\mathit{\beta }}}_1,\dots ,{\hat{\mathbf{\beta }}}_{K-1}$ called discriminant vectors. the classification rule is reduced to projecting a new observation on the space defined by the discriminant vectors and assigning it to the class whose projected mean is closest.

A third formulation that yields the LDA classification rule is optimal scoring, introduced by ⁵, that transform the classification problem into a regression problem by solving

$$\begin{gathered} {\text{minimize}}_{{\mathit{\theta }}_k,{\mathit{\beta }}_k}\left\{{||Y{\mathit{\theta }}_k-\mathit{X}{\mathit{\beta }}_k||}^2\right\} \\ \text{subject to}{\mathit{\theta }}_k^{\text{'}}Y\text{'}Y{\mathit{\theta }}_k=1. \end{gathered}$$ (2.3)

where Y = (y _ik ) denotes a n × K matrix with y _ik = 1 if the i-th observation belongs to the k−th class and y _ik = 0 otherwise. Moreover, θ _K is a scores vector of dimension K and Yθ represents the scores vector of the training data. The solution ${\hat{\mathbf{\beta }}}_k$ to (2.3) are proportional to solution to (2.2) (see ¹⁹).

Among the proposals to extend LDA to the problems of high dimensionality, those that produce more interpretable models are of particular interest. We need that these extensions allow variable selection and classification simultaneously. In linear regression model, lasso is a regularization method that, imposing a penality at the size of the estimated parameters vector, shrunks some coefficients estimates exactly to zero, producing variable selection. Using this approach ⁵ propose SDA by imposing to (2.3) an additional penalization:

$$\begin{gathered} {\text{minimize}}_{{\theta }_k,{\beta }_k}\left\{{||Y{\mathit{\theta }}_k-\mathit{X}{\mathit{\beta }}_k||}^2+\mathit{\gamma }{\mathit{\beta }}_k^{\text{'}}\Omega {\mathit{\beta }}_k+\lambda {||{\mathit{\beta }}_k||}_1\right\} \\ \text{subject to}{\mathit{\theta }}_k^{\text{'}}Y\text{'}Y{\mathit{\theta }}_k=1,{\mathit{\theta }}_k^{\text{'}}Y\text{'}Y{\mathit{\theta }}_i=0,\forall i<k \end{gathered}$$ (2.4)

where Ω is a positive definite matrix and γ and λ are nonnegative tuning parameters.SDA produce sparse discriminant vectors since if λ is large, some entries of β _k will be zero.

Penalized Fisher’s linear discriminant, proposed by ²⁰, is based on imposing an ℓ₁ penalty on equation (2.2):

$$\begin{gathered} {\text{maximize}}_{{\mathit{\beta }}_k}\left\{{\mathit{\beta }}_k^{\text{'}}{\hat{\Sigma }}_b{\mathit{\beta }}_k-\lambda {||{\mathit{\beta }}_k||}_1\right\} \\ \text{subject to} \\ {\mathit{\beta }}_k^{\text{'}}{\hat{\Sigma }}_w{\mathit{\beta }}_k=1\text{and}{\mathit{\beta }}_j^{\text{'}}{\hat{\Sigma }}_w{\mathit{\beta }}_k=0,\forall j<k, \end{gathered}$$ (2.5)

where ${\tilde{\Sigma }}_w$ is some full rank estimate for Σ_w , such as ${\tilde{\Sigma }}_w=diag\left({\widehat{{\sigma }_1}}^2,\cdots ,{\widehat{{\sigma }_p}}^2\right)$, where ${\widehat{{\sigma }_j}}^2$ is the j-th diagonal element of ${\hat{\Sigma }}_w$. This method also produce sparse discriminant vectors. Although there is a correspondence between the critical points of problems (2.4) and (2.5), the solutions are not necessarily the same ²⁰.

3 PLUG-IN LDA

In high dimensionality, to estimate the precision matrix associated to a multivariate Gaussian distribution it usual to assume that it is sparse. The most commonly estimation methods employ an ℓ₁-penalized maximum likelihood (¹⁰; ⁸; ⁴). As an alternative, ¹³ propose a novel distribution free estimator that controls the false positive rate in the selection of the non zero entries of the estimated precision matrix. This proposal begins with an initial estimator ${\hat{\Omega }}_0$- such as debiased graphical lasso (¹²) or debiased ridge estimator (³) and iteratively using a binary search procedure locates the densest estimator near I _p or alternatively the sparsest estimator close to ${\hat{\Omega }}_0$. Considering these sparse estimates and plug-in principle in the Mahalanobis distance we propose a new version of LDA sparse which we will call plug-in LDA.

4 ACCURACY ASSESSMENT MEASURES

In selecting the best classification algorithm a very important task is the choice of an appropriate performance evaluation measures. Accuracy assessment is traditionally conducted using the sample confusion matrix, in which classification results of the validation dataset are compared to “ground truth.” Based on the confusion matrix, a variety of measures can be calculated, including Cohen’s Kappa coefficient ⁸ and the F ₁ score¹⁸

Cohen kappa statistic is a chance-corrected method for assessing agreement. The values of range lie in [−1; 1] with 1 presenting complete agreement and 0 meaning no agreement or independence. A negative statistic implies that the agreement is worse than random.

Table 1:
The confusion matrix for a multi-class classification involving k classes.

	True class	1	2	...	k	Total
Predicted class
1		n 11	n 12	...	n 1k	n 1.
2		n 21	n 22	...	n 2k	n 2.
...		...	...	...	...	...
k		n k1	n k2	...	n kk	n k.
Total		n .1	n .2	...	n .k	N

Using notation similar to Cohen (1960), the kappa coefficient of agreement,(κ), is estimated from confusion matrix:

$\hat{\kappa }=\frac{p_0-p_e}{1-p_e}$ (4.1)

where p ₀ is the proportion of cases correctly classified (i.e. overall accuracy) and p _e is the expected proportion of cases correctly classified by chance (under independency hipotesis):

$p_0=\frac{{\displaystyle \sum _{i=1}^k}{n}_{ii}}{N}$ (4.2)

$p_e=\frac{{\displaystyle \sum _{i=1}^k}{n}_{i.}{n}_{.i}}{N^2}$ (4.3)

Precision and recall (sensitivity) are two important model evaluation measures. For a multiclass problem we can compute the per-class precisio´n and recall as follows:

$recal{l}_i=\frac{n_{ii}}{n_{.i}};1\le i\le k$ (4.4)

$precisio{n}_i=\frac{n_{ii}}{n_{i.}};1\le i\le k$ (4.5)

The F ₁ value is used to combine the precision and recall measurements into a single value. This is practical because it makes it easier to compare the combined performance of precision and completeness between various algorithms. F1 score is a harmonic mean of precision and recall:

${\left(F_1\right)}_i=\frac{2}{{\displaystyle \frac{1}{recal{l}_i}}+{\displaystyle \frac{1}{precisio{n}_i}}};1\le i\le k$ (4.6)

The weighted F ₁ is defined as

$w{F}_1=\frac{{\displaystyle \sum _{i=1}^k}{n}_{.i}F{1}_i}{N}$ (4.7)

To get a more accurate estimate of the measurements, the training data set was randomly divided into three groups, or folds, of equal size (n = p). The first fold is treated as a training set, and remaining folds as validation set. This process results in three estimates of each measure, which are averaged in order to obtain a unique value ¹¹.

5 APPLICATIONS

Experiments were conducted to test the performance of the SDA; penalized - LDA; plug-in-LDA and RF algorithms with one simulated image and one real image.

For SDA, Penalized-LDA and RF we will use the R-packages sda, sparseMatEst and randomForest, respectively.

5.1 Simulated Data

The first experimental dataset (simulated image) consist of 100 × 100 pixels, 200 bands and four classes. Each pixel x _ij ; 1 ≤ i, j ≤ 100 is a p-dimensional random vector, where p represent the number of bands. We assume that ${\mathit{x}}_{ij}~N\left({\mathit{\mu }}_k;\Sigma \right)$ if pixel (i, j) belongs to class k, 1 ≤ k ≤ 4. Figure 1a and 1b represent an ideal class image with four regions, acting as the reference image, and the RGB-composition of simulated random image assuming a normal distribution. The size of the training set is 160 (40 pixels of each class). Table 2 shows the Frobenius norm and kappa estimated coefficient for the five classifiers. The validation is made with the entire image. The Figure 2 and the Table 2 show that, penalized-LDA has the best performance followed by SDA and plug-in-LDA, while RF has the worst performance.

Figure 1:
(a) Ideal class image with four regions. (b) RGB-composition of simulated random image assuming a normal multivariate distribution.

Figure 2:
Classification results simulated image.

Table 2:
Frobenius norm. (kappa estimated coefficient).

training sample	RF	SDA	Pen-LDA	Plug-in LDA
	36.3	11.3	6.3	17.4
	(0.953)	(0.995)	(0.998)	(0.989)

5.2 Real Data

Figure 3 shows an 209 x 167 extract of an EO-1 Hyperion image obtained in http://eo1.usgs.gov. The study area is EO1H2290822012007110P1₁ T . 44 non-informative bands were removed, as well as water absorption bands 120-132, 165-182, 185-187, 221-224. So, 160 bands were established to be useful for further analysis. By prior knowledge of the area, it was decided to classify the image into four classes: water, urban area, high vegetation index area, and low vegetation index area (blue, white, green and brown, respectively). The size of the training sample is 160.

Figure 3:
Real image: false color composite image (R-G-B = band 120-60-2).

Figure 5 show the thematic maps produced by the five classifiers. When executing the LDA function of the R-package MASS, the following error message appears: ”variables are collinear”. This problem, caused by the insufficient size of the training sample to estimate the covariance matrix, produces an unacceptable thematic map (see Figure 5a).

Figure 4:
Google maps extract (a) city area (b) mountain area.

Figure 5:
Classification results real image.

If the algorithms are ranking based on adequacy measures presented in the Table 3, the first one are SDA and penalized-LDA, followed by plug-in-LDA and finally RF. A visual comparison of the thematic maps (see Figure 5c, 5d, and 5e) with the truth of the terrain (Figure 4a and 4b), seems to show that SDA had a better performance, being very similar to penalized-LDA and plug-in-LDA, while RF had the worst performance.

Table 3:
Kappa estimated coeficient and weighted F ₁ score averaged over k = 3 folds.

	LDA	RF	SDA	Pen-LDA	Plug-in LDA
wF1 score	0.43	0.90	0.94	0.94	0.92
$\hat{\kappa }$	0.34	0.87	0.92	0.92	0.90

6 CONCLUSION

This paper proposes a new version of sparse LDA, called plug-in LDA, for the supervised classification of hyperspectral images. It consists of replacing the maximum likelihood estimator of the precision matrix used in Mahalanobis distance with a sparse estimate recently proposed. The results obtained show that plug-in LDA outperforms to RF and has a similar behavior to SDA and Penalized LDA, both in the real image and in the simulated one.

In the context of hiperspectral image classification, extracting totally pure training samples is difficult, since there are often incorrectly labeled pixels, wich can have a strong negative impact on the classification result ¹⁶; ¹⁷. The main advantage of the plug-in strategy is the possibility of replacing in the Mahalanobis’s distance the estimators of both the position and the covariance matrix with robust alternatives which could turn plug-in into a robust classification algorithm.

REFERENCES

• ¹
  J. Bai & S. Shi. Estimating High Dimensional Covariance Matrices and Its Applications. Annals of Economics and finance, 12(2) (2011), 199-215.
• ²
  L. Breiman. Random forests. Machine Learning, 1 (2001), 5-32.
• ³
  P. Buhlmann. Statistical significance in high-dimensional linear models. Bernoulli, 19 (2015), 1212-1242.
• ⁴
  T. Cai & W. Liu. Adaptive thresholding for sparse covariance matrix estimation. Journal of the American Statistical Association, 106 (2011), 672-684.
• ⁵
  L. Clemmensen, T. Hastie, D. Witten & B. Ersbøll. Sparse Discriminant Analysis. Technometrics, 53(4) (2011), 406-413. doi:10.1198/TECH.2011.08118.
  » https://doi.org/10.1198/TECH.2011.08118
• ⁶
  T.N. Do, P. Lenca, S. Lallich & N.K. Pham. Classifying Very-High-Dimensional Data with Random Forests of Oblique Decision Trees. Advances in knowledge discovery and management, 1 (2009).
• ⁷
  J. Fan, Y. Feng & Y. Wu. Network exploration via the adaptive LASSO and SCAD penalties. The Annals of Applied Statistics, 3(2) (2009), 521-541. doi:10.1214/08-AOAS215.
  » https://doi.org/10.1214/08-AOAS215
• ⁸
  O.B.L. Ghaoui & A. d’Aspremont. Model selection through sparse maximum likelihood estimation for multivariate gaussian or binary Data. Journal of Machine Learning Research, 9 (2008), 485-516.
• ⁹
  Y.L. J. Fan & H. Liu. An overview on the estimation of large covariance and precision matrices. The Econometrics Journal, 19(1) (2016), 1-34. doi:https://doi.org/10.1111/ectj.12061
  » https://doi.org/10.1111/ectj.12061
• ¹⁰
  T.H. J. Friedman & R. Tibshirani. Sparse inverse covariance estimation with the graphical lasso. Biostatistics, 9(3) (2007), 432-441.
• ¹¹
  G. James, D. Witten, T. Hastie & R. Tibshirani. “An Introduction to Statistical Learning with Applications in R”. Springer, 2 ed. (2017).
• ¹²
  J. Jankova & S. Van De Geer. Confidence intervals for high-dimensional inverse covariance estimation. Electronic Journal of Statistics, 9 (2015), 1205-1229.
• ¹³
  A. Kashlak. Non-asymptotic error controlled sparse high dimensional precision matrix estimation. arXiv:1903.10988 [stat.ME], (2019).
• ¹⁴
  C. Li. The Application of high-dimensional Data Classification by Random Forest based on Hadoop Cloud Computing Platform. Chemical Engineering Transactions, 51 (2016).
• ¹⁵
  W. Liu & X. Luo. Fast and adaptive sparse precision matrix estimation in high dimensions. Journal of Multivariate Analysis, 135 (2015), 153-162. doi:https://doi.org/10.1016/j.jmva.2014.11.005 URL https://www.sciencedirect.com/science/article/pii/S0047259X14002607
  » https://doi.org/10.1016/j.jmva.2014.11.005» https://www.sciencedirect.com/science/article/pii/S0047259X14002607
• ¹⁶
  M. Popov, S. Alpert, V. Podorvan, M. Topolnytskyi & S. Mieshkov. Method of Hyperspectral Satellite Image Classification under Contaminated Training Samples Based on Dempster-Shafer’s Paradigm. Central European Researchers Journal, 1 (2015), 82-93.
• ¹⁷
  L. Qingming, Y. Haiou & J. Junjun. Label Noise Cleansing with Sparse Graph for Hyperspectral Image Classification. Remote Sensing, 11 (2019), 1116. doi:10.3390/rs11091116.
  » https://doi.org/10.3390/rs11091116
• ¹⁸
  Y. Sasaki. The truth of the F-measure. Teach Tutor Mater, (2007), 1-5.
• ¹⁹
  A.B. T. Hastie & R. Tibshirani. Penalized Discriminant Analysis. The Annals of Statistics, 23(1) (1995), 73-102.
• ²⁰
  D. Witten & R. Tibshirani. Penalized classification using Fisher’s linear discriminant. J R Stat Soc Series B Stat Methodol, 73(54) (2011), 753-772. doi:10.1111/j.1467-9868.2011.00783.x.
  » https://doi.org/10.1111/j.1467-9868.2011.00783.x
• ²¹
  H. Zou. Classification with high dimensional features. WIREs computational Statistics, 11 (2019). doi:https://doi.org/10.1002/wics.1453
  » https://doi.org/10.1002/wics.1453

Publication Dates

• Publication in this collection
  05 Sept 2022
• Date of issue
  Jul-Sep 2022

History

• Received
  30 Sept 2021
• Accepted
  24 Mar 2022