Real world data is often high dimensional, yet concentrates near a lower dimensional manifold, embedded in the high dimensional ambient space. In many applications, finding a low-dimensional representation of the data is necessary to efficiently handle it and the representation usually reveals meaningful structures within the data. In recent years, different manifold learning methods have been developed for high dimensional data analysis, which are based on the geometry of the data, i.e. preserving distances within local neighborhoods of the data. These include kernel PCA [1], ISOMAP [2], locally linear embedding (LLE) [3], Laplacian eigenmaps [4], Hessian eigenmaps [5] and Diffusion maps [6]. Some of these methods are spectral methods, based on the eigenvectors of adjacency or affinity matrices of graphs on the data. These methods are capable of capturing a smooth representation of the data and have been shown to be robust to noise and outliers. The ability to capture the underlying low-dimensional structure of data has made these methods appropriate for dimensionality reduction. Unlike classical dimensionality reduction methods, such as principal component analysis (PCA), these methods are nonlinear, which is essential as real-world data typically does not lie on a hyperplane. In addition, they preserve local structures in the data, while disregarding distances between points that are far apart, which are typically meaningless in high-dimensional data. These approaches are very popular in machine learning, signal processing and data mining applications.

When the data set is very large, or when processing online sequential data, it is impractical to directly compute an embedding for the entire dataset. The computational complexity of calculating the affinity matrix and the eigen-decomposition of the matrix become infeasible in terms of memory and running-time. Since these non-linear techniques do not provide an explicit mapping from the data to the embedding, out-of-sample extension (OOSE) methods are used to extend the embedding to new data points [7–14]. In such cases, the low-dimensional embedding is constructed for a representative sample of the data and is then extended to all remaining, or new points. This is a common approach in image processing applications, for example, especially for high-resolution images [15,16].

Deep learning has gained popularity in the past years, achieving state-of-the-art results in machine learning, computer vision and speech processing applications, handling increasingly large datasets [17–19]. Deep neural nets are capable of learning increasingly abstract representations for the data [20], some of which are robust to small perturbations around training points [21,22]. However, these representations are built globally, without incorporating local geometry or density of the data. Jia et al. [23] recently introduced Laplacian Autoencoders, which impose locality preserving constraints via the weighted affinity matrix. Their goal is to improve pre-training of autoencoders in constructing neural networks. However, their formulation lacks both a parameter to balance between the reconstruction loss and the affinity regularizer, and an embedding to compare to, so there is no way to ensure that the final representation is anywhere near the true manifold eigenfunctions. Because of this, there are no theoretical guarantees on the convergence of the algorithm or on the usefulness of the representations.

In this paper, we propose a new approach to out-of-sample extension, applying a deep neural network to learn the mapping between the data and the embedding. We employ deep learning from a manifold learning perspective, by explicitly incorporating a manifold embedding of the data in the deep learning framework. We address three closely connected problems: OOSE of the embedding to new points, a pre-image solution [24, 25] that can include regularization, and outlier detection on test data. This third goal is important and often neglected in OOSE methods. If new data is an outlier and does not fit the model of the training data, the extension is ill-defined and its new representation will be insufficient.

To accomplish these three goals, we propose to train a neural network-based encoder and decoder, and combine them in an autoencoder. First, we train a multi-layer encoder to approximate the low-dimensional embedding on the data; specifically we use diffusion maps [6] for embedding the data. The encoder performs OOSE of the embedding. We then train a multi-layer decoder whose input is the diffusion map, to learn the inverse mapping between the embedding and the data. The decoder enables to recover the pre-image of the embedding, mapping new points in the diffusion space to the high-dimensional data space. Finally, by stacking the two networks, we obtain an autoencoder, termed the diffusion net. The diffusion net enables to perform outlier detection, indicating when the extension of a given point via the encoder is faulty due to its being an outlier that does not follow the model of the data. The diffusion net also performs denoising of the data, reconstructing a clean version of noisy data. Our approach is both computationally efficient and has low memory costs. Once the diffusion net has been trained, it is unnecessary to retain the training data and embeddings as required in other methods. Thus, harnessing the efficiency of deep learning networks enables to efficiently process large quantities of data.

This paper is organized as follows. Section 2 provides background on manifold learning and deep neural networks. In Sections 3 and 4 we propose a deep learning approach for manifold learning, enabling out-of-sample extension, pre-image computation and outlier detection. A proof on bounding the convergence of a single layer encoder is presented in Section 5. We present experimental results in Section 6. Future research directions are discussed in Section 7.

2.1 Diffusion Maps

Diffusion maps is a popular manifold learning technique, based on the construction of the graph Laplacian of the data set [9]. It has been used successfully in various signal processing, image processing and machine learning applications [10, 16, 26–32]. In this section, we briefly review the diffusion maps construction. For a detailed discussion on diffusion maps, see [6].

Given a high dimensional set , a weighted graph is constructed with the data points as nodes and the weight of the edge connecting two nodes, measure of the similarity between the two data points. The affinity matrix required to be symmetric and non-negative, where a common choice is a radial basis function (RBF) kernel

where 0 is a global scale parameter. A local scale can be set for each point as in [33]. In practice, K can be computed using only the nearest neighbors of every point and K[i, j] is set to zero for that are not among the nearest neighbors of

We apply a normalization of the data to obtain an approximation of the Laplace-Beltrami operator on the data, so the embedding will not rely on the distribution of the points [6, 10]. The kernel is normalized by the degree of each point

A random walk is then created on the normalized data set by:

The row-stochastic matrix viewed as the transition matrix of a Markov chain on the data set X. The eigen-decomposition of P yields a sequence of biorthogonal left and right eigenvectors, respectively, and steps of the Markov chain can be calculated as

A diffusion distance ) between two points is defined by

where is the stationary probability distribution on the graph. This metric is robust to noise, since the distance between two points depends on all possible paths of length t between the points. Due to the spectrum decay of P, the diffusion distance can be approximated using only the first d eigenvectors. Equation (3) implies that a mapping can be defined between the original space and the eigenvectors . Retaining only the first d eigenvectors, the mapping Ψembeds the data set X into the Euclidean space , where the diffusion distance is equal to the Euclidean distance in this new embedding:

Note that is not used in the embedding because it is a constant vector. In this paper we set t = 1, but our approach can be used for estimating the embedding for general t.

2.2 Out-of-sample Function Extension

