Contents

1 Introduction

The Gizzly Peak Caldera is located Colorado’s Mineral Belt approximately 15 miles southeast of Aspen. It is the product of a volcanic eruption $\sim 35$ million years ago that ejected some 600 ${\text{ km}}^3$ of magma. With a volcanic explosive index of 7, it was at least four times larger than Tambora.[3, sec 1.5]

Figure 1: West Red (Ruby Mountain) From Enterprise Peak

Coulter sought to identify acidic conditions by the different weathering states of iron oxide, which forms Hematite (red) in acidic conditions, and Goethite (lite yellow) and Jarosite (dark yellow) as pH increases, Jarosite being the most immediate oxidation product of iron sulfide (pyrite).

Figure 2: West Red Iron Endmembers from Aviris:

Red: Hematite	Green: Goethite	Blue: Jarosite
(high pH)		(low pH)

Note on Partial Unmixing:

”Unmixing is critical in imaging spectrometry since virtually every pixel contains some macroscopic mixture of materials. The theory of, and methods for unmixing of spectroscopic signatures are found in a number of sources. Hapke (1993) provides models for linear and non-linear spectral mixing and a discussion of the criteria for using each approach. In Earth remote sensing, a linear mixing model is typically used. Full unmixing which assumes that the spectra of all possible components are known, is described by van der Meer (2000) and Boardman (1989). Since it is often impossible to identify, a priori, all possible components, partial unmixing is an important tool. Match Filtering (MF) (Harsanyi and Chang 1994), Constrained Energy Minimization (CEM – similar to Match Filtering) (Farrand and Harsanyi 1997; Farrand 2001), and Mixture Tuned Match Filtering (${\text{ MTMF}}^{\text{ tm}}$) (Boardman et al. 1995) are commonly used methods for partial unmixing. ${\text{ MTMF}}^{\text{ tm}}$is probably the most popular unmixing method used for geologic remote sensing.”[3, Coulter sec 1.4.4]

(Emphasis Added.)

2 Image Classification and Statistics

2.1 Image Mean Vector and Covariance Matrix

Image Mean Vector and Covariance Matrix [10, White 2005] If image has $P$ bands and $N$ pixels, the mean vector is

where ${\overrightarrow{f}}_j$ is the jth pixel vector of the image

The image covariance matrix $C_{\text{ PxP}}$ is[2, Chang 2013 (6.6)]

where $F_{\text{ PxN}}$ is the matrix of $N$ pixel vectors each of length $P$

$M_{\text{ PxN}}$ is the matrix of $N$ identical mean vectors ($P$ rows by $N$ columns):

where $1_{\text{ PxN}}$ is an $P\times N$ matrix of ones.

2.2 Principal Component Transformation

Principal Component Transformation [7, Smith 1985] Karhunen-Loeve Transformation[11, White 2005] GRASS imagery i.pca: Let

$F_{\text{ PxN}}$

input-image multi-pixel vector ($P$ bands by $N$ pixels)

$M_{\text{ PxN}}$

mean vector matrix,

$Z_{\text{ PxN}}$

output-image multi-pixel vector,

$A_{\text{ PxP}}$

$P\times P$ matrix whose rows are the eigenvectors of the covariance matrix $C_{\text{ PxP}}$, arranged by decreasing magnitude of eigenvalue, as typically returned by SVD routines.

Magnitudes of ${\lambda }_i$ impose ordering on transformed component vectors $Z_{ik}={\sum <!--FUNCTION\; APPLICATION-->}_{j=1}^P{a}_{ij}{X}_{jk}$. Those with the largest ${\lambda }_i$ s.t. ${\lambda }_i∕{\lambda }_{max}>tol$ are the Principal Components. Tolerance tol should be related to the noise floor.

2.3 Minimum Noise Fraction

Minimum Noise Fraction[9, pg. 38] [4, ] We wish to find a particular coefficient matrix $\{a_{ij};i,j=1,\dots <!--FUNCTION\; APPLICATION-->P\}$ that in some sense maximizes the image S/N, assuming the image pixel vectors $\{{\overrightarrow{X}}_k,k=1,\dots <!--FUNCTION\; APPLICATION-->N\}$ are the sum of uncorrelated signal and noise:

