Deep nets for local manifold learning

The problem of extending a function $f$ defined on a training data $\mathcal{C}$ on an unknown manifold $\mathbb{X}$ to the entire manifold and a tubular neighborhood of this manifold is considered in this paper. For $\mathbb{X}$ embedded in a high dimensional ambient Euclidean space $\mathbb{R}^D$, a deep learning algorithm is developed for finding a local coordinate system for the manifold {\bf without eigen--decomposition}, which reduces the problem to the classical problem of function approximation on a low dimensional cube. Deep nets (or multilayered neural networks) are proposed to accomplish this approximation scheme by using the training data. Our methods do not involve such optimization techniques as back--propagation, while assuring optimal (a priori) error bounds on the output in terms of the number of derivatives of the target function. In addition, these methods are universal, in that they do not require a prior knowledge of the smoothness of the target function, but adjust the accuracy of approximation locally and automatically, depending only upon the local smoothness of the target function. Our ideas are easily extended to solve both the pre--image problem and the out--of--sample extension problem, with a priori bounds on the growth of the function thus extended.


Introduction
Machine learning is an active sub-field of Computer Science on algorithmic development for learning and making predictions based on some given data, with a long list of applications that range from computational finance and advertisement, to information retrieval, to computer vision, to speech and handwriting recognition, and to structural healthcare and medical diagnosis.In terms of function approximation, the data for learning and prediction can be formulated as {(x, f x )}, obtained with an unknown probability distribution.Examples include: the Boston housing problem (of predicting the median price f x of a home based on some vector x of 13 other attributes [60]) and the floor market problem [59,3] (that deals with the indices of the wheat floor pricing in three major markets in the United States).For such problems, the objective is to predict the index f x in the next month, say, based on a vector x of their values over the past few months.Other similar problems include the prediction of blood glucose level f x of a patient based on a vector x of the previous few observed levels [54,51], and the prediction of box office receipts (f x ) on the date of release of a movie in preparation, based on a vector x of the survey results about the movie [57].It is pointed out in [40,38,19] that all the pattern classification problems can also be viewed fruitfully as problems of function approximation.While it is an ongoing research to allow non-numeric input x (e.g., [10]), we restrict our attention in this paper to the consideration of x ∈ R D , for some integer D ≥ 1.
In the following discussion, the first component x is considered as input, while the second component f x is considered the output of the underlying process.The central problem is to estimate the conditional expectation of f x given x.Various statistical techniques and theoretical advances in this direction are well-known (see, for example [61]).In the context of neural and radial-basis-function networks, an explicit formulation of the input/output machines was pointed out in [25,24].More recently, the nature of deep learning as an input/output process is formulated in the same way, as explained in [35,55].To complement the statistical perspective and understand the theoretical capabilities of these processes, it is customary to think of the expected value of f x , given x , as a function f of x.The question of empirical estimation in this context is to carry out the approximation of f given samples {(x, f (x))} x∈C , where C is a finite training data set.In practice, because of the random nature of the data, it may be possible that there are several pairs of the form (x, f x ) in the data for the same values of x.In this case, a statistical scheme, such as some kind of averaging of f x being the simplest one, can be used to obtain a desired value f (x) for the sample of f at x, x ∈ C. From this perspective, the problem of extending f from the traning data set C to x not in C in machine learning is called the generalization problem.
We will illustrate this general line of ideas by using neural networks as an example.To motivate this idea, let us first recall a theorem originating with Kolmogorov and Arnold [37,Chapter 17,Theorem 1.1].According to this theorem, there exist universal Lipschitz continuous functions φ 1 , • • • , φ 2D+1 and universal scalars λ 1 , • • • , λ D ∈ (0, 1), for which every continuous function f : [0, 1] D → R can be written as where g is a continuous function that depends on f .In other words, for a given f , only one function g has to be determind to give the representation formula (1.1) of f .A neural network, used as an input/output machine, consists of an input layer, one or more hidden layers, and an output layer.Each hidden layer consists of a number of neurons arranged according to the network architecture.Each of these neurons has a local memory and performs a simple non-linear computation upon its input.The input layer fans out the input x ∈ R D to the neurons at the first hidden layer.The output layer typically takes a linear combination of the outputs of the neurons at the last hidden layer.The right hand side of (1.1) is a neural network with two hidden layers.The first contains D neurons, where the j-th neuron computes the sum The next hidden layer contains 2D + 1 neurons each evaluating the function g on the output of the j-th neuron in the first hidden layer.The output layers takes the sum of the results as indicated in (1.1).
From a practical point of view, such a network is clearly hard to construct, since only the existence of the functions φ j and g is known, without a numerical procedure for computing these.In the early mathematical development of neural networks during the late 1980s and early 1990s, instead of finding these functions for the representation of a given continuous function f in (1.1), the interest was to study the existence and characterization of universal functions σ : R → R, called activation functions of the neural networks, such that each neuron evaluates the activation function upon an affine transform of its input, and the network is capable of approximating any desired real-valued continuous target function f : K → R arbitrarily closely on K, where K ⊂ R D is any compact set.
For example, a neural network with one hidden layer can be expressed as a function Here, the hidden layer consists of n neurons, each of which has a local memory.The local memory of the k-th neuron contains the weights w k ∈ R D , and the threshold b k ∈ R. Upon receiving the input x ∈ R D from the input later, the k-th neuron evaluates σ(w k • x + b k ) as its output, where σ is a non-linear activation function.The output layer is just one circuit where the coefficients {a k } are stored in a local memory, and the evaluates the linear combination as indicated in (1.2).Training of this network in order to learn a function f on a compact subset K ⊂ R D to an accuracy of ǫ > 0 involves finding the parameters The most popular technique for doing this is the so called back-propagation, which seeks to find these quantities by minimizing an error functional usually with some regularization parameters.We remark that the number n of neurons in the approximant (1.2) must increase, if the tolerance ǫ > 0 in the approximation of the target function f is required to be smaller.From a theoretical perspective, the main attraction of neural networks with one hidden layer is their universal approximation property as formulated in (1.3), which overshadows the properties of their predecessors, namely: the perceptrons [52].In particular, the question of finding sufficient conditions on the actvation function σ that ensure the universal approximation property was investigated in great detail by many authors, with emphasis on the most popular sigmoidal function, defined by the property σ(t) → 1 for t → ∞ and σ(t) → 0 for t → −∞.For example, Funahashi [22] applied some discretization of an integral formula from [28] to prove the universal approximation property for some sigmoidal function σ.A similar theorem was proved by Hornik, Stinchcombe, White [27] by using the Stone-Weierstrass theorem, and another by Cybenko [17] by applying the Hahn-Banach and Riesz Representation theorems.A constructive proof via approximation by ridge functions was given in our paper [11], with algorithm for implementation presented in our follow-up work [12].A complete characterization of which activation functions are allowed to achieve the universal approximation property was given later in [49,36].
However, for neural networks with one hidden layer, one of the severe limitations to applying training algorithms based on optimization, such as back-propagation or those proposed in the book [61] of Vapnik, is that it is neccessary to know the number of neurons in N in advance.Therefore, one major problem in the 1990s, known as the complexity problem, was to estimate the number of neurons required to approximate a function to a desired accuracy.In practice, this gives rise to a trade-off: to achieve a good approximation, one needs to have a large number of neurons, which makes the implementation of the training algorithm harder.
In this regard, nearly a century of research in approximation theory suggests that the higher the order of smoothness of the target function, the smaller the number of neurons should be, needed to achieve the desired accuracy.There are many different definitions of smoothness that give rise to different estimates.For example, under the condition that the Fourier transform of the target function f satisfies Barron [1] proved the existence of a neural network with O(ǫ −2 ) neurons that gives an L 2 ([0, 1] D ) error of O(ǫ).While it is interesting to note that this number of neurons is essentially independent of the dimension D, the constants involved in the O term as well as the number of derivatives needed to ensure the condition on the target function may increase with D. Several authors have subsequently improved upon such results under various conditions on the activation function as well as the target function so as to ensure that the constants depend polynomially on D (e.g., [31,32,44] and references therein).
The most commonly understood definition of smoothness is just the number of derivatives of the target function.It is well-known from the theory of n-widths that if r ≥ 1 is an integer, and the only a priori information assumed on the unknown target function is that it is r-times continuously differentiable function, then a stable and uniform approximation to within ǫ by neural networks must have at least a constant multiple of ǫ −D/r neurons.In [42], we gave an explicit construction for a neural network that achieves the accuracy of ǫ using O(ǫ −D/r ) neurons arranged in a single hidden layer.It follows that this suffers from a curse of dimensionality, in that the number of neurons increases exponentially with the input dimension D. Clearly, if the smoothness r of the function increases linearly with D, as it has to in order to satisfy the condition in [1], then this bound is also "dimension independent".While this is definitely unavoidable for neural networks with one hidden layer, the most natural way out is to achieve local approximation; i.e., given an input x, construct a network with a uniformly bounded number of neurons that approximates the target function with the optimal rate of approximation near the point x, preferably using the values of the function also in a neighborhood of x.Unfortunately, this can never be achieved as we proved in [13].Furthermore, we have proved in [14] that even if we allow each neuron to evaluate its own activation function, this local approximation fails.Therefore the only way out is to use a neural network with more than one hidden layer, called deep net (for deep neural network).Indeed, local approximation can be achieved by a deep net as proved in our papers [40,41].In this regard, it is of interest to point out that an adaptive version of [40,41] was derived in [48] for prediction of time series, yielding as much as 150% improvement upon the state-of-the-art at that time, in the study of the floor market problem.
Of course, the curse of dimensionality is inherent to the problem itself, whether with one or more hidden layers.Thus, while it is possible to construct a deep net to approximate a function at each point arbitrarily closely by using a uniformly bounded number of neurons, the uniform approximation on an entire compact set, such as a cube, would still require an approximation at a number of points in the cube, and this number increasing exponentially with the input dimension.Equivalently, the effective number of neurons for approximation on the entire cube is still exponentially increasing with the input dimension.
In addition to the high dimensionality, another difficulty in solving the function approximation problem is that the data may be not just high dimensional but unstructured and sparse.A relatively recent idea which has been found very useful in applications, in fact, too many to list exhaustively, is to consider the points x as being sampled from an unknown, low dimensional sub-manifold X of the ambient high dimensional space R D .The understanding of the geometry of X is the subject of the bulk of modern research in the area of diffusion geometry.An introduction to this subject can be found in the special issue [9] of Applied and Computational Harmonic Analysis.The basic idea is to construct the so-called diffusion matrix from the data, and use its eigen-decomposition for finding local coordinate charts and other useful aspects of the manifold.The convergence of the eigen-decomposition of the matrices to that of the Laplace-Beltrami and other differential operators on the manifold is discussed, for example, in [2,33,58].It is shown in [29,30] that some of the eigenfunctions on the manifolds yield a local coordinate chart on the manifold.In the context of deep learning, this idea is explored as a function approximation problem in [53], where a deep net is developed in order to learn the coordinate system given by the eigenfunctions.
On the other hand, while much of the research in this direction is focused on understanding the data geometry, the theoretical foundations for the problems of function approximation and harmonic analysis on such data-defined manifold are developed extensively in [38,20,21,46,47,15].The theory is developed more recently for kernel construction on directed graphs and analysis of functions on changing data in our paper [39].However, a drawback of the approach based on data-defined manifolds, known as the out-of-sample extension problem, is that since the diffusion matrix is constructed entirely using the available data, the whole process must be done again if new data become available.A popular idea is then to extend the eigen-functions to the ambient space by using the so called Nyström extension [16].
The objective of this present paper is to describe a deep learning approach to the problem of function approximation, using three groups of networks in the deep net.The lowest layer accomplishes dimensionality reduction by learning the local coordinate charts on the unknown manifold without using any eigen-decomposition.Having found the local coordinate system, the problem is reduced to the classical problem of approximating a function on a cube in a relatively low dimensional Euclidean space.For the next two layers, we may now apply the powerful techniques from approximation theory to approximate the target function f , given the samples on the training data set C. We describe two approaches to construct the basis functions using multi-layered neural networks, and to construct other networks to use these basis functions in the next layer to accomplish the desired function approximation.
We summarize some of the highlights of our paper.
• We give a very simple learning method for learning the local coordinate chart near each point.The subsequent approximation process is then entirely local to each coordinate patch.
• Our method allows us to solve the pre-image problem easily; i.e., to generate a point on the manifold corresponding to a given local coordinate description.
• The learning method itself does not involve any optimization based technique, except probably for reducing the noise in the values of the function itself.
• We provide optimal error bounds on approximation based on the smoothness of the function, while the method itself does not require an a priori knowledge of such smoothness.
• Our methods can solve easily the out-of-sample extension problem.Unlike the Nyström extension process, our method does not require any elaborate construction of kernels defined on the ambient space and commuting with certain differential operators on the unknown manifold.
• Our method is designed to control the growth of the out-of-sample extension in a tubular neighborhood of X, and is local to each coordinate patch.
This paper is organized as follows.In Section 2, we describe the main ideas in our approach.The local coordinate system is described in detail in Section 2.2.Having thus found a local coordinate chart around the input, the problem of function approximation reduces to the classical one.In Section 2.3, we demonstrate how the popular basis functions used in this theory can be implemented using neural networks with one or more hidden layers.The function approximation methods which work with unstructured data without using optimization are described in Section 2.6.In Section 3, we explain how our method can be used to solve both the pre-image problem and the out-of-sample extension problem.