Having calculated a diffusion map Ψ on the data X, various methods have been proposed for extending Ψ to new points. In simple examples, there are analytic mechanisms for creating harmonic extensions when the eigenfunctions of the Laplacian can be derived analytically. However, this is not applicable in the general case. The Nystr¨om extension method is a popular method for general OOSE. Given a new point , the eigenvector is extended to this point as:

Geometric Harmonics [9,10] is an OOSE method which improves upon the Nystr¨om extension method. It treats both the numerical instability due to by extending only the eigenvectors with significant eigenvalues. In addition, it finds an appropriate scale for the kernel in the extension, dependent on the function that is being extended. Rabin and Coifman proposed a Laplacian pyramids-based OOSE method in [11]. The eigenvectors are extended in an iterative multi-scale scheme, where the number of scales is adapted to the complexity of each eigenvector separately. This approach was recently extended in [12] to implicitly incorporate cross validation in the training procedure and avoid over-fitting in the training. Aizenbud, Bermanis and Averbuch [14] introduced an extension method based on a generalized least squares solution for each new test point within its local neighborhood in the training set. This solution is shown to minimize the Mahalanobis distance between the embedding of the training points and the estimated embedding for the new point, in respect to a covariance matrix that incorporates geometric properties of the data and embedding.

The computational complexity of these methods typically depends on the number of training points, since these methods rely on the distances between a new test point and all training points, or its nearest neighbors in the training set (determined by a nearest neighbor search algorithm). Therefore, to perform OOSE, it is necessary to keep all the training points and their corresponding embeddings in memory. If the affinity matrix K in (2) is based on nonEuclidean local metrics, such as [31, 34–37], it is not possible to use nearest neighbors search, since each training point is associated with its own local metric. This increases the complexity of the distance calculation and adds to the memory requirements of the OOSE. The method we propose has no such requirement. After training, our approach enables OOSE without retaining any of the training data or embeddings.

2.3 Deep Neural Networks

Artificial neural networks (ANNs) are networks composed of interconnected computational units termed neurons, which are typically organized in layers. A deep neural network is composed of multiple hidden layers, and is typically designed as a feed-forward network, in which there are no cycles or loops. In our framework, we use multi-layer perceptrons (MLP), which are a popular and important class of neural nets, in which the layers are densely connected. The output of each layer is computed as an affine mapping of the previous layer followed by a non-linear function:

where a is termed the activation, ] denotes the weight associated with the connection between unit is a bias vector, ) is a non-linear function applied element-wise, and . We denote the number of hidden units in layer and the overall number of layers in a network, including the input and output, by L. In our experiments, we used a sigmoid non-linearity in the activation: . Other choices include ) and rectified linear units:

Deep nets have been successfully applied to various tasks such as regression for learning a function over a given dataset, classification, feature learning, etc., achieving state-of-the-art results. The task the network performs is determined by the output layer and the cost function minimized over the network. In supervised learning, the goal is to predict a function or labels on the input data. The cost function of a network consists of a loss function, and a weight regularization term is usually added in order to prevent over-fitting. Given a multi-layer network, we denote the weights and biases of all the layers by . For regression of a multi-dimensional function, , as in our application, the squared error can be used for the loss term:

where is the output of the net for input is the Frobenius norm, parameter. We use an penalty on the weights, which is a very popular choice for regularization. Neural nets are typically trained using variants of stochastic gradient descent (SGD) for calculating the weight and biases in the network that minimize the cost function. The gradient of the loss function relative to the weights is computed efficiently using backpropagation [38], starting from the output layer backwards to the input.

2.3.1 Autoencoders

Deep learning can also be used in an unsupervised manner, such as in training autoencoders. An autoencoder is composed of an encoder function ) and a decoder function ), where the dimension of f is typically smaller that the dimension of the input data. The reconstruction of an input is given by stacking the decoder on the encoder: )), and the autoencoder is trained to minimize a reconstruction error loss L(x, r(x)) over the training points, i.e. trying to approximate the identity function. This is performed by setting the number of output units of the decoder to be the dimension of the input data, and fine-tuning a loss function between the output and the input. The output of the encoder, i.e. the middle layer of the autoencoder, can be seen as a low-dimensional representation of the data. Autoencoders were popularized by Hinton [18]. Basic autoencoders consist of a single hidden layer [21,22], but deep autoencoders have also been proposed for classification, denoising, image retrieval, speech enhancement and more [18,21].

3.1 Problem Setup

We assume the data lies on a smooth, compact, d-dimensional Riemannian manifold M embedded in a high-dimensional space ). Calculating the diffusion map Ψ : for a given training set maps the high-dimensional space to Euclidean space revealing the low-dimensional structure of the data. We address three related problems in this setting: (a) Out-of-sample extension, (b) pre-image solution and (c) outlier detection.

Given new test points , the purpose of an out-of-sample extension method is to extend Ψ to the new points. For the given training points, the extension should be as close as possible to the true embedding: , for small . In addition, the extension ) should preserve the properties of the original embedding, such as preserving local structures in the data.

The second problem we address is calculating the pre-image [24, 25]. Given a point the diffusion space, the pre-image of is a data-point x for which Ψ(. The pre-image has been shown to be closely related to the OOSE problem [25]. Calculating the pre-image is essential to being able to pull back calculations made in the embedding space to the data space.

Finally, we aim to provide a new outlier detection measure, to detects outlier in newly arrived data. In OOSE algorithms, the ability to provide a good embedding for new points depends on how well they follow the model of the training data, i.e. their distance from the training data. For example in kernel-based extension methods such as Nystr¨om, the scale of the Gaussian kernels limits how close points need to be to the training data. Denote by

the nearest-neighbor of a test point within the training set X. If the distance between the two is very large, ), i.e. the OOSE reverts to 1-nearest neighbor prediction. Numerically, if 0, then the affinity to the training set evaluates to zero, and the OOSE is the origin. In related works (e.g. [8–11]), there is typically an implicit assumption that the test data follows the distribution of the training data. However, the embedding obtained for outliers is uninformative and can lead to mistakes in the task the embedding is used for, i.e. classification, signal processing, etc. Therefore, indicating if a point is an outlier within an OOSE framework is important to properly processing the data.

3.2 Our Approach

