Dimension Independent Data Sets Approximation and Applications to Classification

We revisit the classical kernel method of approximation/interpolation theory in a very specific context motivated by the desire to obtain a robust procedure to approximate discrete data sets by (super)level sets of functions that are merely continuous at the data set arguments but are otherwise smooth. Special functions, called data signals, are defined for any given data set and are used to succesfully solve supervised classification problems in a robust way that depends continuously on the data set. The efficacy of the method is illustrated with a series of low dimensional examples and by its application to the standard benchmark high dimensional problem of MNIST digit classification.


Introduction
The problem of finding a function that explains a given set of data is a fundamental problem in mathematics and statistics.If the data are assumed to be the discrete manifestation of a function defined on a continuous (as opposed to discrete) domain of definition, the problem can be viewed as an approximation problem where the data can be leveraged to help identify a sensible approximation to the function.Often one resorts the prior knowledge about the target function to reduce the set of candidates from which to choose a good approximation, if not the best approximation.Within this framework, an extensive mathematical knowledge has been obtained over the past several decades along with a variety of powerful tools (see [1], in particular, for the philosophical approach taken here).In the current world, where data about almost anything you can imagine or wish for is available, one of the most interesting and often challening problems consists in extracting information, knowledge, and structure from high dimensional data.A variety of commonly used approaches belong to a category referred to as Machine Learning.Neural networks in all forms and shapes are particularly widespread due to their success in dealing with a series of challenging problems.Another class is that of Support-Vector Machines (SVM) [2], which were very popular before being somewhat superseded by improved neural networks.While they often use linear classification via hyperplanes, the so-called kernel trick [3,4] makes it possible for them to capture nonlinear decision boundaries albeit by working in a higher dimensional (feature) space to which the data is mapped but that admits an efficient computation of scalar products (via a suitable kernel).A connection between SVMs and neural networks was discovered in [5] and has since been investigated by many more authors.The method proposed here is philosophically in the category of SVMs but distinguishes itself by directly working with the data at hand without embedding into a higher dimensional feature space.This is made possible by directly looking for a good approximation in a large space of functions that allows for nonlinear behavior as opposed to a priori restricting the space of possible approximants by choosing a feature vector that replaces the data and on which the eventual approximant depends linearly upon.It will turn out that the end result of the approach taken here bears similarity with the We first review a classical method of inexact interpolation that yields a continuous function approximating a given data set.We do so by taking a PDE perspective that reveals important features that are exploited in order to obtain the announced stable method of classification even when the distribution of available data is not uniform in space and, possibly, noisy.In high dimension, even large data sets are often sparse due to the so-called curse of dimensionality.Moreover, the data is often supported on lower dimensional manifolds, where even dense data make up a thin slice of the ambient space.In the latter case, it is demonstrated in this paper that, while the data may not be sufficient for the reliable identification of a well-defined global approximant, it can still be used fruitfully (in a global or local fashion) for data analysis purposes.This is mainly due to the fact that the proposed method is capable of connecting the dots (data points) into a manifold by capturing geometric features of the data.
It is widely understood and accepted that a function interpolating discrete data should be at least continuous so as to provide a certain stability of prediction and resilience in the face of noise.In learning problems this is sometimes expressed as local constancy, even if that does not require continuity.Unless the data is known to stem from a very smooth underlying function, but, typically, even in that case, it should also not be exceedingly smooth, as this would lead to some blurring and reduce its ability to capture sharp transitions.It would, moreover, be advantageous for the interpolating function to depend at least continuously on the data set it is constructed from as, in that case, perturbations (due to measurement errors or to other sources) would not have a large impact on the outcome of classification based on the inperpolant.This kind of stability is akin to that discussed in [6] in the context of learning algorithms.
The starting point is a set of data consisting of point/value pairs In order to enforce minimal regularity on the interpolant function we take it from the space where S denotes the Schwartz space of tempered distributions, i.e. the topological dual of the space of smooth, rapidly decreasing functions.Thanks to the general Sobolev inequality it holds that where the containing space is that of bounded and uniformly Hölder continuous functions of exponent 1 2 .In order to obtain stability we may sacrifice some interpolation accuracy by not necessarily requiring the exact validity of Then, for approximate interpolation, u is determined by minimization of the energy functional given by for α > 0 and where the normalizing constant c d will be explicitly given in the next section (right below equation (2.6)).For exact interpolation u is determined by minimization of with constraints (0.2) that shall be summarized as u(X) = Y.Formally, it is set While in many practical problems the dimension d can be very large and the constant c d astronomical, this approach remains viable since the minimizers can be identified by solving a well-posed m × m linear system of equations for any α ≥ 0. The rest of the paper is organized as follows.In the next section we provide a detailed description of the method and obtain some of its basic mathematical properties.In the following section, we discuss a variety of numerical experiments that showcase the viability and efficacy of the method.There are interesting connections between this method and kernel based interpolation

