Support Vector Regression via Mathematica

B. Paláncz
 
 Department of Photogrammetry and Geoinformatics
 Budapest University of Technology  and Economics, 
 H-1111 Budapest, Hungary, Müegyetem rkp. 3.
 
 e-mail:  palancz@epito.bme.hu

Abstract

In this tutorial type paper a Mathematica function for Support Vector Regression has been developed. Summarizing the main definitions and theorems of SVR, the detailed implementation steps of this function are presented and its application is illustrated by solving three 2D function  approximation test problems, employing a stronger regularized universal Fourier  and a wavelet kernel. In addition a real world regression problem, forecasting of the peak of  floodwave,  is also solved. The numeric - symbolic results show how easily and effectively Mathematica can be used for solving  SVR problems.

Introduction

In geoinformatics, it frequently happens, that data collections are irregular and far from substantial for function approximation. This fact leads to overfitting, a poor data generalization, which means a good accuracy for training samples but bad fitting of the test data points[1]. To solve this problem one may use ridge regression or other regularization techniques modifying the object function in some ways [2, 3, 4]. However, in case of finding nonlinear mapping for input/output data of high dimensions, these methods usually result computational problem of  NP complexity.

In the last few years, there have been very significant developments in the theoretical understanding of a relatively new family of algorithms, they present a series of useful features for classification as well as generalization of datasets [5]. Support Vector Machines (SVM) algorithms combine the simplicity and computational efficiency of linear algorithms, such as the perceptron algorithm or ridge regression, with the flexibility of nonlinear systems, like neural networks,  and rigour of statistical approaches,  as regularization methods in multivariate statistics [6]. As a result of the special way they represent functions, these algorithms typically reduce the learning step to a convex optimization problem, that can always be solved in polynomial time, avoiding the problem of local minima typical of neural networks, decision trees and other nonlinear approaches [7].

Their foundation in the principles of  statistical learning theory make them remarkably resistant to overfitting especially in regimes where other methods are effected by the curse of dimensionality. It is for this reason that they have become popular in bioinformatics and text analysis. Another important feature for applications is that they can naturally accept input data that are not in the form of vectors, such as strings, trees and images [8, 9].

In this tutorial type article, it will be shown, that how easy to use this efficient technique for function approximation via Mathematica. The steps of implementation of  support vector regression (SVR) are discussed  and the application of this Mathematica function is illustrated by solving 2D approximation  test problems with different type of kernels.

1. Kernels

1.1 Definition of the kernel

Let X and H denote a linear vector space and a Hilbert space with scalar product ⟨*|*⟩. A continuous symmetric function K: X × X -> ℝ  is said to be a kernel on X if there exists  a map Φ :  X -> H  with

K (x, z) =  ⟨Φ (x), Φ(z)⟩

for which holds

⟨ K x, z⟩ = ⟨ x, K z⟩

for all x, z ∈ X. We call Φ a feature map and H a feature space of K.  For example a possible mapping is the following,

{(Subscript[x, 1]), (Subscript[x, 2])}∈(ℝ^( 2)) -> Φ ((Subscript[x, 1]), (Subscript[x, 2]))={((Subscript[x, 1])^2), ((Subscript[x, 2])^2), (Subscript[x, 1]) (Subscript[x, 2])}∈ H

Note that  in general neither the H nor the function  Φ is unique. However, for the kernels satisfying certain condition given below there exists a unique Hilbert space of  functions, named the Reproducing Kernel Hilbert Space (RKHS).

1.2 Mercer's condition

Let X  be a compact subset of  (ℝ^( n)). Suppose K is a continuous symmetric function for which

(∫_(X × X))K(x, z)f(x) f(z)dx𝕕z

