kernlab – An S4 Package for Kernel Methods in R

kernlab is an extensible package for kernel-based machine learning methods in R 1 . It takes advantage of R ’s new S4 object model and provides a framework for creating and using kernel-based algorithms. The package contains dot product primitives (kernels), implementations of support vector machines and the relevance vector machine, Gaussian processes, a ranking algorithm, kernel PCA, kernel CCA, and a spectral clustering algorithm. Moreover it provides a general purpose quadratic programming solver, and an incomplete Cholesky decomposition method.


Introduction
Machine learning is all about extracting structure from data, but it is often difficult to solve problems like classification, regression and clustering in the space in which the underlying observations have been made.Kernel-based learning methods use an implicit mapping of the input data into a high dimensional feature space defined by a kernel function, i.e., a function returning the inner product ⟨Φ(x), Φ(y)⟩ between the images of two data points x, y in the feature space.The learning then takes place in the feature space, provided the learning algorithm can be entirely rewritten so that the data points only appear inside dot products with other points.This is often referred to as the "kernel trick" (Schölkopf and Smola 2002).More precisely, if a projection Φ : X → H is used, the dot product ⟨Φ(x), Φ(y)⟩ can be represented by a kernel function k k(x, y) = ⟨Φ(x), Φ(y)⟩, (1) which is computationally simpler than explicitly projecting x and y into the feature space H.
One interesting property of kernel-based systems is that, once a valid kernel function has been selected, one can practically work in spaces of any dimension without paying any computational cost, since feature mapping is never effectively performed.In fact, one does not even need to know which features are being used.
Another advantage is the that one can design and use a kernel for a particular problem that could be applied directly to the data without the need for a feature extraction process.This is particularly important in problems where a lot of structure of the data is lost by the feature extraction process (e.g., text processing).The inherent modularity of kernel-based learning methods allows one to use any valid kernel on a kernel-based algorithm.

Software review
The most prominent kernel based learning algorithm is without doubt the support vector machine (SVM), so the existence of many support vector machine packages comes as little surprise.Most of the existing SVM software is written in C or C++, e.g. the award winning libsvm1 (Chang and Lin 2001), SVMlight2 (Joachims 1999), SVMTorch3 , Royal Holloway Support Vector Machines4 , mySVM5 , and M-SVM6 with many packages providing interfaces to MATLAB (such as libsvm), and even some native MATLAB toolboxes7 8 9 .
Putting SVM specific software aside and considering the abundance of other kernel-based algorithms published nowadays, there is little software available implementing a wider range of kernel methods with some exceptions like the Spider10 software which provides a MATLAB interface to various C/C++ SVM libraries and MATLAB implementations of various kernel-based algorithms, Torch11 which also includes more traditional machine learning algorithms, and the occasional MATLAB or C program found on a personal web page where an author includes code from a published paper.

R software
The R package e1071 offers an interface to the award winning libsvm (Chang and Lin 2001), a very efficient SVM implementation.libsvm provides a robust and fast SVM implementation and produces state of the art results on most classification and regression problems (Meyer, Leisch, and Hornik 2003).The R interface provided in e1071 adds all standard R functionality like object orientation and formula interfaces to libsvm.Another SVM related R package which was made recently available is klaR (Roever, Raabe, Luebke, and Ligges 2004) which includes an interface to SVMlight, a popular SVM implementation along with other classification tools like Regularized Discriminant Analysis.
However, most of the libsvm and klaR SVM code is in C++.Therefore, if one would like to extend or enhance the code with e.g.new kernels or different optimizers, one would have to modify the core C++ code.

kernlab
kernlab aims to provide the R user with basic kernel functionality (e.g., like computing a kernel matrix using a particular kernel), along with some utility functions commonly used in kernel-based methods like a quadratic programming solver, and modern kernel-based algorithms based on the functionality that the package provides.Taking advantage of the inherent modularity of kernelbased methods, kernlab aims to allow the user to switch between kernels on an existing algorithm and even create and use own kernel functions for the kernel methods provided in the package.