The method
In order to derive a concrete method it needs to be shown that minimizers of E d can be computed efficiently.We first observe that the functional has a unique minimizer no matter what the given data set is.

Theorem 2.1
The functional E α has a unique minimizer u D,α for any given data set D.
Proof Take α > 0 first.Thanks to the embedding (0.1) the functional E α : H It is clearly also strictly convex as a quadratic functional since the first term is the square of a norm.Strongly lower semi-continuous convex functionals on a Hilbert space are known to be weakly lower-semicontinuous and so is therefore E α .It is also coercive since closed bounded sets in H d+1 2 (R d , R n ) are relatively weakly compact by the Banach-Alaoglu Theorem.A weakly lower-semicontinuous and coercive functional possesses a minimum, which, by strict convexity, is necessarily unique.The case α = 0 uses essentially the same argument where the full function space is replaced by the convex closed subset of functions satisfying u(X) = Y.Now that a unique continuous interpolant u D has been obtained for each given data set D, it needs to be determined in a usable form.The next step consists in deriving the Euler-Lagrange equation for u D .

Theorem 2.2
The minimizer u D,α for α > 0 satisfies the equation (system) in the weak sense, i.e. the equation holds in the space ) Here δ x denotes the Dirac distribution supported at the point x ∈ R d .If α = 0 it holds that in the weak sense for some λ i ∈ R n for i = 1, . . ., m.
Proof Notice that, thanks to (0.1), it holds that δ x ∈ H − d+1 2 (R d , R) for each x ∈ R d .Taking variations in direction of any function ve k for k = 1, . . ., n, where e k is the kth basis element of R n , with arbitrary v ∈ H d+1 2 (R d ) yields the equations , and where u D,k = (u D ) k .The identities amount to the validity of the system (2.2).The equivalence of the latter with (2.1) follows from and the fact that the (pseudo)differential operator (1 − ) d+1 4 is self-adjoint along with the validity of where the latter duality pairing is between the space H d+1 2 (R d , R n ) and its dual H − d+1 2 (R d , R n ).In the case when α = 0, one takes variations in direction of test C ∞ functions with ϕ(X) = to obtain that in the sense of distributions.This means that supp (1 − ) It is know that compactly supported distributions are of finite order and thus it must hold that (1 − ) d+1 2 u is a finite linear combinations of δ x i , i = 1, . . ., m, and derivatives of them.Notice, however, that ∂ l δ x / ∈ H − d+1 2 (R d ) for l = 1, . . ., d and so no derivatives can be present in the linear combination.This yields the claim.
Proof This regularity will readily follow from the reprentation of the solution discussed next.
The minimization of E 0 or a variety of similar functionals has long been recognized to provide an answer to the so-called universal approximation property in the context of learning.It indeed can be shown that The convergence takes place in the space's natural topology as X becomes a finer and finer, not necessarily regular, discrete grid that fills the whole domain.The curse of dimensionality, however, limits the applicability of this approximation procedure as the size of any finite "filling" grid is exponential in the ambient dimension.
A reason the above approach or, more in general, kernel based approximation or interpolation has found widespread use, is its ability to bridge the gap between finite and infinite dimensional spaces.This property amounts in the case at hand in the possibility of computing the continuous variable solution u D,α by solving finite dimensional systems.