Maximize

where ${\lambda }_i$ is generalized eigenvalue of $C_{\text{ PxP}}$ wrt $C_{\text{ PxP}}^N$, and $a_i$ are corresponding generalized eigenvectors. Compare with PCA:

Noise Covariance $C_{\text{ PxP}}^N$ Green suggests $C_{\text{ PxP}}$ be of unit variance and band-to-band uncorrelated.

• Sensor error estimate: (Aviris .rcc file) Call supplied error vector $\{e_j;j=1,\dots <!--FUNCTION\; APPLICATION-->P\}$. Then

  $$C_{lm}^N=e_m^n{\delta }_{lm}(n=1\text{ or }2)$$

  (25)

  is completely uncorrelated. In ideal case all $e_m$ are equal:

  $$\begin{array}{rcll}{C}_{lm}^N& =& eI,\text{ and }& \text{(26)}\\ a_i^T{C}_{\text{ PxP}}{a}_i^T& =& {\lambda }_i{a}_i^T{C}_{\text{ PxP}}^N{a}_i={\lambda }_{i}e\left(a_i^{T}I{a}_i\right)& \text{(27)}\\ & =& {\lambda }_{i}e(a_i^T\cdot a_i)={\lambda }_{i}e={\lambda }_i^{PCA}& \text{(28)}\end{array}$$

  since the eigenvectors $a_i$ are orthonormal. ${\lambda }_i={\lambda }_i^{PCA}$ if the variance $e\equiv 1$.

Homogeneous Area Method [8, sec 2.9.1] If possible find homogenous area of $N_h$ pixels in image:

Local Means and Local Variances [8, sec 2.9.2]

1.

Divide image into small $N_b$ pixel blocks (4x4, 5x5,...)

2.

For each (block, band) get local mean and variance: $$\begin{array}{rcll}{M}_{\text{ L}i}& =& \frac{1}{N_b}\sum _{k=1}^{N_b}{X}_{ik}i=1,2,\dots P\text{ bands}& \text{(34)}\\ {{\sigma }_{\text{ L}i}}^2& =& \frac{1}{N_b-1}\sum _{k=1}^{N_b}{(X_{ik}-M_{\text{ L}i})}^2& \text{(35)}\\ C_{\text{ L}ij}& =& \frac{\left({\delta }_{ij}\right)}{N_k-1}\sum _{k=1}^{N_b}(X_{ik}-M_{\text{ L}i})(X_{jk}-M_{\text{ L}j})& \text{(36)}\end{array}$$

$C$ is the local $N_b\times N_b$covariance matrix.

3.

bin $\left\{{{\sigma }_{\text{ L}i}}^2\right\}$ into classes between band min and max values.

4.

bin with most blocks represents mean noise of image. Hope this bin is same for all bands...

Local Means and Local Variances (con’t)

5.

Suppose ”most popular” bin is a P-cube, each side ranging $[{\sigma }_i^{\ast },{\sigma }_i^{\ast }+△{\sigma }_i^{\ast }]$, contains $N^{\ast }$ points. Then the average value over the bin

$$\overline{C_{\text{ L}ij}}=\frac{1}{N^{\ast }}\sum _{k=1}^{N^{\ast }}{\left(C_{\text{ L}ij}\right)}_k$$

(37)

is desired noise covariance matrix.

6.

Caveat: Assumes image is ”slowly varying enough” that enough of the $N_{\text{ total}}$ blocks are homogeneous, i.e. their covariance is really due to noise, not true features in the image. Blocks and bins must both be ”small enough” – but not too small!

Other methods: ”unsupervised training” derived-endmember classification schemes e.g. LAS’ search ([12, ]) and GRASS’ cluster/maxlik are based upon local covariance minimization.

3 Virtual Dimensionality

We implement the Harsani-Farrand-Chang (HFC) and noise-whitened Harsani-Farrand-Chang (NWHFC) esimators of the image virtual dimensionality VD = $p\le L$ as described in [1, C.-I Chang and Du 2004] and [2, C.-I Chang 2013] Chapter 5.3.6.

3.1 HFC

