Approximation of almost diagonal non-linear maps by lattice Lipschitz operators

Lattice Lipschitz operators define a new class of nonlinear Banach-lattice-valued maps that can be written as diagonal functions with respect to a certain basis. In the $n-$dimensional case, such a map can be represented as a vector of size $n$ of real-valued functions of one variable. In this paper we develop a method to approximate almost diagonal maps by means of lattice Lipschitz operators. The proposed technique is based on the approximation properties and error bounds obtained for these operators, together with a pointwise version of the interpolation of McShane and Whitney extension maps that can be applied to almost diagonal functions. In order to get the desired approximation, it is necessary to previously obtain an approximation to the set of eigenvectors of the original function. We focus on the explicit computation of error formulas and on illustrative examples to present our construction.


Introduction and notation
Diagonalizable operators play a fundamental role in classical operator theory.Indeed, they have special geometric features that allow the mathematical treatment of properties in many fields of science, and can be used to approximate linear operators in infinite-dimensional Banach spaces.However, there is no analog for the case of Lipschitz maps in Euclidean spaces, although there are many situations in which both classes (linear and Lipschitz maps) can be naturally associated.Several authors have recently been interested in the determination and characterization of classes of operators for which the notion of diagonalization could be adapted.To do this, one must begin by giving a definition of eigenvector that can make sense in the specific context in which one wishes to investigate.Nonlinear spectral theory is already a classical subject and has a long development, both from a theoretical and practical point of view ([2, 5, 10]).For example, for the case of multilinear mappings [6,7,8,16], polynomials [14] or Lipschitz maps [9].In particular, we have studied suitable versions of these notions for Lipschitz maps in [3,9].
Our aim in the present paper is to adapt for the case of lattice Lipschitz operators the classical extension/representation of diagonalizable linear operators once a basis of eigenvectors is known.The idea is to use the same representation pattern (but in an approximate form) for the case of functions that are approximable by lattice Lipschitz maps, which is a broad class of functions containing linear operators.As we will explain, our method is based on the use of the order in the Euclidean lattice, rather than linearity.
Thus, the class of lattice Lipschitz operators has recently been introduced to fill the gap in the helpful definition of diagonalizable (linear) operators for the case of Lipschitz [3] maps.Using the results presented in the aforementioned paper, we explain here an approximation method for "almost diagonalizable" Lipschitz operators, i.e., Lipschitz functions that can be approximately computed as lattice Lipschitz maps.In this paper we provide specific tools for the (approximate) representation of general Lipschitz maps by means of their eigenvectors using the vector-valued lattice versions of the classical McShane and Whitney formulas that were obtained in [3].
Although some theoretical statements will also be presented, this paper is written from the point of view of applied mathematics, with the idea of providing an efficient algorithm for the approximation of Lipschitz maps with the requirements, as we have said, of being "almost diagonal".An explanation of the "initial class" of lattice Lipschitz operators, together with some clarifying examples, are given in Section 2. General results on error bounds for such approximation can be found in [3], but we will go further in this direction in Section 3 of the present paper, which is mainly devoted to the approximate calculation of a "sufficient" subset of eigenvectors.As will be explained, our method is based on the McShane-Whitney extension which is based on the order of the studied operator restricted to that subset of eigenvectors, as if it were a lattice Lipschitz operator.Once such a set is obtained, we present in Section 4.1 how to define a suitable order in Euclidean space that allows the use of McShane-Whitney expressions.Sections 4.2 and 4.3 are devoted to give the representation formulas and the final computational algorithms.In Section 5 we show an illustrative example, and in Section 6 we give some final comments.The algorithm used to compute a set of approximate eigenvectors, programmed in R, is given in Appendix Appendix A.
We will use the notation of general topology and mathematical analysis.We will write R for the set of the real numbers endowed with its standard Euclidean distance.If (M, d) and (D, ρ) are metric spaces, we say that a map T : . The Lipschitz constant of T is the infimum of all such constants K.
We will center our attention on Lipschitz-type functions between Euclidean spaces E. Recall that the classical McShane-Whitney Theorem ( [4,15,17]) establishes that, if S is a subset of a metric space (D, d) and T : S → R is a Lipschitz function with Lipschitz constant K, we can always find an extension to D with the same Lipschitz constant.There are two classical ways of computing such an extension, that are provided by the formulas Suitable versions of such formulas, adapted for Lipschitz operators T : E → E, will be the main computational tools in this paper.
Given a Euclidean space E of dimension n and a basis in it, we can always define an associated order in E that is provided by the coordinatewise ordering of the vectors (two vectors x = (x 1 , ..., x n ), y = (y 1 , ..., y n ) ∈ E are ordered, x ≤ y, if x i ≤ y i for i = 1, ..., n).This gives, when the 2-norm of the coordinates is considered, a Banach lattice structure (E, • , ≤).From this point of view, each vector represented by its coordinates in the chosen basis can be considered as a function f, {1, ..., n} w → f (w) = x w .The standard symbols ∨ and ∧ are used for the maximum and the minimum of two (or several) vectors in the lattice, respectively.