is strictly positive, that is for all  f  from  the normed space of the continuous functions on X with infinite norm \(TraditionalForm\`\(\(∥\)\(\(*\)\( ∥ \_∞\)\)\)\) , namely  f ∈ C [X ]. Then for  K  exists an  RKHS.

1.3 Universal kernel

A function f : X -> ℝ is induced by the kernel K if there exists an element w ∈ H such that  f = ⟨w, Φ(*)⟩.  For example,

f ((Subscript[x, 1]), (Subscript[x, 2]))=⟨{(Subscript[w, 1]),  (Subscript[w, 2]),  (Subscript[w, 3])},Φ ((Subscript[x, 1]), (Subscript[x, 2]))⟩={(Subscript[w, 1]),  (Subscript[w, 2]),  (Subscript[w, 3])} . {((Subscript[x, 1])^2), ((Subscript[x, 2])^2), (Subscript[x, 1]) (Subscript[x, 2])} = 
          (Subscript[w, 1]) ((Subscript[x, 1])^2) +  (Subscript[w, 2]) ((Subscript[x, 2])^2) +  (Subscript[w, 3]) (Subscript[x, 1]) 
              (Subscript[x, 2])

A continuous kernel is called universal if the set of all induced functions is dense in the space of continuous function on X,  C [X ], which means that for all  g ∈ C [X ] and ϵ > 0 there exists a function f induced by K with \(TraditionalForm\`\((\(∥\)\(f\  - \ g\)\(∥\))\)\_∞\) ≤ ϵ.

If  K can be represented by a Taylor series, a necessary and sufficient condition for universality of K can be expressed with the help of the nonvanishing  Taylor coefficients. A kernel is universal if and only if its RKHS  is dense in C [X ] [10].

1.4 Some typical kernels

Here, we  mention  some typical kernels and display them  in case of X⊂(ℝ^( 1)).

1.4.1 The Gaussian RBF universal kernel

The Gaussian Radial Basis Function universal kernel with β>0 and all compact X⊂(ℝ^( n)),

K (x, z) = exp(-β \(TraditionalForm\`\((\(∥\)\(x\  - z\)\(∥\))\)\^2\))

with β = 0.1

β = 0.1; 
    K[x_, z_] := Exp[-β (x - z) (x - z)];

Plot3D[K[x, z], {x, -10, 10}, {z, -10, 10}, 
        ViewPoint -> {3.849,  6.909,  4.295}, 
        PlotPoints -> {30, 30}];

1.4.2  Polynomial kernel

Polynomial kernel of degree d > 0,

K (x, z) = \(TraditionalForm\`\((c\  + \ ⟨\ x, \ z⟩)\)\^d\)

with

c  = 1; d = 4; K[x_, z_] := (((c + x z))^d);

Plot3D[K[x, z], {x, 0, 10}, {z, -10, 10}, 
        ViewPoint -> {-7.654,  -1.091,  4.608}, 
        PlotPoints -> {30, 30}];

1.4.3  Vovk's real infinity polynomial universal kernel

Vovk's real infinity polynomial universal kernel, for d > 0 and all compact X⊂{ x ∈(ℝ^( n)) : \(\(∥\)\(x\)\(∥\)\(\ \)\(<\)\(\ \)\(1\)\)},

K (x, z) = \(TraditionalForm\`\((1\  - \ ⟨\ x, \ z⟩)\)\^\(-d\)\)

with

d = 5; K[x_, z_] := (((1 - x z))^(-d));

Plot3D[K[x, z], {x, -0.95, 0.95}, {z, -0.95, 0.95}, 
        ViewPoint -> {-5.906,  4.688,  4.913}, 
        PlotPoints -> {30, 30}];

1.4.4  Stronger regularized universal Fourier kernel

The stronger regularized universal Fourier kernel with 0 < q < 1 and all compact X ⊂ (([0,  2 π ))^( n))

K(x, z)=\(∏\+\(i = 1\)\%n\( 1 - q\^2\)\/\(2\ \((1 - 2\ q\ cos[x\_i - z\_i] + q\^2)\)\)\)

with

q  =  0.5;  
    K[x_, z_] := ((1 - (q^2))/(2 (1 - 2 q Cos[x - z] + (q^2))));

Plot3D[K[x, z], {x, 0, 2 π}, {z, 0, 2 π}, 
        ViewPoint -> {3.849,  6.909,  4.295}, 
        PlotPoints -> {30, 30}];

1.4.5  Wavelet kernel

Wavelet kernel with  a ∈ (ℝ^( 1)) and all compact  X⊂(ℝ^( n)),

K(x, z)=(∏_(i = 1)^n)(cos[
              1.75 (((Subscript[x, i]) - (Subscript[z, i]))/a)] exp[-(( ((((Subscript[x, i]) - 
                            (Subscript[z, i])))^2))/(2 (a^2)))])

with

a = 4;  
    K[x_, z_] := 
      Cos[1.75 ((x - z)/a)] Exp[-((((x - z))^2)/(2 (a^2)))];

Plot3D[K[x, z], {x, -10, 10}, {z, -10, 10}, 
        ViewPoint -> {4.172,  5.346,  5.917}, 
        PlotPoints -> {30, 30}];

1.5 Making kernels from kernels

The following proposition can be viewed as shown that kernels satisfy a number of closure properties, allowing us to create more complicated kernels from simple building blocks. 

Let (Subscript[K, 1])and (Subscript[K, 2]) be kernels over X × X, x, z ∈ X ⊆ (ℝ^n), a ∈ (SuperPlus[ℝ]) and f a real- valued function on X, in addition exist an other function Φ : X -> (ℝ^m), with (Subscript[K, 3])over (ℝ^m)× (ℝ^m),then the following functions are kernels.

K (x, z) = (Subscript[(K ), 1])(x, z) +(Subscript[K, 2]) (x, z)

K (x, z) = a (Subscript[K, 1]) (x, z)

K (x, z) = (Subscript[(K ), 1])(x, z) (Subscript[K, 2]) (x, z)

K (x, z) = f(x) f(z)

K (x, z) = (Subscript[K, 3])(Φ(x),Φ(z))

In addition, if p(x) a polynomial with positive coefficient, then

K (x, z) = p((Subscript[K, 1]) (x, z) )

K(x, z) = exp(K(x, z))

2. Support vector machine for regression

2.1 ϵ - Insensitive loss function

The problem of regression is that of finding a function which approximates mapping from an input domain to the real numbers based on a training sample. We refer to the difference between the hypothesis output and its training value as the residual of the output., an indication of the accuracy of the fit at this point. We must decide how to measure the importance of this accuracy, as small residuals may be inevitable while we wish to avoid large ones. The loss function determines this measure. Each choice of loss function will result in a different overall strategy for performing regression. For example least square regression uses the sum of the squares of the residuals.
Although several different approaches are possible, we will provide an analysis for generalization of regression by introducing a threshold test accuracy θ, beyond which we consider a mistake to have been made. We therefore aim to provide a bound on the probability that a randomly drawn test point will have accuracy less than θ. If we access the training set performance using the same θ, we are effectively using the real - valued regressors as classifiers and the worst case lower bounds apply. What we must do in order to make use of dimension free bounds is to allow a margin in the regression accuracy that corresponds to the margin of a classifier. We will use the symbol γ to denote the margin, which measures the amount by which the training and test set accuracy can differ. It should be emphasized that we are therefore using a different loss function during training and testing, where γ measures the discrepancy between the two losses, implying that training point counts as mistake if its accuracy less than θ - γ.  One way of visualizing this method of assessing performance is to consider a band of size ± (θ - γ) around the hypothesis function  any training points lying outside this band are considered to be training mistakes. Test points count as mistakes only if they lie outside the wider band ± θ.

The linear ϵ-insensitive loss function (L^ϵ)(x, y, f) is defined by

(L^ϵ)(x, y, f)=\(\((\(|\)\(y\  - \ f \((x)\)\)\(|\))\)\_ϵ\  = \ 
          max\ \((0, \(\(|\)\(y\  - \ 
                  f \((x)\)\)\(|\)\(\ \)\(-\ ϵ\)\))\)\)

where f is a real - valued function on a domain X, x ∈ X and y ∈ ℝ. Similarly the quadratic ϵ-insensitive loss is given by

((Subscript[L, 2])^ϵ)(x, y, f)=\(\((\(|\)\(y\  - \ f \((x)\)\)\(|\))\)\_ϵ\^2\)

2.2 Support Vector Regression

SVR uses an admissible kernel, which satisfies the Mercer 's condition to map the data in input space to a high-dimensional feature space in which we can process a regression problem in linear form. Let x ∈ (ℝ^n) and y ∈ ℝ, where (ℝ^n) represents input space. By some nonlinear mapping Φ, the vector x is mapped into a feature space in which a linear regressor function is defined,

y = f (x,  w) = ⟨w, Φ(x)⟩+ b

We seek to estimate this  f  function based on independent uniformly distributed data {{(Subscript[x, 1]), (Subscript[y, 1])},...,{(Subscript[x, m]), (Subscript[y, m])}}, by finding w which minimizing the quadratic ϵ - insensitive losses, with ϵ = θ - γ, namely the following function should be minimize [ 11]

c (∑_(i = 1)^m) ((Subscript[L, 2])^ϵ)((Subscript[x, i]), (Subscript[y, i]), 
                f)+ \(\(1\/2\) \((\(∥\)\(w\)\(∥\))\)\^2\)-> min

where w is weight vector and c is a constant parameter.
Considering dual representation of a linear regressor in (19),  f (x) can be expressed as [12],

\(f \((x)\) = ∑\+\(i = 1\)\%m β\_i\ \(y\_i\) ⟨Φ \((x\_i)\), Φ \((x)\)⟩ + b\)

which means that the regressor can be expressed as a linear combination of the training points. Consequently using an admissible kernel, a kernel satisfying the Mercer' s condition, we get

f (x) = (∑_(i = 1)^m) (Subscript[β, i]) (Subscript[y, i]) 
              K ((Subscript[x, i]), x) + 
          b = (∑_(i = 1)^m) (Subscript[α, i]) K ((Subscript[x, i]), x) + 
          b

By using Lagrange multiplier techniques, the minimization problem of (20) leads to the following dual optimization problem [ 12],

\(maximize\ \ \ W \((α)\)\  = ∑\+\(i = 1\)\%m y\_i\ α\_i - \ ϵ\ ∑\+\(i = 1\)\%m ∣ α\_i ∣ \(-\(1\/2\)\) \(∑\+\(i, j = 1\)\%m\( α\_i\) α\_j\ \((K \((x\_i, 
                      x\_j)\) + \(1\/c\) δ\_ij)\)\)\)
subject to (∑_(i = 1)^m) (Subscript[α, i])  =  
      0

Let

\(f \((x)\) = ∑\+\(i = 1\)\%m\( α\_i\^*\)\ K \((x\_i, 
              x)\) + \(b\^*\)\)

where \(TraditionalForm\`\(α\^*\)\) is the solution of the quadratic optimization problem and  \(TraditionalForm\`\(b\^*\)\) is chosen so that f((Subscript[x, i])) = (Subscript[y, i])- ϵ - \(TraditionalForm\`\(α\_i\^*\)\/c\)  for any \(TraditionalForm\`\(α\_i\^*\)\) > 0.  

For samples are inside the ϵ - tube, {(Subscript[x, i]): ∣ f((Subscript[x, i])) - (Subscript[y, i])∣ < ϵ, the corresponding \(TraditionalForm\`\(α\_i\^*\)\)is zero. It means, we do not need these samples to describe the weight vector w. Consequently

\(f \((x)\) = ∑\+\(i ∈ \ 𝒮𝒱\)\(α\_i\^*\)\ K \((x\_i, x)\) + \(b\^*\)\)

where

𝒮𝒱={i : ∣f((Subscript[x, i]))-(Subscript[y, i])∣ ≥ \(ϵ\ \ OR\ \ \(α\_i\^*\)\)> 0}

These (Subscript[x, i]) sample vectors, {(Subscript[x, i]) : i ϵ 𝒮𝒱 } ,  that come with nonvanishing coefficients (SuperStar[(Subscript[α, i])]) are called support vectors.

2.3 The main features of SVR

We summarize briefly the main features of the SVR,

- application of nonlinear mapping of data from input space into a feature space, in which linear regressor (machine or learning machine) is used. The dual representation of linear regressor leads to the employment of kernel function.

-  quadratic ϵ - insensitive loss function is used as objective function with regularization term for estimation of the weight vector in the kernel representation of the regressor function f, ensuring ϵ accuracy, and leading the a part of samples, called support vector, only which influence the weight vector.

3. Implementation of SVR in Mathematica

3.1 Steps of implementation

The dual optimization problem can be solved conveniently using Mathematica. In this section, the steps of  the implementation of  SVR algorithm are shown by solving a function approximation problem [ 13 ]. The function to be approximated is,

z[{x_, y_}] := ((x^2) - (y^2)) Sin[0.5 x]

Let us display it

p1 = 
        Plot3D[z[{x, y}], {x, -10, 10}, {y, -10, 10}, 
          PlotPoints -> {30, 30}];

A training set can be created in form {{(Subscript[x, i]), (Subscript[y, i])}, z ((Subscript[x, i]), (Subscript[y, i]))}

data = 
        Table[{{x, y}, z[{x, y}]}, {x, -10, 10, 2.5}, {y, -10, 10, 
            2.5}];

Let us separate the dependent and independent variables,

{xym, zm} = Transpose[Flatten[data, 1]];

Wavelet kernel with parameter a = 4,  in case of  dimension  n = 2

n = 2; a = 4. ;

K[u_, 
        v_] := (∏_(i = 1)^n)(Cos[
              1.75 ((u[[i]] - v[[i]])/a)] Exp[-(( (((u[[i]] - v[[i]]))^2))/(2 (a^2)))])

The number of the data pairs in the training set, m  is

m = Length[xym]

Create the objective function W (α) to be maximized, with parameters

ϵ = 0.0025; c = 200. ;

First, we prepare an matrix M

M = (Table[
              N[K[xym[[i]], xym[[j]]]], {i, 1, m}, {j, 1, 
                m}] + (1/c) IdentityMatrix[m]);

then the objective function can be expressed as,

W = (∑_(i = 1)^m) (Subscript[α, i]) 
              zm[[i]] - ϵ (∑_(i = 1)^m) 
                Abs[(Subscript[α, i])] - 
          (1/2) (∑_(i = 1)^m)(∑_(j = 1)^m)((Subscript[α, i]) (Subscript[α, j]) M[[i, j]]);

The constrains for the unknown variables are

g = 
        Apply[And, 
          Join[Table[-c < (Subscript[α, i]) ≤ c, {i, 1, 
                m}], {((∑_(i = 1)^m) (Subscript[α, i]) == 
                  0)}]];

However the maximization problem (23) is a convex quadratic problem, from practical reasons to maximize the objective function the built in function NMaximize is applied. NMaximize implements several algorithms for finding constrained global optima. The methods are flexible enough to cope with functions that are not differentiable or continuous, and are not easily trapped by local optima. Possible settings for the Method option include "NelderMead", "DifferentialEvolution", "SimulatedAnnealing" and "RandomSearch". 
Here we use  DifferentialEvolution, which is a genetic algorithm that maintains a population of specimens, (Subscript[x, 1]), …, (Subscript[x, n]), represented as vectors of real numbers (“genes”). Every iteration, each (Subscript[x, i]) chooses random integers a, b, and c and constructs the mate (Subscript[y, i]) = (Subscript[x, i]) + γ((Subscript[x, a]) + ((Subscript[x, b]) - (Subscript[x, c]))), where γ is the value of ScalingFactor. Then (Subscript[x, i]) is mated with (Subscript[y, i]) according to the value of CrossProbability, giving us the child (Subscript[z, i]). At this point (Subscript[x, i]) competes against (Subscript[z, i]) for the position of (Subscript[x, i]) in the population. The default value of SearchPoints is Automatic, which is Min[10*d, 50], where d is the number of variables.

We need the list of unknown variables α,

vars = Table[(Subscript[α, i]), {i, 1, m}];

Then the solution of the maximization problem is,

sol = 
      NMaximize[{W, g}, vars, 
        Method -> DifferentialEvolution]

The consistency of this solution can be checked by computing values of  b  for every data points. Theoretically, these values should be same for any data points, however, this is only approximately true.

bdata = 
      Table[(zm[[j]] - (∑_(i = 1)^m) (Subscript[α, i]) 
                  K[xym[[j]], 
                    xym[[i]]] - ϵ - ((Subscript[α, j])/c)) /. 
          sol[[2]], {j, 1, m}]

The value of b can be chosen as the average of these values

b = Apply[Plus, bdata]/m

Then the regressor function is,

f[{x_, 
            y_}] := ((∑_(i = 1)^m) (Subscript[α, i]) 
                  K[{x, y}, xym[[i]]] + b) /. 
          sol[[2]];

The result in symbolic form is,

Short[f[{x, y}], 5]

This function can be compared with the original one,

p2 = 
        Plot3D[f[{x, y}], {x, -10, 10}, {y, -10, 10}, 
          PlotPoints -> {30, 30}, 
          DisplayFunction -> Identity];

Show[GraphicsArray[{p1, p2}], 
        DisplayFunction -> $DisplayFunction];

These figures indicate a very good approximation.

3.2 Mathematica module for SVR

These steps can be collected in a module, where the vector xm contains the input vectors and the  vector ym contains the corresponding scalar output values,

SupportVectorRegression[{xm_, ym_}, K_, ϵ_, c_] := 
        Module[
{m, n, M, i, j, W, g, vars, sol, 
            b}, 
m = Length[ym]; 
          n = Length[xm[[1]]]; 
M = 
            Table[K[xm[[i]], xm[[j]]], {i, 1, m}, {j, 1, 
                  m}] + (1/c) 
                IdentityMatrix[
                  m]; 
W = (∑_(i = 1)^m) (Subscript[α, i]) ym[[i]] - ϵ (∑_(i = 1)^m) 
                    Abs[(Subscript[α, i])] - 
              (1/2) (∑_(i = 1)^m)(∑_(j = 1)^m)((Subscript[α, i]) (Subscript[α, j]) M[[i, 
                            j]]); 
g = 
            Apply[And, 
              Join[Table[-c < (Subscript[α, i]) ≤ c, {i, 1, 
                    m}], {((∑_(i = 1)^m) (Subscript[α, i]) == 
                      0)}]]; 
vars = 
            Table[(Subscript[α, i]), {i, 1, 
                m}]; 
sol = NMaximize[{W, g}, vars, 
                Method -> 
                  DifferentialEvolution][[2]]; 
b = (1/m) Apply[Plus, 
                Table[(ym[[j]] - (∑_(i = 1)^m) (Subscript[α, i]) 
                            K[xm[[i]], 
                              xm[[j]]] - ϵ - ((Subscript[α, j])/c)) /. 
                    sol, {j, 1, 
                    m}]]; 
{((∑_(i = 1)^m) (Subscript[α, i]) K[xm[[i]], Table[(Subscript[x, j]), {j, 1, n}]] + b) /. sol, 
            vars /. sol}];

The module outputs are the regression function in symbolic form and the values of α 's. To test our module,  let us solve another problem, [14]. The following function should be approximated in  [-5, 5] ×  [-5, 5],

z[{x_, y_}] := Sin[(Sqrt[(x^2) + (y^2)])]/(Sqrt[(x^2) + (y^2)])

using a stronger regularized universal Fourier kernel with parameter, q = 0.25

q = 0.25;

K[u_, 
        v_] := (∏_(i = 1)^n)(( 1 - (q^2))/(2 (1 - 2 q Cos[u[[i]] - v[[i]]] + (q^2))))

Because this kernel is valid for x, y ∈  X ⊂ (([0,  2 π ))^( n)), therefore the variables of  x, y should be scaled,

z[{x_, y_}] := 
      Sin[(Sqrt[(((((5  x)/π) - 5))^2) + (((((5  y)/π) - 5))^2)])]/(Sqrt[(((((5  x)/π) - 5))^2) + (((((5  y)/π) - 5))^2)])

Let us display this function,

p1 = 
        Plot3D[z[{x, y}], {x, 0, 2 π}, {y, 0, 2 π}, 
          PlotPoints -> {30, 30}];

Generating data,

data = 
        Table[{{x, y}, z[{x, y}]}, {x, 0, 2 π, ((2  π)/9)}, {y, 0, 
              2 π, ((2  π)/9)}] // N;

{xym, zm} = Transpose[Flatten[data, 1]];

then calling the SVR module,

F = 
        SupportVectorRegression[{xym, zm}, K, ϵ, c];

Symbolic form of the result,

Short[F[[1]], 18]

The surface approximation can be qualified  by comparing it with the original one,

p2 = 
        Plot3D[F[[1]], {(Subscript[x, 1]), 0, 2 π}, {(Subscript[x, 2]), 0, 2 π}, 
          PlotPoints -> {30, 30}, 
          DisplayFunction -> Identity];

Show[GraphicsArray[{p1, p2}], 
        DisplayFunction -> $DisplayFunction];

The third test problem is a well-known benchmark, the approximation of the following function from noisy training samples [15],

z[{x_, y_}] := 
        1.9 (1.35 + (E^x) Sin[
                  13 (((x - 0.6) )^2)] (E^(-y)) Sin[7 y]);

Let us display it

p1 = 
        Plot3D[z[{x, y}], {x, 0, 1}, {y, 0, 1}, PlotPoints -> {30, 30}, 
          Shading -> False, 
          ViewPoint -> {-1.097,  -1.601,  0.963}];

Noise of  normal distribution N(μ, σ), with μ = 0 and σ = 0.15 was added to the function.

<<Statistics`ContinuousDistributions`

zn[{x_, y_}] := 
      z[{x, y}] + Random[NormalDistribution[0. , 0.15]]

A training set can be created in form {{(Subscript[x, i]), (Subscript[y, i]), z ((Subscript[x, i]), (Subscript[y, i]))}

data = 
        Flatten[Table[
            N[{i, j, zn[{i, j}]}], {i, 0, 1, (1/9)}, {j, 0, 1, (1/9)}], 
          1];

<< Graphics`Graphics3D`

p2 = 
        ScatterPlot3D[data, BoxRatios -> {1, 1, 1}, 
          PlotStyle -> PointSize[0.019], 
          ViewPoint -> {-1.014,  -1.580,  1.237}, 
          DisplayFunction -> Identity];

Show[{p1, p2}, 
        DisplayFunction -> $DisplayFunction];

Let us separate the dependent and independent variables,

{xym, zm} = 
        Transpose[
          Map[{{#[[1]], #[[2]]}, #[[3]]} &, data]];

Wavelet kernel with parameter a = 0.23,  in case of  dimension  n = 2

n = 2; a = 0.23;

K[u_, 
        v_] := (∏_(i = 1)^n)(Cos[
              1.75 ((u[[i]] - v[[i]])/a)] Exp[-(( (((u[[i]] - v[[i]]))^2))/(2 (a^2)))])

The number of the data pairs in the training set, m  is

m = Length[xym]

Create the objective function W (α) to be maximized, with parameters

ϵ = 0.0025; c = 200. ;

F = 
        SupportVectorRegression[{xym, zm}, K, ϵ, c];

Symbolic form of the result,

Short[F[[1]], 18]

p3 = 
        Plot3D[F[[1]], {(Subscript[x, 1]), 0, 1}, {(Subscript[x, 2]), 0, 1}, 
          PlotPoints -> {30, 30}, 
          ViewPoint -> {-1.014,  -1.580,  1.237}, 
          DisplayFunction -> Identity];

Show[GraphicsArray[{p1, p3}], 
        DisplayFunction -> $DisplayFunction];

These figures indicate a  good approximation.

4. Engineering application

4.1 Problem definition

We should like to develop a statistical model for forecasting the peak of the  floodwave of  the Danube at Budapest. The following 24 data triplets registered from the year of 1896 up to 1955 available in  [16],

({{x, y, z}, {rainfall

           
          in mm, water - 
            level at the beginning of raining 

          in cm, peak of water - level 

          in cm}})

The measured data are,

data = {{58, 405, 590}, {52, 450, 660}, {133, 350, 780}, {179, 285, 
            770},  {98, 330, 710}, {72, 400, 640}, {72, 550, 670}, {43, 480, 
            520}, {62, 450, 660}, {67, 610, 690}, {64, 380, 500}, {33, 460, 
            460}, {57, 425, 610}, {62, 560, 710}, {54, 420, 620}, {48, 620, 
            660}, {86, 390, 620}, {74, 350, 590}, {95, 570, 
            740}, 
{44, 710, 730}, {77, 580, 720}, {46, 
            700, 640}, {123, 560, 805}, {62, 430, 673}};

The model z = f (x, y) should generalizes these data and ensures forecasting z from measured x and y values. Now we have a real regression problem, because the data are not on a smooth surface,

<< Graphics`Graphics3D`

\(\(p0 = 
        ScatterPlot3D[data, BoxRatios -> {1, 1, 1}, 
          PlotStyle -> PointSize[0.03], 
          AxesLabel -> {"Rainfall ", "Starting Level            ", "Peak Level"}];\)\)

A general model means that the set of input/output relationships, derived from the training set, apply equally well to new sets of data from the same problem not included in the training set. The main goal is thus the generalization to new data of the relationships learned on the training set.
In order to test the generality of our model to be developed, we divide the data into a training set and a test set. The number of data is

ndata = Length[data]

Let consider every fourth data as test data

dataV = Table[data[[i]], {i, 1, ndata, 4}]

The remaining elements are for training

dataT = Complement[data, dataV]

Let us display the training and test set

\(\(p1 = 
        ScatterPlot3D[dataT, BoxRatios -> {1, 1, 1}, 
          PlotStyle -> {RGBColor[0, 0, 1], PointSize[0.03]}, 
          AxesLabel -> {"Rainfall ", "Starting Level            ", "Peak Level"}, DisplayFunction -> Identity];\)\)

\(\(p2 = 
        ScatterPlot3D[dataV, BoxRatios -> {1, 1, 1}, 
          PlotStyle -> {RGBColor[1, 0, 0], PointSize[0.03]}, 
          AxesLabel -> {"Rainfall ", "Starting Level            ", "Peak Level"}, DisplayFunction -> Identity];\)\)

Show[GraphicsArray[{p1, p2}], 
        DisplayFunction -> $DisplayFunction];

Preparation of   these sets for learning and testing

xym = Map[Drop[#, -1] &, data];

zm = Map[Take[#, {3, 3}] &, data] // Flatten;

xymT = Map[Drop[#, -1] &, dataT];

zmT = Map[Take[#, {3, 3}] &, dataT] // Flatten;

4.2 Solution with SVR

Employing wavelet kernel with parameter a = 0.3

n = 2; a = 0.3;

K[u_, 
        v_] := (∏_(i = 1)^n)(Cos[
              1.75 ((u[[i]] - v[[i]])/a)] Exp[-(( (((u[[i]] - v[[i]]))^2))/(2 (a^2)))])

and with the following parameters for  SVR

ϵ = 0.025; c = 200. ;

F = 
        SupportVectorRegression[{xymT, zmT}, K, ϵ, c];

The regressor function is,

Short[F[[1]], 15]

f[{x_, y_}] = F[[1]] /. {(Subscript[x, 1]) -> x, (Subscript[x, 2]) -> y};

Let us display the relative error of the approximation on the whole data set

<< Graphics`Graphics`

BarChart[Abs[zm - Map[f[#] &, xym]]/zm 100, 
        Ticks -> {{1, 5, 9, 13, 17, 21, 25}, {0. , 10, 20}}, 
        PlotRange -> {0, 20},  AspectRatio -> 0.5];

The figure shows, that on most of the test data, which were not included in the training, the error is considerably high, while on the elements of the training set the error is negligible. The maximum error is

Max[Abs[zm - Map[f[#] &, xym]]/zm 100]

and its standard deviation  is

StandardDeviation[Abs[zm - Map[f[#] &, xym]]/zm 100]

Let us see the values of α 's,

F[[2]]

As we know, the input vector (Subscript[xymT, i])  is a support vector, if the corresponding  (Subscript[α, i])  ≠  
        0.

\(We\ shall\ consider\ α\_i = 0, \ 
    if\ its\ absolute\ value\ is\ smaller\ than\ \(\(10\^\(-3\)\)\(.\)\(\ \)\)\)

Then the support vectors can be  selected in the following way,

supportvectors = 
      Extract[xymT, 
        Position[F[[2]], _?(Abs[#] > (10^(-3)) &)]]

Now, all of the training vector are support vector. It is easy to visualize them,

<< Graphics`PlotField`

ListPlotVectorField[Map[{{0, 0}, #} &, supportvectors], 
        AspectRatio -> 0.5];

In order to eliminate the approximation error on the test set, to improve the generalization of our model, the parameter ϵ should be increased. Now let us carry out the computation with ϵ = 10.

ϵ = 10. ;

F = 
        SupportVectorRegression[{xymT, zmT}, K, ϵ, c];

The regressor function is,

Short[F[[1]], 15]

f[{x_, y_}] = F[[1]] /. {(Subscript[x, 1]) -> x, (Subscript[x, 2]) -> y};

Let us display again the relative error of the approximation on the whole data set

BarChart[Abs[zm - Map[f[#] &, xym]]/zm 100, 
        Ticks -> {{1, 5, 9, 13, 17, 21, 25}, {0. , 10, 20, 30}}, 
        PlotRange -> {0, 20},  AspectRatio -> 0.5];

This figure shows, that although the maximum of the approximation error  decreased just a little,  the error distribution is tending to a more uniform one, which indicates higher statistical confidence of the model.

Max[Abs[zm - Map[f[#] &, xym]]/zm 100]

and its standard deviation  is

StandardDeviation[Abs[zm - Map[f[#] &, xym]]/zm 100]

Let us see the values of α 's,

F[[2]]

Now, we have two small α 's indicating less support vector as before

supportvectors = 
        Extract[xymT, 
          Position[F[[2]], _?(Abs[#] > (10^(-3)) &)]];

ListPlotVectorField[Map[{{0, 0}, #} &, supportvectors], 
        AspectRatio -> 0.5];

5. Conclusions

Support vector regression method has been implemented in Mathematica code as a simple function and tested on three different problems. The solutions of them, especially the third one, a real world civil engineering problem,  forecasting of the peak of a floodwave, clearly demonstrated the generalization ability of SVR as well as its robustness.
Additional computations were also carried out with polynomial regression with different orders (2, 3  and 4) and RBF neural network with different number of neurons (3, 5 and 10). In case of these methods the maximum of the relative error was always considerably higher than that of SVR and the phenomena of overfitting always took place with increasing number of  model parameters.

References

1. Tzafestas S.G., Dalianis P.J. and Anthopoulos G. On the overtraining phenomenon of back propagation neural networks, Mathematics and Computers in Simulation,  Vol. 40. pp. 507 - 521, 1996.

2. Neumaier A. Solving ill - conditioned and singular linear systems: a tutorial on regularization, WWW:http://solo.cma.univie.ac.at/~neum/

3. Friess T.T. and Harrison R.F. A kernel based Adaline or function approximation, Intelligent Data Analysis, Vol.3. pp. 307 - 313, 1999.

4. Burger M and Neubauer A. Analysis of Tikhonov regularization for function approximation by neural network, Neural Networks, Vol.16. pp.79 - 90, 2003.

5. Berthold M. and Hand D.J. (Eds.) "Intelligent Data Analysis, An Introduction", Springer, 2003.
 
6. Burgers C.J.C. A tutorial on support vector machines for pattern recognition, Data Mining and Knowledge Discovery 2, pp. 121 - 167, 1998.

7. Hearst M. A. Suport Vector Machines, IEEE Intelligent Systems, pp. 18 -28, July/August, 1998.

8. Kim K. I., Jung K., Park S.H. and Kim H.J. Support Vector Machines for Texture Classification, IEEE Trans. Pattern Analysis and Machine Intelligence, Vol. 24. No.11. pp. 1542 - 1550, Nov. 2002.
 
9. Genov R. and Cauwenberghs G. Keltron: Support Vector "Machine" in Silicon, IEEE Trans. Neural Networks, Vol. 14. No.5. pp.1426 -1433, Sept. 2003.

10. Steinwart I. On the Optimal Parameter Choice for ν - Support Vector Machines, IEEE Trans. Pattern Analysis and Machine Intelligence, Vol. 25.

11. Schölkopf B. and Smola A.J. "A Tutorial on Support Vector Regression ", NeuroCOLT2 Technical Report Series, NC2 - TR - 1998 - 030, October, 1998.

12.Cristianini N. and Shawe-Taylor J. "An introduction to Support Vector Machines and other kernel - based learning methods", Cambridge, University Press, 2003.

13. Zhang L, Zhou W. and Jiao L. Wavelet Support Vector Machine, IEEE Trans. Systems, Man and Cybernetics - Part B: Cybernetics, Vol. 4. No.1. pp. 34 -39, Febr. 2004  

14. Hong X., Sharkey P. M. and Warwick. A Robust Nonlinear Identification Algorithm Using PRESS Statistic and Forward Regression, IEEE Trans. Neural Networks, Vol. 14. No. 2., pp.454 - 458,  March. 2003.

15. Cherkassy V., Gehring  D. and Mulier F. Comparison of adaptive methods for function estimation from samples, IEEE Trans. Neural Networks 7, pp. 969-984, 1996.

16. Szesztay, K. Einige Methoden der Vorhersage der Abflussverhältnisse, Vizgazdálkodási Tudományos  Kutató intézet. 1961.