Let $R_{\text{ LxN}}\equiv [r_1,r_2,\dots <!--FUNCTION\; APPLICATION-->r_N]$ be the input image matrix of $N$ pixel (sample) vectors $r_n$ each of length $L$ spectral bands: $r_i={(r_{i1},r_{i2},\dots <!--FUNCTION\; APPLICATION-->r_{i,L})}^T;1\le i\le N$. The image spectral correlation matrix $R_{\text{ LxL}}$ and covariance matrix $K_{\text{ LxL}}$ are then

Let $\{{\widehat{\lambda }}_1\ge {\widehat{\lambda }}_2\ge \cdots \ge {\widehat{\lambda }}_L\}$ and $\{{\lambda }_1\ge {\lambda }_2\ge \cdots \ge {\lambda }_L\}$ denote the (ordered) correlation and covariance eigenvalues. By assuming that signal sources are nonrandom unknown positive constants and noise is white with zero mean, we may expect

where ${\sigma }_{nl}^2$ is the noise variance in the $l$th spectral channel. Formulate the VD determination as a binary hypothesis problem:

The null hypothesis $H_0$ and the alternative hypothesis $H_1$ represent the case that the correlation eigenvalue is equal to its corresponding covariance eigenvalue (no signal), and the case that the corrleation eigenvalue is greater than its corresponding covariance eigenvalue. When $H_1$ is true (i.e. $H_0$ fails), it implies there is an endmember contributing to the correlation eigenvalue in addition to noise, since the noise energy represented by the eigenvalue of $R_{\text{ LxL}}$ in that particular component is the same as the one represented by the eigenvalue of $K_{\text{ LxL}}$ in its corresponding component.

If we assume the noise in each spectral dimension has zero mean and varance ${\sigma }_{nl}^2$, then ${\widehat{\lambda }}_l={\mu }_l^2+{\sigma }_{nl}^2$ and ${\lambda }_l={\sigma }_{nl}^2$ where ${\mu }_l$ is the sample mean in the $l$th spectral dimension and ${\sigma }_{nl}^2$ is the channel noise. Despite the fact the $\widehat{\lambda }$ and $\lambda$ are unknown constants, we can model each pair of eigenvalues $(\widehat{\lambda },\lambda )$ under hypothesis $H_0$ and $H_1$ as random variables with asymptotic probability densities given by

respectively, where ${\mu }_l$ is an unknown constant for each $l$ and the variance ${\sigma }_{zl}^2$ is given by

“When the total number of samples $N$ is sufficiently large, $\text{ Var}\left[{\widehat{\lambda }}_l\right]\simeq 2{\widehat{\lambda }}_l^2∕N$ and $\text{ Var}\left[{\lambda }_l\right]\simeq 2{\lambda }_l^2∕N$. Therefore, the noise variance ${\sigma }_{zl}^2$ in (42) can be estimated and approximated using (49).” (Although its not obvious exactly how. Play statisitics on collections of subimages, perhaps? In his MATLAB code, Chang approximates ${\sigma }_{zl}^2\simeq (2∕N)({\widehat{\lambda }}_l^2+{\lambda }_l^2)$, so we do as well.)

Chang and Du ([1, eq. 8]) use Schwartz’ inequality to bound

From (46), (47), and (49), we define the false alarm probability and detection probability as

We then choose a given false-alarm probability $P_F$ (for example 0.001) and invert (54) to give

A case of ${\widehat{\lambda }}_l-{\lambda }_l>{\tau }_l$ fails the null hypothesis test, and indicates there is signal energy assumed to contribute to the eigenvalue ${\widehat{\lambda }}_l$ in the $l$th spectral dimension. Do note the threshold ${\tau }_l$ is different for each spectral dimension $l$.

3.2 NWHFC

As noted by Chang ([2, Section 5.3.6]), the signature variance ${\sigma }_{sl}^2$ is generally very small and can be influenced by interband covariances. Interband covariance may be reduced by “noise whitening”: preprocessing by multiplying the input image by a noise covariance matrix estimated by a technique developed by Roger and Arnold, which may be performed prior to the HFC method. This is discussed far too sparsely in ([1, Chang and Du]) and hardly at all in ([2, Chang 2013]), so we’ll merely reference the former and translate Matlab code from the latter (Section A.1).

Noise Estimation :