S4 objects
kernlab uses R's new object model described in "Programming with Data" (Chambers 1998) which is known as the S4 class system and is implemented in the methods package.
In contrast with the older S3 model for objects in R, classes, slots, and methods relationships must be declared explicitly when using the S4 system.The number and types of slots in an instance of a class have to be established at the time the class is defined.The objects from the class are validated against this definition and have to comply to it at any time.S4 also requires formal declarations of methods, unlike the informal system of using function names to identify a certain method in S3.An S4 method is declared by a call to setMethod along with the name and a "signature" of the arguments.The signature is used to identify the classes of one or more arguments of the method.Generic functions can be declared using the setGeneric function.Although such formal declarations require package authors to be more disciplined than when using the informal S3 classes, they provide assurance that each object in a class has the required slots and that the names and classes of data in the slots are consistent.
An example of a class used in kernlab is shown below.Typically, in a return object we want to include information on the result of the method along with additional information and parameters.Usually kernlab's classes include slots for the kernel function used and the results and additional useful information.

Namespace
Namespaces were introduced in R 1.7.0 and provide a means for packages to control the way global variables and methods are being made available.Due to the number of assignment and accessor function involved, a namespace is used to control the methods which are being made visible outside the package.Since S4 methods are being used, the kernlab namespace also imports methods and variables from the methods package.

Data
The kernlab package also includes data set which will be used to illustrate the methods included in the package.The spam data set (Hastie, Tibshirani, and Friedman 2001) set collected at Hewlett-Packard Labs contains data on 2788 and 1813 e-mails classified as non-spam and spam, respectively.The 57 variables of each data vector indicate the frequency of certain words and characters in the e-mail.
Another data set included in kernlab, the income data set (Hastie et al. 2001), is taken by a marketing survey in the San Francisco Bay concerning the income of shopping mall customers.It consists of 14 demographic attributes (nominal and ordinal variables) including the income and 8993 observations.The ticdata data set (van der Putten, de Ruiter, and van Someren 2000) was used in the 2000 Coil Challenge and contains information on customers of an insurance company.The data consists of 86 variables and includes product usage data and socio-demographic data derived from zip area codes.The data was collected to answer the following question: Can you predict who would be interested in buying a caravan insurance policy and give an explanation why?
The promotergene is a data set of E. Coli promoter gene sequences (DNA) with 106 observations and 58 variables available at the UCI Machine Learning repository.Promoters have a region where a protein (RNA polymerase) must make contact and the helical DNA sequence must have a valid conformation so that the two pieces of the contact region spatially align.The data contains DNA sequences of promoters and non-promoters.The spirals data set was created by the mlbench.spiralsfunction in the mlbench package (Leisch and Dimitriadou 2001).This two-dimensional data set with 300 data points consists of two spirals where Gaussian noise is added to each data point.

Kernels
A kernel function k calculates the inner product of two vectors x, x ′ in a given feature mapping Φ : X → H.
The package includes implementations of the following kernels: • the linear vanilladot kernel implements the simplest of all kernel functions which is useful specially when dealing with large sparse data vectors x as is usually the case in text categorization.
• the Gaussian radial basis function rbfdot which is a general purpose kernel and is typically used when no further prior knowledge is available about the data.
• the polynomial kernel polydot which is used in classification of images.
• the hyperbolic tangent kernel tanhdot which is mainly used as a proxy for neural networks.
• the Bessel function of the first kind kernel besseldot is a general purpose kernel and is typically used when no further prior knowledge is available and mainly popular in the Gaussian process community.
• the Laplace radial basis kernel laplacedot which is a general purpose kernel and is typically used when no further prior knowledge is available.
• the ANOVA radial basis kernel anovadot performs well in multidimensional regression problems where x k is the kth component of x.