Lattice Lipschitz operators
As has been shown in [3], the operators belonging to the special class of Lipschitz maps, that are called lattice Lipschitz operators, have a relevant property which motivates the method we propose in this paper.
Property.Each lattice Lipschitz operator T can be written as the Mc-Shane extension (equivalently, the Whitney extension) of the restriction T | R(B) of the operator itself to the set of the rays R(B) generated by a basis B of eigenvalues of T.
Although we cannot expect all Lipschitz maps to behave like lattice Lipschitz operators, we can use this specific class as an approximation family.The idea is to find suitable extensions for each (approximated) eigenvector basis we can obtain for E, and mix them in a common structure to find good approximations to such a Lipschitz operator.
On the other hand, the concept of diagonalizable linear operator on finite dimensional Euclidean spaces can be extended to the setting of Lipschitz maps by means of the notion of lattice Lipschitz operator.Essentially, these are Lipschitz maps that can be written as combination of real valued Lipschitz functions that works independently in the directions of a given basis of the space.The main characterizations and results regarding this class or maps have been extensively studied in [3].In what follows we recall the main technical definitions and results, as well as some illustrative examples.
Let E be a finite dimensional linear space R n .Let B be a fixed basis of E. We will consider the order ≤ provided by the pointwise order of the coordinates of the vectors of E in the basis B. As usual, using the coordinate representation given by B, every vector in x = (α 1 , ..., α n ) ∈ E can be understood as a function x : {1, ..., n} → R, x(w) = α w , w ∈ {1, ..., n} = Ω.
In this setting, we recall the definition of our main tool.
Definition 2.1.(Definition 1 in [3]) Let E 0 be a subset of E. A Lipschitz operator T : E 0 → E is lattice Lipschitz (with respect to the order ≤ associated to the basis B) if there is a function The pointwise infimum of the functions K in the inequality above is called the associate function to T.
As we will see below, lattice Lipschitz operators can be identified with another class of functions that allows a geometric description.A Lipschitz operator T : We call the functions f i the coordinate functions of T with respect to the basis B.
Example 2.2.Let us give an example of diagonal map.Consider the function Both the coordinates in the domain (x, y) and in the range are assumed to be with respect to the canonical basis of R 2 .To find the eigenvectors and the eigenvalue functions, just notice that we have to find the vectors (x, y) such that S(x, y) and (x, y) are linearly dependent.This means that that is, y(x 2 + y 2 ) = 2x 2 y, what leads to the solutions y = 0 with all x ∈ R, or y 2 = x 2 , and so, the set of eigenvectors is Consider the basis for R 2 given by B = {(1, 1), (1, −1)}.Then we can consider the change or coordinates provided by (x, y) = z(1, 1) + v(1, −1), what gives and so In this new basis, we can write the map S as Consequently, S is a diagonal map defined by the Lipschitz eigenvalue functions (that is, the functions that provide the eigenvalue associated to a given eigenvector) , that are also Lipschitz in the domain.
An easy computation taking into account that the function is defined on B R 2 gives that the associate function is However, note that S is not diagonal with respect to the basis D = {(1, 0), (1, 1)}.Indeed, to be diagonal would mean that for the change of coordinates provided by (x, y) = z(1, 0) + t(1, 1), we should have a diagonal representation.But note that Since (z 2 + 2t 2 + 2zt, 2zt + 2t 2 ) cannot be written as s 1 (z)(1, 0) + s 2 (t)(1, 1), we have that S is not diagonal with respect to D.
Next result establishes the identification among diagonal and lattice Lipschitz operators.
Theorem 2.3.(Theorem 3 in [3]) A Lipschitz operator T : E → E is lattice Lipschitz (with respect to a certain order ≤) with associate function K : Ω → R if and only if T is diagonal with respect to the basis B (associated to the order ≤) with Lipschitz constant functions K(i), i = 1, ..., n.
Note that, following Theorem 2.3 and taking into account the last part of Example 2.2, we know that an operator can be diagonal (and therefore lattice Lipschitz) with respect to a given eigenvector basis B, and not be diagonal with respect to another eigenvector basis D. In this case, using the theorem we obtain that the operator is not lattice Lipschitz with respect to the order defined by D.
Although we do not use it in this paper, we should note that Theorem 2.3 allows a "local" version, i.e., it can be adapted for operators that are only defined on subsets of E satisfying certain properties.In its generality, this result opens the door to the approximation of Lipschitz operators by means of lattice Lipschitz extensions.Since having an eigenvector basis will play a central role in this, we will fix in the next section a method to deal with the (approximate) computation of an eigenvector basis in the case of Lipschitz operators.
Example 2.4.The previous result can be directly applied to characterize lattice Lipschitz operators.Let us explain with an easy example how to do it.Take the function G : R 2 → R 2 given by the formula where (x, y) are the coordinates with respect to the canonical basis, as well as the coordinates in the range.
1) Let us compute the conditions for being eigenvectors of the map.The requirement (that has to be computed pointwise) is given by the equation We directly get that this relation is satisfied for all vectors as (x, 0), x ∈ R, and what gives x = y.Therefore, all vectors as (x, 0) and (x, x), x ∈ R, are eigenvectors.
2) Consider the basis B = {(1, 0), (1, 1)}, and take the coordinates of the vectors in this basis to be (z, t).The equation gives the change of coordinates Therefore, the formula for (z, t) ∈ R 2 , gives a diagonal representation for the function G.
3) The functions g 1 (z) = z and g 2 (t) = |t|/(1+|t|) are real valued Lipschitz functions with constant 1.It is obvious for g 1 .For g 2 , just note that for all t 1 , t 2 ∈ R, what means that the Lipschitz constant is less or equal to one, and doing t 1 → ∞ and t 2 = 0 we see that this constant is 1.The associate function is then given by K(1) = 1 and K(2) = 1.
Summing up, all the vectors as z(1, 0) and t(1, 1) are eigenvectors, with eigenvalues e 1 z(1, 0) = 1 and e 2 t(1, 1) = |t|/t(1 + |t|).Using Theorem 2.3, we can say that G is a lattice Lipschitz operator when the order is the one inherited by the basis B.