The noise covariance matrix $K_{\text{ noise}}$ is the LxL diagonal matrix whose elements are the inverse of the diagonal elements of $K_{\text{ LxL}}^{-1}$. For noise whitening we want the square root of its inverse, denoted $K_n^{-(1∕2)}$, also diagonal:

$Y$is the noise-whitened image matrix. The NWHFC VD estimate is obtained simply by running the HFC algorithm on $Y$rather than $R$.

3.3 NSP

The Noise Subspace Projection method ([1, ]) estimates VD solely from the spectral covariance matrix $K_{\text{ LxL}}$. Let the matrix $K_n^{-(1∕2)}$ be as above and define the noise-whitened covariance matrix

As result, the noise variance of each band in the whitened $\overline{K}$ is reduced to unity. Let ${\left\{u_l\right\}}_{l=1}^L$ be a set of eigenvectors generated by $\overline{K}$, which can then be expressed as

where ${\left\{u_l\right\}}_{l=1}^{\text{ VD}}$ and ${\left\{u_l\right\}}_{l=\text{ VD}+1}^L$ span signal subspace and noise subspace respectively. The variances of the noise components of (67) have been whitened and normalized to unity, so ${\overline{\lambda }}_l=1;l=\text{ VD}+1,\dots <!--FUNCTION\; APPLICATION-->,L$ and $\overline{\lambda }>1$ otherwise. The problem of VD estimation can then be formulated as a binary hypothesis testing problem:

where

where ${\mu }_l$ is an unknown constant for each $l$ and the variance ${\sigma }_{yl}^2$ is given by

which can be further reduced under hypothesis $H_0$ to

Finally, find the Neyman-Pearson detector ${\delta }_{\text{ NP}}$ for (68) to determine VD. The false-alarm probability $P_F$ is the probability that we think we’ve detected a signal when none was actually present. In terms of eigenvalues ${\overline{\lambda }}_l$ this the probability of measuring a value ${\overline{\lambda }}_l>1$ when the “true” value was its minimum noise-only value ${\overline{\lambda }}_l=1$. For a desired $P_F$ we seek a threshold ${\tau }_l$ such that

with solution

The virtual dimensionality VD is then the number of eigenvalues ${\overline{\lambda }}_l$ whose value exceeds the threshold $1+{\tau }_l$: ${\overline{\lambda }}_l\ge 1+{\tau }_l;l=1,2,\dots <!--FUNCTION\; APPLICATION-->,\text{ VD}$.

3.4 OSP

Orthogonal Subspace Projection ([2, 5.4.1]). OSP is a general technique that also works well as part of VCA (4).

Assume that there are $p$ signal sources with L bands, $\{s_1,s_2,\dots <!--FUNCTION\; APPLICATION-->s_p\}$, present in the data to be processed and every data sample vector $r_i$ can be expressed by a linear mixture of these p signal sources as

where $S_p=[s_1{s}_2\cdots s_p]$ is a signal matrix made up of the p signal sources, $\{s_1,s_2,\dots <!--FUNCTION\; APPLICATION-->s_p\}$, and $n_i$ can be interpreted as the noise vector or model error vector.

Let $P_p=S_p{\left(S_p^T{S}_p\right)}^{-1}{S}_p^T$ be the p-signal projection matrix formed by the signal matrix $S_p$. It can then be used to map all data sample vectors $\{r_1,r_2,\dots <!--FUNCTION\; APPLICATION-->,r_N\}$ into the spectral space linearly spanned by the $p$ signal sources $\{s_1,s_2,\dots <!--FUNCTION\; APPLICATION-->s_p\}$. In other words, every data sample vector $r_i$ can be expressed by a linear mixture of the $p$ signal sources,$\{s_1,s_2,\dots <!--FUNCTION\; APPLICATION-->s_p\}$, specified by (83) via the $p$-signal projection matrix $P_p$ where the noise $n_i$ in (83) is included to account for the linear mixture model error:

Since each sample vector produces a different model error and residual, the sample mean vector is used to represent an averaged model error and residual. In this case, from (83) the sample mean vector $\mu$ can be expressed by