Kernel utility methods
The package also includes methods for computing commonly used kernel expressions (e.g., the Gram matrix).These methods are written in such a way that they take functions (i.e., kernels) and matrices (i.e., vectors of patterns) as arguments.These can be either the kernel functions already included in kernlab or any other function implementing a valid dot product (taking two vector arguments and returning a scalar).In case one of the already implemented kernels is used, the function calls a vectorized implementation of the corresponding function.Moreover, in the case of symmetric matrices (e.g., the dot product matrix of a Support Vector Machine) they only require one argument rather than having to pass the same matrix twice (for rows and columns).
The computations for the kernels already available in the package are vectorized whenever possible which guarantees good performance and acceptable memory requirements.Users can define their own kernel by creating a function which takes two vectors as arguments (the data points) and returns a scalar (the dot product).This function can then be based as an argument to the kernel utility methods.For a user defined kernel the dispatch mechanism calls a generic method implementation which calculates the expression by passing the kernel function through a pair of for loops.The kernel methods included are: kernelMatrix This is the most commonly used function.It computes k(x, x ′ ), i.e., it computes the matrix K where K ij = k(x i , x j ) and x is a row vector.In particular, K <-kernelMatrix(kernel, x) where the x i are the columns of X and kernelFast This method is different to kernelMatrix for rbfdot, besseldot, and the laplacedot kernel, which are all RBF kernels.It is identical to kernelMatrix, except that it also requires the squared norm of the first argument as additional input.It is mainly used in kernel algorithms, where columns of the kernel matrix are computed per invocation.In these cases, evaluating the norm of each column-entry as it is done on a kernelMatrix invocation on an RBF kernel, over and over again would cause significant computational overhead.Its invocation is via Here a is a vector containing the squared norms of x1.
kernelMult is a convenient way of computing kernel expansions.It returns the vector f = (f (x 1 ), . . ., f (x m )) where The need for such a function arises from the fact that K may sometimes be larger than the memory available.Therefore, it is convenient to compute K only in stripes and discard the latter after the corresponding part of Kα has been computed.The parameter blocksize determines the number of rows in the stripes.In particular, f <-kernelMult(kernel, x, alpha) kernelPol is a method very similar to kernelMatrix with the only difference that rather than computing . This means that K <-kernelPol(kernel, x, y) computes the matrix K ij = y i y j k(x i , x j ) where the x i are the columns of x and y i are elements of the vector y.Moreover, . Both x1 and x2 may be matrices and y1 and y2 vectors.
An example using these functions :

Kernel methods
Providing a solid base for creating kernel-based methods is part of what we are trying to achieve with this package, the other being to provide a wider range of kernel-based methods in R. In the rest of the paper we present the kernel-based methods available in kernlab.All the methods in kernlab can be used with any of the kernels included in the package as well as with any valid user-defined kernel.User defined kernel functions can be passed to existing kernel-methods in the kernel argument.