Estimates of the eigenvectors of a Lipschitz operator: the exponential distribution Monte-Carlo model
Following the arguments provided in the previous sections, we need to find sets of vectors that are approximately eigenvectors of any Lipschitz map we want to analyze.In this section we explain how to perform a statistical procedure to obtain such sets.As we have seen, to know a basis of eigenvectors with respect to which a given Lipschitz map is diagonal allows an easy representation of a lattice Lipschitz operator.However, we plan to establish an approximation procedure for a broad class of Lipschitz maps, so we cannot expect them to be lattice Lipschitz.We intend to use the same formulas that provide representations of lattice Lipschitz maps, so we need to determine-at least, approximately-the set of eigenvectors of a given Lipschitz map in order to fix a convenient basis of the space, if possible.
To find a set of (approximated) eigenvectors we mix geometric and stochastic arguments.Suppose that E is a real Euclidean space, and recall that its dual space can be directly identified with the original space E and the dual action is the scalar product.Take a Lipschitz function T : E → E. Following the same ideas used in [9, Section 2] we define the diagonal value λ T (x) of T (x) as the real number that satisfies We will simply write λ(x) instead of λ T (x) if T is fixed in the context.Note that if x is an eigenvector of T with associated eigenvalue β ∈ R, the real number λ(x) defined as above gives the eigenvalue, that is, The main property of the diagonal value is that it represents the optimal diagonal approximation to T (x).Furthermore, for the case of Euclidean spaces, the optimal diagonal approximation coincides with the projection of T (x) on the unit vector in the direction of x.We explicitly write this in Proposition 3.1, including the proof-a straightforward consequence of the Euclidean geometry-for the aim of completeness.
Let us define before the diagonal projection error, ε T (x)(α), as the function of α ∈ R that represents the size of the difference between T (x) and α • x.After normalization, the formula for this quantity is given by We will write ε T (x) for the minimal value of the diagonal projection error and, as in the case of λ, ε(x) if T has been already fixed.
Proposition 3.1.Let X be a Euclidean space, T : X → X a Lipschitz operator and x ∈ X.Then the diagonal value λ T (x) = λ(x) is the real number that minimizes the diagonal projection error.That is, the diagonal error is given by Proof.It is given by a direct computation.Note that the solution of the equation we get the result.
Corollary 3.2.For every point x ∈ E, the error ε(x) is given by the formula Proof.It is a consequence of Proposition 3.1 and a direct calculation involving the definition of λ(x).Indeed, we have that Using this result, the idea is to perform a Monte Carlo procedure to obtain a large enough set of eigenvectors of the Lipschitz operator T. It must be said that most of the effort for the design of Monte Carlo methods for operator diagonalization comes from quantum physics and stochastic analysis [11,13,18,19].In general, these approaches do not provide a fundamental framework to support a general methodology, since they focus on the actual computation of quantities with some physical or mathematical meaning.Consequently, we propose in what follow our own procedure based on the functions defined just above.
We fix a suitable uniform value to be accepted for the error committed when the operator is approximated by its diagonal value.We will use such a set for the computation of the McShane and Whitney representation formulas, as a substitute of the exact set of eigenvalues that could not be computable.An example of such situation coincides with the case in which the Lipschitz operator T is the addition of a linear diagonalizable operator L plus a perturbation P with small norm, T = L + P.
Therefore, we will base our method on a sampling procedure supported by probabilistic arguments using the normal distribution and following the next steps.
(1) We start by fixing a bounded set in which we will search for our approximate eigenvectors; if no additional information is known, we choose a product P of intervals in R n centered in 0; in case we know previously that the eigenvectors are located in a particular set M, we use it instead.
(2) We use the uniform distribution to sample a starting set of vectors S 0 .
In case we have some previous knowledge on the set, we can introduce Bayesian procedures to fix a more accurated probability distribution Ψ for doing the sampling.
(3) We compute the functional ε(x)-where ε(•) is given by the error formula provided in Proposition 3.1-for all x ∈ S 0 .Now, we consider the set S 1 of all the points x that are most similar to an eigenvector in terms of having small ε(x).This can be done by selecting the points with smallest ε(x)-for example, 10%. of the points of S 0 -.
(4) Now, for every s ∈ S 1 , we start an iterative process.Sample a fixed number of points around s using a gaussian distribution centered on s and with variance depending of the error ε(x) formula where τ is a fitting parameter.Select the point with the smallest ε(•) and repeat this step a fixed number of times-changing the value of τ if necessary-.
Taking into account that the only property we know about the function T is that it is Lipschitz, this distribution allows to center the sampling near the points in which the diagonal error is small.By the Lipschitz inequality T (x) − T (a) ≤ K x − a , T (x) and T (a) are controled by the distance between x and a, so using the proposed distribution maximizes the probability of getting approximated eigenvectors in the sampling.
The computations involved in the algorithm whose scheme has been presented above are easy to perform using R. To show concrete situations, we will show some numerical examples of functions T : R 2 → R 2 , since they allow a graphical representation.The algorithm used can be found in the Appendix Appendix A of this work.The calculations shown below have followed the next rule.We consider the domain subset [−5, 5] × [−5, 5].We start with N = 500 initial points in the domain, which are obtained randomly.We choose the N 0 = 100 best with respect to the error value.Using the distribution ψ written in step (3) of the algorithm with τ = 5, for each of these points we generate N 1 = 10 points around, from which the best one is selected.We repeat this step 10 times.Let us show the graphical representation of the results with three different values of the parameter r.
• r = 0.The sine term in the first coordinate is eliminated.This makes the example simpler, with an easy to understand representation of a suitable subset of approximate eigenvectors.
• r = 3.In this case, the sine term in the first coordinate causes an important perturbation, producing a more dispersed eigenvector structure.• r = −10.The sine term causes a stronger perturbation (of negative sign).The consequence is that the set of eigenvectors no longer follows (even approximately) clear lines.These examples show that, although we cannot expect a diagonal distribution of eigenvectors for non lattice Lipschitz maps, we can sometimes treat general Lipschitz operators as perturbations of lattice Lipschitz maps, in the case where some diagonal distribution of a set of approximate eigenvectors is still preserved.Since no clearly defined axes are given for such a set, we have to provide a technique for the generation of two straight lines (or n straight lines in the general case) that can play this role.How to do so will be shown in the next section.

