REDUCING THE DIMENSIONALITY OF DATA: LOCALLY LINEAR EMBEDDING OF SLOAN GALAXY SPECTRA

Over the last decade, the PCA or K–L classification has been applied to a wide range of spectral classification problems. From nonparametric classification (Connolly et al. 1995) to extracting star formation properties (Ferreras et al. 2006) to identifying supernovae within SDSS spectra (Madgwick et al. 2003b; Krughoff et al. 2009), the many applications of PCA have demonstrated the importance of dimensionality reduction in classification. This technique seeks to decompose a data set into a linear combination of a small number of principal components, such that each principal component describes the maximum possible variance:

$$\begin{equation} S(\lambda) = \sum _i a_i e_i(\lambda), \end{equation} \tag{ 1 }$$

where S(λ) is the spectrum to be projected, a_i are the expansion coefficients, and e_i(λ) are the orthogonal eigenspectra.

Schematically, PCA is equivalent to finding a best-fit linear subspace to the entire data set, such that the variance of the data projected into this subspace is maximized. Yip et al. (2004a) present a robust PCA approach to spectral classification of the SDSS spectra. In it, they find that the information within the SDSS main galaxy sample can be almost completely characterized by the first dozen eigenmodes. This represents a reduction in data size by a factor of nearly 400 over the complete spectra. Furthermore, the coefficients of the first three eigencomponents correlate with physical properties of the galaxies such as star formation rate and post-starburst activity. In fact, PCA decomposition can lead to a single parameter description of galaxy spectra: the mixing angle ϕ between the first and second eigencoefficients. This mixing angle has been shown to correlate well with the galaxy spectral type (Connolly et al. 1995; Madgwick et al. 2003a) and can be used for an approximate separation between quiescent and emission-line galaxies.

For continuum emission, PCA has a proven record in representing the variation in the spectral properties of galaxies. It does not, however, perform well when reconstructing high-frequency structure within a spectrum (i.e., the distribution of emission lines, lines ratios, and line widths). This is due to two effects: principal components preserve the total variance of the system, i.e., the overall color of the spectrum, and principal components express linear relations between components. High-frequency, or local, features do not contribute significantly to the total variance and are, therefore, not represented in the primary eigencomponents. Variations in line widths and, often, line ratios are also inherently nonlinear in nature (due to the impact of dust or variations in the galaxy mass) and are not compactly represented as a linear combination of orthogonal components. This is shown in Yip et al. (2004a), where the majority of galaxies can be represented by three eigenspectra but emission-line galaxies require up to eleven components to express their variance. PCA performs poorly when distinguishing between emission-line galaxies, e.g., separating broad-line QSOs from narrow-line QSOs or star-forming galaxies.

Line-ratio diagrams, sometimes known as BPT plots (Baldwin et al. 1981) or Osterbrock diagrams (Osterbrock & De Robertis 1985), are one answer to the problems associated with PCA classification. Whereas PCA classifies on the basis of the total variance within a spectrum, which weights the continuum emission more than individual emission lines, line-ratio diagrams are sensitive to emission-line characteristics while ignoring continuum properties. This makes them appropriate for diagnostics that distinguish between galaxies with narrow emission lines. Line-ratio diagrams, however, do not account for continuum emission that might provide additional information on the physical state of a galaxy. They are ill-suited to classify nonemitting galaxies and galaxies with a low signal-to-noise ratio (S/N) for individual lines, and fail completely for galaxies with certain types of emission, e.g., broad-line QSOs (often defined as those with emission line FWHM larger than ∼1200 km s⁻¹). Line-ratio diagrams, then, are useful for the classification of only a small subset of observed objects.

Clearly, PCA and line-ratio diagrams serve complimentary purposes. PCA takes into account broad, low-frequency features, and can distinguish galaxies based on their continuum properties. Line-ratio diagrams, by design, depend only on the detailed features of emission lines, and therefore distinguish spectra based on their narrow, high-frequency features. If both types of information could be combined into a joint classification scheme, one that treats the spectrum as a whole rather than trying to isolate individual features, more insight may be gained into the physical state of a given object. This would be particularly helpful for automated classification within spectroscopic surveys, where a single, coherent technique could be used to identify objects that could then be analyzed further by specialized pipelines. The requirements of such a general classification are: it must be sensitive to both low and high frequency spectral information, it must be able to account for nonlinear relationships between the spectral properties, it should be robust to outliers within the data, and it should be able to express the properties of a galaxy in terms of a small number of physically motivated parameters. In this paper, we consider LLE as a solution to these questions.