Support vector machine
Support vector machines (Vapnik 1998) have gained prominence in the field of machine learning and pattern classification and regression.The solutions to classification and regression problems sought by kernel-based algorithms such as the SVM are linear functions in the feature space: for some weight vector w ∈ F .The kernel trick can be exploited in this whenever the weight vector w can be expressed as a linear combination of the training points, w = n i=1 α i Φ(x i ), implying that f can be written as A very important issue that arises is that of choosing a kernel k for a given learning task.Intuitively, we wish to choose a kernel that induces the "right" metric in the space.Support Vector Machines choose a function f that is linear in the feature space by optimizing some criterion over the sample.
In the case of the 2-norm Soft Margin classification the optimization problem takes the form: Based on similar methodology, SVMs deal with the problem of novelty detection (or one class classification) and regression.
kernlab's implementation of support vector machines, ksvm, is based on the optimizers found in bsvm12 (Hsu and Lin 2002b) and libsvm (Chang and Lin 2001) which includes a very efficient version of the Sequential Minimization Optimization (SMO).SMO decomposes the SVM Quadratic Problem (QP) without using any numerical QP optimization steps.Instead, it chooses to solve the smallest possible optimization problem involving two elements of α i because they must obey one linear equality constraint.At every step, SMO chooses two α i to jointly optimize and finds the optimal values for these α i analytically, thus avoiding numerical QP optimization, and updates the SVM to reflect the new optimal values.The SVM implementations available in ksvm include the C-SVM classification algorithm along with the ν-SVM classification formulation which is equivalent to the former but has a more natural (ν) model parameter taking values in [0, 1] and is proportional to the fraction of support vectors found in the data set and the training error.
For classification problems which include more than two classes (multi-class) a one-against-one or pairwise classification method (Knerr, Personnaz, and Dreyfus 1990;Kreßel 1999) is used.This method constructs k 2 classifiers where each one is trained on data from two classes.Prediction is done by voting where each classifier gives a prediction and the class which is predicted more often wins ("Max Wins").This method has been shown to produce robust results when used with SVMs (Hsu and Lin 2002a).Furthermore the ksvm implementation provides the ability to produce class probabilities as output instead of class labels.This is done by an improved implementation (Lin, Lin, and Weng 2001) of Platt's posteriori probabilities (Platt 2000) where a sigmoid function is fitted on the decision values f of the binary SVM classifiers, A and B are estimated by minimizing the negative log-likelihood function.To extend the class probabilities to the multi-class case, each binary classifiers class probability output is combined by the couple method which implements methods for combing class probabilities proposed in (Wu, Lin, and Weng 2003).
Another approach for multIn order to create a similar probability output for regression, following Lin and Weng (2004), we suppose that the SVM is trained on data from the model where f (x i ) is the underlying function and δ i is independent and identical distributed random noise.Given a test data x the distribution of y given x and allows one to draw probabilistic inferences about y e.g. one can construct a predictive interval Φ = Φ(x) such that y ∈ Φ with a certain probability.If f is the estimated (predicted) function of the SVM on new data then η = η(x) = y − f (x) is the prediction error and y ∈ Φ is equivalent to η ∈ Φ. Empirical observation shows that the distribution of the residuals η can be modeled both by a Gaussian and a Laplacian distribution with zero mean.In this implementation the Laplacian with zero mean is used : Assuming that η are independent the scale parameter σ is estimated by maximizing the likelihood.The data for the estimation is produced by a three-fold cross-validation.For the Laplace distribution the maximum likelihood estimate is : i-class classification supported by the ksvm function is the one proposed in Crammer and Singer (2000).This algorithm works by solving a single optimization problem including the data from all classes: where where the decision function is argmax m=1,...,k ⟨x i , w n ⟩ This optimization problem is solved by a decomposition method proposed in Hsu and Lin (2002b) where optimal working sets are found (that is, sets of α i values which have a high probability of being non-zero).The QP sub-problems are then solved by a modified version of the TRON13 (Lin and More 1999) optimization software.
One-class classification or novelty detection (Schölkopf, Platt, Shawe-Taylor, Smola, and Williamson 1999;Tax and Duin 1999), where essentially an SVM detects outliers in a data set, is another algorithm supported by ksvm.SVM novelty detection works by creating a spherical decision boundary around a set of data points by a set of support vectors describing the spheres boundary.The ν parameter is used to control the volume of the sphere and consequently the number of outliers found.Again, the value of ν represents the fraction of outliers found.Furthermore, ϵ-SVM (Vapnik 1995) and ν-SVM (Schölkopf, Smola, Williamson, and Bartlett 2000) regression are also available.The problem of model selection is partially addressed by an empirical observation for the popular Gaussian RBF kernel (Caputo, Sim, Furesjo, and Smola 2002), where the optimal values of the hyper-parameter of sigma are shown to lie in between the 0.1 and 0.9 quantile of the ∥x − x ′ ∥ statistics.

Relevance vector machine
The relevance vector machine (Tipping 2001) is a probabilistic sparse kernel model identical in > set.seed(123)> x <-rbind(matrix(rnorm(120),,2),matrix(rnorm(120,mean=3),,2)) > y <-matrix(c(rep(1,60),rep(-1,60))) > svp <-ksvm(x,y,type="C-svc") > plot(svp,data=x) where α n are the model "weights" and K(• , • ) is a kernel function.It adopts a Bayesian approach to learning, by introducing a prior over the weights α governed by a set of hyper-parameters β, one associated with each weight, whose most probable values are iteratively estimated for the data.Sparsity is achieved because in practice the posterior distribution in many of the weights is sharply peaked around zero.Furthermore, unlike the SVM classifier, the non-zero weights in the RVM are not associated with examples close to the decision boundary, but rather appear to represent "prototypical" examples.These examples are termed relevance vectors.
kernlab currently has an implementation of the RVM based on a type II maximum likelihood method which can be used for regression.The functions returns an S4 object containing the model parameters along with indexes for the relevance vectors and the kernel function and hyperparameters used.