General algorithm and examples
Let us show in this section how the extension/representation formulas (McShane and Whitney versions) that work for the case of lattice Lipschitz operators can be adapted for any Lipschitz operator T : S ⊂ E → E in order to obtain an approximate functional expression for T (Definition 2, Proposition 1 and Proposition 2 in [3]).In the case of a lattice Lipschitz operator L, if we fix S to be the union of the rays defined by a certain basis of eigenvectors of L for E, both of these formulas give exact representations of L (Theorem 4 in [3]).That is, the operator L coincides with both (L| S ) M and (L| S ) W .In case we consider S an arbitrary subset of E, the McShane and Whitney formulas give approximations, for which the error expressions are known.These adapted formulas are where each w denotes the index on the corresponding element in the basis B, and the functions K(w) is the pointwise evaluation of the Lipschitz constant for each coordinate.Therefore, the use of these formulas explicitly requires an order in space.Indeed, the expression |s − v|(w) appearing in them is calculated using the order provided by a basis.In the next subsection we propose a method for defining such an order for the case of general Lipschitz maps.

The definition of the order for the approximate representation of a Lipschitz operator
As we are designing an approximation method, the procedure to find a good basis has to be related to the mass distribution of the set of approximate eigenvalues.The main idea consists in defining a partition of the set of approximate eigenvalues into n sets that are intended to describe the mass distribution.The centres of mass of the subsets of the partition give the n−dimensional basis necessary for the definition of the lattice order.Any clustering method could provide a technique for doing so.In this section we explain a method based on observing the mass distribution of the set of approximate eigenvectors obtained, from Principal Component Analysis of the point cloud defined by this set.Fix an operator T : R n → R n .The proposed method follows the next steps, that give different solutions depending on the symmetry of the set of approximated eigenvectors that are obtained.
(1) The direct case: if the approximated eigenvectors are distributed around a set of n vectors that are linearly independent, we take them as the adequate basis B.
(2) Otherwise, we compute the PC (Principal Components) of the cloud of approximated eigenvectors of T, that has been obtained.This technique is widely used for determined the main trends that can be detected in a point cloud (see for example [1,12]).We define the new (orthonormal) basis B using the PC.This is a candidate for being a good basis in case the symmetry of the sets that define the distribution of mass of the eigenvectors coincide with the directions of the vectors provided by the PC.The main problem is that this method always gives an orthogonal basis, which, as we have seen, is not necessarily the best way to describe the distribution of the eigenvector distribution, even if we have a lattice Lipschitz map.
(3) Suppose now that the cloud of approximate eigenvectors is not oriented following the direction of any set of n vectors, but can be found in a particular region of space.In this case we consider the octants (or hyperoctants) defined by the orthogonal basis provided by the PCA.We will choose a new basis C defined by the vectors crossing these octants along their axes of symmetry.That is, take σ to be any of the elements of [−1, 1] n .We consider the vectors, expressed in the orthogonal basis provided by the PC, as For example, we get the first vector of C to be (1, 1, 1, ..., 1)/n 1/2 , and we choose other n − 1 vectors as the ones above to complete a basis.
In the spirit of setting a concrete procedure for this article, we will follow the rule explained above for the definition of the order in the lattice in the next section.As we have said, this is not the only way that can be proposed to obtain such an order.In general, any rule for defining a suitable basis should depend on the symmetry of the problem.