Theorem 2.4
The minimizer u D,α , α > 0, is completely determined by its values u X,α = u D,α X = u D,α (x i ) i=1,...,m , on the set of arguments X = x i i = 1, . . ., m .The latter can be determined by solving the well-posed linear system where the matrix M D ∈ R m×m is given by for any x ∈ R d .For α = 0, it holds that Proof The Fourier transform of the Laplace kernel is known.Indeed where For the purposes of this paper the parameter ε is set to be 1.The right-hand side of the above identity is the symbol of the pseudo-differential operator where the second equality sign follows from (2.6) and well-known properties of the Fourier transform.It only remains to evaluate this identity on the arguments X to obtain the finite linear system.When α = 0, representation (2.3) entails that and therefore that as claimed.

Proposition 2.5
The matrix M is invertible and it holds that positive kernel as follows from its Fourier transform and the well-known characterization of positivity.Continuity of the inversion function inv Remark 2.6 The proof of Lemma 2.3 is now obvious since the explicit representation of u D,α reveals that singularities are only found on the set X.

Remark 2.7
The convergence u X,α → u X,0 as α → 0 implies convergence u D,α to u D,0 in the topology of H d+1 2 (R d , R n ) and hence, uniform convergence as well, thanks the embedding (0.1).
Remark 2.8 Since the continuous minimization problem (0.5) has a unique solution, the linear system (2.4) is assured to be solvable and, in fact, well conditioned also in parameter ranges of interest (but of course not in general).It also follows that its solution u X , as well as its extension u D to R d , depend continuously on the data D since the forcing term in (2.1) depends continuously on D. The latter follows from the linear dependence on Y and the fact that Dirac distributions depend continuously on the location of their support in the topology of H − d+1 2 (R d ), again a known consequence of (0.1).
Remark 2.9 Depending on the data set D, the values u X will be close or not so close to the prescribed values y i i = 1, . . ., m .Thus, if u D is considered an interpolant of the data, it will not be exact, but only approximately capture the data.In many applications, some of which are considered in next section, this is a small price to be paid for the gain in robustness that the approach guarantees.
Remark 2.10 Using directly that1 This shows that u D (x i ) = y i − αλ i and, in particular, reiterates the point about the convergence as α → 0. In practice, it is more convenient to work with this system in order to compute u D .
While the parameter α ≥ 0 plays an important role, it will be dropped from the notation from now on.The understanding is that its value can be inferred from the context and that, whatever its value is, it is kept fixed.In this paper we are particularly focussed on the case n = 1 and on the trivial value set Y = 1, where y i = 1 for i = 1, . . ., m.

Definition 2.11
If Y = 1, we say that u D is the (continuous) signal generated by the data X.We sometimes denote it by u X or u X,1 .
The signal is the inexact interpolation of the characteristic function of the data set 1 .It can be strong, if u X (x i ) 1, i = 1, . . ., m.This is the case, as was observed above, when X is a fine and locally filling discretization of the ambient space, such as when approximating a set of positive measure by a set of discrete points.More often, however, the signal will be weak in the sense that u X (x i ) is significantly less than 1 for i = 1, . . ., m.In this paper we contend that the usefulness of the signal u X does not only stem from its approximation or interpolation properties, but also (and perhaps mainly) from the fact that most of its level sets are very smooth due to Proposition 2.3 and, in fact, deliver smooth manifold approximations of X that can effectively be employed as stable decision boundaries in supervised classification problems.Superlevel sets of the signal are reliable approximations with positive measure of any discrete data set that prove robust against noise.They, in a sense, connect the dots and capture the shape of the data.Thus the main philosophical difference between the traditional view point, that considers the data as the manifestation of a function that needs to be reconstructed, and the view point taken in this paper is that here the data set itself is approximated by the mostly smooth (super)level sets of a function that may not even fit the data well at all.We are indeed more interested in the level surfaces generated by the data's signal than we are in its values.It will be shown that data signals, whether they are strong or weak, can be succesfully exploited in this sense.The practical experiments run in the next section will make use the following proposition.