Main ideas and results
The purpose of this paper is to develop a deep learning algorithm to learn a function f : X → R, where X is a d dimensional compact Riemannian sub-manifold of a Euclidean space R D , with d ≪ D, given training data of the form {(x j , f (x j ))} M j=1 , x j ∈ X.It is important to note that X itself is not known in advance; the points x j are known only as D-dimensional vectors, presumed to lie on X.In Sub-section 2.1, we explain our main idea briefly.In Sub-section 2.2, we derive a simple construction of the local coordinate chart for X.In Sub-section 2.3, we describe the construction of a neural network with one or more hidden layers to implement two of the basis functions used commonly in function approximation.While the well known classical approximation algorithms require a specific placement of the training data, one has no control on the location of the data in the current problem.In Section 2.6, we give algorithms suitable for the purpose of solving this problem.

Outline of the main idea
Our approach is the following.
1. X is a finite union of local coordinate neighborhoods, and x belongs to one of them, say U. We find a local coordinate system for this neighborhood in terms of Euclidean distances on R D , say Φ : U → [−1, 1] d , where d is the dimension of the manifold.Let y = Φ(x), and with a relabeling for notational convenience, {x j } K j=1 be the points in U, y j = Φ(x j ).This way, we have reduced the problem to approximating g = f •Φ : [−1, 1] d → R at y, given the values {(y j , g(y j ))} K j=1 , where g(y j ) = f (x j ).We note that {y j } is a subset of the unit cube of low dimensional Euclidean space, representing a local coordinate patch on X.Thus, the problem of approximation of f on this patch is reduced that of approximation of g, a well studied classical approximation problem.
2. We will summarize the solution to this problem using neural networks with one or more hidden layers, e.g., an implementation of multivariate tensor product spline approximation using multi-layerd neural network.
Thus, the layers of our deep learning networks will have three main layers.
1.The bottom layer receives the input x, figures out which of the points x j are in the coordinate neighborhood of x, and computes the local coordinates y, y j .
2. The next several layers compute the local basis functions necessary for the approximation, for example, the B-splines and their translates using the multi-layered neural network as in [40].
3. The last layer receives the data {(y j , g(y j ))} K j=1 , and computes the approximation described in Step 2 above.