Weakening the lattice Lipschitz inequality
Although the definitions of the lattice versions of the McShane and Whitney extensions provide accurate results for diagonal Lipschitz maps, we cannot expect such good behavior for operators that are not exactly lattice Lipschitz.If the mapping T is not diagonalizable or the set of axes is not exactly determined, the assumptions of Theorem 2.3 are not satisfied, so that T may not be a lattice Lipschitz operator.This makes that the associated function K(w) may be too large, and therefore also the error of the approximation.
The solution we present for this problem is to change the condition (1) to where 0 ≤ α ≤ 1 and || • || is a norm on R n .Writing it in terms of the order of the space E, as where 1 denotes the constant function one in Ω.
Observe that if α = 0, the condition (3) is the same as (1) and, if α = 1, it is equivalent to every coordinate function T i of T being a real Lipschitz function.In the case that T = L + P is the addition of a lattice Lipschitz operator L plus a perturbation P , which we assume to have a small Lipschitz constant, T satisfies where K is the associated function of L and C the Lipschitz constant of C. Note that, to control the function K of (1), α can be smaller the smaller is the perturbation P .If a function defined on a subset S of R n , T : S → R n satisfies (3), we can also define the McShane and Whitney extensions as Following the arguments presented in [3], we find that the error bounds become in this case The idea now is to use the previous method to approximate a function, using for doing so a small α > 0.