For notational conventions, refer to Appendix A. We initially consider a series of points in a high-dimensional space x_i (in the case of SDSS, each spectrum is represented by point in a D_in = 4000 dimensional space) such that the set of these points, X, is given by $\mathbf {X} = [\mathbf {x_0},\mathbf {x_1},\mathbf {x_2},\ldots, \mathbf {x_N}],\ \mathbf {x_i} \in \mathbb {R}^{D_{{\rm in}}}$. We map this to a lower dimensional space $\mathbf {Y} = [\mathbf {y_0},\mathbf {y_1},\mathbf {y_2},\ldots,\mathbf {y_N}],\ \mathbf {y_i} \in \mathbb {R}^{D_{\rm out}}$, with D_in>D_out. For each point x_i, let n⁽ⁱ⁾ = [n⁽ⁱ⁾₁, n⁽ⁱ⁾₂, ..., n⁽ⁱ⁾_K]^T be the indices of the nearest neighbors, such that $\mathbf {x}_{n_j^{(i)}}$ is the jth nearest neighbor of x_i. Also, let w⁽ⁱ⁾ = [w⁽ⁱ⁾₁, w⁽ⁱ⁾₂, ..., w⁽ⁱ⁾_K]^T represent the reconstruction weights associated with the K nearest neighbors of x_i.¹

The key assumption is that each point x_i lies near a locally linear, low-dimensional manifold, such that a tangent hyperplane is locally a good fit. If this is the case, a point can be accurately represented as a linear combination of its K nearest neighbors, with K ⩾ D_out. The error in reconstruction can be thought of as the distance from this manifold to the point in question. A convenient form of this error is the reconstruction cost function

$$\begin{equation} \mathcal {E}_1^{(i)}(\mathbf {w}^{(i)}) = \Bigg | \mathbf {x_i} - \sum _{j=1}^K w_j^{(i)}\mathbf {x}_{ n_j^{(i)} }\Bigg |^2. \end{equation} \tag{ 2 }$$

By minimizing $\mathcal {E}_1^{(i)}(\mathbf {w}^{(i)})$ subject to the constraint

$$\begin{equation} \sum _j w_j^{(i)} = 1, \end{equation} \tag{ 3 }$$

we find the optimized local mapping. This local mapping has two important properties: first, it is invariant under scaling, translation, and rotation; second, it encodes the relationship between points of the local neighborhood in a way that is agnostic to the dimensionality of the parameter space; i.e., it equally well encodes the properties of the neighborhood in observational parameter space, as well as the properties of the neighborhood as projected onto the local tangent hyperplane. These two facts lie at the core of the LLE algorithm.

Once the weights w⁽ⁱ⁾ are determined for each point x_i, the same weights are used to determine the projected vectors y_i that minimize the global cost function:

$$\begin{equation} \mathcal {E}_2(\mathbf {Y}) = \sum _{i=1}^N \Bigg | \mathbf {y_i} - \sum _{j=1}^K w_j^{(i)}\mathbf {y}_{n_j^{(i)}}\Bigg |. \end{equation} \tag{ 4 }$$

Note the symmetry inherent in Equations (2) and (4). Intuitively, the fact that the input data X and the projected data Y optimize these linear forms shows that the local neighborhoods of X and Y have similar properties.

Computationally, each of these steps can be implemented efficiently using optimized linear algebra methods (see, for example, Roweis & Saul 2000; de Ridder & Duin 2002). We discuss algorithmic considerations in Appendices B–D.

We provide a simple and classical example of the performance of the algorithm in Figure 1. The first plot shows a particularly simple nonlinear test set: the "S-curve," with N = 2000 data points in three dimensions uniformly selected from the two-dimensional bounded manifold described by

$$\begin{equation} \begin{array}{ll} \left. \begin{array}{l@{}} \dsty x = \sin (\theta)\\ \dsty z = \frac{\theta }{|\theta |}(\cos (\theta)-1))\\ \end{array} \right\rbrace & -1.5\pi \le \theta \le 1.5\pi \\ 0<y<5. & \\ \end{array} \end{equation} \tag{ 5 }$$

A linear dimensionality reduction such as PCA would not recover the inherently nonlinear shape of the two-dimensional manifold: it would find three principal axes within the data. An ideal nonlinear reduction would, in effect, unwrap this manifold and map it to a flat plane that preserves local clustering, i.e., the colors would remain unmixed.