Proposition 2.12 Let
the notation means that the union is disjoint) be a given data set consisting of N ∈ N subsets (or classes).For Proof The signal u D 0 and the signals u D 1 , . . ., u D N relative to X 0 , are solutions of the linear equations where 2 and where The claim therefore follows from the fact that Remark 2.13 Notice that signals do not, however, behave additively, in the sense that in general, even when If one is given a labeled data set X 0 , where N is the number of labels and X l , l = 1, . . ., N are the subsets consisting of the data corresponding to label l, i.e.L(x) = l for x ∈ X l , then one obtains a classification algorithm by computing the relative signals u X l for l = 1, . . ., N and assigning any unlabeled datum z ∈ R d to the class that exhibits the strongest relative signal, that is, (2.7) Remark 2.14 An important aspect of the proposed approach (and of kernel based methods as well) is that the fundamental solution of the (pseudo)differential operator 1 2 used to obtain the finite linear system is essentially dimension independent.It depends on it only through the Euclidean distance function, which is a minimal ingredient that can hardly be avoided.It would of course be extremly difficult to work directly with the (pseudo)differential operator or the energy functional in high dimension.An application to the well-known MNIST classification problem will be discussed in the next section using an approach based on the signal generated by the training data.In that case one has that d = 784.
Remark 2.15 Just as for exact interpolation and due to the fact there are no constraints like e.g.boundary conditions, the method is completely local.This means that an approximant can be computed based on a subset of the original data set that is confined or restricted to a subregion of interest.This feature will be exploited in the MNIST classification problem.

Remark 2.16
It is well-known that the parameter α > 0 has a regularizing effect that can be used to deal with noisy data when performing interpolation.It turns out that it also helps smooth out the level sets of u D .This is also illustrated in the next section.
Remark 2.17 It is sometimes convenient to modify the decay rate of the exponential "basis" functions, especially if the data undergoes some initial normalization.This can be done without significant consequences other than a slight modification of the objective functional or, equivalently, of the corresponding differential operator.Indeed, for γ > 0 and using the well-known scaling and translation properties of the Fourier transform F d , it holds that Here we use the notation u = F d (u) for the Fourier transform of the function u as well as • x and • ξ as place holders for the independent variables x and ξ in order to distinguish a function from its values.Furthermore τ y denotes translation, that is, and σ γ scaling, or σ γ u (•) = u(γ •).Replacing u by exp −2π| • | and y = x i for any i = 1, . . ., m it is seen that Thus the use of the modified exponentials merely corresponds to an inconsequential modification of the Euler-Lagrange equation (or its generating functional).

Related approaches
In this section we highlight some important connections between the point of view just described and well established frameworks.

Reproducing kernel Hilbert spaces (RKHS)
When α = 0, it can be seen that H d+1 2 is a RKHS with kernel K given by the fundamental solution G of (1 − ) d+1 2 via K (x, y) = G(x − y).Indeed, it holds by construction that As δ x 1 , . . ., δ x m are linearly independent vectors when the points are distinct and (1 − ) 2 is an isomorphism between H d+1 2 and H − d+1 2 , the functions K (•, x 1 ), . . ., K (•, x m ) are linearly independent and thus H d+1 2 is a fully interpolating RKHS according to the definition given in [7].It follows, in particular, that equation (2.4) is the equation that determines the minimal norm interpolant for the data (X, Y).In applications it is more common to start with a positive definite kernel (see [7] for a definition).While this may make little difference from a purely pragmatic point of view, it somewhat obfuscates the exact nature of the corresponding RKHS space and its norm.It is our contention that the norm can be chosen in a way as to shape the features of the associated kernel and interpolation.The discriminating power of the method arguably owes more to the choice of norm (and, hence, kernel) than to any adjustable parameter that may be present in the kernel function.The commonly used polynomial, exponential, and sigmoid kernels are all smooth in contrast to the kernel chosen here which has a carefully chosen regularity that determines the properties of the corresponding interpolant.The RKHS space is chosen as large as possible compatible with the continuity requirement discussed earlier.