Local coordinate learning
We assume that 1 ≤ d ≤ D are integers, X is a d dimensional smooth, compact, connected, Riemannian submanifold of a Euclidean space R D , with geodesic distance ρ.
Before we discuss our own construction of a local coordinate chart on X, we wish to motivate the work by describing a result from [30].Let {λ 2 k } ∞ k=0 be the sequence of eigenvalues of the (negative of the) Laplace-Beltrami operator on X, and for each k ≥ 0, φ k be the eigenfunction corresponding to the eigenvalue λ 2 k .We define a formal "heat kernel" by The following result is a paraphrasing of the heat triangulation theorem proved in [30, Theorem 2.2.7] under weaker assumptions on X.
Let B ⊂ X be the geodesic ball of radius c 4 R, centered at x * 0 , and Since the paper [30] deals with a very general manifold, the mapping Φ jms is not claimed to be a diffeomorphism, although it is obviously one-one on B.
We note that even in the simple case of a Euclidean sphere, an explicit expression for the heat kernel is not known.In practice, the heat kernel has to be approximated using appropriate Gaussian networks [33].In this section, we aim to obtain a local coordinate chart that is computed directly in terms of Euclidean distances on R D , and depends upon d + 2 trainable parameters.The construction of this chart constitutes the first hidden layer of our deep learning process.As explained in the introduction, once this chart is in place, the question of function extension on the manifold reduces locally to the well studied problem of function extension on a d dimensional unit cube.
To describe our constructions,we first develop some notation.
In this section, it is convenient to use the notation , which we will use in the rest of the sections.If 1 ≤ d ≤ D is an integer, and x ∈ R d , x d denotes the Euclidean norm of x.If x ∈ R D , we will write π There exists δ * > 0 with the following properties.The manifold is covered by finitely many geodesic balls such that for the center x * 0 ∈ X of any of these balls, there exists a diffeomorphism, namely, the exponential coordinate map u = (u 1 , • • • , u D ) from B d (0, δ * ) to the geodesic ball around x * 0 = u(0) [18, p. 65].If J is the Jacobian matrix for u, given by J i,j (y) = D i u j (y), y ∈ B d (0, δ * ), then (2.4) Further, there exists κ > 0 (independent of x * ) such that In turn, this leads to , be chosen with the following properties: and, with the matrix function U defined by we have J(0)U (0)y d ≥ γ > 0, y d = 1. (2.10) Any set {x * ℓ } with these properties will be called coordinate stars around x * .We note that the matrix J(0)U (0) has columns given by π c (x * − x * j ), j = 1, • • • , d, and hence, can be computed without reference to the map u.Let Then U is a neighborhood of x * 0 and Φ maps U diffeomorphically onto [−1, 1] d .We oberve that X is a union of finitely many neighborhoods of the form U, so that any x ∈ X belongs to at least one such neighborhood.Moreover, Φ(x) can be computed entirely in terms of the description of x in terms of its D-dimensional coordinates.Remark 2.3 Since the mapping Φ in Remark 2.1 is a quadratic polynomial in x, it can be implemented as a neural network with a single hidden layer using the activation function given in (2.22) as described in Sub-Section 2.4.
Example 2.1 Let 0 < a < 1, and M = M a be the circular helix defined by u(s) = (cos as, sin as, 1 − a 2 s) T .
Clearly, M is a one dimensional manifold, and s is the arclength parameter, measured from (1, 0, 0) T .The curvature at any point is a 2 .For any point z 0 ∈ M , U z0 = M , with the diffeomorphism given by u(tan −1 (πt/2)), t ∈ (−1, 1).An interesting fact is that u(s + 2π) − u(s) 3 = 2π √ 1 − a 2 can be made arbitrarily small by choosing a close to 1, even though the geodesic distance between u(s + 2π) and u(s) is 2π.Let s 0 ∈ R, s 0 + π/4 ≤ s 1 ≤ s 0 + 3π/8, and Therefore, using the well known estimates we obtain that and Hence, for any points u(t 1 ), u(t 2 ) ∈ U , we have We note that the neighborhood U where this estimate holds and the constants are independent of the curvature.
The remainder of this section is devoted to the proof of Theorem 2.2.

