Paul Honeine
email: paul.honeine@utt.fr

Cédric Zineb Noumir
email: zineb.noumir@utt.fr

Cédric Richard
email: cedric.richard@unice.fr

Multiclass classification machines with the complexity of a single binary classifier

Keywords: Multiclass classification, machine learning, SVM, one-versus-all, least-squares classification

come

Multiclass classification machines with the complexity of a single binary classifier

Introduction

The multiclass classification problem has been widely investigated in machine learning and data mining, in many fields including signal and image processing. Based on a set of observations from several classes, the aim is to learn a decision rule that accurately classifies new observations. Most classification techniques were initially derived for binary task, i.e., solving two-class classification problems. This is the case of many kernel-based classifiers, including the large-margin technique in support vector machines (SVM) [START_REF] Vapnik | Statistical Learning Theory[END_REF] and the least-squares machines (LSM) 1 [START_REF] Suykens | Least squares support vector machine classifiers[END_REF][START_REF] Rifkin | Regularized least squares classification[END_REF]. These classification techniques have been thoroughly investigated, and their performance well studied and understood. One seeks to generalize these approaches for multiclass tasks [START_REF] El-Yaniv | Better multiclass classification via a margin-optimized single binary problem[END_REF][START_REF] Liu | A novel and quick SVM-based multi-class classifier[END_REF], which find numerous applications in signal processing [START_REF] Galluccio | Graph based k-means clustering[END_REF][START_REF] Zeng | Kernel-based nonlinear discriminant analysis using minimum squared errors criterion for multiclass and undersampled problems[END_REF][START_REF] Wahlberg | Methods for alignment of multi-class signal sets[END_REF][START_REF] Aladjem | Multiclass discriminant mappings[END_REF]. Multiclass machines can be roughly divided into two categories. The first category attempts to incorporate simultaneously all the classification constraints within a single-machine scheme, in a similar way to binary classifiers (for instance, by maximizing the margin between classes). However, such a scheme significantly increases the problem complexity and thus requires advanced optimization techniques [START_REF] Weston | Support Vector Machines for Multi-Class Pattern Recognition[END_REF][START_REF] Crammer | On the algorithmic implementation of multiclass kernel-based vector machines[END_REF][START_REF] Bredensteiner | Multicategory classification by support vector machines[END_REF].

The second category consists of a divide-to-conquer strategy, such as the one-versus-all2 (OvA) [START_REF] Rifkin | In defense of one-vs-all classification[END_REF], the one-versus-one (OvO) [START_REF] Fürnkranz | Round robin classification[END_REF], and the decision directed-acyclic-graph (DAG) [START_REF] Platt | Large margin DAGs for multiclass classification[END_REF] schemes. In this case, the multiclass problem is decomposed into multiple binary tasks. All these subproblems are solved separately, and their results combined to infer the multiclass solution. For a survey, see [START_REF] Crammer | On the algorithmic implementation of multiclass kernel-based vector machines[END_REF][START_REF] Dietterich | Solving multiclass learning problems via error-correcting output codes[END_REF][START_REF] Ou | Multi-class pattern classification using neural networks[END_REF][START_REF] Allwein | Reducing multiclass to binary: a unifying approach for margin classifiers[END_REF].

In [START_REF] Rifkin | In defense of one-vs-all classification[END_REF], Rifkin and Klautau conducted a comparative study of multiclass machines, including those presented in all the aforementioned papers. Since these results are well-known in the literature, we give here a small review: multiclass schemes such as OvO, OvA, and DAG, have essentially the same accuracy as single-machine schemes. The LSM performs just as well as SVM, as illustrated in many studies [START_REF] Rifkin | Regularized least squares classification[END_REF][START_REF] Suykens | Least Squares Support Vector Machines[END_REF] including [START_REF] Van Gestel | Benchmarking least squares support vector machine classifiers[END_REF] with binary and multiclass classification. Therefore, a simple scheme such as OvA (for SVM or LSM) is preferable to a more complicated single-machine or error-correcting coding scheme. It is often argued that the main drawback of an OvA scheme is its computational complexity. For an ℓ class classification task, it requires to train ℓ binary classifiers, each one involving all the training data. This should be compared with the OvO and DAG schemes, which both need to train ℓ(ℓ -1)/2 binary classifiers, each individual subproblem being much smaller. See [START_REF] Fürnkranz | Round robin classification[END_REF] for more details.

In this paper, we study the multiclass classification problem by providing algorithms with the complexity of a single binary classifier. To this end, we derive a framework for multiclass problems by considering two possible forms of the solution, which are equivalent in the binary case [START_REF] Kimeldorf | Some results on tchebycheffian spline functions[END_REF][START_REF] Schölkopf | A generalized representer theorem[END_REF]. The first form includes explicitly the class labels in the coefficients of the decision function to be estimated. The second form includes them implicitly. This leads to two classes of machines.

On the one hand, when labels are explicitly included in the coefficients of the expansion, we explore a new trend of deriving low computational algorithms based on learning vector-valued functions. Initially introduced for multi-task learning [START_REF] Micchelli | On learning vector-valued functions[END_REF][START_REF] Evgeniou | Learning multiple tasks with kernel methods[END_REF], this idea has been naturally applied for multiclass classification problems, by using a SVM-like formulation [START_REF] Szedmak | Multiclass learning at one-class complexity[END_REF][START_REF] Szedmak | Learning via linear operators: Maximum margin regression; multiclass and multiview learning at one-class complexity[END_REF]. More recently, we derived a least-squares formulation in [START_REF] Noumir | Multi-class least squares classification at binary-classification complexity[END_REF][START_REF]Classification multi-classes au prix d'un classifieur binaire[END_REF]. In this paper, we provide a framework for deriving several algorithms based on binary classification algorithms, including LSSVM, RLS, and SVM.

On the other hand, we consider a model where the labels do not appear explicitly. We show that the nature of the LSM with the OvA scheme provides an algorithm, with essentially the same computational complexity as a single binary classifier. The proposed method is designated by oneLSM in this paper. We conduct a study on the model coefficients and we establish connections between the labels in multiclass problems. A single-machine formulation is proposed, and connections with other related work are presented [START_REF] Cortes | A general regression framework for learning string-to-string mappings[END_REF].

We also conduct a comprehensive study on the class coding in multiclass classification. Several techniques have been proposed in the literature, as extensions of the binary problem. In this paper, we show that they are equivalent, including the ±1 coding [START_REF] Dietterich | Solving multiclass learning problems via error-correcting output codes[END_REF][START_REF] Allwein | Reducing multiclass to binary: a unifying approach for margin classifiers[END_REF], the standard basis coding (or indicators) inspired from artificial neural networks (ANN) [START_REF] Bishop | Neural networks for pattern recognition[END_REF], the alignment-based coding [START_REF] Guermeur | SVM multiclasses, théorie et applications[END_REF][START_REF] Noumir | Multi-class least squares classification at binary-classification complexity[END_REF], the consistency-based coding used with the inductive principle in multiclass SVM [START_REF] Lee | Multicategory support vector machines, theory, and application to the classification of microarray data and satellite radiance data[END_REF], and the minimum-correlation coding proposed in [START_REF] Szedmak | Learning via linear operators: Maximum margin regression; multiclass and multiview learning at one-class complexity[END_REF]. We study these label codings and investigate their equivalence. Experimental results corroborate these findings. A nonlinear extension is also proposed by applying a kernel on the labels.

The rest of the paper is organized as follows. Sec. 2 provides a succinct presentation of the classification problem in kernel-based machine learning. We describe the proposed framework in Sec. 3, and provide analogy with the binary classification case. This leads to two classes of machines, the vector-output machines defined in Sec. 4, and the oneLSM given in Sec. 5 where we revisit the least-squares machines. In Sec. 6, we study the label coding for multiclass classification. Sec. 7 is devoted for further discussions. Sec. 8 gives comparative results in terms of accuracy and computational time, as well as a statistical test. Conclusions and further directions are given in Sec. 9.

Notations:

In this paper, we use the superscript index to designate an entry related to a subproblem, such as f (k) (•) and α (k) associated to the k-th subproblem classification. The subscript index is related to a given sample, such as y i and α i corresponding to the i-th sample x i .

From binary to multiclass classification: a primer

Binary classification with kernel machines

In a two-class classification problem, one seeks a decision rule that predicts well the class membership of a given sample, based on a set of training data with available class membership. Such decision rule, for any new sample x, takes the form

if f (x) < 0, then y = -1 if f (x) > 0, then y = +1 (1)
namely, comparing the sign of f (x) to determine if x belongs to the (-1)-class or to the (+1)-class 3 . This function f (•) is estimated using a training set of labeled data, {(x 1 , y 1 ), (x 2 , y 2 ), . . . , (x n , y n )}, with label y i = ±1 depending on which class sample x i belongs to. This decision rule is equivalent to

arg max y=±1 yf (x), (2)
as well as arg min y=±1 |f (x) -y| 2 .

A straightforward solution for the binary classification problem is obtained when one seeks a function from the space H of linear functions,

f (x) = w ⊤ x. (3)
Consider the least-squares classification problem (also known as Ridge Regression) defined by min

f ∈H n j=1 |w ⊤ x j -y j | 2 + γ w 2 , ( 4
)
where γ is a tunable parameter that controls the trade-off between errors on the training data and regularity of the solution. Taking the derivative of the above function with respect to w, namely 2 n j=1 (w ⊤ x j -y j ) x j + 2γ w, and setting it to zero at the optimum, we observe that w can be written in the form4

w ⊤ = n i=1 α i x ⊤ i , (5)
with α i = (y iw ⊤ x i )/γ. Substituting the above form of w in the previous one, we obtain the linear system of n equations with n unknowns

n j=1 α j x ⊤ j x i + γ α i = y i , for i = 1, 2, . . . , n.
In matrix form, we have w = X α with X = [x 1 x 2 • • • x n ] and α a column vector with entries α j for j = 1, 2, . . . , n. Therefore, the above linear system can be written as

(K + γI) α = y,
where K is a n-by-n matrix with entries x ⊤ i x j , for i, j = 1, 2, . . . , n, i.e., K = X ⊤ X, y is a column vector with entries y j , for j = 1, 2, . . . , n, and I is the identity matrix of appropriate size (n-by-n here).

The resulting coefficients are considered in the decision rule for any x, with f (x) = α ⊤ κ x , where κ x is a column vector whose j-th entry is x ⊤ j x. The decision rule is given by comparing the evaluation of this function f (x) to a threshold according to [START_REF] Vapnik | Statistical Learning Theory[END_REF]. This threshold (called bias) is often set to 0, as advised in many studies [START_REF] Poggio | b[END_REF].

There exists another formulation of the Ridge Regression problem, when the labels y i = ±1 are applied.

In this case, the least-squares problem (4) is equivalent to