Kernel support vector machines (K-SVM)
There is also a connection to SVM that use the kernel trick in order to maintain the idea of linear separation but to apply it to a feature vector generated by a special nonlinear transformation of the data set X.As the method reduces to computations involving only the scalar product of feature vectors, that can be easily computed via a kernel, this procedure does not render the method impractical (the feature vectors typically live in a possibly much higher dimension than the original data).In this paper, the motivation is rather to allow all nonlinear functions of a large function space to compete in order to interpolate the "data manifold" via their level sets.The choice of space is made in order to ensure the possibility of sharp (but continuous) transitions while maintaining the smoothness of (almot all) level sets.Notice that the choice of kernel in K-SVM is a straighforward trick that, however, conceals potential difficulties due to the increased dimensionality not necessarily present in the original data.We refer to [8,Section 12.3.4]for a more detailed discussion and only point out here that a choice of kernel may not necessarily be aligned with the structure of the data and make it hard for the method to identify the "data manifold" in the presence of noise.

Ridge regression
Another related method admitting solutions that can be described by kernels is that of ridge regression.The connection appears when α > 0. In kernel ridge regression [8] a data interpolant is looked for in the form where λ ∈ R n and h = (h 1 , . . ., h n ) is a vector of basis functions, by minimizing a measure of the error with a regularizing term given as a multiple of |λ|2 .In that case, if one defines the K n (x, y) = n j=1 h j (x)h j (y), a solution can be found in the form u = m i=1 λ i K (•, x i ) leading to a system very much like (2.4).Formally letting n tend to ∞, one would obtain a problem in infinite dimensions if the functions h j were chosen to be the eigenfunctions of some kernel (where available).The penalizing norm would, however, converge to an L 2 -type norm 2 and, moreover, a discrete set of eigenfunctions may not be available in general, as it is the case for the Laplace kernel used here.It is interesting to observe that the method described in this paper, like RKHS methods (α = 0), finds an infinite dimensional problem that necessarily has a solution which can be obtained by solving a finite dimensional one.
In the general category of data-smoothing models one can also find those developed and studied by Wahba and her school, see [9] for a nice account, where, starting with the one dimensional case, the functional is used to generate piecewise cubic polynomials as interpolants on an interval [a, b] ⊂ R.
The method was also extended to include other interpolants by using other differential operators and made practical for moderate dimensions by a clever (and not immediately obvious) use of tensor products.The approach taken in this paper is similar in spirit but is truly ambient dimension independent.

Numerical experiments
In this section a series of experiments are performed to illustrate the effectiveness of the method described in the previous section.First we consider two dimensional problems to highlight important aspects and in order to motivate and justify the use of the method in a high dimensional context.The section then concludes with an application to the classification of the MNIST data set.

Stability of signals' level sets
Working in the context of approximating measurable functions, simple functions play an important role as they are the building block of any measurable function.While the approximation property is well-known we present a few examples to illustrate the efficacy of the use of the data signal's level sets for classification purposes.First we will consider situations where the approximation is good, then examples when it is rather poor.It will be shown that, in all cases, i.e., regardless of how good the approximation is, the signal's level sets still provide very useful information.This observation is crucial since it opens the door to applications to high dimensional data, where it is inconceivable that the data arguments X provide a fine grid of even a small portion of the ambient space.

Characteristic functions of sets of positive measure
Take the three subsets of R 2 and consider the associated characteristic function χ S j for j = 1, 2, 3.The first data set consists of the values of these functions on a regular grid, that is, It follows that if a characteristic function has to be recovered or inferred from a data set, thresholding based on the interpolant u D m,j is an effective strategy and the decision boundary u D m,j = .5max(u D m,j ) is a solid choice across a range of values of the regularization parameter.Figure 2 depicts the same experiments using the denser data sets D 32,j for j = 1, 2, 3.
In Figs. 3 and 4, it is shown how the method performs in the presence of data corruption.In Fig. 3 2% of the data is misclassified, whereas the misclassification rate in Fig. 4 is 5%.By this we mean that a mistake is made, with the given probability, when a value is assigned to an argument by evaluating the corresponding characteristic function.These examples clearly demonstrate the usefulness of the regularizing parameter which leads to data signals whose decision level sets are more stable in the presence of classification errors.