General algorithm
Using the formulas explained, we apply the complete procedure following the steps below.
(1) Fix an operator T : R n → R n that one wants to reconstruct from its eigenvectors.If the set of eigenvectors is known, sample some points in it, calculate the value of T and go to the next step.Otherwise, the eigenvector set is approximated by a Monte Carlo procedure based on diagonal error minimization as is explained in Section 3.Call S the set of -approximate-eigenvectors.
(2) For fixing the order on R n , use the method exposed on Section 4.1 based on PCA -other procedures could also be possible-.
(3) Fix a small α > 0 (for example α = 0.1) and compute the best K(w) possible by using the formula and the McShane and Whitney lattice-type formulas associated to the order provided in the previous step.
(4) An interpolation of these formulas provide a(n approximate) representation of the original operator T. We control the error commited by these formulas by using the error associated to the Lipschitz inequality and the formulas in (4), together with the diagonal error of the approximation of the set of eigenvectors.

A numerical example
Let us show how the method explained in Section 4.2 works in a concrete numerical example.Let E = R 2 and consider the function f : R 2 → R 2 defined as the function (x, y) → (x + y, x − y) (which is a diagonal map, and a lattice Lipschitz operator), with a small perturbation: First of all, we approximate the set of eigenvectors with the Monte Carlo method explained above.Let S 0 be a sample 250 points using a uniform distribution on the set P = [−5, 5] × [ −5, 5].Compute the diagonal error ε(•) at each point of S 0 , and select the 50 of them with smallest error (20% of the points on S 0 ).Write S 1 for this set.Now, for each element x on S 1 , we start an iterative process to find a set of approximated eigenvectors as explained before.After 5 iterations, the result (the set S 2 ) is plotted in Figure 4.In the next step we choose the axes that will define the lattice structure of E. In this case, we apply Principal Component Analysis (PCA) [1].The resulting new axes can be seen in Figure 5.These axes "look good" because they have a distribution similar to that of the true eigenvectors.
The final step is to compute the McShane and Whitney formulas (with α = 0.1 and the Euclidean norm || • ||), using the points of S 2 and the order given by the orthonormal basis provided by the PCA.The best function K is in this case K(1) = 2.24, K(2) = 3.26.
Observe that if the norm modification is not considered on the lattice Lipschitz inequality (α = 0), the best K possible is much larger: K(1) = 16.2,K(2) = 220.4,which would cause a worst approximation.The approximation result computed as the mean value of the McShane and Whitney formulas, can be seen in Figure 6.
In order to compare our approximation with the original f , we compute the error using a Monte-Carlo procedure, and we obtain ≈ 0.65.
The pointwise error is bounded; using the formulas (4) we obtain for each component, i = 1, 2.