min f ∈H n j=1 |y j f (x j ) -1| 2 + γ f 2 . ( 6
)
By considering the space of linear functions, of the form f (x) = w ⊤ x, we obtain the following solution

w ⊤ = n i=1 β i y i x ⊤ i . (7)
The connection between (7) and ( 5) is obvious, with α i = β i y i . Both criteria, f (x j ) -y j and y j f (x j ), are considered in [START_REF] Suykens | Least Squares Support Vector Machines[END_REF] for regression, as opposed to [START_REF] Suykens | Least squares support vector machine classifiers[END_REF] where only the latter is used. It is worth noting that the equivalence between both criteria is true only when the (±1)-label coding is considered in classification.

Finally, thanks to the kernel trick, one can easily extend this linear model into a nonlinear one, by substituting some kernel function κ(x i , x j ) for x ⊤ i x j . In fact, inner products are involved, both in the estimation of the optimal coefficients α (with K) and in the decision function f (x) = α ⊤ κ x (with κ x ). The most used kernels are the polynomial kernel, of the form κ(x i , x j ) = (x ⊤ i x j + c) p , and the Gaussian kernel, with κ(x i , x j ) = exp(x ix j 2 /2σ 2 ), where σ is the bandwidth parameter.

Multiclass classification

Consider an ℓ-class classification task. Using a divide-to-conquer strategy, the problem can be solved by working on a collection of binary classification subproblems, then combining each of the resulting real-valued decision functions. Thus, each function f (k) (•) is defined from a binary classification subproblem by using the same training data, x 1 , x 2 , . . . , x n , and assigning its proper labels, y

2 , . . . , y

n . In a one-versus-all (OvA) scheme, ℓ binary classifiers are trained, where each subproblem confronts a class with all the rest. In this case, we have y (k) j ∈ {-1; +1} depending on the k-th class membership of x j . This requires estimating ℓ × n coefficients, α (k) j , for j = 1, 2, . . . , n and k = 1, 2, . . . , ℓ. In a one-versus-one (OvO) scheme, ℓ(ℓ -1)/2 binary classifiers are trained by taking data from each pair of classes. For this purpose, y (k) j ∈ {-1; 0; 1} where the value zero corresponds to samples that do not belong to the classes under examination. Although ℓ(ℓ -1)/2 binary classifiers are trained, the individual problems are much smaller. Note that each sample is used by n -1 decision functions. It is well-known that these schemes (as well as directed-acyclic-graph, complete codes, ...) perform essentially identically (see Sec. 1 for references). Therefore, the simplest scheme is preferable to more complex ones, making OvA and OvO good candidates.

Still, one needs to solve several binary classification subproblems.

Any multiclass machine, defined by a set of ℓ functions f (k) (•), is completely described by the coefficients to be determined, either β (k) j or α (k) j for j = 1, 2, . . . , n and k = 1, 2, . . . , ℓ. This yields an estimation problem with n × ℓ unknown parameters. In this paper, we show that, without sacrificing performance, one can naturally reduce significantly the number of unknowns, by imposing some relation between the functions f (k) (•)'s. This yields multiclass classification machines with essentially the same computational complexity as a single binary classifier.

The proposed framework

Let f (1) (•), f (2) (•), . . . , f (ℓ) (•) be the functions to be learnt 5 , in the same spirit as the OvA scheme while not limiting ourselves to it. In what follows, we regroup the ℓ functions in the vector-of-functions6

f (•) = [f (1) (•) f (2) (•) • • • f (ℓ) (•)] ⊤
For the k-th classification subproblem, defined by the function f (k) (•), we assign the n-entry label vector y (k) , namely the (±1)-label coding as defined in Table 1 (See Sec. 6 for a study of several codings of the labels). Let Y be the ℓ-by-n matrix with each row corresponding to the labels associated to a classification subproblem, namely Y ⊤ = [y (1) y (2) • • • y (ℓ) ]. In other words, each column of Y contains y i , the label associated to x i in all the classification subproblems. Thus, the i-th entry of y (k) is the k-th entry of y i .

We explore the differences, between (4) and [START_REF] Galluccio | Graph based k-means clustering[END_REF], and between the forms (5) and [START_REF] Zeng | Kernel-based nonlinear discriminant analysis using minimum squared errors criterion for multiclass and undersampled problems[END_REF], now written in a multiclass classification formulation. By coupling these diversities, we investigate several classification methods. The following table summarized these differences, in terms of the optimized criterion, the parameters to be estimated, β i 's or α i 's, as well as the sparsity of the solution:

Labelbook

[

y (k) ] i = [y i ] k y ⊤ i y j (±1)-label +1 if x i belongs to class k; -1 otherwise ℓ if y i = y j ; ℓ -4 otherwise Standard basis (or indicators) 1 if x i belongs to class k; 0 otherwise 1 if y i = y j ; 0 otherwise Alignment-based ℓ-1 ℓ if x i belongs to class k; -1 √ ℓ(ℓ-1) otherwise 1 if y i = y j ; -1 ℓ-1 otherwise Consistency-based 1 if x i belongs to class k; -1 ℓ-1 otherwise ℓ ℓ-1 if y i = y j ; -ℓ (ℓ-1) 2
otherwise Table 2: Analogy between the binary and the multiclass formulations.

binary formulation multiclass formulation models

f (•) f (x) = [f (1) (•) • • • f (ℓ) (•)] ⊤ w W = [w (1) • • • w (ℓ) ] f (x) = w ⊤ x f (x) = W ⊤ x w = Xα W = X[α (1) • • • α (ℓ) ] w ⊤ = i β i y i x ⊤ i W ⊤ = i β i y i x ⊤ i criterion y j f (x j ) y j f (x j ) |f (x j ) -y j | 2 f (x j ) -y j 2 w 2 W 2 F decision max y=±1 yf (x) max y∈Y y ⊤ f (x) min y=±1 |f (x) -y| 2 min y∈Y f(x) -y 2
vector-output machines (see Section 4) OvA (see Section 5)

RLS LSSVM SVM LSM criterion f (x j ) -y j y j f (x j ) f (x j ) -y j f (x j ) -y j β i 's or α i 's β i β i β i α i full or sparse full full sparse full
Finally, once the model parameters estimated, the decision rule is given by analogy with (2), for a given x :

arg max y y ⊤ f (x) (8)
The decision rule is studied in detail in Sec. 7.1.

The analogy between the binary formulation and the multiclass formulation is given in Table 2, where the vector w is used in the former and the matrix W in the latter. Thus, the Euclidean vector norm w is replaced by the matrix norm W F , where • F denotes Frobenius norm defined as

W 2 F = i,j [W ] i,j 2 = trace(W ⊤ W ),
We associate these methods with several label codes, beyond the conventional (±1)-label. See Table 1.

The alignment-based label coding can be viewed as a natural extension of the binary ±1 label coding. This is illustrated in Fig. 1. We make a thorough study of the label coding in Sec. 6. Independent of the choice of the label coding, we derive next algorithms for multiclass classification.

Vector-output machines for multiclass classification

In this section, we show that one can easily derive multiclass classification algorithms, in the same spirit as binary classification algorithms and specifically RLS, SVM and LSSVM.

Following [START_REF] Szedmak | Learning via linear operators: Maximum margin regression; multiclass and multiview learning at one-class complexity[END_REF] and more recently [START_REF] Noumir | Multi-class least squares classification at binary-classification complexity[END_REF][START_REF]Classification multi-classes au prix d'un classifieur binaire[END_REF], we propose the following expression for f (x)

f (x) = W ⊤ x, (9)
with W a d-by-ℓ matrix defined, by analogy with [START_REF] Zeng | Kernel-based nonlinear discriminant analysis using minimum squared errors criterion for multiclass and undersampled problems[END_REF], as

W ⊤ = n j=1 β j y j x ⊤ j .
In this case, we have

f (x) = n j=1 β j y j x ⊤ j x, (10)
which means that the vector f (x) belongs to the span of the training label vectors y 1 , y 2 , . . . , y n .

An interesting property of this model is that the functions f (•)'s in f (•) share the same coefficients, namely

β j = β (k) j for all k = 1, 2, . . . , ℓ -1 0 +1 y (2) y (3)
y (1) 0

Figure 1: Illustration of the class-membership encoding. Left: In the binary case, y = ±1 are unit-norm and defined by a segment centered at the origin. Right: In the multiclass case, the vector labels are unit-norm vertices and having the origin as centroid, as defined in Table 1 (here ℓ = 3) with

y (1) = √ 2 √ 3 -1 √ 6 -1 √ 6 , y (2) = -1 √ 6 √ 2 √ 3 -1 √ 6 , y (3) = -1 √ 6 -1 √ 6 √ 2 √ 3
.

In [START_REF] Weston | Support Vector Machines for Multi-Class Pattern Recognition[END_REF], these coefficients are weighted by the label associated to x j in the k-th class membership, i.e, [y (k) ] j = ±1, which leads to the only difference between the functions in f (•). This allows us to have only n unknowns, in the same way as a single binary classifier. Next, we show how to estimate these coefficients, by adapting several classical classification algorithms to operate for vector output:

• multiclass vo-RLS

• multiclass vo-LSSVM

• multiclass vo-SVM which provides a sparse solution to the latter It is worth noting that both RLS and LSSVM are equivalent in the scalar-output case, when y j = ±1.

Finally, once the model parameters β i 's estimated, the decision rule is defined by combining ( 8) with [START_REF] Weston | Support Vector Machines for Multi-Class Pattern Recognition[END_REF], for any given x, with arg max

y n j=1 β j y ⊤ y j x ⊤ j x. (11)

Multiclass vo-RLS

In the least-squares sense, we consider the following optimization problem min

f n j=1 f (x j ) -y j 2 + γ β 2 , ( 12
)
where

β = [β 1 β 2 • • • β n ] ⊤ .
By substituting the expression of f (x) in this optimization problem, we get the following cost function to be minimized with respect to the β i 's:

n i=1 y ⊤ i y i + n i=1 n j,k=1 β j β k x ⊤ i x j x ⊤ i x k y ⊤ j y k + γ n i=1 β 2 i -2 n i=1 y ⊤ i n j=1 β j x ⊤ i x j y j .
This optimization problem can be written in matrix form as,

min β Y ⊤ Y -2 d β + β ⊤ Gβ + γ β ⊤ β, ( 13
)
where G is a matrix whose (j, k)-th entry is

[G] j,k = y ⊤ j y k x ⊤ j x k 2 ,
and d is a vector whose j-th entry is

[d] j = n i=1 y ⊤ i y j x ⊤ i x j .
By taking the gradient of the objective function in [START_REF] Rifkin | In defense of one-vs-all classification[END_REF] with respect to β, namely -d + Gβ + γ β, and setting it to zero, we obtain the final solution

(G + γ I) β = d. (14)
One can also consider other regularization terms in the above optimization problem. For instance, we consider the following optimization problem min