We propose a solution to these problems based on deep neural nets, training a neural network-based encoder and decoder, and combining them in a deep multi-layer autoencoder. Thus, instead of training a general autoencoder in which the middle layer can be seen as providing a low-dimensional representation of the data, we incorporate the given embedding on the data in the training procedure. The encoder and the decoder are both MLPs composed of L layers, whose output layer is a regression layer. In our approach, the encoder learns a mapping from the manifold to the diffusion space, by minimizing the squared loss between the output of the encoder and Ψ(x). In addition, we impose a new constraint to preserve properties of the embedding. The decoder learns the inverse mapping from the low-dimensional diffusion space back to the high-dimensional space of the data, solving the pre-image problem. By stacking the encoder and decoder we obtain an autoencoder, whose inner-most layer computes the diffusion embedding of the data. The autoencoder can be used for both outlier detection and denoising. Our framework is presented in Algorithms 1, 2 and 3.

3.3 Encoder

The encoder learns the mapping between the high-dimensional data space and the embedding space. The encoder is designed as an MLP, minimizing the loss between the diffusion map Ψ and the output of the net, i.e. ) in (6), where denotes the output matrix of the encoder for all training points, . To facilitate the output of the encoder approximating the diffusion map, we add a new constraint to the training objective of the encoder. Since the

coordinates of the diffusion map are eigenvectors of the random walk matrix on the data,P in (2), we add a new cost term on the output of the encoder being an eigenvector of this matrix. The additional term, termed the eigenvector (EV) constraint, is

where -th row of the output matrix, -th eigenvalue of the affinity matrix, and is an optimization cost parameter.

This new term imposes smooth locality on the training points, and maintains the local geometry of the points in the encoder function. In contrast to general constraints in deep learning, this constraint “mixes” output values of different data points together, via the randomwalk matrix. The EV constraint enforces that an output value for a given data point can be reconstructed as a weighted average of the outputs of the local neighborhood of the data point. The locality of the neighborhood depends on the locality of the kernel used to define the affinity matrix (1). This constraint serves to minimize the loss error on the output, such that the output is not a general regression of the desired function, but an eigenvector of the matrix, thus maintaining the geometric properties of the embedding. Note that for a high value of this constraint forces the output to the origin, which is a trivial solution to minimizing this constraint.

Although we use diffusion maps for manifold learning, other kernel methods may be used and imposed with the EV constraint. This only requires replacing the matrix P with the matrix that is decomposed in other approaches, for example the normalized affinity matrix

as in spectral clustering [39] or the unnormalized graph Laplacian matrix Laplacian Eigenmaps [4]. Adding the EV constraint to the cost function of the encoder, the total cost is:

where denotes the set of the weights and biases of all the layers of the encoder. To incorporate this constraint in the gradient calculation in back-propagation, the gradient of this constraint in regards to the output layer is:

For new data points, the encoder performs out-of-sample extension, based on the mapping learned from the training points. This extension is computationally efficient as it relies on a few affine transforms and non-linear element-wise functions. In addition, once trained, it does not depend at all on the training data, making it efficient also in terms of memory.

3.4 Decoder

The decoder has the same architecture as the encoder, only in reverse. The input layer is . There are no new constraints in the decoder cost, and it consists of a loss term and weight regularization term as in (6):

where where denotes the set of the weights and biases of all the layers of the decoder, and is the output for input Ψ(). Note that the decoder weights are not tied to those of the encoder , therefore the number of units in the hidden layers can be different than for the encoder. Enforcing tied weights in the autoencoder will be explored in future work.

The decoder solves the pre-image problem and has several interesting applications. As the decoder learns the non-linear inverse mapping between the data and the diffusion space, it enables to output points in the data space given new points in the diffusion space. This enables data visualization, performed by covering the diffusion space with new points, and inputting these points to the decoder. This also enables to perform data augmentation for increasing the training set in machine learning applications. Another benefit of the decoder is the ability to perform calculations in the embedding space and then pull back to the data-space, for example, calculating centroids in a clustering application, or interpolation of points in the embedding space. In our current formulation we do not impose constraints on the decoder, however, in future work we will explore including a PDE regularizer in the cost function, such as a harmonic constraint. This will enable recovering interesting surfaces, such as “minimal surfaces” from the embedding.

3.5 Autoencoder

Having trained the encoder and the decoder, the two networks are stacked together to obtain an autoencoder. This architecture is displayed in Figure 1. The network is composed from two stacks, one for the encoder and one for the decoder. In most of our experiments each stack has 2 hidden layers, and in some we use simpler stacks of a single hidden layer.

One application of the autoencoder is to use the autoencoder reconstruction error to detect outliers in the data. Denoting the output of the autoencoder by r, the training data provides an average reconstruction error:

Figure 1: Autoencoder. Left: encoder with two hidden layers, Right: decoder with two hidden layers. The output of the encoder is a prediction of the diffusion map. The output of the decoder is a reconstruction of the data.

This enables setting a threshold on new test points to determine whether they fit the model of the data learned by the autoencoder. If the reconstruction error for a new point than this threshold, the point is considered an outlier and its OOSE can be disregarded. This condition is given by

where C is a constant determined by the user and the training data. Performing outlier detection within the framework of OOSE was previously addressed in [14], where a Mahalanobis distance in the embedding space was used for anomaly detection. Our approach differs in that our measure depends on the reconstruction of the point in the data space and not in the embedding space.

A second application of the autoencoder is denoising. When applied to noisy data, the diffusion map recovers a smooth manifold. Noise in the data, which relates to the relaxation time in the diffusion process, is typically manifested in eigenvectors relating to small eigenvalues [26]. These are are usually disregarded when choosing the number of dimensions to retain in the embedding. Thus, noisy points which relate to the same original clean point are embedded at identical coordinates, and this mapping is learned by the encoder. The decoder, which learns the inverse mapping from the embedding to the data, can only recover a smooth version of the data. Variations in the recovered data are such that they are along the low-dimensional manifold, whereas variations due to noise, that are not seen by the manifold, are suppressed. Thus, inputting new noisy points into the autoencoder, outputs a clean version of the points. Vincent et al. [21] provide an interpretation of their denoising autoencoder based on the manifold assumption, however they do not explicitly incorporate the manifold in their training procedure. In addition, they are intentionally corrupting the data by adding noise to the input and then training the network to recover the clean data. We do not assume the data is clean, but that the embedding provides a smooth representation of the data.

This section provides implementation details regarding the optimization of the diffusion net and the computational complexity of our approach.

4.1 Pre-training