The second plot in Figure 1 shows the two-dimensional LLE projection applied to these data, with K = 15. It is clear that LLE, in this case, successfully recovers the desired two-dimensional embedding. Points of similar colors (which are clustered within the higher dimensional space) are clustered in the projected space. As discussed in Section 5.2, the effectiveness of this projection will depend on the size of the region that we consider local (i.e., the number of nearest neighbors) as well as the assumed dimensionality of the embedded manifold.

The above test data densely sample a simply connected manifold. In the case of sparsely sampled, or nonconnected manifolds (e.g., due to missing data and survey detection limits), LLE is less successful at recovering the lower dimensional projection. This issue can be addressed using various modifications to the algorithm (see, for example, Donoho & Grimes 2003; Chang & Yeung 2006; Zhang & Wang 2007), some of which are implemented in the publicly available source code for this work (see Appendix E).

The distance of a galaxy along any one of the branches within Figure 2 correlates well with its physical characteristics (as reflected in their spectra). Figures 3–7 show progressions along various regions, along with the mean spectrum of each labeled region. When calculating the mean, we exclude excessively noisy objects, identified as those with flux F(λ) satisfying ∫|F(λ)|dλ>1.5∫F(λ)dλ. Because the spectra are preprocessed with a normalization, the few spectra satisfying this inequality would otherwise contribute disproportionately to the mean spectrum in each region.

Zoom In Zoom Out Reset image size

Figure 3 shows the progression of spectra along the emission galaxy track. To do this, we calculate the mean spectrum for galaxies contained within the regions associated with the bounding boxes E1–E5. As shown later, this sequence tracks not only increasing line strengths, but also changing line ratios (see Figure 12 later).

Figure 4 shows the progression of spectra along the broad-line QSO track. We find that, as we move from Q1 to Q5, the spectral type of the QSO evolves from a Seyfert 1.9 to a classic broad line QSO. As would be expected for such a transition, the N ii emission line strength decreases, while the Hβ line strength increases as the spectra become progressively dominated by the accretion disk emission. Associated with these emission-line properties we find that the continuum slope of the QSO spectra is bluer for the Q5 class than for Q1. This not only supports the classification of the spectra from Seyfert 1.9 to broad-line active galactic nucleus (AGN), it also demonstrates how LLE can use jointly the information contained within the continuum and line emission.

Zoom In Zoom Out Reset image size

Figure 5 shows the progression of spectra along the narrow-line QSO track. As with the emission-line galaxy track, the narrow QSO track traces increasing emission from N1 to N5. N1–N3 have narrow Hβ and Hα, with line ratios that indicate power-law ionization, and are spectrally consistent with LINERs. N3–N5 have increasingly broad Hα, which suggests classification as Seyfert 1.9–2.0. Line-ratio diagnostics confirm this (see Figure 12 later).

Zoom In Zoom Out Reset image size

In Figure 6, we consider the progression of spectra across the different spectral groupings (from broad-line QSOs to galaxies with narrow emission lines) as opposed to along an individual track. The X3 sources are, as before, broad-line QSOs with strong Hβ emission. As we transition to X2, the width of the emission line decreases until X2 has the spectral characteristics of a Seyfert type 1.5. By X1, the spectral properties are consistent with an H ii region or narrow emission line galaxy (for wavelengths <5500 Å) while showing evidence for an AGN component in the broad-line Hα. This one parameter sequence from broad-line QSO to H ii region spectra (including regions where the populations are mixed) is described very naturally with this LLE projection.

Zoom In Zoom Out Reset image size

The last subsample we consider are the quiescent galaxies in Figure 7. The transition from G1 to G5 tracks the same spectral classification from quiescent to active star formation as found using the PCA analysis of Yip et al. (2004a). In terms of spectral properties, galaxies classified as G1 correspond to the E/S0 spectral type from Kennicutt (1998) and galaxies in G5 correspond to the Sbc/Sc Kennicutt spectrum. Along this progression, the emission-line strengths for Hα and N ii increase and the strength of the Mg iii absorption line decreases. It is also noted that the depth of the break at 4000 Å decreases with increasing index number. Comparisons with the results of Yip et al. (2004a) indicate that this trajectory maps directly to the star formation rate within the galaxy (see Figure 11 later).

Zoom In Zoom Out Reset image size