Gaussian processes
Gaussian processes (Williams and Rasmussen 1995) are based on the "prior" assumption that adjacent observations should convey information about each other.In particular, it is assumed that the observed variables are normal, and that the coupling between them takes place by means of the covariance matrix of a normal distribution.Using the kernel matrix as the covariance matrix is a convenient way of extending Bayesian modeling of linear estimators to nonlinear situations.Furthermore it represents the counterpart of the "kernel trick" in methods minimizing the regularized risk.
For regression estimation we assume that rather than observing t(x i ) we observe y i = t(x i ) + ξ i where ξ i is assumed to be independent Gaussian distributed noise with zero mean.The posterior distribution is given by and maximizing ln p(α | y) for α to obtain the maximum a posteriori approximation yields Knowing α allows for prediction of y at a new location x through y = K(x, x i )α.In similar fashion Gaussian processes can be used for classification.
gausspr is the function in kernlab implementing Gaussian processes for classification and regression.

Ranking
The success of Google has vividly demonstrated the value of a good ranking algorithm in real world problems.kernlab includes a ranking algorithm based on work published in (Zhou, Weston, Gretton, Bousquet, and Schölkopf 2003).This algorithm exploits the geometric structure of the data in contrast to the more naive approach which uses the Euclidean distances or inner products of the data.Since real world data are usually highly structured, this algorithm should perform better than a simpler approach based on a Euclidean distance measure.First, a weighted network is defined on the data and an authoritative score is assigned to every point.The query points act as source nodes that continually pump their scores to the remaining points via the weighted network, and the remaining points further spread the score to their neighbors.The spreading process is repeated until convergence and the points are ranked according to the scores they received.Suppose we are given a set of data points X = x 1 , . . ., x s , x s+1 , . . ., x m in R n where the first s points are the query points and the rest are the points to be ranked.The algorithm works by connecting the two nearest points iteratively until a connected graph G = (X, E) is obtained where E is the set of edges.The affinity matrix K defined e.g. by K ij = exp(−σ∥x i − x j ∥ 2 ) if there is an edge e(i, j) ∈ E and 0 for the rest and diagonal elements.The matrix is normalized as > data(spirals) > ran <-spirals[rowSums(abs(spirals) < 0.55) == 2,] > ranked <-ranking(ran, 54, kernel = "rbfdot", kpar = list(sigma = 100), edgegraph = TRUE) > ranked [54, 2] <-max(ranked[-54, 2]) > c<-1:86 > op <-par(mfrow = c(1, 2),pty="s") > plot(ran) > plot (ran, cex=c[ranked[,3]]/40) Figure 3: The points on the left are ranked according to their similarity to the upper most left point.Points with a higher rank appear bigger.Instead of ranking the points on simple Euclidean distance the structure of the data is recognized and all points on the upper structure are given a higher rank although further away in distance than points in the lower structure.
is iterated until convergence, where α is a parameter in [0, 1).The points are then ranked according to their final scores f i (t f ).
kernlab includes an S4 method implementing the ranking algorithm.The algorithm can be used both with an edge-graph where the structure of the data is taken into account, and without which is equivalent to ranking the data by their distance in the projected space.