Conclusions
The notion of lattice Lipschitz operator in finite-dimensional normed spaces has been introduced to provide a suitable set of Lipschitz-type operators that can be used for the design of approximation algorithms.Since any lattice Lipschitz operator always allows a diagonal representation, the family of functions to which this approximation method is applied is composed of nonlinear functions with the property that they allow an "almost diagonal" representation.
This makes it necessary to find an approximation method to find the "almost eigenvalues" of the objective function, and, in a second step, to determine a good set of lattice Lipschitz maps that can be used as an approximation family for it.We propose a concrete algorithm, and show with an example how it works, taking into account the measure of the error made when using the approximation, whose formulas have also been obtained in the paper.

Figure 1 :
Figure 1: Left: first 500 random points with the 100 chosen points for r = 0. Right: final set of approximated eigenvectors.

Figure 2 :
Figure 2: Left: first 500 random points with the 100 chosen points for r = 5.Right: final set of approximated eigenvectors.

Figure 3 :
Figure 3: Left: first 500 random points with the 100 chosen points for r = −10.Right: final set of approximated eigenvectors.

Figure 4 :
Figure4: On the left, first approach to the set of eigenvectors, S 0 in gray and S 1 in black.The set S 2 obtained after 5 iterations can be seen on the right.

Figure 5 :
Figure 5: The new axes resulting from PCA; in red, the first principal component and, in blue, the second one.