From these initial comparisons, it is clear that by projecting onto a three-dimensional space using LLE, we can differentiate between line-emitting galaxies and those dominated by continuum emission. We can also distinguish between sources as a function of their emission-line characteristics. As such, LLE encapsulates both the line emission and continuum emission properties in its classification scheme. That LLE can express the variation in continuum properties of galaxies when projecting from 1000 to 3 dimensions is not altogether surprising as the PCA analyses of Connolly et al. (1995), Yip et al. (2004a), and others, have demonstrated that quiescent galaxies occupy only a small region of the available parameter space and that the dominant direction in this space is governed by the star formation rate. For emission-line galaxies, and in particular QSOs, standard PCA approaches require between 10 and 50 components to express the variation in line strengths as a function of galaxy type. This is due to the inherently nonlinear nature of the reconstruction of emission lines.

The fact that LLE can be used to separate both broad and narrow emission lines into separate groups or families, that the transition between narrow and broad line spectra is smooth, and that the position of a galaxy within each of these groups is dependent on the individual line strengths demonstrates the ability of LLE to concisely represent inherently nonlinear systems in an intuitive manner.

LLE is useful in that it maps a large-dimensional data set to a smaller dimensional parameter space. Within this projected space, many familiar classification schemes can be applied (see Hastie et al. 2001 for a thorough introduction to the topic). In the case of SDSS spectra, LLE succeeds in mapping physically similar spectra to the same region of parameter space. More importantly, physically different spectra are mapped to distinct, well-separated regions, leading to the possibility of a very intuitive classification scheme.

Because of this, we can quickly locate objects that are misclassified by both the SDSS pipeline and the line-diagnostic plots. Figure 8 displays a few examples of these. M1 is classified by the SDSS pipeline as a broad-line QSO, but falls on the emission-line galaxy track. This misclassification is likely due to contamination of the Hα line from the neighboring Nitrogen features, leading to an erroneous large line width. M2 is not recognized as an emitting object, because the Hα feature is missing. In fact, the entire narrow track on which M2 lies consists of objects that look like emission galaxies/type-2 Seyferts, with the Hα line missing due to saturation of the lines, coincident sky lines, and bad pixels. LLE isolates these points on their own track of outliers. M3 and M4 fall on the broad-line QSO track, but are not recognized as broad emitters by the SDSS pipeline, probably due to high background noise. The LLE projection correctly groups these with other strong emitters. Because the LLE projection is entirely neighborhood-based, any given point will have similar properties to nearby points in the projection space.

Zoom In Zoom Out Reset image size

Next, we compare the classification based on LLE with the classification by more traditional methods, e.g., emission line widths, line strengths, and line ratios. This requires a choice of classification technique for points in the LLE projection space. For illustrative purposes, we create a simple classification scheme and define two well-separated regions that correspond to physically different objects. The "broad-line QSO region" is defined as all points on the broad-line track with LLE component three greater than 0.02. The "emission galaxy region" is defined as all points on the emission galaxy track with LLE component one greater than 0.015 (see Figure 2).

Out of the 624 galaxies in the broad-line QSO region, 31 were classified by SDSS as having narrow line emission, and 11 were classified as quiescent. Through a visual inspection, 29 of the 31 and 1 of the 11 were found to show broad emission features, with the rest of the spectra too noisy to accurately classify. The SDSS pipeline, therefore, misclassifies about 5% of their broad-line sources as narrow line.

Out of the 675 galaxies in the emission galaxy region, 45 have line ratios on the narrow-line QSO side of the Kewley et al. (2001) cutoff, two have Hα emission broader than 1200 km s⁻¹, and four have prominent emission features not detected by the SDSS pipeline.

Combining the above tallies, we can estimate how successfully LLE clusters physically similar spectra. From a subsample that by design contains a high fraction of abnormal spectra (see Section 5.3), less than one percent of the spectra were physically dissimilar from the average population in each region. These errors are primarily due to very noisy spectra. For other cases, e.g., the 45 spectra where classifications based on LLE and Kewley et al. (2001) line diagnostics disagree, it is not immediately clear which classification is correct, if any. The majority of these sources are transitional, falling near the dividing boundary of the line-ratio diagram. This dividing line is an arbitrary analytic separation between two regions of a certain subset of the parameter space. We note, however, that the LLE-based classification is without any a priori training. We do not preselect sources based on known properties nor fine-tune the parametrization of the classification, as opposed to the optimized SDSS pipelines.

Figure 9 shows a sample of some significant outliers identified with LLE. O1–O3 are outliers due to their noise. O4 and O6 appear to be sky spectra that were misclassified based on emission features. O5 appears to be a QSO with its Hα line missing, possibly due to atmospheric absorption. O7 has an artificial feature likely due to faulty sky-noise subtraction. In all of these cases, the spectrum's status as an outlier based on LLE reflects its unusual spectral characteristics.