Online learning with kernels
The onlearn function in kernlab implements the online kernel algorithms for classification, novelty detection and regression described in (Kivinen, Smola, and Williamson 2004).In batch learning, it is typically assumed that all the examples are immediately available and are drawn independently from some distribution P .One natural measure of quality for some f in that case is the expected Since usually P is unknown a standard approach is to instead minimize the empirical risk Minimizing R emp [f ] may lead to overfitting (complex functions that fit well on the training data but do not generalize to unseen data).One way to avoid this is to penalize complex functions by kernlab -An S4 Package for Kernel Methods in R instead minimizing the regularized risk.
where λ > 0 and ∥f ∥ H = ⟨f, f ⟩ 1 2 H does indeed measure the complexity of f in a sensible way.The constant λ needs to be chosen appropriately for each problem.Since in online learning one is interested in dealing with one example at the time the definition of an instantaneous regularized risk on a single example is needed The implemented algorithms are classical stochastic gradient descent algorithms performing gradient descent on the instantaneous risk.The general form of the update rule is : where f i ∈ H and ∂ f < is short hand for ∂ ∂f (the gradient with respect to f ) and η t > 0 is the learning rate.Due to the learning taking place in a reproducing kernel Hilbert space H the kernel k used has the property ⟨f, k(x, • )⟩ H = f (x) and therefore where l ′ (z, y) The onlearn function implements the online learning algorithm for regression, classification and novelty detection.The online nature of the algorithm requires a different approach to the use of the function.An object is used to store the state of the algorithm at each iteration t this object is passed to the function as an argument and is returned at each iteration t + 1 containing the model parameter state at this step.An empty object of class onlearn is initialized using the inlearn function.

Spectral clustering
and normalized as Then the top k eigenvectors (where k is the number of clusters to be found) of the affinity matrix are used to form an n × k matrix Y where each column is normalized again to unit length.Treating each row of this matrix as a data point, kmeans is finally used to cluster the points.
kernlab includes an S4 method called specc implementing this algorithm which can be used through an formula interface or a matrix interface.The S4 object returned by the method extends the class "vector" and contains the assigned cluster for each point along with information on the centers size and within-cluster sum of squares for each cluster.In case a Gaussian RBF kernel is being used a model selection process can be used to determine the optimal value of the σ hyperparameter.For a good value of σ the values of Y tend to cluster tightly and it turns out that the within cluster sum of squares is a good indicator for the "quality" of the sigma parameter found.We then iterate through the sigma values to find an optimal value for σ.

Kernel principal components analysis
Principal component analysis (PCA) is a powerful technique for extracting structure from possibly high-dimensional datasets.PCA is an orthogonal transformation of the coordinate system in which we describe the data.The new coordinates by which we represent the data are called principal components.Kernel PCA (Schölkopf, Smola, and Müller 1998) performs a nonlinear transformation of the coordinate system by finding principal components which are nonlinearly related to the input variables.Given a set of centered observations x k , k = 1, . . ., M , x k ∈ R N , PCA diagonalizes the covariance matrix C = 1 M M j=1 x j x T j by solving the eigenvalue problem λv = Cv.The same computation can be done in a dot product space F which is related to the input space by a possibly nonlinear map Φ : R N → F , x → X.Assuming that we deal with centered data and use the covariance matrix in F , the kernel principal components are then computed by taking the eigenvectors of the centered kernel matrix K ij = ⟨Φ(x j ), Φ(x j )⟩.
kpca, the the function implementing KPCA in kernlab, can be used both with a formula and a matrix interface, and returns an S4 object of class kpca containing the principal components the corresponding eigenvalues along with the projection of the training data on the new coordinate system.Furthermore, the predict function can be used to embed new data points into the new coordinate system.

Kernel feature analysis
Whilst KPCA leads to very good results there are nevertheless some issues to be addressed.First the computational complexity of the standard version of KPCA, the algorithm scales O(m 3 ) and secondly the resulting feature extractors are given as a dense expansion in terms of the of the training patterns.Sparse solutions are often achieved in supervised learning settings by using an l 1 penalty on the expansion coefficients.An algorithm can be derived using the same approach in feature extraction requiring only n basis functions to compute the first n feature.Kernel feature analysis (Smola, Mangasarian, and Schölkopf 2000) is computationally simple and scales approximately one order of magnitude better on large data sets than standard KPCA.Choosing Ω This setting leads to the first "principal vector" in the l 1 context Subsequent "principal vectors" can be defined by enforcing optimality with respect to the remaining orthogonal subspaces.Due to the l 1 constrain the solution has the favorable property of being sparse in terms of the coefficients α i .The function kfa in kernlab implements Kernel Feature Analysis by using a projection pursuit technique on a sample of the data.Results are then returned in an S4 object.