Sample data
Finally we demonstrate that the method offers a degree robustness when the data arguments are randomly sampled from a uniform distribution supported on the box B. The resulting decision boundaries of half maximal value are depicted in Fig. 5 along with the sampled argument data sets X m .The sampling rate clearly affects the smoothness of the level sets, a deterioration that is to some extent counteracted by the regularization.
Fig. 1 Level lines of the signal u Dm,j for m = 16, j = 1, 2, 3, and α = 0.1, 1.0, 2.0.Depicted are the data set (blue dots correspond to the value 1 while magenta dots to the value 0) and three level lines corresponding to levels at 20%, 50%, and 80% of the signal's maximal value, respectively.Darker lines correspond to higher levels.The parameter α grows from left to right.The boundary of the set S j appears as a dashed black line Fig. 2 Level lines of the signal u Dm,j for m = 32, j = 1, 2, 3, and α = 0, 1, 1, 2. Depicted are the data set (blue dots correspond to the value 1 while magenta dots to the value 0) and three level lines corresponding to levels at 20%, 50%, and 80% of the signal's maximal value, respectively.Darker lines correspond to higher levels.The parameter α grows from left to right.The boundary of the set S j appears as a dashed black line

Classification
Continuity of the interpolant and its level sets (almost all of them actually smooth) are obtained at the cost of approximate interpolation.Such an approximation can still be accurate when the argument data set covers the function's domain of definition uniformly and the value set is accurate, but the real advantage of this method is its applicability to incomplete data and/or noisy data sets.This point is further reinforced with the next series of experiments, where the data build a lower dimensional manifold of the ambient space and its signal is weak, that is, information about the underlying function is limited to sets of zero measure, or the data is not deterministic (in its argument set) but only has a probability distribution for its location.
Consider the data sets D k , k = 1, 2, consisting of points belonging to two distinct classes: points along a circle and points on the union of two segments with two different densities as depicted in Fig. 6.For k = 1, 2, denote the circular data sets by and the union of segments data sets by In the spirit of the previous examples, we create two pairs of data sets and compute the associated signals Figure 7 shows the level sets for different values of the regularization parameter.The region in their interior (i.e. the one containing the corresponding data set) can be considered as a smooth fattening of the data set to a set of positive measure.It can be obtained for any data set regardless of the intensity of its signal.
Next we turn our interest to the question of classification: given a point z ∈ R 2 that needs to be classified, we use the decision algorithm defined by (2.7).This gives which, in this case, yields the level lines (hypersurfaces in higher dimension) as the decision boundary.This is illustrated in Fig. 8 for the classification problem of the two class pairs (C k , S k ) for k = 1, 2 introduced above.
If the data pairs (C k , S k ), k = 1, 2, are considered the ground truth, then the above decision boundary is arguably optimal.If, on the other hand, it is known that the actual sets are the continuous circle and the union of two segments, then the data are only a sampling of these sets.In this case, the decision boundary may be biased by the relative oversampling of the one set compared to the other.This is evident when one compares the decision boundaries of Fig. 8.In concrete situations, if information about the dimensionality of the ground truth is known, this effect can be mitigated by using comparable sampling rates for the different classes (see next section for an example of this procedure) or by normalizing the data fidelity term to read 1 2m m i=1 u(x i ) − y i 2 .We conclude this section with a classification problem for data X 0 split into three classes X l , l = 1, 2, 3, each consisting of a set of points which are normally distributed with mean p l ∈ R 2 and different covariance matrices.The data set and the corresponding decision regions computed based on (2.7) are depicted in Fig. 9.In these experiments α = 1.