Zoom In Zoom Out Reset image size

As with any classification scheme, it is important that we can account for missing and noisy data in a robust manner. As LLE is a locally weighted projection, it is expected that outliers (even outliers for a single pixel) can impact the definition of the local neighborhood and might bias any classification scheme.

In practice it is difficult to treat this problem completely as it essentially replaces the least squares analysis of each neighborhood with a total least squares/errors in variables analysis for each of the N points. These cannot be solved in a closed form, and require a computationally expensive iterative process (see van Huffel & Vandewalle 1991 for a detailed exploration of the total least squares problem).

We can, however, follow the methodology of Chang & Yeung (2006) and address the problem of outliers using a robust LLE (RLLE) technique. In RLLE, a reliability score for each data point is determined based on its efficacy in predicting the location of the points in its neighborhood. Outliers will have two general properties: first, they will not be members of large number of local neighborhoods; second, within the neighborhoods in which they appear, they will be far from the best-fit hyperplane. With this in mind, reliability scores are calculated using the following method:

• 1.  
  For each point, an iteratively reweighted least squares analysis (Holland & Welsch 1977) is performed to determine optimal weights for a weighted PCA reconstruction of the point from its neighbors. Schematically, this is akin to finding the best-fit hyperplane, weighting points based on their distance from that hyperplane, and repeating until the process converges. As a result of this process, each point in a local neighborhood is assigned a local weight that quantifies its contribution to the local tangent space.
• 2.  
  Each point may be a member of many neighborhoods. The reliability score for a point is determined by taking the sum of its local weights from each neighborhood of which it is a member.

Thus, a higher score will be given to points that are part of many neighborhoods, and points that contribute more to the reconstruction in each neighborhood. With this in mind, outliers can be identified based on their low reliability score.

In Figure 10, we compare the performance of LLE and RLLE on the S-curve data set from Section 3.2. In this case, we assign 15% of the points within the data set to be randomly selected from the full parameter space. The lower left panel shows the application of the standard LLE approach. As is clear from this figure, outliers can change the local weights of the nearest neighbors and, thereby, distort the resulting lower dimensional manifold; causing the points of different colors to be heavily mixed. Effectively, the random noise within the point distribution fills in enough of the three-dimensional space that LLE cannot extract the lower dimensional manifold. For RLLE, as long as the number of random points is small relative to the number of nearest neighbors, we can isolate and exclude the contaminating points from the local weights. Figure 10 shows that we can recover the underlying two-dimensional manifold in the presence of noise and outliers through this approach.

Zoom In Zoom Out Reset image size

There are three free parameters that must be chosen for LLE: K, the number of nearest neighbors, D_out, the output dimension, and r, the regularization parameter. K and D_out will be discussed here, while r, which depends on details of the implementation, is addressed in Appendix B.

Often in PCA, instead of specifying an output dimension D_out, one chooses to specify a percentage of the sample variance that should be preserved. In LLE, dimensionality can be determined in an analogous way for each local neighborhood (de Ridder & Duin 2002). In the computation of an LLE projection, a variant of the covariance matrix is computed for each local neighborhood. This matrix differs from the covariance matrix used in PCA by only a translation and scaling, neither of which affects the relative magnitudes of the eigenvalues. Taking the lead from PCA, we can specify a percentage of the local variance that we would like to conserve, and thereby find the minimum number of dimensions required locally, simply by doing an eigenanalysis of the local covariance matrix. Taking the mean of the N local dimensionality determinations gives an estimate of the optimal value of D_out for the global projection. For our SDSS analysis, the first three components contain an average of 75% of the local variance. The first ten components contain 92% of the local variance, significantly less than the 99.9% of total variance within the first ten PCA eigenspectra of Yip et al. (2004a). This discrepancy is likely due to noise and intrinsic variation at the local level, which contributes more to the small, local neighborhoods that LLE probes than to the overall variance that PCA probes. Also, our sampling technique (Section 5.3) chooses sources based on dissimilarity, leading to a more intrinsic variation in the subsample.