Kernel canonical correlation analysis
one wishes to find those values of w 1 and w 2 which maximize the correlation between y 1 and y 2 .Similar to the KPCA algorithm, CCA can be extended and used in a dot product space F which is related to the input space by a possibly nonlinear map Φ : R N → F , x → X as Following (Kuss and Graepel 2003), the kernlab implementation of a KCCA projects the data vectors on a new coordinate system using KPCA and uses linear CCA to retrieve the correlation coefficients.The kcca method in kernlab returns an S4 object containing the correlation coefficients for each data set and the corresponding correlation along with the kernel used.

Interior point code quadratic optimizer
In many kernel based algorithms, learning implies the minimization of some risk function.Typically we have to deal with quadratic or general convex problems for support vector machines of the type minimize f (x) subject to c i (x) ≤ 0 for all i ∈ [n]. (37) f and c i are convex functions and n ∈ N. kernlab provides the S4 method ipop implementing an optimizer of the interior point family (Vanderbei 1999) which solves the quadratic programming problem minimize This optimizer can be used in regression, classification, and novelty detection in SVMs.

Incomplete cholesky decomposition
When dealing with kernel based algorithms, calculating a full kernel matrix should be avoided since it is already a O(N 2 ) operation.Fortunately, the fact that kernel matrices are positive semidefinite is a strong constraint and good approximations can be found with small computational cost.The Cholesky decomposition factorizes a positive semidefinite N × N matrix K as K = ZZ T , where Z is an upper triangular N × N matrix.Exploiting the fact that kernel matrices are usually of low rank, an incomplete Cholesky decomposition (Wright 1999) finds a matrix Z of size N × M where M ≪ N such that the norm of K − Z ZT is smaller than a given tolerance θ.The main difference of incomplete Cholesky decomposition to the standard Cholesky decomposition is that pivots which are below a certain threshold are simply skipped.If L is the number of skipped pivots, we obtain a Z with only M = N − L columns.The algorithm works by picking a column from K to be added by maximizing a lower bound on the reduction of the error of the approximation.kernlab has an implementation of an incomplete Cholesky factorization called inc.chol which computes the decomposed matrix Z from the original data for any given kernel without the need to compute a full kernel matrix beforehand.This has the advantage that no full kernel matrix has to be stored in memory.

Conclusions
In this paper we described kernlab, a flexible and extensible kernel methods package for R with existing modern kernel algorithms along with tools for constructing new kernel based algorithms.It provides a unified framework for using and creating kernel-based algorithms in R while using all of R's modern facilities, like S4 classes and namespaces.Our aim for the future is to extend the package and add more kernel-based methods as well as kernel relevant tools.Sources and binaries for the latest version of kernlab are available at CRAN 14 under the GNU Public License.
A shorter version of this introduction to the R package kernlab is published as Karatzoglou, Smola, Hornik, and Zeileis (2004) in the Journal of Statistical Software.

Figure 1 :
Figure 1: A contour plot of the SVM decision values for a toy binary classification problem using the plot function

Figure 2 :
Figure 2: Relevance vector regression on data points created by the sinc(x) function, relevance vectors are shown circled.

Figure 4 :
Figure 4: Clustering the two spirals data set with specc

Figure 6 :
Figure 6: Projection of the spam data on two features using an RBF kernel The notion of a kernel is obviously central in the making of any kernel-based algorithm and consequently also in any software package containing kernel-based methods.Kernels in kernlab are S4 objects of class kernel extending the function class with one additional slot containing a list with the kernel hyper-parameters.Package kernlab includes 7 different kernel classes which all contain the class kernel and are used to implement the existing kernels.These classes are used in the function dispatch mechanism of the kernel utility functions described below.Existing kernel functions are initialized by "creator" functions.All kernel functions take two feature vectors as parameters and return the scalar dot product of the vectors.An example of the functionality of a kernel in kernlab: The sigest function uses a sample of the training set to estimate the quantiles and returns a vector containing the values of the quantiles.Pretty much any value within this interval leads to good performance.An example for the ksvm function is shown below.