In general neural nets, a greedy layer-wise pre-training procedure was first developed by Hinton, Osindero and Teh [17] for Restricted Boltzman Machines, and later extended to autoencoders [19]. This procedure was shown to improve the optimization of the neural network. In this approach, each layer is first trained in an unsupervised manner as an autoencoder whose input is the output of the previous layer. This enables to “initialize” the parameters the network before performing backpropagation, termed fine-tuning, over the complete network for the supervised learning task. It was shown empirically that this procedure yields improved results compared to initializing with random weights which can get stuck in poor local minimum solutions.

Recently, deep networks trained using very large quantities of labeled data have achieved successful results, with pre-training having very little impact, if at all. It appears that employing large amounts of data in the training cancels the advantage gained by pre-training [20]. In our case, the amount of the data used to train the network is not large, therefore pre-training is beneficial.

We pre-train every hidden layer as a denoising autoencoder [21] with a sparsity term. This term encourages the average activation of every hidden unit to be small so that the hidden representation of the data is sparse [40]. The loss function is given by

where . The input point is corrupted by setting a random subset of entries to zero and ˆis the reconstruction of by the autoencoder. The third term is the sparsity loss where where is the number of hidden units, is a sparsity parameter, typically a small value close to zero, ˆis the average activation of hidden unit is the Kullback-Leibler divergence between the probability mass functions of Bernoulli random variables with parameters . This constraint imposes that each hidden units “responds” only to specific input patterns in the data.

4.2 Optimiztion

The cost function in neural networks is highly non-convex. However, convex optimization methods are used, as typically a local minimum yields good results. We perform training with Limited memory BFGS (L-BFGS) [40, 41] with line search procedure, a quasi-Newton method implemented in Matlab’s minFunc package, and the gradients are computed using standard back-propagation [38]. L-BFGS is a suitable approach when the dataset is not too large, and even out-performs SGD, which is a more popular approach in deep learning, in certain settings [40].

The weights and biases in all layers are initialized with random values from a normal distribution with zero mean and small variance. The weight regularization parameter was set to be small, i.e. order of 10. In the experimental results in Section 6 we examine the performance of the new constraint (7) for different values of

4.3 Mini-batch training

Mini-batch training is typically used with large datasets in deep learning and has been shown to improve performance. In this method, a random subset of samples is chosen for each iteration (or number of iterations) and the gradients are computed and averaged only for this small subset. Since the constraint we added (7) entails multiplying the output of the encoder by the normalized random walk matrix, it is essentially “mixing” between the outputs of different training points. This is problematic in mini-batch training, since a point’s local neighborhood which most influences the eigenvector constraint will not necessarily be included in the mini-batch. If extracting only a subset of points, a new matrix needs to be calculated for the mini-batch by extracting the relevant rows and columns, and renormalizing the new matrix. In this case, this constraint is reduced to the output being only an approximation of the eigenvector of the full random walk matrix. However, since the size of training data used in manifold learning is usually not large due to computational complexity constraints, it is unnecessary to train using mini-batches. In this our setting differs from typical large-scale applications in deep learning such as Imagenet, CIFAR-100, etc.

4.4 Computational complexity

Training is an iterative sequential process, performed offline. Once the encoder has been trained, out-of-sample extension on new data is calculated with complexity and . In our experiments, we used an encoder with two hidden layers so the complexity was of order This is opposed to other OOSE methods in which the complexity depends on m, the number of training points, where typically, example, the complexity of OOSE using Nystr¨om is O((d+n)m).

An advantage of our method is that there is no need to retain the training points and embedding, once the network has been trained. Only the weight matrices and bias vectors of all the layers of the network are necessary. Thus, our approach requires memory on the order of ). Other methods, on the other hand, require retaining all training points and embeddings. For Nystr¨om and Geometric Harmonics [9] this results in memory on order of O((d + n)m). The memory cost of the PCA-based approach in [14] is higher, requiring an additional ) to save covariance matrices for the embeddings of the training points. If using a non-Euclidean metric as in [34, 36], retaining the covariance matrices of the training points leads to an additional memory cost on the order of ). Thus, our method has a large advantage in applications and systems in which the memory and run-time are limited.

In this section, we provide a theoretical bound on the error rate for approximating eigenfunctions of the Laplacian using a single layer network of sigmoid hidden units. The full derivation is given in A. Suppose M is a smooth compact d-dimensional submanifold of is equipped with a metric which is locally bi-Lipschitz with respect to the Euclidean metric, meaning

Theorem 5.1. Let M be a smooth compact submanifold of equipped with a metric which is locally bi-Lipschitz with respect to the Euclidean metric. Let be a Euclidean ball of radius r such that be an eigenfunction of the Laplacian of M with eigenvalue be an extension of the eigenfunction to the ambient dimension via

Then there is a linear combination sigmoidal units such that

Corollary 5.2. Under the same conditions as Theorem 5.1, let be eigenfunctions of the Laplacian and be the extension of . Then there exists a single hidden layer network with Kd sigmoidal units and output

Note that Theorem 5.1 guarantees existence of a solution, yet since the optimization problem is non convex, it does not guarantee convergence.

In this section, we present experimental results for several toy problems and real image data. We demonstrate the performance of the encoder and decoder separately and then joining them in an autoencoder. We evaluate the effect of the eigenvector constraint and demonstrate how it improves the performance of our system, especially in noisy scenarios. Finally, the autoencoder is successfully used in outlier detection in images. This demonstrates that our solution not only performs OOSE, but also means by which to verify that new data agrees with the model inferred from the training data.

6.1 Encoder

Our first example is a closed 3-D curve parametrized by:

We examine the effect of adding the eigenvector constraint (7) to the encoder cost function (8) for this toy problem. We address both the effect of the architecture of the network, i.e. number of layers and units, and the effect of noise in the data. We add white Gaussian noise standard deviation (std.) 05 to the data:

and sample 2000 points from this noisy curve. The 2D diffusion map of these points is a smooth circle (Fig. 2).

Figure 2: (a) 3D closed curve. (b) First two coordinates of the diffusion map. Points colored by

Since there is a great deal of randomness in the optimization, we trained encoders for 10 realizations of the data, and averaged the MSE over the realizations. The encoder MSE for a single realization of the data is calculated as

Figure 3: (a) MSE of the encoder for the training data, comparing 1 hidden layer (red ’x’) vs. 2 hidden layers (blue circle) for varying values. (b) MSE of the encoder for the training data, comparing various