f n j=1 f (x j ) -y j 2 + γ R(f ), (15)
where the regularization term is given by

R(f ) = n i,j=1 β i β j y ⊤ i y j x ⊤ i x j . (16)
From above, it is easy to see that we obtain the final solution

(G + γ H) β = d, ( 17
)
where H is a matrix whose (i, j)-th entry is

[H] i,j = y ⊤ i y j x ⊤ i x j .
In this paper, we study the performance of both optimization problems. The vo-RLS problem defined by ( 12)-( 14) is denoted by vo-RLS(β), while the use of the regularization term [START_REF] Dietterich | Solving multiclass learning problems via error-correcting output codes[END_REF], establishing the problem ( 12)-( 14), is denoted by vo-RLS(f).

We conclude that a multiclass regularized least squares classifier can be obtained by solving a single system of n linear equations with n unknowns. This requires the inversion of a n-by-n matrix. The computational complexity of this algorithm is cubic in the number of training data, n, but essentially independent of the number of classes ℓ. This is made possible here thanks to the fact that the label vectors only appear in terms of inner products, with y ⊤ i y j in H, G and d. This reflects an analogy with the kernel trick in kernel machines, where data are involved only in terms of inner products, namely x ⊤ i x j . See Sec. 7.2 for a discussion on applying nonlinear kernels on the labels.

Multiclass vo-LSSVM

Another way to tackle the problem, is to solve an optimization problem with equality constraints, such as min

W ,ξ 1 2 W 2 F + γ 1 2 n i=1 ξ 2 i , subject to y ⊤ i (W ⊤ x i ) = 1 -ξ i , for i = 1, 2, . . . , n.
The corresponding Lagrangian is defined by

L(W, ξ, β) = 1 2 W 2 F + γ 1 2 n i=1 ξ 2 i - n i=1 β i (y ⊤ i (W ⊤ x i ) -1 + ξ i ),
where β i 's are the Lagrangian multipliers. The optimality conditions are given by

                 ∂L ∂W = 0 ⇒ W ⊤ = i β i y i x ⊤ i ∂L ∂ξi = 0 ⇒ β i = γ ξ i ∂L ∂βi = 0 ⇒ y ⊤ i (W ⊤ x i ) -1 + ξ i = 0
This leads to the following linear system

H + γ -1 I β = 1 n
where 1 n is a n-entry column vector of ones. Solving this linear system, i.e., identifying β, requires the inversion of a n-by-n matrix. The computational cost of such inversion is about O(n 3 ). This should be compared with Suykens and Vanderwalle's multiclass LSSVM algorithm in [START_REF] Suykens | Multiclass least squares support vector machines[END_REF] which requires the inversion of an (ℓn + ℓ)-by-(ℓn + ℓ) matrix, with a computational complexity that scales cubically with the number of classes, namely O((ℓn) 3 ).

Multiclass vo-SVM

In the same spirit as the support vector machines [START_REF] Vapnik | The nature of statistical learning theory[END_REF], we consider the following optimization problem: min

W ,ξ 1 2 W 2 F + γ n i=1 ξ i (18)
subject to

y ⊤ i (W ⊤ x i ) ≥ 1 -ξ i , for i = 1, 2, . . . , n, and
ξ i ≥ 0, for i = 1, 2, . . . , n.
One can solve this constrained optimization problem using Lagrangian, L(W , ξ, β, δ) with β and δ the Lagrangian multipliers, namely min

W ,b,ξ max β,δ 1 2 trace(W ⊤ W ) + γ n i=1 ξ i - n i=1 δ i ξ i - n i=1 β i y ⊤ i (W ⊤ x i ) -1 + ξ i .
The minimum of the Lagrangian with respect to W and ξ is given by

       ∂L ∂W = 0 ⇒ W ⊤ = i β i y i x ⊤ i ∂L ∂ξi = 0 ⇒ β i + δ i = γ
This leads to the following dual form max

β n i=1 β i - 1 2 n i,j=1 β i β j y ⊤ i y j x ⊤ i x j (19)
subject to

0 ≤ β 1 , β 2 , . . . , β n ≤ γ. (20)
In matrix form, we can write max

β 1 ⊤ n β - 1 2 β ⊤ Hβ subject to 0 n ≤ β ≤ γ1 n ,
where 0 n is a n-entry column vector of zeros. This is a quadratic programming problem, which can be solved using any off-the-shelf optimization technique. The is essentially similar to the SVM binary-classifier, thus the same optimization routines can be used for both binary and multiclass classification tasks 7 .

Model with a bias

In this paper, we study the bias-free model [START_REF] Aladjem | Multiclass discriminant mappings[END_REF]. One may also consider an offset in the model, with

f (x) = W ⊤ x + b,
where b is the bias parameter. It is easy to see that our approach extends naturally to this model, where the constraint y ⊤ i (W ⊤ x i ) = 1 -ξ i is substituted with

y ⊤ i (W ⊤ x i + b) = 1 -ξ i ,
in the vo-LSSVM and vo-SVM algorithms. In their dual forms, we get the following additional constraint:

Y β = 0 ℓ . (21)
This constraint can be easily incorporated within the proposed algorithms, as illustrated here.

Consider the vo-LSSVM, this leads to the following augmented linear system

0 Y Y ⊤ H + γ -1 I b β = 0 ℓ 1 n .
where 0 is the square matrix of zeros of appropriate size (ℓ-by-ℓ here). Solving this linear system requires the inversion of a (n + ℓ)-by-(n + ℓ) matrix. Since the number of classes is significantly smaller than the number of available data, the computational cost of such inversion is about O(n 3 ). This still outperforms the multiclass LSSVM algorithm proposed by Suykens and Vanderwalle's in [START_REF] Suykens | Multiclass least squares support vector machines[END_REF] with O((ℓn) 3 ).

In the case of the vo-SVM, we get max

β 1 ⊤ n β - 1 2 β ⊤ Hβ
7 For instance, one can use the same matlab function quadprog for both problems.

subject to

Y β = 0 ℓ and 0 n ≤ β ≤ γ1 n .
This quadratic programming problem is similar to the SVM binary-classifier. The only difference is the linear constraint, here n i=1 β i y i = 0 ℓ as opposed to n i=1 β i y i = 0 for the binary case. This is essentially the multiclass SVM algorithm proposed Szedmak et al. in [START_REF] Szedmak | Learning via linear operators: Maximum margin regression; multiclass and multiview learning at one-class complexity[END_REF] (see also [START_REF] Szedmak | Multiclass learning at one-class complexity[END_REF]).

Discussion: biased versus unbiased model

The use of a biased or an unbiased model is still an open question in machine learning, in the binary classification case [START_REF] Poggio | b[END_REF] as well as in the multiclass case [START_REF] Abril | A note on the bias in svms for multiclassification[END_REF]. For instance in SVM, and by analogy with the binary case, one can define the optimal value of the bias by averaging, over all support vectors x i , the expression

b i = y i - n j=1 β j y j x ⊤ j x i .
In [38, page 203], it is advised not to use an optimal value, but to change it in order to adjust the number of false positives and false negatives. Many authors disallow the bias term completely, i.e., b = 0 ℓ , see for instance [START_REF] Ertekin | Nonconvex online support vector machines[END_REF][START_REF] Poggio | b[END_REF]. It is worth noting that in SVM implementations that disallow the bias term, the linear constraint n i=1 β i y i = 0 ℓ is removed from the SVM problem. For the least-squares approaches, many studies motivate the use of the unbiased version, see for instance [START_REF] Poggio | b[END_REF][START_REF] Rifkin | Everything old is new again: A fresh look at historical approaches in machines learning[END_REF].

In practice, it turns out that the unbiased model often outperforms the biased one. This is illustrated in Sec. 8, e.g., by comparing our vo-SVM with the one proposed by Szedmak et al. in [START_REF] Szedmak | Learning via linear operators: Maximum margin regression; multiclass and multiview learning at one-class complexity[END_REF]. A statistical test is also conducted on the performance of the latter, showing that it highly depends on the choice of the label coding. Still the classification accuracy remains poor in general, as shown in our experimentations.

The multiclass least-squares machines : oneLSM

In this section, we revisit the least-squares solution with

w ⊤ = n i=1 α i x ⊤ i ,
and show how our framework is a natural choice for the multiclass problem. To this end, we recall the least-squares problem for the binary case, as defined in Sec. 2. The solution is given as (K + γI) α = y, which leads to f (x) = y ⊤ (K + γI) -1 κ x .

Multiclass LSM in one shot

Consider the OvA scheme for an ℓ-class classification task. This corresponds to estimate ℓ decision functions, each taking the form

f (k) (x) = n i=1 α (k) i x ⊤ i x, (22)
for k = 1, 2, . . . , ℓ. The ℓ vectors of unknowns are denoted by α (k) , with entries α (k) i for i = 1, 2, . . . , n. It is obvious that all these functions share the same input data, x 1 , x 2 , . . . , x n . The only difference between the binary classifiers is the assigned labels in the subproblems. For the k-th binary classification subproblem, we assign the label vector y (k) as given in Table 1. See Sec. 6 for a study of the label coding and its impact on the solution.

Then the multiclass problem is defined by the expressions α (1) = (K + γI) -1 y (1) ,

α (2) = (K + γI) -1 y (2) , . . . ( 23
)
α (ℓ) = (K + γI) -1 y (ℓ) .
It is obvious that one only needs to inverse a single matrix, (K + γI), for all the binary classification subproblems. Moreover, the above ℓ equations (each implicitly involving n equalities) can be easily written in matrix form with

[α (1) α (2) • • • α (ℓ) ] = (K + γI) -1 Y ⊤ ,
and consequently we obtain

f (x) = Y (K + γI) -1 κ x . (24)
It is obvious that this approach, denoted oneLSM in this paper, is fundamentally an OvA scheme of LSM, without naively computing ℓ binary classifiers. Finally, the label of any given x is determined by

arg max y y ⊤ Y (K + γI) -1 κ x . (25)
By studying the computational requirements, the matrix inversion is the most demanding operation. The proposed trick allows us to compute a single matrix inversion of an n-by-n matrix, rather than ℓ inversions of the same matrix for the naive implementation.

Relationship between the parameters

While expression [START_REF] Schölkopf | A generalized representer theorem[END_REF] shows that there are n × ℓ unknowns, the α

i 's, there exists a connection between these parameters. The following proposition can be easily derived from (23) : Proposition 1. In the oneLSM multiclass problem (as well as in LSM with the OvA scheme), the coefficients of binary classifiers are connected to each others with the following relation :

α (k) ⊤ y (k ′ ) = α (k ′ ) ⊤ y (k) ,
for all pairs k, k ′ = 1, 2, . . . , ℓ.

This result leads to other expressions, such as [α (1) k) . These results give very useful insights, and can be illustrated on specific types of labels, such as the standard basis (or indicators) given in Table 1, namely [y (k) ] j = 1 if x j belongs to class k and 0 otherwise. In this case, we have