Figure 6 :
Figure 6: Approximation f (in blue) of f (in orange)-left: first component, right: the second one-.
Appendix A. R code l i b r a r y ( g g p l o t 2 ) l i b r a r y ( t i d y v e r s e ) ev lambda <− f u n c t i o n ( v , f ){ l a = i f e l s e ( a l l ( v==c ( 0 , 0 ) ) , 0 , sum ( f ( v ) * v ) / sum ( v ˆ2 ) ) r e t u r n ( l a ) } e v e r r o r <− f u n c t i o n ( v , f ){ e v l a <− ev lambda ( v , f ) i f ( a l l ( v==c ( 0 , 0 ) ) ) { r e t u r n ( 0 ) } e l s e { r e t u r n ( s q r t ( sum( ( f ( v)− e v l a * v ) ˆ2 ) ) / s q r t ( sum ( v ˆ2 ) ) ) } } a p r o x e i g e n R 2 <− f u n c t i o n (FUN, range = 1 , N, N0, N1 , tau , s t e p s = 1){ # FUN: t h e o r i g i n a l f u n c t i o n # range : work on t he r e c t a n g e [− range , range ] ˆ2 # N: number o f p o i n t s t o s t a r t t h e f i r s t time # N0 : number o f p o i n t s ( o f th e N) t o s e l e c t # N1 : number o f p o i n t s t o ge n e r a t e # o f each o f N0 t o s e l e c t th e b e s t # tau : t he f a c t o r o f t he v a r i a n c e f o r th e g a u s s i a n # s t e p s : s t e p s t o r e p e a t # output : a l i s t with a l l t h e i n f o r m a t i o n r e s <− l i s t ( S= l i s t ( ) ) # STEP 0 d f <− data .frame ( x = r u n i f (N, min=−range , max=range ) , y = r u n i f (N, min=−range , max=range ) ) %>% mutate ( e v e r r=map2 dbl ( x , y , f u n c t i o n ( x , y ){ e v e r r o r ( c ( x , y ) , FUN) } ) ) %>% a r r a n g e ( e v e r r ) r e s $ S 0 0 <− d f r e s $ p l o t s t e p 0 <− g g p l o t ( ) + geom point ( data=df , mapping=a e s ( x , y ) , c o l o r ="gray " ) + geom point ( data=d f [ 1 : N0 , ] , mapping=a e s ( x , y ) , c o l o r ="b l a c k " ) + c o o r d c a r t e s i a n ( xlim=c(−range , range ) , ylim=c(−range , range ) ) d f <− d f [ 1 : N0 , ] r e s $ S 0 <− d f # STEPS >= 1 f o r ( s i n 1 : s t e p s ){ f o r ( i i n 1 : N0){ x = c ( df $ x [ i ] , rnorm ( N1 , mean=df$ x [ i ] , sd=tau * d f $ e v e r r [ i ] ) ) y = c ( df $ y [ i ] , rnorm ( N1 , mean=df$ y [ i ] , sd=tau * d f $ e v e r r [ i ] ) ) e v e r r <− map2 dbl ( x , y , f u n c t i o n ( x , y ){ e v e r r o r ( c ( x , y ) ,FUN) } ) i nd <− e v e r r %>% which .min () i f ( x [ i nd ]> range ){ x [ i n d ]= range } i f ( x [ i nd]<−range ){ x [ i n d]=−range } i f ( y [ i nd ]> range ){ y [ i n d ]= range } i f ( y [ i nd]<−range ){ y [ i n d]=−range } d f [ i , ] <− c ( x [ in d ] ,y [ in d ] , e v e r r [ in d ] ) %>% data .frame ( ) %>% t ( ) } r e s $ S [ [ s ] ] <− d f } r e s $ e i g e n v e c <− d f %>% s e l e c t ( x , y ) %>% mutate ( lambda = map2 dbl ( x , y , ˜ev lambda ( c ( .x , .y ) , FUN) ) ) r e s $ m e a n e r r <− mean ( d f [ , ' e v e r r ' ] ) r e s $ p l o t s t e p e n d <− g g p l o t ( ) + geom point ( data=df , mapping=a e s ( x , y ) , c o l o r ="b l a c k " ) + c o o r d c a r t e s i a n ( xlim=c(−range , range ) , ylim=c(−range , range ) ) r e t u r n ( r e s ) } ### EXAMPLE −−−−−− REF −−−−−− R0 <− f u n c t i o n ( v ) { x = v [ 1 ] y = v [ 2 ] r e t u r n ( c (8 * x , 4 * x ˆ2 + 4 * x * y +y ˆ2 − 2 * x − 1/5 * s q r t ( abs ( x+y ) ) ) ) } ev <− a p r o x e i g e n R 2 (FUN = R0 , range = 5 , N = 5 0 0 , N0 = 1 0 0 , N1 = 1 0 , tau = 5 , s t e p s = 10 ) e v $ p l o t s t e p e n d Funding: The first author was supported by a contract of the Programa de Ayudas de Investigación y Desarrollo (PAID-01-21), Universitat Politècnica de València.