Local basis functions
Having found a local coordinate map Φ on a neighborhood U of x on X, the problem of extending f from There is, of course, 100+ years of research on this subject.We restrict ourselves to two examples, which can be implemented using neural networks with one or more hidden layers.One of the most popular activation function in the deep learning literature (e.g., [35]) is the rectified linear unit function Since this function is not continuously differentiable, there are some technical difficulties to use common algorithms like back-propagation with these activation functions.Although we do not need back-propagation in our theory, we prefer to deal with a rectified quadratic unit function defined for t ∈ R by which is continuously differentiable on R. Our theory will work in general with any activation function of order k ≥ 2; i.e., with a function σ that satisfies but for the sake of clarity of exposition, we will use only the activation function σ defined in (2.22).

Polynomials
The most basic class of classical approximants is the set of all polynomials.For n > 0, we denote the class of all algebraic polynomials of coordinatewise degree at most n in d variables by Π d n .(It is convenient to use the same notation also when n is not an integer; in this case, Π d n is just Π d ⌊n⌋ .The basic implementation of polynomials is given in [11, Proof of Theorem 3.1], where an explicit construction is given for finding the weights {w k }, the thresholds {b k } and the coefficients {a k } used in (2.24) below.
Here, the weights {w k } and the thresholds {b k } are independent of P and the coefficients {a k } are linear functionals on Π d N .
We observe that so that the expression on the right hand side of (2.24) can be expressed as a neural network with log 2 N hidden layers.
We note that a neural network with one hidden layer is given in [42], but using a C ∞ activation function σ; e.g., σ(t) = (1 + e −t ) −1 .This uses the fact that for w, x ∈ R d and b ∈ R, such that none of the derivatives σ (j) (b), A finite difference scheme to implement this differentiation yields a neural network with one hidden layer, containing exactly dim(Π d n ) neurons, and should be stable for C ∞ functions.If stability is a greater concern, then one may use other numerical differentiation schemes to implement this formula, e.g., spectral methods [26].

B-splines
For t ∈ R, and integer m ≥ 1, let A tensor product cardinal B-spline at y ∈ [−1, 1] d is defined by It is explained in [40,41] that the quantity N m (y) can be computed using a neural network with a sigmoidal function of order m − 1 consisting of finitely many neurons arranged in multiple hidden layers (the number of neurons and layers depending on m and d alone).Thus, if m is a power of 2, then each of the terms (y j − k j ) m + can be implemented as an iterated power of (y j − k j ) 2 + (cf.(2.26).)The product of d such expressions can be implemented using either Theorem 2.3 as a network with mulitple hidden layer and utilizing the rectified quadratic unit function as the activation function, or a discretization of the formula (2.27) using a C ∞ sigmoidal function as explained in Sub-section 2.4.

Function approximation
In this section, if y ∈ R d , y ∞ is the ℓ ∞ norm of y.

Spline based approximation
In [7, Section 4.5], [8], a quasi-interpolatory spline function is defined by where λ * m are compactly supported linear functionals, designed specifically to ensure that Q m (P ) = P for every polynomial P of coordinatewise degree at most m − 1 in d variables.With Q m,h (f )(y) = Q m (f (h(•)))(y/h), h > 0, one has the approximation bound for small h: (2.30) The linear functionals λ * m are based on finitely many samples of f at the grid points in a compact subset of Z d .In our context, the data for approximating f is not in this form.Therefore, we may use the following algorithm given in [50], where we assume that λ * m is scaled so as to be supported on [−1, 1] d .Given: and δ(C) be sufficiently small.
2. Choose C 0 ⊆ C, so that each subcube has exactly one point of C 0 .
Substituting γ in place of λ * m in the definition of Q m (f ) yields the desired spline approximation and it is proved in [50] that the estimate (2.30) holds with Qm replacing Q m .

Polynomial quasi-interpolation
A standard method for polynomial approximation is to consider a filtered projection defined in (2.38) below.
The Chebyshev polynomials (of first kind) are defined recursively for t ∈ R and integer m ≥ 0 by In terms of monomials, the Chebyshev polynomials are given by T m j (y j ). (2.36) We choose a smooth low pass filter h; i.e., an even function h : R → [0, 1] such that h(u) = 1 if |u| ≤ 1/2, and h(u) = 0 if |u| ≥ 1, and abuse the notation as usual to define With this filter, we define the kernel (2.37) Then the filtered projection operator is defined by It is well known that if f is any continuous function on [−1, 1] d , then {V n (f )} converges uniformly to f at the near optimal rate of approximation.For example, if f has partial derivatives up to order r in each variable, then analogously to (2.30), but for large n rather than small h, Theoretically, the question then is to compute V n (f ) using the data C as in Sub-section 2.6.1.The procedure we describe below from [43,45,51] also describes the choice of the parameter n depending upon the data.
and δ • (C) be sufficiently small.We also fix an integer n > 0.
Objective: To find real numbers {w ξ } ξ∈C such that the functional ) , for all d-variate polynomials P of coordinatewise degree ≤ n − 1. Steps: 2. Choose C 0 ⊆ C, so that each subcube has exactly one point of C 0 .
3. Solve the following (underdetermined) system of equations for the unknowns w ξ , ξ ∈ C 0 . ξ∈C It is proved in the papers cited above that with the discretized operator one obtains the near best rates of approximation.In particular, if f has partial derivatives up to order r in each variable, then (2.39) holds with V n,C (f ) replacing V n (f ).In practice, we choose n to be the largest integer such that either the condition number of the system of equations in (2.41) is "reasonable" or else by checking the resulting errors in (2.41) [34].The formula (2.42) can be re-written in the form where