The MNIST data set
The final application is to the standard machine learning example and toy problem of digit classification for the MNIST data set.
Remark 2.18 The purpose of this example is to illustrate how the choice of regularizer allows for the method to be used for data that live in high dimension.There is no effort to optimize the parameter choices nor to compete with state-of-the-art algorithms.The idea is rather to show how a theoretically justified and derived method with a small number of parameters can perform an acceptable job on a high dimensional problem.The method has its limitations as it assumes that the data classes (almost) lie on separable manifolds.This may be approximately true for MNIST but is certainly not true for other data sets.
The data set consists of 28×28 grayscale images of hand-written digits stored as vectors in [0, 255] 784 .Here it is considered that the ambient space is simply R 784 .The argument data is normalized to have unit Euclidean norm, that is, each original vector x is replaced by x/|x|.In this way the maximal Euclidean distance between any two data points is 2 √ 2. The data set is split into a training set containing 60,000 data points x and their corresponding label d(x) indicating which digit is represented, and a testing data set of size 10,000.The labels of the testing data are known but need to be inferred from any knowledge that can be gleaned from the training set.This an example of when, due to the so-called curse of dimensionality, the data does not have any hope to fill the ambient space uniformly and thus, even if one assumed the existence of an underlying function d : R 784 → {0, 1, . . ., 9}, the data would never be sufficient to accurately approximate it.It has to be said of course, that the testing data mostly does not stray away significantly from the training data and the different digits in the latter build thin subsets of the ambient space.This fact is typically captured by saying that the data lives in some lower dimensional manifold(s).We know from the previous section and from the two dimensional experiments, however, that the (training) data still generates a significant, if not strong, signal.The classification method described in the sequel exploits this signal and does not require any kind of training based on the minimization of non-convex functionals, as is often the case for machine learning algorithms based on neural nets.It is, in fact, based on the solution of low dimensional, well-posed linear systems as is about to be explained.First, in order to strengthen the signal somewhat, the training set is expanded to include rotations by ±10 • and horizontal/vertical translations by ±2 pixels of each image.Then, given a test image z, the closest 5 training images of each digit class are determined X = x i j j = 1, . . ., 50 , where d(x i j ) = j/5 .The idea is now to use system (2.4)This approach yields an accuracy of 98.56%3 on the test set.In Table 1 we record the detailed outcome of the classification, performed with α = 1.5.
Recall that this method is stable and depends continuously on the data set and hence delivers a robustness that methods with higher classification rates typically do not.Moreover, unlike neural networks, it does not require any training but uses the training set directly in a fully transparent way.

Fig. 1
some level lines of the interpolant u D m,j are shown for the three functions χ S j , j = 1, 2, 3, and for different values of the regularizing parameter α > 0 and m = 16.While the size of the data set clearly correlates with the "accuracy" of the interpolation, the approximating function does an excellent job at generating meaningful and smooth level sets.Their smoothness is affected mainly by the parameter α and the their proximity to the level sets corresponding to the highest values.

Fig. 3 Fig. 4 Fig. 5
Fig. 3 Decision boundary of the signal u Dm,j for m = 16, j = 2, 3, and a = 0, 1, 1, 2 with 2% data corruption rate.Depicted are the data set (blue dots correspond to the value 1 while magenta dots to the value 0) and the level line corresponding to 50% of the signal's maximal value.The parameter α grows from left to right.The boundary of the set S j appears as a dashed black line

Fig. 6
Fig. 6 From left to right, the data set pairs (C k , S k ), k = 1, 2. They can be thought as different samplings of the same pair of "continuous" sets

Fig. 7 2 Fig. 8 2 Fig. 9
Fig. 7 The 50% of maximum level lines for the signals u D C and u D S of the data sets based on the two data classes C k and S k for k = 1, 2. The first two rows correspond to k = 1 and the second two to k = 2. Withing each row, from left to right, the regularization parameter is α = .1,1, 2 in order to produce approximate interpolants u d := u D d of the characteristic functions of each digit class given by the data setD d = x i j , δ d,d(x i j ) j = 1, . . ., 50 ,where d = 0, . . ., 9, andδ d, d = ⎧ ⎨ ⎩ 1, d = d, 0, d = d.Finally the approximative characteristic functions u d will compete to determine the digit d(z) to be associated with the test image z via (2.7), in this case d(z) = argmax d u d (z).