where ${\overline{\alpha }}_p=(1∕N){\sum <!--FUNCTION\; APPLICATION-->}_{i=1}^N{\alpha }_i$ and $\overline{n}=(1∕N){\sum <!--FUNCTION\; APPLICATION-->}_{i=1}^N{n}_i$. There are two ways to find ${\text{ VD}}^{\text{ OSP}}$.

without involving noise covariance matrix estimation. Theoretically, the value of OSP($p$) in (89) increases as the value of $p$ increases. For any given error threshold $𝜖$, VD can be determined by a stopping rule by $𝜖$. The value of $p$ determining the OSP($p$) in (89) or (90) is denoted by ${\text{ VD}}^{\text{ OSP}}$. Two criteria are developed to detect the abrupt change of the OSP($p$) value. One that is based on the gradient, denoted by “$\nabla <!--FUNCTION\; APPLICATION-->$”, is defined by

The other is based on the difference and is denoted by minus “-”:

For the difference criterion, the sample mean vector $\mu$ is normalized before orthogonal projection so that the values of threshold $𝜖$ for these two criteria are comparable for analysis. The threshold $𝜖$ in (91) and (92) is generally selected according to a sudden drop or a clear gap between two consecutive values of $p$ in plots of the gradient in (91) and difference (92) versus the value of $p$.