α (2) • • • α (ℓ) ] ⊤ y (k) = [y (1) y (2) • • • y (ℓ) ] ⊤ α (
xj ∈C (k ′ ) α (k) j = xi∈C (k) α (k ′ ) i , for any k, k ′ = 1, 2, . . . , ℓ,
where C (k) denotes the k-th class. For a given binary classifier k, we consider the ℓ -1 above equations obtained for C (k ′ ) of the remaining classes, which sums to

k ′ =k xj ∈C (k ′ ) α (k) j = xi∈C (k) k ′ =k α (k ′ )
i , for a given k = 1, . . . , ℓ.

The above left-hand side can be written as xi ∈C (k) α (k) j . Therefore, the above equation can be read as follows: For a given binary classifier k-against-all the rest, the contributions (in terms of α's) in classifier k of the data belonging to all the remaining classes is equal to the contribution in all the remaining classifiers of the data belonging to the k-th class.

To the best of our knowledge, these connections between the parameters and the labels were never stated or proved before. These results are to be compared with the equality constraints [START_REF] Kimeldorf | Some results on tchebycheffian spline functions[END_REF] in the SVM and LSSVM machines, as given in Sec. 4.4. In the well known binary classification case with (±1)-label coding, we have

n i=1 β i y i = 0, thus xj ∈C β j = xi ∈C β i .

Multiclass oneLSM as a single-machine problem

The proposed approach can be derived using a single-machine problem as illustrated in next proposition.

To this end, we write [START_REF] Evgeniou | Learning multiple tasks with kernel methods[END_REF] as

f (x) = W ⊤ x, where W = [w (1) w (2) • • • w (ℓ)
] with w (k) associated to the k-th binary classifier.

Proposition 2. The multiclass LSM problem is equivalent to the optimization problem

min W n i=1 W ⊤ x i -y i 2 + γ W 2 F , whose solution W = X(X ⊤ X + γI) -1 Y ⊤ defines the decision function f (x) = W ⊤ x.
Proof. To see this, consider the derivative with respect to W ⊤ of the above cost function, namely

n i=1 W ⊤ x ⊤ i x i -y i x ⊤ i + γW ⊤ ,
and set it to zero:

W ⊤ (XX ⊤ + γI) = Y X ⊤ .
This leads to W ⊤ = Y (X ⊤ X +γI) -1 X ⊤ , where the matrix identity X ⊤ (XX ⊤ +γI) -1 = (X ⊤ X +γI) -1 X ⊤ is considered.

Labels

For the sake of simplicity, we will use label to name a label vector, and labelbook Y to denote the set of ℓ label vectors, each one being associated to a class. Thus, for any x i we associate y i ∈ Y, with y i = y j if and only if x i and x j belong to different classes. The choice of a labelbook for a given classification task is studied in this section.

Several labelbooks can be used as extensions of the binary case to the multiclass case. Expressions of the most used ones are given in Table 1: • A straightforward generalization of the binary case to the multiclass case is to keep all the entries with -1 except a single +1 to identify the class. This labelbook is often suggested in the literature. See for instance [START_REF] Dietterich | Solving multiclass learning problems via error-correcting output codes[END_REF][START_REF] Allwein | Reducing multiclass to binary: a unifying approach for margin classifiers[END_REF].

[+1 -1 -1] ⊤ (±1)-label [1 0 0] ⊤ Standard basis +1 -1 2 -1 2 ⊤ Consistency-based √ 2 √ 3 -1 √ 6 -1 √ 6 ⊤ Alignment-based
• The standard basis (also known as indicators) is used in artificial neural networks [START_REF] Bishop | Neural networks for pattern recognition[END_REF][START_REF] Huang | Extreme learning machine for regression and multiclass classification[END_REF] (e.g., Boltzmann networks and extreme learning machines), and is called dummy variables by statisticians.

In this case, the labelbook is defined by the columns of the identity matrix.

• The alignment-based labelbook is optimal in the sense of the kernel-target alignment criterion [START_REF] Guermeur | SVM multiclasses, théorie et applications[END_REF],

and has been systematically applied in many multiclass machines [START_REF] Noumir | Multi-class least squares classification at binary-classification complexity[END_REF][START_REF] Szedmak | Learning via linear operators: Maximum margin regression; multiclass and multiview learning at one-class complexity[END_REF].

• In [START_REF] Lee | Multicategory support vector machines, theory, and application to the classification of microarray data and satellite radiance data[END_REF], the authors study the consistency of the inductive principle in multiclass SVM. This desirable property is satisfied by the consistency-based labelbook.

See Fig. 2 for a three-dimensional illustration, i.e., a 3-class classification problem. For any of these labelbooks, we see that all the entries of a given label y i are identical except for a single entry, which corresponds to the class membership. In what follows, we call such category of labelbooks the one-per-class labelbooks (short for one entry that identifies the class). This is similar to the notion of class-symmetric matrix as defined in [START_REF] Rifkin | In defense of one-vs-all classification[END_REF], where the authors consider only {-1, +1} entries. Still, in our definition, the values may be arbitrary. Table 1 presents several well-known labelbooks. We also present the corresponding inner products, since this is the essential information of any labelbook, as illustrated in the parameters estimation and the decision rule, in [START_REF] Crammer | On the algorithmic implementation of multiclass kernel-based vector machines[END_REF] and ( 25), namely arg max

y y ⊤ Y (K + γI) -1 κ x .
where only inner products between labels is required.

Each one-per-class labelbook defines ℓ distinct (vector) labels of ℓ entries. Thus, the ℓ vectors lie in a ℓ-dimensional space. One may also consider a labelbook of ℓ vectors in an ℓ -1 dimensional space. For instance, when three classes are in competition, this corresponds to three vectors spanning an equilateral triangle. In [START_REF] Szedmak | Learning via linear operators: Maximum margin regression; multiclass and multiview learning at one-class complexity[END_REF], the authors generalize this result to any number of classes, by imposing that correlation between any distinct vectors is the same, and minimized. The resulting vectors are given as the solution of an eigen-problem.

Proposition 3. The minimum-correlation labelbook is equivalent to the alignment labelbook.

Proof. From [26, Proposition 2], it is shown that the minimum-correlation labelbook is defined by ℓ vectors spanning a ℓ -1 dimensional subspace and the inner-products of distinct vectors equal to -1/(ℓ -1). It turns out that this is exactly the inner-products of the alignment-based labelbook (see third column in Table 1). Since both labelbooks have the same y ⊤ i y j , they are equivalent for both vector-output machines and the oneLSM machines.

The consistency-based and the alignment-based labelbooks are identical up to the multiplicative constant (ℓ -1)/ℓ (see Fig. 2). Thus, they give comparable results in multiclass classification, since the scale factor will be absorbed in optimization. In what follows, we study the one-per-class labelbooks (e.g., those given in Table 1) for the oneLSM machines.

Labels in oneLSM machines

Next, show that the one-per-class labelbooks are equivalent for the oneLSM optimization problem ( 24)- [START_REF] Szedmak | Multiclass learning at one-class complexity[END_REF]. But before, we show the link between these labelbooks.

Lemma 4. The one-per-class labelbooks can be generated from each others using the linear transformation

ay i + b1, (26)
from some arbitrary a and b, where y i spanning some given labelbook and 1 is the column vector with ℓ ones.

The proof is straightforward, by generating any one-per-class labelbook from the standard basis labelbook, and vice versa. Therefore, we can define any of the one-per-class labelbook with elements y i such that

[y i ] k = a + b if x i belongs to class k; b otherwise
where the standard basis is implicitly used.

Theorem 5. The one-per-class labelbooks are equivalent for the oneLSM.

Proof. Consider the transformation according to [START_REF] Szedmak | Learning via linear operators: Maximum margin regression; multiclass and multiview learning at one-class complexity[END_REF], when all the labels, y ⊤ and Y , are transformed into ay + b1 and aY + b11 ⊤ , respectively. Next, we show that the decision rule ( 25) is invariant to such transformation:

arg max y∈Y (ay + b1) ⊤ (aY + b11 ⊤ )(K + γI) -1 κ x = arg max y∈Y a 2 y ⊤ Y + ab y ⊤ 11 ⊤ (K + γI) -1 κ x = arg max y∈Y a 2 y ⊤ Y (K + γI) -1 κ x = arg max y∈Y y ⊤ Y (K + γI) -1 κ x
where the first equality follows from removing terms independent of y, and the second equality is due to the fact that y ⊤ 1 is constant for all y ∈ Y in a given labelbook.

Corollary 6. The one-per-class labelbooks are equivalent to the minimum-correlation labelbook.

This corollary follows directly from combining Proposition 3 and Theorem 5.

Beyond establishing the equivalence between several labelbooks, Theorem 5 also provides the equivalence between the decision rule ( 27) and the maximum-value rule used in the classical OvA scheme. The proof is straightforward for the standard basis labelbook, since arg max y∈Y y ⊤ f (x) yields arg max 1≤k≤ℓ f (k) (x).

The extension of this result to other labelbooks is due to Theorem 5.

Discussions

The decision rule

Several decision rules can be implemented. In classical OvA scheme, the multiclass decision for any new sample x is obtained by selecting the class whose corresponding classifier f k (x) has maximum value, i.e., the winner-takes-all strategy. This can be easily applied in our approach, by simply inspecting the resulting vector-of-functions f (x). Other decision rules are the classical Hamming distance, or the less-common Bayesian distance measure which gives an estimation of the posterior probability; see [START_REF] Dietterich | Solving multiclass learning problems via error-correcting output codes[END_REF].

By analogy with (2), we consider the following decision rule for any given observation x:

arg max y∈Y y ⊤ f (x), (27)
where the identified label determines the membership class. By substituting f (x) with its expression in [START_REF] Evgeniou | Learning multiple tasks with kernel methods[END_REF], it is easy to see that the labels are only involved in terms of inner products, with y ⊤ Y . This decision rule is motivated in [START_REF] Szedmak | Learning via linear operators: Maximum margin regression; multiclass and multiview learning at one-class complexity[END_REF] for maximum margin in vector-output SVM. Since by construction one may set the norm of the label to be constant for a given labelbook8 , this decision criterion is equivalent to

arg min y∈Y f (x) -y 2 .
It is worth noting that the proposed decision rule allows the use of any labelbook. However, the decision rule [START_REF] Noumir | Multi-class least squares classification at binary-classification complexity[END_REF] shows a dominant influence of the label coding on the result. In the next subsection, we show that the considered inner-product form provides a non-linear kernel-like expansion.

Kernels on the labels

It is well known that linear classification techniques can be easily extended for nonlinear classification, i.e., where data are not linearly separable, thanks to the kernel trick. Moreover, the same rule can be applied for the labels. This is true since the proposed framework leads to algorithms that depend on the labels in terms of inner products y ⊤ i y j , as illustrated in this paper with vo-RLS, vo-LSSVM, vo-SVM, and oneLSM (see for instance [START_REF] Szedmak | Multiclass learning at one-class complexity[END_REF]). Thus one may also propose a kernel function for the labels, namely

κ y (y i , y j ) = Ψ(y i ) ⊤ Ψ(y j ),
where the nonlinear function Ψ(•) maps the output space to some higher dimensional space. Nevertheless, the labels are by construction separable, as illustrated in Fig. 1. Therefore, we do not think that a kernel function for labels is necessary in multiclass classification. Moreover, the choice of the optimal kernel as well as the optimal value of its parameter(s) is an open issue. We have conducted some preliminary experiments on applying a Gaussian kernel to the labels, and found no significant difference in performance when compared to the linear kernel on the labels.

Output kernels may be interesting as they allow us to include prior knowledge, e.g., incorporating the probability priors on the classes to deal with cases when the sizes of the classes are highly unbalanced by applying uneven weights. Another application of the label kernels provides an approach to tackle with the complex label structures, e.g. hierarchical relations. For example in [START_REF] Tsochantaridis | Large margin methods for structured and interdependent output variables[END_REF], the authors propose a SVM-like algorithm by applying joint kernels between inputs and outputs. Joint kernels can be easily implemented within the framework proposed in paper, as illustrated by the different forms of vo-RLS, vo-LSSVM, SVM, and oneLSM.

In [START_REF] Cortes | A general regression framework for learning string-to-string mappings[END_REF], the authors study transduction learning by solving a problem similar to the one given in Proposition 2 but with label kernels, i.e., by including implicitly some nonlinear transformation Ψ(•) on the labels.

Predicting structured outputs, such as in graphs, the resulting problem requires to find the inverse transformation of Ψ(•). This is an ill-posed problem, where function Ψ(•) is not explicited thanks to the kernel trick. This inverse problem, known as the pre-image problem, is nonconvex and nonlinear. See [START_REF] Honeine | Preimage problem in kernel-based machine learning[END_REF] for a definition of the pre-image problem with techniques to solve this problem.

Experimentations

To illustrate the relevance of the proposed approach, we have tested all the algorithms on several wellknown datasets (available at the UCI Repository): iris, wine, glass, vowel, and yeast. We also considered two datasets with large size and/or large dimensions and/or large number of classes: letters and USPS handwritten digit data [START_REF] Hull | A database for handwritten text recognition research[END_REF]. See Table 3 for some statistics, including the number of samples n, the dimension d, and the number of classes ℓ. In order to provide a study that can be comparable with previous work, we considered a configuration given in [START_REF] Hsu | A comparison of methods for multi-class support vector machines[END_REF], summarized as follows. Training data were normalized into the range [-1, 1] with the appropriate scaling and bias factors, which were applied to the unlabelled test data. We also applied the same kernel function, which is the Gaussian kernel. To estimate the classification error, a ten-fold cross-validation was used, with parameters optimized by grid search over γ ∈ {2 -4 ; 2 -3 ; • • • ; 2 3 ; 2 4 } and σ ∈ {2 -4 ; 2 -3 ; • • • ; 2 3 ; 2 4 }. This configuration is very similar to the one given in [START_REF] Rifkin | In defense of one-vs-all classification[END_REF], where the authors opted for optimizing separately both parameters, thus sub-optimal as opposed to our joint optimization over both γ and σ.

In Table 4, we give the accuracy rate for each machine, where the values were averaged over 10 Monte Carlo simulations. For a comparative study to other machines, the same partitions for cross-validation were applied to conventional OvA and OvO schemes with LSM machines, as well as vector-output SVM as defined in [START_REF] Szedmak | Learning via linear operators: Maximum margin regression; multiclass and multiview learning at one-class complexity[END_REF]. We borrowed the last column of the Table 4 from [START_REF] Huang | Extreme learning machine for regression and multiclass classification[END_REF], where ANN extreme learning machines were applied, the optimal parameters being estimated on a wider grid search {2 -24 ; 2 -23 ; • • • ; 2 24 ; 2 25 }. We also borrowed 5 columns of the Table 4 from [START_REF] Hsu | A comparison of methods for multi-class support vector machines[END_REF], with essentially the same configuration while requiring more sophisticated algorithms to train these SVM machines. The proposed single matrix inversion scheme is as competitive as these machines in terms of accuracy.

The computation time was estimated over the ten-fold cross-validation with the grid search for the optimal parameters identification, as illustrated in Table 5 where values were averaged over 10 Monte Carlo simulations. The computational time was estimated on a Matlab 64-bit running on a Macbook Pro laptop10 .

As illustrated with the average running time, it is obvious that the proposed approach highly boosts the speed of the multiclass machine.

We also studied the proposed approach on two large scale datasets: letters (from the UCI Repository) and USPS handwritten digit data [START_REF] Hull | A database for handwritten text recognition research[END_REF]. The USPS dataset consists of a large number of scans of the ℓ = 10 digits, in a 16 × 16 grey level images. We considered the same settings as in [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF], with pixel intensities normalized within [-1, 1]. As stated in [START_REF] Huang | Extreme learning machine for regression and multiclass classification[END_REF], both the letters and the USPS datasets require high-performance computers for classification. In order to give a comparative study with different multiclass classification methods, we considered several reduced datasets, with n = 1 000 and n = 2 000 samples, where the USPS images were resampled into d = 8×8 images. The classification error was estimated using a ten-fold cross-validation, with parameters optimized by the same aforementioned grid search. Table 6 shows the estimated classification error and the computational time. It is easy to see that the classical schemes, applying OvA or OVO, are inappropriate for such large scale datasets. The multiclass SVM algorithm proposed by Szedmak et al. in [START_REF] Szedmak | Learning via linear operators: Maximum margin regression; multiclass and multiview learning at one-class complexity[END_REF] is more suitable for this task, still it is outperformed by all the methods proposed in this paper.

Statistical test

Several statistical tests have been derived in the literature to address the problem of comparing two classifiers. In [START_REF] Dietterich | Approximate statistical tests for comparing supervised classification learning algorithms[END_REF], five statistical tests are analyzed, including a new one, the 5 × 2 cv t-test based on 5 iterations of 2-fold cross-validation. The latter test has low type I error and is more powerful than other tests including the McNemar's test and the classical two-tailed t test. This study is extended further in [START_REF] Alpaydin | Combined 5 x 2 cv f test for comparing supervised classification learning algorithms[END_REF],

where a combined 5 × 2 cv F -test is proposed.

The combined 5×2 cv F -test for comparing any two classification algorithms is described as follows. Five replications of two-fold cross-validations are performed, where each replication partitions the data into two halves: one half is used in training 11 and the other one in testing, and vice-versa. Let p (j) i be the difference in the error rates of the two classifiers in fold j of the i-th iteration, for j = 1, 2 and i = 1, 2, . . . , 5. By

denoting p i = 1 2 (p (1)
i + p

i ) the mean, and

s 2 i = (p (1)
ip i ) 2 + (p

i -p i ) 2 the variance, then We analyzed the influence of the label on the performance, and compared all the machines proposed in this paper. We considered the application of the above F -test, and found that there is no significative difference (at 95% level). We also compared the multiclass SVM as derived by Szedmak et al., where a biased model is considered [START_REF] Szedmak | Learning via linear operators: Maximum margin regression; multiclass and multiview learning at one-class complexity[END_REF]. We found that the performance of their method highly depends on the label coding, a property confirmed by the above F -test. Table 4 shows that their method has, in almost all cases, worst performance than all the methods proposed in this paper.

Conclusion

In this paper, our main thesis was that simple multiclass classification machines can be designed with the same complexity as a single binary classifier. Several algorithms were developed for this purpose. We also provided a theoretical study on the coding of the labels, and showed that several well-known labelbooks are equivalent. Experiments conducted on well-known datasets show that the resulting machines are faster than classical ones, and performs just as well as other machines in the literature.

As for future work, we plan to extend the proposed approach to the one-versus-one scheme. While this is less straightforward, one can still take advantage of classical linear algebra to reduce its computational complexity. Also, it would be desirable to derive an optimal labelbook for a given classification task, for instance by optimizing a kernel function applied on the labels. We also plan to consider an optimization on symmetric positive-definite matrices manifold, in the same spirit of [START_REF] Laanaya | Learning general gaussian kernel hyperparameters of svms using optimization on symmetric positive-definite matrices manifold[END_REF]. The sparsity of the solution is also of great interest, in the same spirit as the framework of compressed sensing, with ℓ 1 or even ℓ 0 norms as opposed to the ℓ 2 norm given in this paper.

Figure 2 :

 2 Figure 2: Illustration of the (vector) labels in ℓ = 3 dimensions, for well-known labelbooks.

s 2 i

 2 is approximately F distributed with 10 and 5 degrees of freedom. The null hypothesis that the two algorithms have the same performances is rejected, at statistical significance level of 0.95, if f > 4.74.

Table 1 :

 1 Expressions of well-known labelbooks.

Table 3 :

 3 Databases under investigation, where n denotes the number of available samples with d attributes (dimensions) and ℓ classes.

		n	ℓ	d
	iris	150	3	4
	wine	178	3	13
	glass	214	6	13
	vowel	528 11	10
	yeast	1 484 10	8
	letters up to 2 000 26	16
	USPS	1 000 10 up to 64

Table 4 :

 4 A comparative study of the classification error of several vector-output algorithms, using the same configuration as in Table4(the errors are given in the format mean| worst best ).

	iris	wine	glass	vowel	yeast

Table 5 :

 5 A comparative study of the estimated computation time for several multiclass algorithms, in hh:mm:ss.

		least-squares machines	vector-output machines
		OvA	OvO oneLSM vo-LSSVM vo-RLS vo-SVM
	iris	1:42	1:44	38	43	1:06	1:00
	wine	3:41	3:39	1:14	1:19	2:31	9:32
	glass	7:10	18:22	1:22	1:27	2:45	12:12
	vowel 1:48:28 8:56:00	10:32	11:48	40:51	45:35
	yeast 8:26:45 37:36:24	54:29	1:07:09	3:21:06 2:40:40
				see Sec. 5		see Sec. 4

Table 6 :

 6 A comparative study of the classification error and computation time (in hh:mm:ss) between several multiclass algorithms for the letters (ℓ = 26 classes) and USPS (ℓ = 10 classes) datasets.

			USPS (ℓ = 10)	letters (ℓ = 26)
		n =	1 000	2 000	2 000
	this paper	d = oneLSM vo-LSSVM vo-SVM	8 × 8 6.50 (7:42) 8.41 (7:54) 8.80 (10:02)	8 × 8 4.59 (46:32) 6.19 (46:02) 6.54 (41:23)	16 10.69 (25:27) 14.08 (29:08) 14.08 (27:15)
		LS OvA	6.50 (1:06:36) (7:04:15) 4.59	10.69 (10:29:03)
		LS OvO	6.50 (4:56:37) (31:06:11) 4.59	10.69 (133:10:22)
		[26]	16.61 (9:42)	12.69 (41:52)	15.04 (28:23)

Short of one class versus all the rest, to be more specific.

One may also consider the Fisher coding, with y = -(n/n -) or y = n/n + , where n -and n + are the number of samples in each class.

This is the Representer Theorem[START_REF] Kimeldorf | Some results on tchebycheffian spline functions[END_REF], well known in regularized optimization problems, which states that the solution is a linear combination of the input training data. It also applies to the kernel-based formulation[START_REF] Schölkopf | A generalized representer theorem[END_REF], where H is a reproducing kernel Hilbert space associated to some kernel κ.

In principle, ℓ functions can encode up to 2 ℓ different classes. However, the design of an optimal coding matrix for such task for the maximum number of classes requires prior information on the samples, often unavailable in practice. See for instance[START_REF] Voloshynovskiy | Information-theoretic multiclass classification based on binary classifiers[END_REF] for an information-theoretic point of view. This study is out of scope of this paper, where we consider ℓ functions to encode ℓ classes.

To be more precise, f (•) is an ℓ-tuple of the discriminant functions.

The squared norm of any y ∈ Y is y 2 = y ⊤ y, for instance ℓ for the (±1)-label. See the last column in Table1for other labelbooks.

In[START_REF] Huang | Extreme learning machine for regression and multiclass classification[END_REF], the Vowel dataset contains 528 + 462 instances for training+test, as opposed to the dataset studied in several studies[START_REF] Hsu | A comparison of methods for multi-class support vector machines[END_REF][START_REF] Rifkin | In defense of one-vs-all classification[END_REF], including this paper.[START_REF] Aladjem | Multiclass discriminant mappings[END_REF] In[START_REF] Huang | Extreme learning machine for regression and multiclass classification[END_REF], the Vowel dataset contains 528 + 462 instances for training+test, as opposed to the dataset studied in several studies[START_REF] Hsu | A comparison of methods for multi-class support vector machines[END_REF][START_REF] Rifkin | In defense of one-vs-all classification[END_REF], including this paper.

Mackbook Pro with a 2.53 GHz Intel Core 2 Duo processor and 4 GB RAM.

As recommended in[START_REF] Cantú-Paz | An empirical comparison of combinations of evolutionary algorithms and neural networks for classification problems[END_REF], the training set is also used for estimating the optimal parameters. To this end, we used a ten-fold cross-validation on the training set, with parameters optimized by grid search over γ ∈ {2 -4 ; 2 -3 ; • • • ; 2 3 ; 2[START_REF] El-Yaniv | Better multiclass classification via a margin-optimized single binary problem[END_REF] } and σ ∈ {2 -4 ; 2 -3 ; • • • ; 2 3 ; 2 4 }. The testing set is not considered at all by the algorithms in the training stage.

✩ This work was partly supported by ANR-08-SECU-013-02 VigiRes'Eau.


Maya Kallas
email: maya.kallas@utt.fr

Paul Honeine
email: paul.honeine@utt.fr

Clovis Francis
email: cfrancis@ul.edu.lb

Hassan Amoud
email: hassan.amoud@gmail.com

Kernel autoregressive models

Keywords: kernel machines, autoregressive model, time series prediction, Yule-Walker equations, pre-image problem

come

Introduction

The autoregressive (AR), or linear predictive, model is pervasive in science and technology, with an essential role in the analysis of time series in applications ranging from financial forecasting, to meteorological analysis, to speech processing. For instance to maintain a phone conversation, every cell phone estimates a linear model every 20 milliseconds [START_REF] Dutoit | Applied Signal Processing[END_REF]. The underlying mathematics that govern the AR model are the Yule-Walker equations. The scientific community has made an ever-growing investment to master these equations for the linear prediction [START_REF] Chevalier | Autoregressive models with time-dependent log area ratios[END_REF]. The Yule-Walker equations are the building block of the linear AR model, connecting its parameters to the covariance function of the process. The model parameters are therefore estimated from the covariances of the time series. Forecasting can be considered by applying the resulting predictive model. However, the linearity assumption is often insufficient to explain nonlinear phenomena.

A first attempt to derive a nonlinear Yule-Walker like procedure for a specific nonlinear, high-order, model is given in [START_REF] Cusani | Parameter identification of frequency-selective noisy fast-fading rayleigh digital channels via nonlinear yule-walker-like equations[END_REF]. Nevertheless, up to our knowledge, there is no work that combine the power of the Yule-Walker equations with the proliferating kernel-based methods.

Kernel machines are essentially based on a nonlinear transformation of the data, by using a mapping function from the input space to some feature space, prior to applying a linear procedure in the latter space [START_REF] Shawe-Taylor | Kernel Methods for Pattern Analysis[END_REF]. Nevertheless, it is not necessary to explicitly define the nonlinear transformation, but implicitly by considering a (positive semi-definite) kernel function. The use of kernel machines has received considerable attention since Vapnik's Support Vector Machines (SVM) [START_REF] Vapnik | Statistical Learning Theory[END_REF]. Many nonlinear techniques have been derived, such as the kernel principal component analysis, kernel Fisher discriminant analysis, and SVM novelty detection, only to name a few. In the same spirit, some kernel-based methods were considered for the analysis and prediction of time series data [START_REF] Richard | Online prediction of time series data with kernels[END_REF], including the SVM regression and kernel Kalman filter [START_REF] Ralaivola | Time series filtering, smoothing and learning using the kernel kalman filter[END_REF].

In this paper, we derive nonlinear prediction models by taking full advantage of the Yule-Walker equations. This leads to the estimation of the model parameters by using lagged expected kernels. It is worth noting that the concept of expected kernels has shown its efficiency in recent research [START_REF] Anderson | Channel-robust classifiers[END_REF][START_REF] Anderson | Expected kernel for missing features in support vector machines[END_REF]. Two models are under investigation in this paper.

The first model is based on the underlying concept behind kernel machines, namely mapping data from an input space to a feature space. By operating an AR model on the images of the samples, prediction is defined in the feature space. To provide the predicted sample, one needs to get back to the input space, namely to the space of samples. This is the pre-image problem, with a solution (sample in the input space) that has an image as close as possible to the predicted feature (in the feature space). By following recent developments in the resolution of this ill-posed problem [START_REF] Honeine | The pre-image problem in kernel-based machine learning[END_REF], we derive an iterative technique to provide a fine-tuned solution to this problem. We propose to bypass the pre-image problem, by deriving another model. In the second model, we propose an hybrid formulation, as a tradeoff considering the computational time and the precision, compared to the iterative, fine-tuned, model.

The rest of the paper is organized as follows: In the next section, we introduce the linear AR model and present the Yule-Walker equations for estimating the model parameters, and give the main idea behind kernel machines in Section 3. The first model is derived in 4, by applying the AR model on the images of the samples, and solving the pre-image problem to interpret the prediction in the input space. Section 5 provides pre-image-free techniques, by deriving an AR model on the kernel values. Finally, section 7 illustrates the efficiency of the proposed models on several time series data, and provides a comparative study with well-known prediction methods.

The Yule-Walker equations of the linear autoregressive model

The linear AR model defines each sample as a linear combination of previous samples. Let x 1 , x 2 , . . . , x n be a time series, the p-order AR model is described by

x i = p j=1 α j x i-j + ε i , (1)
for i = p + 1, . . . , n, and where ε i is the unfitness error, often assumed white Gaussian with zero mean. Figure 1 illustrates the concept of the AR model. The parameters α 1 , α 2 , . . . , α p are directly connected with the covariance function of the process. One can therefore determine these parameters from the autocorrelation function. This is the essence of the Yule-Walker equations, as illustrated here.

[Figure 1 about here.]

Let the data be centered, thus let µ be the expectation of x i , namely,

µ = IE[x i ],
where IE[•] is the expectation1 . If we apply the expectation on each side of (1), we get that (1 -

p j=1 α j )µ = IE[ε i ].
For any positive lag τ , we can evaluate the autocorrelation function of each time series. Let r be the empirical counterpart of the autocorrelation function of the time series, then r(τ ) = p j=1 α j r(τ -j), for any lag τ ≥ 1. Since the autocorrelation function is even, i.e., r(-τ ) = r(τ ), we obtain the matrix form of the Yule Walker equations

r = R α, where r = [r(1) r(2) • • • r(p)] , α = [α 1 α 2 • • • α p ] ,
and

R =             r(0) r (
            .
Assuming that the p×p symmetric matrix R is invertible, the coefficients α are estimated by α = R -1 r. Once the coefficients are estimated, the AR model can be applied to predict future samples, with

x k = p j=1 α j x k-j .
While this technique is easy to implement, it is not adapted for nonlinear systems.

Next, we derive Yule-Walker-like equations for nonlinear models, within the framework of kernel machines. But before, we prepare the ground by briefly describing the main idea behind the kernel machines.

Kernel machines

A kernel is a symmetric and continuous function defined by κ : X × X → IR, where

X is an input space. If the kernel verifies i,j α i α j κ(x i , x j ) ≥ 0 for all α i , α j ∈ IR and all x i , x j ∈ X , then the kernel is positive semi-definite. The Moore-Aronszajn theorem [START_REF] Aronszajn | Theory of reproducing kernels[END_REF] states that each positive semi-definite kernel defines a unique (up to an isometry) feature space and vice-versa. This feature space H is obtained using a mapping function κ(x i , x j ) = Φ(x i ), Φ(x j ) , for any x i , x j ∈ X , where • , • denotes the corresponding inner product in H.

Table 1 summarizes some of the most commonly used kernel functions, grouped into two classes: projective kernels written in terms of dot product x i • x j and radial kernels in terms of Euclidean distance x i -x j . It is worth noting that some kernels induce infinite-dimensional feature spaces, such as the Gaussian kernel.

[Table 1 about here.]

Autoregressive Φ-model

To derive a nonlinear autoregressive model, a straightforward approach consists of applying a nonlinear transformation on the samples, prior to the autoregressive model. This is illustrated in Figure 2, where the samples are mapped from the input space to a feature space using a nonlinear transformation. An AR model is applied to predict a feature in the feature space. Finally, to predict a future sample, one needs to get back to the input space, by solving the so-called pre-image problem. These two stages are derived next.

[Figure 2 about here.]

Autoregressive model in feature space

Let us start by evaluating the autoregressive model in the feature space. To this end, let κ(•, •) be a (reproducing) kernel, defining a nonlinear map Φ(•) from the input space X to some feature space H. Thus, each sample x i is mapped to its corresponding Φ(x i ).

Therefore, to a time series, given by the samples x 1 , x 2 , . . . , x n ∈ X , corresponds the samples Φ(x 1 ), Φ(x 2 ), . . . , Φ(x n ) in H. By applying the AR model ( 1) for the latter, we can write

Φ(x i ) = p j=1 α j Φ(x i-j ) + ε Φ i , (2)
where these terms are given in the feature space, including ε Φ i ∈ H which represents the unfitness of the model. While the samples x i 's are assumed zero-mean in the input space, this is not the case for the Φ(x i )'s in the feature space.

Let µ Φ denote the expectation of the latter, namely

µ Φ = IE[Φ(x i )].
On the one hand, by applying the expectation on both sides of (2), we get (1 -

p j=1 α j )µ Φ = IE[ε Φ i ],
where the process is assumed stationary. On the other hand, we have

Φ(x i ) -µ Φ = p j=1 α j Φ(x i-j ) + ε Φ i -µ Φ = p j=1 α j Φ(x i-j ) -µ Φ + ε Φ i -1 - p j=1 α j µ Φ .
By combining these results, and by taking the inner product (in the feature space) of both sides of the above equation with (Φ(x i-τ ) -µ Φ ), for some positive lag τ , we get

Φ(x i ) -µ Φ , Φ(x i-τ ) -µ Φ = ε Φ i -IE[ε Φ i ], Φ(x i-τ ) -µ Φ + p j=1 α j Φ(x i-j ) -µ Φ , Φ(x i-τ ) -µ Φ . (3)
By analogy with the linear AR case, we assume that the noise ε Φ i and Φ(x i-τ ) are uncorrelated for every positive lag τ . Therefore, taking the expectations of expression

(3) and assuming a stationarity sequence, we get for any τ ≥ 1:

IE[κ c (x i , x i-τ )] = p j=1 α j IE[κ c (x i-j , x i-τ )], (4)
where

κ c (•, •) is the centered version of the kernel κ(•, •), defined with κ c (x i , x j ) = Φ(x i ) -µ Φ , Φ(x j ) -µ Φ .
Finally, we considered all the lag values, and write the expression (4) in matrix form,

r κc = R κc α,
where

r κc = IE[κ c (x i , x i-1 )] IE[κ c (x i , x i-2 )] • • • IE[κ c (x i , x i-p )] ,
and R κc is the matrix described by the expected kernels with

R κc =             IE[κ c (x i , x i )] IE[κ c (x i , x i-1 )] • • • IE[κ c (x i , x i-p+1 )] IE[κ c (x i , x i-1 )] IE[κ c (x i , x i )] • • • IE[κ c (x i , x i-p+2 )] . . . . . . . . . . . . IE[κ c (x i , x i-p+1 )] IE[κ c (x i , x i-p+2 )] • • • IE[κ c (x i , x i )]            
.

The vector of coefficients α is obtained by inverting the matrix R κc , with

α = R -1 κc r κc .
In practice, the expectations are estimated over a set of n available samples (See [START_REF] Cortes | Algorithms for learning kernels based on centered alignment[END_REF] for a recent review). The centered version of the kernel is evaluated using

κ c (x i , x j ) = κ(x i , x j ) - 1 n n k=1 κ(x i , x k ) - 1 n n k=1 κ(x j , x k ) + 1 n 2 n k,k =1 κ(x k , x k ).
Once the model parameters are determined on a set of n available samples, one can predict by

ψ i = p j=1 α j Φ(x i-j ), (5)
where the predicted ψ i lies in the feature space. Still, one is more interested in the predicted sample in the original input space. Thus, we need to map back ψ i from the feature space to the input space, namely to the so-called pre-image.

The exact pre-image may not exist in general, and even if it exists, it may not be unique. This is the ill-posed pre-image problem, where one identifies the best x * i in the input space whose image Φ(x * i ) is as close as possible to ψ i . See [START_REF] Honeine | The pre-image problem in kernel-based machine learning[END_REF] for a recent review. We propose to solve the optimization problem by minimizing the distance between elements in the feature space, in the same spirit as in [START_REF] Kallas | Nonlinear feature extraction using kernel principal component analysis with non-negative pre-image[END_REF], namely

x * i = arg min x ψ i -Φ(x) 2 .
We propose to use the iterative fixed-point method. By injecting the autoregressive model (5) into the above expression, we obtain the following optimization problem for any x i :

x * i = arg min x p j=1 α j Φ(x i-j ) -Φ(x) 2 .
This optimization problem can be written as

x * i = arg min x J i (x),
where J i (x) is the cost function defined by

J i (x) = - p j=1 α j κ(x i-j , x) + 1 2 κ(x, x), (6)
where the term independent of x has been removed.

Consider the case of the Radial Basis Functions, with kernels of the form

κ(x j , x j ) = f ( x j -x j 2 ). ( 7
)
The gradient of the function κ(x i-j , x) with respect to x is given by

∇ x κ(x i-j , x) = 2 (x i-j -x) f (1) ( x i-j -x 2 ),
where f (z) is the first derivative of f (•) w.r.t z. By combining this expression with the gradient of the cost function J i (x), we get

∇ x J i (x) = 2 p j=1 α j (x i-j -x) f (1) ( x i-j -x 2 ).
Such expression simplifies further for several kernel functions, such as the Gaussian kernel

with f (ζ) = exp -1 2σ 2 ζ , thus f (1) (ζ) = -1 2σ 2 f (ζ).
Let us now write the pre-image using a linear combination of the available data, that is x * i = p j=1 δ * j x i-j . Many pre-image techniques from the literature [START_REF] Kwok | The pre-image problem in kernel methods[END_REF][START_REF] Honeine | Solving the pre-image problem in kernel machines: a direct method[END_REF][START_REF] Zheng | Penalized preimage learning in kernel principal component analysis[END_REF] have validated this statement. We prove this statement for the radial and projective kernels.

Theorem 1. Any pre-image x * i can be written as a linear combination of the available data, namely

x * i = p j=1 δ * j x i-j
for some weights δ * j .

Proof. First, we study the class of radial kernels, defined by expression [START_REF] Ralaivola | Time series filtering, smoothing and learning using the kernel kalman filter[END_REF]. In such case, the term ∂κ(x, x)/∂x vanishes. The gradient at the optimum can be written as

p j=1 α j ∂κ(x i-j , x * i ) ∂x * i = 0, 10
with the left-hand-side given as

p j=1 α j ∂κ(x i-j , x * i ) ∂x * i = p j=1 α j ∂f ( x i-j -x * i 2 ) ∂( x i-j -x * i 2 ) 2(x * i -x i-j ).
As a consequence, the final result can be expressed as

x * i = p j=1 α j f (1) ( x i-j -x * i 2 ) p k=1 α k f (1) ( x i-k -x * i 2 ) x i-j , thus of the form x * i = p j=1 δ * j x i-j .
We now study the projective kernels, of the form g(x i • x j ). In this case, at the optimum, the gradient is written as

p j=1 α j ∂κ(x i-j , x * i ) ∂x * i = 1 2 ∂κ(x * i , x * i ) ∂x * i .
We evaluate the right-hand-side and the left-hand-side respectively using the form of the kernel, and we combine both expressions to get

x * i = p j=1 α j g (1) (x i-j • x * i ) g (1) (x * i • x * i ) x i-j , of the form x * i = n i=1 δ * j x i-j .

Autoregressive kernel-based models

The Φ-model proposed above requires to solve the ill-posed pre-image problem, for each predicted sample. In this section, we propose to overcome this drawback. To this end, we investigate the estimation of a function defined on X p , where the desired output is x i and the input x i represents the p previous sample with

x i = [x i-1 x i-2 • • • x i-p ] .
Having X p as the input space, we consider kernels on X p ×X p . Taking into consideration this formulation, we propose two different nonlinear models, and derive the corresponding Yule-Walker equations.

It is worth noting that, for the Gaussian kernel, this approach corresponds to considering the kernel

κ(x i-j , x j ) = exp - 1 2σ 2 x i-j -x i 2 ,
where the distance between two vectors x i-j and x i is given as

x i-j -x i 2 = p k=1 x i-j-k -x i-k 2 .
It is easy to see that this leads to the following expression for the Gaussian kernel

κ(x i-j , x i ) = p k=1 κ x i-j-k , x i-k ,
thus connecting the multivariate kernel to the univariate one.

A hybrid autoregressive model

As opposed to the above method, where we apply the autoregressive model to the Φ-images in the feature space, we consider here an autoregressive model on the kernel values. The proposed model is defined by predicting the sample x * i . We investigate the estimation of a function ϕ(•), with ϕ : X p → X . The proposed model is defined by

ϕ(x i ) = p j=1 β j κ(x i-j , x i ) x i-j + ε i , (8)
such that ϕ(x i ) = x * i is the value of the predicted sample in the input space and x i represents the p previous samples with

x i = [x i-1 x i-2 • • • x i-p ] . See Section 5.2 for
the motivation behind this model.

Let µ x be the expectation of ϕ(x i ), and µ xj be the expectation of κ(x i-j , x i ) x i-j , for any j = 1, 2, . . . , p, namely

µ x = IE[ϕ(x i )] µ xj = IE[κ(x i , x i-j ) x i-j ].
we get

IE[ ϕ(x i ) -µ x , ϕ(x i-τ ) -µ x ] = p j=1 β j IE[ κ(x i , x i-j ) x i-j -µ xj , ϕ(x i-τ ) -µ x ]. (9)
for any lag τ = 1, 2, . . . , p. By considering all the lag values, we get the matrix form

r κ = R κ x β
, where r κ is defined by

r κ = IE[ ϕ(x i ) -µ x , ϕ(x i-1 ) -µ x ] • • • IE[ ϕ(x i ) -µ x , ϕ(x i-p ) -µ x ] (10)
and R κ x is the matrix described by

Rκ x =             IE[ κ(xi, xi-1)xi-1 -µx1, ϕ(xi-1) -µx ] . . . IE[ κ(xi, xi-p)xi-p -µxp, ϕ(xi-1) -µx ] IE[ κ(xi, xi-1)xi-1 -µx1, ϕ(xi-2) -µx ] . . . IE[ κ(xi, xi-p)xi-p -µxp, ϕ(xi-2) -µx ] . . . . . . . . . IE[ κ(xi, xi-1)xi-1 -µx1, ϕ(xi-p) -µx ] . . . IE[ κ(xi, xi-p)xi-p -µxp, ϕ(xi-p) -µx ]            
The vector of coefficients β is obtained by inverting the matrix R κ , with β = R -1 κ x r κ .

Once the β parameters are evaluated, the hybrid autoregressive model can predict future samples by applying x * i = ϕ(x i ), namely

x * i = p j=1 β j κ(x i-j , x i ) x i-j . (11)
With this proposed method, we evaluate directly the predicted value without the search for a function in the feature space, where the evaluation is in the input space.

To this end, we lose the precision we obtain in the feature space, since we evaluate in the input space. In this case, we do not obtain a more precise value for the prediction.

However, we gain in time when applying the hybrid technique since we do not need to solve the pre-image problem in order to obtain the predicted value.

Connecting the hybrid model with the aforementioned model

The main motivation behind the hybrid model is its connection with the aforementioned model, as illustrated here. We consider the Φ-model, as derived in Section 4.

Theorem 1 states that the predicted sample at instant i takes the form

x * i = p j=1 δ * j x i-j .
By restricting the solution to δ * j = β j κ(x i-j , x i ), we get the hybrid model defined in [START_REF] Aronszajn | Theory of reproducing kernels[END_REF]. More precisely, by taking the result of the Theorem 1 for the Gaussian kernel, the coefficients β j will be β j = α j

1 p k=1 α k f (1) ( x i-k -x * i 2 ) .

Autoregressive model on the kernel values

One may also consider the autoregressive model to the Φ-images in the feature space.

Therefore, we consider here an autoregressive model on the kernel values. The proposed model is defined by predicting the sample x * i , with

ϕ(x i ) = p j=1 β j κ(x i-j , x i ) + ε i , (12)
where x i represents the p previous sample with

x i = [x i-1 x i-2 • • • x i-p ] , and ϕ(•)
is the function defining the prediction of the future sample such as ϕ(x i ) = x i .

Let µ j the expectation of κ(x i-j , x i ), for each j = 1, 2, . . . , p, namely

µ j = IE[κ(x i , x i-j )],
and

µ x = IE[ϕ(x i )]
. By following the developments given in Section 4, and by analogy with expression (4), we get, for any lag τ = 1, 2, . . . , p

r κ = R κ β,
where r κ defined by [START_REF] Honeine | The pre-image problem in kernel-based machine learning[END_REF] and R κ is the matrix described by the expected kernels with

R κ =             IE[ κ(x i , x i-1 ) -µ 1 , ϕ(x i-1 ) -µ x ] . . . IE[ κ(x i , x i-p ) -µ p , ϕ(x i-1 ) -µ x ] IE[ κ(x i , x i-1 ) -µ 1 , ϕ(x i-2 ) -µ x ] . . . IE[ κ(x i , x i-p ) -µ p , ϕ(x i-2 ) -µ x ] . . . . . . . . . IE[ κ(x i , x i-1 ) -µ 1 , ϕ(x i-p ) -µ x ] . . . IE[ κ(x i , x i-p ) -µ p , ϕ(x i-p ) -µ x ]            
.

The vector of coefficients β is obtained by inverting the matrix R κ , with

β = R -1 κ r κ .
This model provides a coarse representation of the kernel values and does not require any pre-image solving technique. In practice, this AR model performs poorly. For this reason, it will be discarded from the experimental part of this paper.

Study on the complexity for each of the aforementioned techniques

In this section, we study the complexity of the AR methods described in this paper.

We consider a set of n samples for the training stage and another set of n samples for the testing stage. For both models, we have to evaluate the centered version of the version, the computational complexity needed to center such kernel remains low, and it can be determined as in iterative/online methods with low computational cost, see Annexe A. 3 in [START_REF] Kim | Kernel Hebbian algorithm for iterative kernel principal component analysis[END_REF] for more details.

Let us first start with the AR model in the feature space. In order to evaluate the coefficients, we have to invert a matrix R κ x and multiply it by the vector r κ . The dimension of the matrix is p × p, where p is usually a small. The complexity of the inversion is p 3 . Once the parameters are estimated, and the coefficients are evaluated, we have to estimate the pre-image using an iterative technique.

We now consider the hybrid model. In order to evaluate the coefficients, we have to invert a matrix R κ x and multiply it by the vector r κ . The dimension of the matrix is p × p, where p is usually a small. The complexity of the inversion is p3 . Once the parameters are estimated, and the coefficients are evaluated. We can directly predict future samples by applying the equation [START_REF] Aronszajn | Theory of reproducing kernels[END_REF]. To this end, we have to evaluate the kernel and multiply it by a vector of the p previous sample. Therefore, the complexity here is n × p.

Experiments

In this section, we study the relevance of the different models proposed in this paper.

We compare them with the classical linear AR model and nonlinear models such as the multilayer perceptron and the support vector regression. To provide a benchmark study, we considered the same experimental settings as in [START_REF] Ralaivola | Time series filtering, smoothing and learning using the kernel kalman filter[END_REF]. Four time series are under investigation: two chaotic time series2 [START_REF] Mukherjee | Nonlinear prediction of chaotic time series using support vector machines[END_REF], an electrocardiogram (ECG) signal 3 and an electromyogram (EMG) signal having a neuropathy problem 4 .

The Mackey-Glass univariate time series models the blood cells production evolution.

It is defined by the delay differential equation

dx(t) dt = -0.1 x(t) + 0.2 x(t -τ ) 1 + x(t -τ ) 10
which, for values of τ greater than 16.8, shows some highly nonlinear chaotic behavior.

In our case, we set τ = 30, and denote the time series by M G 30 .

equations                dx(t) dt = -ax + ay dy(t) dt = -xz + rx -y dz(t) dt = xy -bz
For our simulations, we set a = 10, r = 28 and b = 8/3. This multivariate system is decomposed into the estimation of three models, for x, y, and z, separately.

The ECG signal is ECG recordings of subjects referred to the Arrhythmia Laboratory where subjects included in this database were found to have had no significant arrhythmias. The ECG signals are considered to be stationary, for they are taken in a 1 minute of time.

The EMG signal is EMG recording taken from a clinical test used to assess function of muscles and the nerves that control them. The EMG studies are used to help in the diagnosis and management of disorders such as the muscular dystrophies and neuropathies. Nerve conduction studies that measure how well and how fast the nerves conduct impulses are often performed in conjunction with EMG studies.

Each time series is decomposed into two parts. The first n = 300 samples are used in the learning stage, for estimating the optimal value of the order p from p ∈ [1, 2, . . . , 5],

the coefficients in the AR expansion, and the bandwidth σ of the Gaussian kernel. The next n = 300 samples are used to evaluate the relevance of the resulting model, by considering the mean square error (MSE) estimated with

ε err = 1 n 2n i=n+1 x * i -x i 2 ,
where x * i is the predicted value at instant i, and x i is the true value of the time series at the same time.

Table 2 gives the mean square error and the computational time measured on an Intel Core 2, with a speed of 2.40 GHz and a random access memory of 1.00 GB, for each of the three proposed nonlinear AR models, as well as the linear model. For the preimage technique, the stopping criterion is given by a lower bound on the tolerance, set to 10 -6 , while limiting the maximum number of iterations to 50. It is obvious that the AR Φ-model with pre-image technique presents the best MSE, however such fine-tuning requires significant computational resources. The hybrid model presents a good tradeoff between accuracy and computational complexity, and outperforms the linear AR model.

The results obtained for the three time series are illustrated in Figure 3.

The experimental settings are similar to the ones given in [START_REF] Ralaivola | Time series filtering, smoothing and learning using the kernel kalman filter[END_REF] for Lorenz and M G 30 time series. This allows us to provide a comparative study with different nonlinear techniques provided in the aforementioned paper, such as multilayer perceptron, support vector regressor, and a kernel Kalman filter. Table 3 shows the MSE evaluated on each time series using different nonlinear prediction technique. As we can see, the kernel AR with pre-image presents the best MSE for the M G 30 , and the hybrid AR model presents the best MSE for the Lorenz attractor time series. An adequate model is contributed to each time series depending on the nature of the data.

[Table 2 about here.]

[Table 3 about here.]

[Figure 3 about here.]

Conclusion and future work

In this paper, we proposed to derive Yule-Walker equations for nonlinear autoregressive models. To this end, we combined the simplicity of the AR model with the efficiency of the kernel machines in machine learning. Two models were derived, with a third one that can be used. The first model, the AR Φ-model was based on the mapping defined by a kernel, but requires a pre-image technique to get back to the input space. To circumvent the pre-image problem, we considered an hybrid model. We derived Yule-Walker equations for all these AR models, and showed the ease of estimating the model parameters. The relevance of the proposed techniques were illustrated on three time series, and compared to several well-known nonlinear techniques. As we have concluded, the kernel AR with the pre-image presents the best MSE, however it requires significant computational resources. The hybrid model represents a good compromise between efficiency and computational cost, and outrun the linear AR model.

As for future work, we are considering the estimation of the optimal order, by adapting several criteria such as the Akaike Information Criterion, the Bayesian Information

x i x i-1 x i-(p-1) x i-p α 1 α p α p-1

Figure 1 :Figure 2 :Figure 3 : 2 )

 1232 Figure 1: Illustration of the AR model, where the x i is defined by a linear combination of the p previous samples x i-k 's, with weight parameters α 1 , α 2 , . . . , αp.

  Illustration of the AR model, where the x i is defined by a linear combination of the p previous samples x i-k 's, with weight parameters α 1 , α 2 , . . . , α p . 23 2 Illustration of the AR Φ-model: the samples are mapped from X to H, where the AR model is applied. Once ψ i is estimated, one needs to map it back to the input space X , in order to provide the prediction x * i . . . . . 24 3 Visualization of the results obtained when predicting the four time series using the linear AR model (first column), the AR Φ-model with pre-image (second column), and the hybrid AR model (last column). The predicted time series is given in red dash-point (-•) lines, while the blue lines define the original data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

	1
	Criterion and the partial autocorrelation function. Also, we intend to apply other tech-
	niques to estimate the model parameters, such as the Levinson-Durbin method. Finally,
	we will extend the use of the kernel machines for the Autoregressive Moving Average
	model.

Table 1 :

 1 Commonly used reproducing kernels in machine learning, with parameters c, σ > 0, and p ∈ IN +

			M G 30	Lorenz	ECG	EMG-Neuro
	Linear AR model	Time (s)	0.0107	0.0539	0.0126	0.0133
	(with (1))	MSE	0.0655	0.2907	0.0332	0.1397
	AR Φ-model (pre-image) Time (s)	5.3201	10.5241	7.8921	9.0315
	(with (4))	MSE	1 • 10 -5	0.1498	6 • 10 -4	1 • 10 -4
	Hybrid AR model	Time (s)	0.0861	0.6091	0.1834	0.1280
	(with (9))	MSE	0.0623	0.0213	0.0290	0.0061

Table 2 :

 2 Estimated computational time and mean square error (MSE) between the predicted values and the original ones.

In this paper, all expectations are taken on the index i.

The time series are available at http://www.bme.ogi.edu/∼ericwan/data.html.

The electrocardiogram signal is taken from the MIT-BIH Normal Sinus Rhythm Database, and is available at http://physionet.org/physiobank/database/nsrdb/

The electromyogram signal is taken from the Examples of Electromyograms, and is available at http://http://physionet.org/physiobank/database/emgdb/

List of Tables