Figure 3(a) plots the MSE for various values of for an encoder with 1 hidden layer and with 2 hidden layers. In Fig. 4, we plot several results from this experiment. The top row displays the output of the encoder with one hidden layer and the bottom row with two hidden layers. Column (a) is the result without the added constraint. Columns (b)-(d) are for respectively.

A first conclusion from these simulations is that adding depth to the network improves the results. This follows known observations in deep learning, that functions that can be compactly represented by networks with L layers, can require an exponential number of units in a network with 1 layers [19]. Using one hidden layer, the reconstruction of the embedding is noisy and the added constraint is not enough. Using two hidden layers, the results are significantly improved and the eigenvector constraint serves to smooth out the noise in the reconstruction. In terms of the MSE, adding a layer decreases the error by a factor of 10. When the cost parameter is too large as in Fig. 4(d), this constraint collapses the output of the encoder to the origin, which is the trivial solution to minimizing this constraint. The optimal results were achieved for 100. Not imposing the EV constraint (twice as high and results in a noisy approximation.

To examine the effect of noise on the approximation, we ran simulations with varying values of . The encoder was trained with 2 hidden layers, with 20 hidden units in each layer. We averaged the training MSE over 20 realizations of the data, where for each realization several encoders were trained, one for each value of in the EV constraint, displayed in Fig. 3(b). The best results are consistently obtained by no noise (is meaningful for estimating the diffusion embedding (values). If is high, the solution collapses to zero, resulting in a high error. For low well so there is some trade-off between the SNR and how high

To examine the dependence of the mean error on the number of units in the single layer encoder, we trained the network with a varying number of hidden units bounded by , so that empirically we are achieving a tighter bound than the bound guaranteed in Corollary 5.2.

Regarding the parameters used in the simulations, the number of units in the hidden layers in the simulations was The cost parameter of the weight regularization term was set as . For the sparsity constraint in pre-training the autoencoder (11), we tried various values of 1 and they do not affect the errors in any

Figure 4: Top row: Encoder with 1 hidden layer, Bottom row: Encoder with 2 hidden layers. Blue plot is original diffusion map, red plot is output of the encoder. The columns examine the effect of the eigenvector constraint for increasing values of is too large, the eigenvector constraint dominates the cost function and the output collapses to the trivial solution that all points equal zero. For a one hidden layer encoder, the eigenvector constraint is not enough to smooth out the noise and obtain a good embedding. For a two layer net, increasing within a reasonable range smooths the noise and yields an improved output compared to not including this constraint.

significant manner. Therefore, we set

6.2 Decoder

The decoder solves the pre-image problem. It can be used to extend data from the diffusion embedding space to the data space, thus creating new points that lie on the manifold. This enables a better visualization of the data belonging to the manifold, in the data space.

Our first example is based on the clean version of the 3D curve, given in (15). We draw 2000 random samples from this curve. The diffusion map is calculated with 1 for the scale in the kernel (1). The first two diffusion coordinates are a circle, as in the previous example. We train a decoder using these 2000 samples and their embedding. We then cover the diffusion space with points enclosed within a circle extended beyond the radius of the embedding and use the decoder to predict the data in the 3D data space. Figure 6 top-left and bottom-left display the points in the embedding space colored by angle and radius, respectively. The remaining three columns display different views of the predicted points in the 3D data space, where the colors correspond to the color of the points in the diffusion space. We can see the points are restricted to a 2D surface enclosed within the 3D curve. Note that the embedding of the origin is handled smoothly, with no singularities in the data space. Also note that points in the diffusion space which are located beyond the original circle, are decoded in the data space along the boundary of the surface. Thus, the range of extension is limited and follows the geometry of the original data.

The next example is of more complex data, whose embedding is also a 2D circle. Given an

Figure 5: Embedding estimation error vs. number of hidden units in 1-layer encoder (solid line). The empirical error is bounded by (dashed line), which is a tighter bound than the theoretical bound guaranteed in Corollary 5.2.

image of a noisy periodic function (std. of the noise is 25), shown in Fig. 7(a), we extract all overlapping patches sized 5 5 pixels, and construct a random walk on the patches, to obtain the diffusion map Ψ. Figure 7(b) displays examples of 49 patches extracted from the image. The first two coordinates of the diffusion embedding are a circle shown in Fig. 7(c). Training the a decoder with two hidden layers for the diffusion map and the patches, yields a function that “produces” image patches. The input dimension is The decoder enables visualization of the data. Figure 8(a) displays the patches obtained by inserting several test points from the diffusion space into the decoder, where the position of the patches in the circle corresponds to the location of the points entered into the decoder. This example shows that the radius in the diffusion embedding represents the amplitude of the periodic function in the image patches. At the origin of the diffusion space, you get a smooth patch. Note the patches are clean as opposed to the training patches.

Figure 8(b) shows how this can be used for image manipulation. A rotation of by 0.45 is applied to the diffusion map as

The values of Ψ are inputted into the decoder. Then the output values are used to reconstruct the image where each value is assigned to its original pixel location. The resulting image is indeed a shift of the original image in Figure 7 and the amplitude has been decreased.

6.3 Autoencoder

We train an autoencoder for the noisy 3D data described in Sec. 6.1. First, we train a decoder from the embedding to the noisy data, using either one or two hidden layers, corresponding to what was used in the encoder. We then stack the decoder on top of the trained encoders, and obtain an autoencoder. We sample 1000 new test points from the noisy curve, and calculate their reconstruction via the autoencoders. These training and test phases were repeated for 10 realizations of the data. For each realization of the data, we train several encoders for different values of and one decoder. Thus, the MSE of the autoencoder depends on and not on the optimization of the decoder, i.e. a lower cost is due to the encoder since all autoencoders share the same decoder. We then average the reconstruction MSE over all realizations, for each separately. The reconstruction error is given by

Figure 6: Left column: the diffusion space is covered in points within a circle. The black points are the diffusion coordinates for the training data and are the input points to the decoder. The remaining three columns are different views of the output of the decoder, where the black points are the original training datapoints. Top row: points colored by angle in the diffusion space. Bottom row: points colored by radius.

Note that we are comparing the output of the autoencoder to the , although the autoencoder is trained using the noisy data . This demonstrates the denoising capabilities of the network. Figure 9(a) summarizes the results for various values of for autoencoder with 1 hidden layer, and with 2 hidden layers. Adding a layer decreases the reconstruction MSE by a factor close to 10. We can see that in all simulations, the training and test reconstruction errors in the data space are of the same order, so that the algorithm does not over-fit the training data. In addition, using 2 layers in the network improves both encoder and autoencoder errors. In addition, as previously shown, if is too large, for example , this corrupts both the embedding and the data. The best overall results were achieved for

To examine the effect of noise on the diffusion net, we trained the autoencoder for varying values of . The network was trained with 2 hidden layers in the encoder, and 2 hidden layers in the decoder, with 20 hidden units in each layer. Following the experimental results of the encoder, we set error over 20 realizations of the data, displayed in Fig. 9(b). The reconstruction error is between the reconstructed output and the original clean data, demonstrating the denoising capabilities of the diffusion net for increasing noise std. As in the previous example, the training and test errors are very similar, implying that the training does not over-fit the data, both in the encoder and in the decoder, even in increasingly noisy scenarios.

6.4 Outlier Detection

We now apply the autoencoder to real image data and demonstrate that the autoencoder performs outlier detection. As stated in Sec. 3, test data will not necessarily follow the model of the data used to calculate the embedding. OOSE applied to new data that significantly differs from the training data will assign embedding values that do not distinguish it from the data. Therefore, it is important to be able to determine when the embedding is unreliable.

Figure 7: (a) Image of noisy periodic function. (b) Example of 49 training patches extracted from the image. (c) Diffusion map of training image.

Figure 8: (a) Decoding image patches. Various points along and within the circle are inputted to the decoder. We display the obtained patches at the locations in the diffusion space that were used an input to the decoder. This demonstrates that the radius of the embedding corresponds to the amplitude of the periodic function in the image space. (b) Image manipulation: the training diffusion map is rotated by 180 degrees and the values multiplied by 0.45. Inputting these points into the decoder and reconstructing the image from output, shows the period of the image has indeed shifted and the amplitude decreased.

Figures 10(a) and 11(a) display two images of patterned semiconductor wafers acquired by a scanning electron microscope, sized 200 200 pixels. Both wafer images have a defect near the center of the image. For each image separately, we randomly extract 2500 patches from the image, sized 8 8 pixels. This training set is used to calculate a diffusion map, reducing the data from dimensionality n = 64 to d = 2 dimensions. The training set is used to train an autoencoder, as in Algorithm 3, and the average training reconstruction error (10) is calculated. We then input all overlapping image patches from the image into the autoencoder and calculate the reconstruction error of each patch. Figures 10(b) and 11(b) display for all pixels in the image, revealing that this approach easily separates the defects from the background.

This is a result of the diffusion map capturing the main structure of the data, i.e. the pattern of the wafer, as represented in the training data. Patches which differ from the training set, as in the case of the defects, are not represented in the diffusion map. Thus, when applying the autoencoder, these patches are not properly reconstructed by the mappings learned from the data space to the diffusion space and back. This result obtained by the autoencoder indicates that, for these patches, the embedding provided by the encoder does not properly represent them.

Figure 9: (a) MSE of data reconstruction via the autoencoder, comparing 1 hidden layer (red ’x’) vs. 2 hidden layers (blue circle) for varying values. Dashed line is training data, solid line is test data. (b) MSE of data reconstruction vs the std. of the noise. Reconstruction is performed via the autoencoder for dashed line is test error.

Figure 10: (a) SEM image of a semiconductor wafer with defect near the middle of the image. (b) Reconstruction error of the image relative to the average training error: reveals the wafer defect.

6.5 OOSE Error vs. Variation in In-sample Embedding

The purpose of an out-of-sample extension algorithm in the manifold learning setting is to provide an extension of the embedding to new points, such that this extension on the new points is close to the embedding of these points, if the embedding was calculated over all the points. Due to the discrete nature of the data, there is no “true” embedding; rather the value of the eigenvectors for a given training point depends on the other points in the set and the scale used in the affinity matrix. If the eigenvalues of the decomposed matrix have a geometric multiplicity, than the eigenvectors belonging to the same eigenvalue span a subspace, and for the same data they are similar up to a rotation. In our experiments in Sec. 6.1, the data is closed curve, so that its continuous equivalent is the heat equation with Neumann boundary condition. Thus continuous solutions to this PDE are the trigonometric functions sin(discrete 3D data, the eigenvectors approximate these continuous periodic eigenfunctions. The first two eigenvectors belong to the same eigenvalue and form a 2D circle in the embedding space.

To demonstrate the performance of our OOSE approach, we conducted an experiment similar to that proposed in [7]. We took m = 2000 training points along the curve and calculated the

Figure 11: (a) SEM image of a semiconductor wafer with defect near the middle of the image. (b) Reconstruction error of the image relative to the average training error: reveals the wafer defect.

diffusion embedding Ψ for these points. We then added points to the training points, and calculated am embedding points. Since the embedding in both cases is a circle, we calculated the rotation between both embeddings using the shared points, i.e. the training points. This was calculated by [32]

The SVD of S is . Then the rotation from Ψ to Ψ is given by

We then calculated the error

We calculated this error over 10 realization of the data and for 4 This error is the variation in the in-sample embedding due to increasing the number of points for which the embedding is calculated. Figure 12 compares this error to the error achieved by our best encoder, trained using 100, the out-of-sample extension error of our encoder is lower than the variation in the in-sample embeddings. Since, OOSE will typically be performed for much higher than 100, this demonstrates that our method provides a good out-of-sample extension.

We have presented a new framework employing deep learning for manifold learning. We proposed designing an encoder and decoder that learn the mapping between a given high-dimensional dataset and low-dimensional embedding, and vice-versa. To this end, we proposed a new constraint in training the encoder, which preserves the locality of the points in the embedding. We demonstrated empirically that this constraint improves the approximation of the embedding. Our encoder enables very efficient out-of-sample extension of the non-linear embeddings to new points, with low memory costs. The decoder provides a solution to the pre-image problem, enabling data visualization and augmentation. Finally, stacking the two networks together as a deep autoencoder enables both denoising and outlier detection of the data, as seen via the embedding. Calculating the reconstruction error of the autoencoder for

Figure 12: Dashed plot is , solid line is DN encoder training error.

new points allows to evaluate whether the OOSE provided by the encoder properly represents these new points. We presented experimental results in noisy scenarios for simulated and real data, demonstrating the properties of the proposed architecture.

Our main focus has been on the encoder for performing out-of-sample extension for data whose distribution follows the distribution of the training data. However, in different applications, such as sequential signal processing, the nature of the data can change over time. In manifold learning, the embedding is usually calculated once for training points, and does not adapt over time for new points, as opposed to online dictionary learning in sparse representations, for example. This could lead to the embedding not providing a good representation of the data as it evolves, and requires re-calculating the embedding again and again. In future work, we propose to develop a method based on online fine-tuning of the autoencoder that will adapt the embedding to new points which do not fit the model of the training data. Instead of performing “regular” fine-tuning of the autoencoder, constraints can be added that will maintain the middle layer as an approximation of the embedding, as we proposed with the encoder in this work. In this case, we will fine-tune with both the test and training data, where the training data regularizes the autoencoder so that its middle layer remains an approximation of the embedding for the training points. By fine-tuning the network so that it reconstructs the new test data, the middle layer should recover a new embedding for the test data. This adaptive approach will be explored in future work.

A second direction is to further explore the decoder and how including different regularizations affects the solution of the pre-image problem. Including a harmonic constraint for example should enable recovering a minimal surface as the example shown in Sec. 6.2. The error rate we provided on the encoder does not apply to the decoder as it requires the function that is being approximated to be band-limited, which does not hold for the decoder in a general case. In future work on the decoder, we intend to provide a theoretical analysis of the decoder, and to expand our theoretical results to multi-layer nets. Computing the pre-image is important in different applications in which interpolating the data by averaging in the high-dimensional data space is meaningless, such as the possibility of performing image texture synthesis. We will analyze datasets in which the high-dimensional data is more complicated and examine how this affects the required complexity of the decoder architecture.

A third research direction is to examine improving deep learning applications. Our network enables to determine the number of nodes needed to learn the geometry of the data and can be used to infer the maximal number of nodes needed to model the complexity of the system for an unconstrained neural net. In addition, we intend to explore whether incorporating our encoder into a deep network will improve deep neural networks. This is motivated by previous works that have shown that implicitly incorporating the manifold assumption in the construction of deep networks improves classification results. Therefore, we expect that explicitly including the embedding in the network via the encoder should be beneficial.

This research was supported by the Israel Science Foundation (grant no. 1130/11). Alexander Cloninger is supported by NSF Award No. DMS-1402254. The authors would like to thank Ronald Coifman, Ronen Talmon and Roy Lederman for helpful discussions and suggestions.

Proof. The proof of the theorem requires two key results in the literature, relying on theorems by Barron [42] and Coifman and Lafon [9].

Theorem A.1. (Theorem 1 from [42]) Let be a function with bounded first moment of its Fourier transform

Let be a Euclidean ball around zero with radius r, which we assume contains our data. Then for every there exists a linear combination sigmoidal units such that

where is the Lebesuge measure on

Now it suffices to show that the extension function from (12) has a bounded first moment of its Fourier transform. To show this, we rely on a second result.

Theorem A.2. (Proposition 11 from [9]) Let M be a submanifold in eigenfunction of its Laplacian with eigenvalue be an approximation level. Let f be an extension function as in (12). Then there exists a band limited function band

We need several more intermediate claims before addressing the result.

Claim A.3. The function

Claim A.4. Let b be as in Theorem A.2, and be its Fourier Transform. Then and

Proof of Claim A.4. Clearly

Because b is band limited with band , meaning supp() is contained inside a ball of radius

where ) is the volume of a ball of radius 1 in

Now we prove Theorem 5.1. By setting in Theorem A.2, we combine the above results to show

where

[1] B. Sch¨olkopf, A. Smola, and K.-R. M¨uller, “Nonlinear component analysis as a kernel eigenvalue problem,” Neural computation, vol. 10, no. 5, pp. 1299–1319, 1998.

[2] J. B. Tenenbaum, V. de Silva, and J. C. Langford, “A global geometric framework for nonlinear dimension- ality reduction,” Science, vol. 290, no. 5500, pp. 2319–2323, Dec. 2000.

[3] S. T. Roweis and L. K. Saul, “Nonlinear dimensionality reduction by locally linear embedding,” Science, vol. 290, pp. 2323–2326, 2000.

[4] M. Belkin and P. Niyogi, “Laplacian eigenmaps for dimensionality reduction and data representation,” Neural Computation, vol. 15, no. 6, pp. 1373–1396, 2003.

[5] D. L. Donoho and C. Grimes, “Hessian eigenmaps: New locally linear embedding techniques for high- dimensional data,” PNAS, vol. 100, pp. 5591–5596, 2003.

[6] R. R. Coifman and S. Lafon, “Diffusion maps,” Appl. Comput. Harmon. Anal., vol. 21, no. 1, pp. 5–30, July 2006.

[7] Y. Bengio, J.-F. Paiement, P. Vincent, O. Delalleau, N. L. Roux, and M. Ouimet, “Out-of-Sample Extensions for LLE, Isomap, MDS, Eigenmaps, and Spectral Clustering,” in Advances in Neural Information Processing Systems 16, S. Thrun, L. Saul, and B. Sch¨olkopf, Eds. MIT Press, 2004, pp. 177–184.

[8] C. Fowlkes, S. Belongie, F. Chung, and J. Malik, “Spectral grouping using the Nystr¨om method,” IEEE Trans. Pattern Analysis and Machine Intelligence, vol. 26, no. 2, pp. 214–225, Jan. 2004.

[9] R. R. Coifman and S. Lafon, “Geometric harmonics: a novel tool for multiscale out-of-sample extension of empirical functions,” Appl. Comput. Harmon. Anal., vol. 21, pp. 31–52, 2006.

[10] S. Lafon, Y. Keller, and R. R. Coifman, “Data fusion and multicue data matching by diffusion maps.” IEEE Trans. Pattern Analysis and Machine Intelligence, vol. 28, no. 11, pp. 1784–97, Nov. 2006.

[11] N. Rabin and R. R. Coifman, “Heterogeneous datasets representation and learning using diffusion maps and Laplacian pyramids,” in Proc. 12th SIAM International Conference on Data Mining, 2012.

[12] A. Fern´andez-Pascual, N. Rabin, and J. R. Dorronsoro, “Auto-adaptative Laplacian pyramids for high-dimensional data analysis,” CoRR, vol. abs/1311.6594, 2013. [Online]. Available: http: //arxiv.org/abs/1311.6594

[13] S. Mousazadeh and I. Cohen, “Out-of-sample extension of band-limited functions on homogeneous manifolds using diffusion maps,” Signal Processing, vol. 108, pp. 521–529, 2015.

[14] Y. Aizenbud, A. Bermanis, and A. Averbuch, “PCA-based out-of-sample extension for dimensionality re- duction,” 2013, submitted for publication.

[15] J. He, L. Zhang, Q. Wang, and Z. Li, “Using diffusion geometric coordinates for hyperspectral imagery representation,” IEEE Geosci. Remote Sens. Letters, vol. 6, no. 4, pp. 767–771, Oct. 2009.

[16] Z. Farbman, R. Fattal, and D. Lischinski, “Diffusion maps for edge-aware image editing,” ACM Trans. Graph., vol. 29, no. 6, pp. 145:1–145:10, Dec. 2010.

[17] G. Hinton, S. Osindero, and Y. Teh, “A fast learning algorithm for deep belief nets,” Neural Computation, vol. 18, no. 7, pp. 1527–1554, July 2006.

[18] G. E. Hinton and R. R. Salakhutdinov, “Reducing the dimensionality of data with neural networks,” Science, vol. 313, no. 5786, pp. 504–507, 2006.

[19] Y. Bengio, “Learning Deep Architectures for AI,” Found. Trends Mach. Learn., vol. 2, no. 1, pp. 1–127, Jan. 2009. [Online]. Available: http://dx.doi.org/10.1561/2200000006

[20] Y. Bengio, A. Courville, and P. Vincent, “Representation learning: A review and new perspectives,” Pattern Analysis and Machine Intelligence, IEEE Transactions on, vol. 35, no. 8, pp. 1798–1828, Aug 2013.

[21] P. Vincent, H. Larochelle, Y. Bengio, and P.-A. Manzagol, “Extracting and composing robust features with denoising autoencoders,” in Proceedings of the 25th International Conference on Machine learning (ICML-08). ACM, 2008, pp. 1096–1103.

[22] S. Rifai, P. Vincent, X. Muller, X. Glorot, and Y. Bengio, “Contractive auto-encoders: Explicit invariance during feature extraction,” in Proceedings of the 28th International Conference on Machine Learning (ICML-11), 2011, pp. 833–840.

[23] K. Jia, L. Sun, S. Gao, Z. Song, and B. E. Shi, “Laplacian Auto-Encoders: An explicit learning of nonlinear data manifold,” Neurocomputing, vol. 160, pp. 250–260, 2015.

[24] J. T.-Y. Kwok and I. W.-H. Tsang, “The pre-image problem in kernel methods,” Neural Networks, IEEE Transactions on, vol. 15, no. 6, pp. 1517–1525, 2004.

[25] P. Arias, G. Randall, and G. Sapiro, “Connecting the out-of-sample and pre-image problems in kernel methods,” in IEEE, 2007, pp. 1–8.

[26] A. Singer, Y. Shkolnisky, and B. Nadler, “Diffusion interpretation of nonlocal neighborhood filters for signal denoising,” SIAM Journal Imaging Sciences, vol. 2, no. 1, pp. 118–139, Jan. 2009.

[27] R. Talmon, I. Cohen, and S. Gannot, “Single-channel transient interference suppression with diffusion maps,” IEEE Trans. Audio, Speech Lang. Process., vol. 21, no. 1, pp. 130–142, Apr. 2012.

[28] G. David and A. Averbuch, “Hierarchical data organization, clustering and denoising via localized diffusion folders,” Applied and Computational Harmonic Analysis, vol. 33, no. 1, pp. 1 – 23, 2012.

[29] S. Gepshtein and Y. Keller, “Image completion by diffusion maps and spectral relaxation,” IEEE Trans. Image Process., vol. 22, no. 8, pp. 2839–2846, 2013.

[30] G. Mishne and I. Cohen, “Multiscale anomaly detection using diffusion maps,” IEEE J. Sel. Topics Signal Process., vol. 7, pp. 111 – 123, Feb. 2013.

[31] A. Haddad, D. Kushnir, and R. R. Coifman, “Texture separation via a reference set,” Appl. Comput. Harmon. Anal., vol. 36, no. 2, pp. 335 – 347, Mar. 2014.

[32] R. R. Coifman and M. J. Hirn, “Diffusion maps for changing data,” Appl. Comput. Harmon. Anal., vol. 36, no. 1, pp. 79 – 107, 2014.

[33] L. Zelnik-Manor and P. Perona, “Self-tuning spectral clustering,” in NIPS 17, 2005, pp. 1601–1608.

[34] A. Singer and R. R. Coifman, “Non-linear independent component analysis with diffusion maps,” Appl. Comput. Harmon. Anal., vol. 25, no. 2, pp. 226 – 239, 2008.

[35] R. Talmon, I. Cohen, S. Gannot, and R. R. Coifman, “Supervised Graph-based Processing for Sequential Transient Interference Suppression,” IEEE Trans. Audio, Speech Lang. Process., vol. 20, no. 9, pp. 2528– 2538, 2012.

[36] R. Talmon and R. R. Coifman, “Empirical intrinsic geometry for nonlinear modeling and time series filter- ing,” Proceedings of the National Academy of Sciences, 2013.

[37] G. Mishne, R. Talmon, and I. Cohen, “Graph-based supervised automatic target detection,” IEEE Trans. Geosci. Remote Sens., vol. 53, no. 5, pp. 2738–2754, May 2015.

[38] R. Rojas, Neural networks: a systematic introduction. Springer Science & Business Media, 1996.

[39] Y. Weiss, “Segmentation using eigenvectors: a unifying view,” in Computer vision, 1999. The proceedings of the seventh IEEE international conference on, vol. 2. IEEE, 1999, pp. 975–982.

[40] Q. V. Le, J. Ngiam, A. Coates, A. Lahiri, B. Prochnow, , and A. Y. Ng, “On optimization methods for deep learning,” in Proceedings of the 28th International Conference on Machine Learning (ICML-11), 2011, pp. 265–272.

[41] S. J. Wright and J. Nocedal, Numerical optimization. Springer New York, 1999, vol. 2.

[42] A. Barron, “Universal approximation bounds for superpositions of a sigmoidal function,” Information Theory, IEEE Transactions on, vol. 39, no. 3, pp. 930–945, May 1993.