It should be noted the above ${\text{ VD}}^{\text{ OSP}}$ definitions involve two key parameters: one is the error threshold $𝜖$,and the other is the algorithm used to produce the $p$-signal matrix $S_p$, which is not specified in (91) and (92). We shall use “VCA” for “algorithm” at least initially, but there are other choices including target-specific algorithms such as ATGP, as well as generic SVD. See ([2, 5.4.1] for further discussion.

4 Vertex Component Analysis

Vertex Component Analysis is an “unsupervised training” derived-endmember classification scheme. We use the notation of [6, Nascimento and Dias], which is slightly different than that used in previous sections.

4.1 VCA Algorithm

Assume linear mixing, and let

$N$ number of pixels in image

$L$ number of spectral bands recorded at each pixel

$p$ number of endmembers, $p\le L$, and usually $p\ll L$ for hyperspectral images.

$m_j$ is an endmember spectral vector of length $L$, $1\le j\le p$

${\gamma }_i$ is a scale factor modeling illumination variability due to surface topography at pixel $i$

${\alpha }_i={[{\alpha }_{i1},{\alpha }_{i2},...{\alpha }_{ip}]}^T$ is the abundance vector containing the fractions of each endmember at pixel $i$. Positivity: ${\alpha }_{ik}\ge 0$ and $1_p^T{\alpha }_i=1$, that is ${\sum <!--FUNCTION\; APPLICATION-->}_{k=1}^p{\alpha }_{ik}=1$.

$x_i$ is the “true” noise-free spectral vector of length $L$ at pixel $i\le N$.

$M_{\text{ Lxp}}=[m_1,m_2,\dots <!--FUNCTION\; APPLICATION-->m_p]$ is an Lxp mixing matrix that maps an abundance vector $\alpha$ to a “true” spectral vector $x$.

Then the recorded spectral vector at pixel $i$ may be given by

Our goal is to find a abundance vectors $\{{\alpha }_i,i=1,\dots <!--FUNCTION\; APPLICATION-->N\}$ corresponding to some endmember set $\{m_j,j=1,\dots <!--FUNCTION\; APPLICATION-->p\}$. An appropriate endmember set $\left\{m\right\}$ is to be determined as part of the VCA algorithm. Endmember identification – matching one or more of the VCA generated endmembers to actual specimen spectral samples such as from USGS spectral library or field ground-truth sampling – may be done in subsequent processing steps.

Since the set $\{\alpha \in {\Re }^p:1^T\alpha =1,\alpha \succcurlyeq 0\}$ is a simplex, the set $S_z=\{z\in {\Re }^L:z=M\alpha ,1^T\alpha =1,\alpha \succcurlyeq 0;\gamma \ge 0\}$ is also a simplex. However, even assuming $n=0$, the observed vector set belongs to convex cone $C_p=\{r\in {\Re }^L:r=M\gamma \alpha ,1^T\alpha =1,\alpha \succcurlyeq 0;\gamma \ge 0\}$ owing to different scale factors $\gamma$.

But the projective projection of the convex cone $C_p$ onto a properly chosen hyperplane is a simplex with vertices corresponding to the vertices of simplex $S_x$. The simplex $S_p=\{y\in {\Re }^L:y=r∕(r^{T}u),r\in C_p\}$ is the projective projection of the convex cone $C_p$ onto the hyperplane $r^{T}u=1$, where the choice of $u$ assures there are no observed vectors $r$ parallel to the hyperplane ( orthogonal to its defining vector $u$): $r_i^{T}u∕\left|u\right|<\left|r_i\right|,i=1,...N$.

Geometry Note : I’m going to suggest there’s just a shade more going on here. Suppose we neglect noise and normalize all the pixel values such that $\left|r_i\right|\equiv 1$, which is easy enough to do. Consider a three band image, such as an rgb color image. the $\left\{r_i\right\}$ are now neatly mapped onto the surface of the first (all axis non-negative) octant of the rgb unit sphere – a pixel’s intensity of any color cannot be negative. For convenience of visualization, assume we may choose $u=(1,1,1)∕\sqrt{3}$ which, if there were pure pixels lying on each unit axis ie (1,0,0), (0,1,0), (0,0,1) assures the simplex of points projected onto the $u$ plane is a pure equilateral triangle whose verteces correspond to the “1” values on the rgb unit axis.

But real photos rarely if ever span the complete rgb bandwidth. Their normalizations will usually map to some sub-patch of that first unit octant. Their projections onto the $u$ plane will still form a simplex, but will it have more vertices than a triangle? The VCA algorithm asserts “no”, and starts by finding the vertices of that projected simplex on the $u$ plane.

4.1.1 SNR

VCA algorithm accuracy is dependent on image SNR. Several ways to estimate:

• Classic definition: SNR = Power(signal)/Power(noise). Select homogeneous blocks from image, areas where the image is expected to be flat so that variations are due to noise. Calculation of data SNR using a (Mean)/(Standard Deviation) method for a homogeneous target. $a_i=(1∕N_b){\sum <!--FUNCTION\; APPLICATION-->}_{k=1}^{N_b}{r}_{ik}$ is the average intensity of the i-th spectral band in the block with $N_b$ pixels. ${\sigma }_i=\sqrt{(1∕N_b){\sum <!--FUNCTION\; APPLICATION-->}_{k=1}^{N_b}{(r_{ik}-a_i)}^2}$ is the standard deviation. Then $SN{R}_i=a_i∕{\sigma }_i$. We may manually select homogeneous regions from the image, such as still water (lakes, ponds) or snowfields, or we may use GRASS i.cluster/i.maklik to (perhaps) find such regions automatically.
• VCA does an image SVD or PCA as an initial dimensionality reduction step. We might try looking at the ratio of the Power in the SVD (and/or PCA) components with large singular values – the ones of interest – with the Power in the components with small singular values. Compare this with the ratio of the sum of the large eigenvalues with the sum of the small eigenvalues – or sums of squares of each, to see if there’s an easy way to compute an “effective SNR” that is usable to select between SVD and PCA in the first step of the VCA algorithm.
• Compare the above with sensor data e.g. AVRIS .rcc files Aviris has very high SNR, over 500 on a good day[5, Sec. 3.3] Satellites e.g. Hyperion might have a tenth this, SNR is an important consideration and worth some synthetic image modelling. Nascimento and Dias use $10^{1.5}L$ for their ${\text{ SNR}}_{\text{ th}}$ and SNR estimator. We note with an SNR greater than 500 (on a good day), AVIRIS is considered extemely good.[5, Kruse JPL 2002] But with 224 bands it will spectacularly fail this particular test, as $10^{1.5}\times 224=31.62\times 224=7083$. We speculate this test is insufficent for any hyperspectral image, and something better related to the number of significant principal components (usually d) might be in order. Nascimento and Dias have quite a bit on the topics of estimating SNR and the virtual dimensionality (VD) of an image. See [6, eqs 10 - 13] and C.-I Chang and Q. Du, Estimation of number of spectrally distinct signal sources in hyperspectral imagery, IEEE Trans. Geosci. Remote Sens., vol. 42, no. 3, pp. 608-619, Mar. 2004. [1, C.-I Chang and Du].

4.1.2 Algorithm 1: Vertex Component Analyis (VCA)

Inputs:

$R_{\text{ LxN}}\equiv [r_1,r_2,\dots <!--FUNCTION\; APPLICATION-->r_N]$ is the image matrix of $N$ pixel vectors $r_n$ each of length $L$ spectral bands.

$p$, desired dimensionality of problem, the number of endmembers. $p$ may be estimated as the number of “large” eigenvalues of the image covariance matrix. $p\le L$.

${\text{ SNR}}_{\text{ th}}$ is the fixed threshold value $10^{1.5}L$ of the SNR measured wrt the signal subspace. Different values might be used for different SNR estimators. See above.

Outputs:

$\widehat{M}$ is the $L\times p$ estimated mixing matrix returned by VCA.

Notations:

${\left[\widehat{M}\right]}_{:,j}$ is the jth column of $\widehat{M}$

${\left[\widehat{M}\right]}_{:,i:k}$ are the ith through kth columns of $\widehat{M}$

${\left[X\right]}_{:,\left[index\right]}$ is the matrix formed from the columns of $X$ listed in index vector $\left[index\right]$

$A^{+}$ denotes the pseudo-inverse of pxp matrix $A$ as computed by SVD. $A^{+}$ is what you think it is. Write the SVD as $A=U\mathrm{\Sigma }{V}^T$ where $\mathrm{\Sigma }$ is the pxp diagonal matrix whose diagonal elements are the (ordered) singular values: $\mathrm{\Sigma }=\text{ diag}\{w_1,w_2,\dots <!--FUNCTION\; APPLICATION-->w_p\}$ Then the pseudo-inverse may be written $A^{+}=V{\mathrm{\Sigma }}^{+}{U}^T$ where ${\mathrm{\Sigma }}^{+}=\text{ diag}\{w_1^{+},w_2^{+},\dots <!--FUNCTION\; APPLICATION-->w_p^{+}\}$ and the pseudo-inverse of any scalar $w^{+}=0$ if $w\simeq 0$ and $w^{+}=w^{-1}$ otherwise. Thus $A^{+}=A^{-1}$ if $A$ is non-singular. If $P=A{A}^{+}$ and $Q=A^{+}A$, then

• $P$ is the orthogonal projector onto the range of $A$
• $Q$ is the orthogonal projector onto the range of $A^T$
• $(I-P)$ is the orthogonal projector onto the kernel (null-space) of $A^T$
• $(I-Q)$ is the orthogonal projector onto the kernel (null-space) of $A$
• $P=A{A}^{+}=U\mathrm{\Sigma }{V}^{T}V{\mathrm{\Sigma }}^{+}{U}^T=U\mathrm{\Sigma }{\mathrm{\Sigma }}^{+}{U}^T$ and $\mathrm{\Sigma }{\mathrm{\Sigma }}^{+}=\text{ diag}\{1,1,\dots <!--FUNCTION\; APPLICATION-->0,\dots <!--FUNCTION\; APPLICATION-->\}$ has d non-zero diagonal entries, $d\le p$. $P$ is a full pxp symetric matrix, none of whose components contain values from lower right p-d block of $U$. If p=d then $P$ and $Q$ are the identity matrix.

See Moore-Penrose_pseudoinverse In particular “ the pseudoinverse for matrices related to $A$ can be computed by applying the Sherman-Morrison-Woodbury formula to update the inverse of the correlation matrix, which may need less work. In particular, if the related matrix differs from the original one by only a changed, added or deleted row or column, incremental algorithms[12][13] exist that exploit the relationship.” – a fact that might be quite useful here.

$U_d$ is an Lxd matrix whose columns are the first d eigenvectors of either $R^{T}R$ (which is LxL) or the LxL spectral covariance matrix $C$. $U_d$ is the left eigenmatrix computed by SVD whose columns are ordered (left to right) in order of decreasing singular value. Its first d columns are the most significant (principal) relative to the rest. We usually select d such that the singular values ${\omega }_i:\{i>d\}$ are in some sense “too small to be significant”.

1: Compute SNR and Virtual Dimensionality $p$
2: if $(\text{ SNR}>{\text{ SNR}}_{\text{ th}})$ then