Choosing the number of nearest neighbors is more difficult. In fact, K need not be constant across all neighborhoods. It could be specified based on local properties of the manifold or the density of sources at a given point. The literature on LLE offers no consensus on how to choose the optimal value of K. Too large a value, and the manifold is oversampled, destroying the assumption of local linearity, and leading to an inherently larger dimensional space. Too low, and the manifold is undersampled, which leads to a loss of information on the relationships between neighborhoods and more sensitivity to the presence of outliers and noise. In this work, we take an empirical approach and try to estimate a series of heuristics that can be applied to the spectral data. We take a range of values for K (10 < K < 30) and determine our ability to maximally distinguish between galaxy populations in the resulting LLE classifications. In the case of galaxy spectra, the K is chosen that maximizes the opening angle between the QSO and emission-line galaxy sequences, measured between any two of the three leading projection dimensions. We find that a value of K = 23 provides the maximal discrimination but that values of K from 20 to 27 are comparable.

It should be noted that the final LLE projection is unique only up to a global translation and rotation. Because of this, the exact nature of the projection is highly dependent on the subset of data chosen, as well as the above-mentioned input parameters. What is invariant, however, is the relationship between each point and its neighbors. In particular, an outlying point in one projection will be an outlier in another projection. This is what makes LLE useful: it allows one to easily visualize the relationships between points in a high-dimensional parameter space.

The LLE algorithm, even when optimized, is computationally expensive. The main bottleneck is the neighbor search, which, for a single point, is $\mathcal {O}[N D_{{\rm in}}]$. Using brute force for all N points is, at worst, $\mathcal {O}[N^2 D_{{\rm in}}]$. While more advanced tree-based algorithms can significantly improve on this, it remains a computational hurdle, especially for high-dimensional data sets (see Berchtold et al. 1998 for a review of fast nearest-neighbor algorithms).

The second bottleneck is the computation of the global projection, which involves determining the bottom D_out + 1 eigenvectors of an [N × N] sparse matrix (see Appendix C). By using iterative techniques such as Arnoldi/Lanczos decomposition, it is possible to determine these eigenvectors without performing a full matrix diagonalization (see Wright & Trefethen 2001 for one well-tested optimized algorithm).

Even fully optimized, it is not feasible to learn the LLE projection from the full SDSS spectral sample, which would amount to hundreds of thousands of points in a 4000-dimensional data space. What is more, a random subsample of these data will be too sparsely populated in some regions of space (e.g., for strong line emitters that account for <0.01% of a randomly sample of spectra from the SDSS), and too densely populated in others (e.g., around L* galaxies). In this case, LLE will not sufficiently probe the entire structure of the manifold. For this reason, a strategy is needed to define a subsample that best spans the entire parameter space, without overpopulating the most probable neighborhoods.

To do this, we follow de Ridder & Duin (2002) and appeal to the fundamental assumption of local linearity on the manifold. For a given point x_i, the first step of LLE is to compute the weights required to construct it from its K neighbors. The procedure outlined in Section 5.2 for determining the intrinsic dimensionality can be used to find the local reconstruction error: i.e., the percentage of total variance not captured in the D_out-dimensional projection. Points for which this local projection error is small add very little information to the overall projection. They can be discarded without much effect. Points for which this reconstruction error is large are not well approximated by a linear combination of neighbors, and, therefore, contain information not present in the local best-fit hyperplane. If they are thrown out, the projected manifold will change significantly.

With this in mind, the following strategy was used to reduce the ∼170,000 unique spectra to a manageable 9000:

• 1.  
  Divide the sample into groups of 2000 spectra each.
• 2.  
  For each group:

  • (a)  
    find the eigenvalues of the local covariance matrix of each point and create a cumulative distribution based on the local reconstruction error; and
  • (b)  
    randomly select 20% of the points based on this cumulative distribution, such that those with the largest reconstruction error are preferentially selected.
• 3.  
  Merge every five subsamples into a new sample of 2000 spectra.
• 4.  
  Repeat this procedure until the total subsample is of the desired size.

Once the projection is determined for a given training sample, it is straightforward to determine the projection of new data onto the training surface. For each new point, the K nearest neighbors are determined within the training data. Weights are determined by minimizing the form of Equation (2). These weights are then used within Equation (4) to determine the projection of the point. Note that this second step amounts to nothing more than a weighted sum of the projected neighbors. This is a key point in the application of LLE: once a projection is defined for a representative training sample, the computational cost of the classification of a set of test-points is largely due to a K-nearest-neighbor search within the training set. This can be done optimally in $\mathcal {O}[D_{{\rm in}}\log N]$.

Reconstructing unknown data from projected coordinates requires the reverse of the above procedure. Weights are determined by minimizing the corollary of Equation (2) within the projected space, and these weights are used to reconstruct the point in the original data space.