Extensions
Since the starting point of diffusion geometry is to consider eigen-decomposition of a diffusion matrix, which is constructed using the available data, the entire computation needs to be redone if a new data becomes available.Since the manifold X is only an abstract model, it is not even clear that the new data will belong to this manifold.This gives rise to two related questions.One is to find new points on the manifold; the so called pre-image problem [53] and the other is the out-of-sample extension problem; i.e., extend the target function to points not necessarily on X.In this section, we make some comments on how to use the theory described in the previous sections for solving these problems.

Pre-image problem
In Remark 2.1, we have given an onto diffeomorphism Φ : U → [−1, 1] d where U ⊂ X.The pre-image problem is the following.Given a point y ∈ [−1, 1] d , find x ∈ X such that Φ(x) = y.This amounts to approximating the D-output function Φ −1 on [−1, 1] d , given its values at the known points {Φ(ξ) : ξ ∈ C}, and can therefore be solved using any of the techniques described in the previous sections.

Out of sample extension
One well known strategy for function extension outside the manifold is the following.One starts with a compact, positive-semi-definite symmetric kernel K : R D × R D and considers the eigen-decomposition of K restricted to X × X; thus, for example, if µ * is the volume measure on X, one finds numbers λ k ≥ 0 and orthonormal functions φ k on X such that A function on X can then be expanded in terms of the orthonormal system of functions {φ k }.The extension to R D is achieved by treating (3.1) as a definition of φ k on R D (which can be done since K is defined on R D × R D ), and the expansion of the original function f where the basis functions are now interpreted as extended by (3.1) as the desired extension.This leads to a variety of theoretical problems related to the judicious construction of kernels defined on the whole space whose eigenfunctions are meaningful as functions on X (e.g., kernels that commute with the Laplace-Beltrami operator) so as to allow such a construction.In the end, it is not clear how well this extension will behave outside of X.
Our construction gives an alternative method for extending a function on X to a tubular neighborhood of X, which we feel is more appropriate for most applications rather than trying to extend the function to the entire ambient space.Toward this goal, we first explain the local coordinate learning phase for tubular neighborhood of X.
Thus, there is no loss of generality in assuming that X is already a s dimensional submanifold of R D , defined by v, and with geodesic distance ρ 1 .This has several consequences.Even if one overestimates the dimension of the original manifold to be s rather than d, the resulting "distance respecting" coordinate system will also be "distance respecting" for the original manifold, except for the presumably small error resulting from the overestimate.If we have no information about d s), we may take s = D.This would answer the question regarding points off the manifold, as well as take noise into account.However, then the advantage of dimension reduction is lost.Also, all the constants will depend upon D rather than s (or d).
Having defined a local coordinate system for the tubular neighborhood of X in this way, we can then construct the local basis functions on this neighborhood as in Section 2.3.However, since the original data C is not dense on the local coordinate patch in the tubular neighborhood, one cannot use the ideas in Section 2.6.We would like instead to keep some control on the growth of the extension operator.For this purpose, we propose to use the minimal Sobolev norm (MSN) interpolant introduced in [5] and used very fruitfully in solutions of partial differential equations [4] and image segmentation [6].
Thus, using the procedures explained in these papers, we consider a differential operator ∆ depending upon the application, and find an integer N and coefficients c * k , |k| ∞ ≤ N , so as the s-variate polynomial over all s-variate polynomials P of coordinatewise degree ≤ N , subject to the conditions P (Φ(x j )) = f (x j ), x j ∈ C. (3.6) The polynomial P * then defines an extension of f to the tubular neighborhood of the local coordinate patch of X.
) 2 −(2j−2ℓ+1) 2 )t 2j+1 .(2.35)For y ∈ R d and multi-integer m = (m 1 , • • • , m d ) ≥ 0, the tensor product Chebyshev polynomial is defined by T m (y) = d j=1 .44)Therefore, rather than evaluating the Chebyshev polynomials as defined in (2.36),(2.35),onecanuse a multilayered network to evaluate V n,C (f ) in a more stable manner as follows.The first layer computes the coefficients h(k/n) f (C, k) using the available data.The output of this layer is input to a recurrent network to execute a multi-variate version of the well known Clenshaw algorithm[23, pp.78-80].