Abstract
Background
Protein fold recognition is a key step in protein threedimensional (3D) structure discovery. There are multiple fold discriminatory data sources which use physicochemical and structural properties as well as further data sources derived from local sequence alignments. This raises the issue of finding the most efficient method for combining these different informative data sources and exploring their relative significance for protein fold classification. Kernel methods have been extensively used for biological data analysis. They can incorporate separate fold discriminatory features into kernel matrices which encode the similarity between samples in their respective data sources.
Results
In this paper we consider the problem of integrating multiple data sources using a kernelbased approach. We propose a novel informationtheoretic approach based on a KullbackLeibler (KL) divergence between the output kernel matrix and the input kernel matrix so as to integrate heterogeneous data sources. One of the most appealing properties of this approach is that it can easily cope with multiclass classification and multitask learning by an appropriate choice of the output kernel matrix. Based on the position of the output and input kernel matrices in the KLdivergence objective, there are two formulations which we respectively refer to as MKLdivdc and MKLdivconv. We propose to efficiently solve MKLdivdc by a difference of convex (DC) programming method and MKLdivconv by a projected gradient descent algorithm. The effectiveness of the proposed approaches is evaluated on a benchmark dataset for protein fold recognition and a yeast protein function prediction problem.
Conclusion
Our proposed methods MKLdivdc and MKLdivconv are able to achieve stateoftheart performance on the SCOP PDB40D benchmark dataset for protein fold prediction and provide useful insights into the relative significance of informative data sources. In particular, MKLdivdc further improves the fold discrimination accuracy to 75.19% which is a more than 5% improvement over competitive Bayesian probabilistic and SVM marginbased kernel learning methods. Furthermore, we report a competitive performance on the yeast protein function prediction problem.
Background
A huge number of protein coding sequences have been generated by genome sequencing projects. In contrast, there is a much slower increase in the number of known threedimensional (3D) protein structures. Determination of a protein's 3D structure is a formidable challenge if there is no sequence similarity to proteins of known structure and thus protein structure prediction remains a core problem within computational biology.
Computational prediction of protein structure has achieved significant successes [1,2]. Focusing on the fold prediction problem of immediate interest to this paper, one computational method known as the taxonomic approach [3,4], presumes the number of folds is restricted and focuses on structural predictions in the context of a particular classification of 3D folds. Proteins are in a common fold if they share the same major secondary structures in the same arrangement and the same topological connections [5,6]. In the taxonomic method for protein fold classification, there are several fold discriminatory data sources or groups of attributes available such as amino acid composition, predicted secondary structure, and selected structural and physicochemical properties of the constituent amino acids. Previous methods for integrating these heterogeneous data sources include simply merging them together or combining trained classifiers over individual data sources [3,4,7,8]. However, how to integrate fold discriminatory data sources systematically and efficiently, without resorting to ad hoc ensemble learning, still remains a challenging problem.
Kernel methods [9,10] have been successfully used for data fusion in biological applications. Kernel matrices encode the similarity between data objects within a given input space and these data objects can include graphs and sequence strings in addition to realvalued or integer data. Thus the problem of data integration is transformed into the problem of learning the most appropriate combination of candidate kernel matrices, representing these heterogeneous data sources. The typical framework is to learn a linear combination of candidate kernels. This is often termed multiple kernel learning (MKL) in Machine Learning, and nonparametric group lasso in Statistics. Recent trends in kernel learning are usually based on the margin maximization criterion used by Support Vector Machines (SVMs) or variants [11]. The popularity of SVM marginbased kernel learning stems from its efficient optimization formulations [1114] and sound theoretical foundation [11,15,16]. Other data integration methods include the COSSO estimate for additive models [17], kernel discriminant analysis [18], multilabel multiple kernel learning [19,20] and Bayesian probabilistic models [21,22]. These methods, in general, can combine multiple data sources to enhance biological inference [21,23] and provide insights into the significance of the different data sources used.
Following a different approach, in this paper we propose an alternative criterion for kernel matrix learning and data integration, which we will call MKLdiv. Specifically, we propose an informationtheoretic approach to learn a linear combination of kernel matrices, encoding information from different data sources, through the use of a KullbackLeibler divergence [2428] between two zeromean Gaussian distributions defined by the input matrix and output matrix. The potential advantage of this approach is that, by choosing different output matrices, the method can be easily extended to different learning tasks, such as multiclass classification and multitask learning. These are common tasks in biological data analysis.
To illustrate the method, we will focus on learning a linear combination of candidate kernel matrices (heterogeneous data sources) using the KLdivergence criterion with a main application to the protein fold prediction problem. There are two different formulations based on the relative position of the input kernel matrix and the output kernel matrix in the KLdivergence objective. For the first formulation, although this approach involves a matrix determinant term which is not convex in general, we elegantly reformulate the learning task as a difference of convex problem, which can be efficiently solved by a sequence of convex optimizations. Hence we refer to it as MKLdivdc. The second KLdivergence formulation for kernel integration, called MKLdivconv, is convex and can be solved by a projected gradient descent algorithm. Experimental results show that these formulations lead to stateoftheart prediction performance. In particular, MKLdivdc outperforms the best reported performance on the important task of protein fold recognition, for the benchmark dataset used.
Methods
In the following we first revisit kernel learning approaches based on SVMs [11] and kernel discriminant analysis [18]. Then, we introduce our novel informationtheoretic approach for data integration based on a KLdivergence criterion. Finally we discuss how to solve the optimization task efficiently. For brevity, we use the conventional notation ℕ_{n }= {1, 2, ..., n} for any n ∈ ℕ.
Background and Related Work
Kernel methods are extensively used for biological data analysis. A symmetric function K : X × X → ℝ is called a kernel function if it is positive semidefinite, by which we mean that, for any n ∈ ℕ and {x_{i }∈ X: i ∈ ℕ_{n}}, the Gram matrix is positive semidefinite. According to [29], its corresponding reproducing kernel Hilbert space (RKHS), usually denoted by ℋ_{K}, can be defined to be the completion of the linear span of the set of functions {K_{x}(·) := K(x, ·): x ∈ X} with inner product satisfying, for any x ∈ X and g ∈ ℋ_{K}, the reproducing property ⟨K_{x}, g⟩_{K }= g(x). By Mercer's theorem, there exists a high dimensional (possible infinite) Hilbert feature space ℱ with inner product ⟨·, ·⟩_{ℱ }and a feature map ϕ: X → ℱ such that K(x, t) = ⟨ϕ (x), ϕ (t)⟩_{ℱ}, ∀ x, t ∈ X. Intuitively, the kernel function K implicitly maps the data space X into a high dimensional space ℱ, see [9,10] for more details.
Within the context of protein fold recognition, we have m different fold discriminatory data sources where samples across each data source can be represented by for ℓ ∈ ℕ_{m }and the outputs are denoted by y = {y_{i }: i ∈ ℕ_{n}}. For kernel methods, for any ℓ ∈ ℕ_{m}, each ℓth data source can be encoded into a candidate kernel matrix denoted by . Depending on the different data sources used, the candidate kernel function K_{ℓ }should be specified a priori by the practitioner. The composite kernel matrix is given by where {λ_{ℓ}: ℓ ∈ ℕ_{m}} are realvalued kernel weights and typically they are restricted to be nonnegative. In this context, the problem of data integration is consequently reduced to the problem of learning a convex combination of candidate kernel matrices: more precisely learning the kernel weights λ. Different optimization criteria over the candidate kernels arise from the particular kernel learning algorithm used. Cristianini et al. [30] proposed a kernel learning approach which uses the cosine of the angle between the two bidimensional vectors K_{λ }and K_{y }representing the Gram matrices. This is achieved by maximizing the kernel alignment:
The above kernel learning formulation can be solved by a semidefinite programming (SDP) approach (see Section 4.7 of [11]). However, an SDP formulation is computationally intensive.
Another widely used criterion for kernel learning is based on the margin concept in SVMs and variants. Denoting the simplex set as Δ = {λ = (λ_{1}, λ_{2}, ..., λ_{m}): }, Lanckriet et al [11] proposed the following formulation for kernel learning:
where 1_{n }is a column vector of all ones, C is a tradeoff parameter, and t = (t_{1}, t_{2}, ..., t_{n}) denotes the binary outputs with t_{i }∈ {1, 1} being the class label for ith instance. This task was reformulated as a quadratically constrained quadratic programming (QCQP) problem and later improved by Sonnenburg et al. [14] who reformulated it as a semiinfinite linear programming (SILP) task. Moreover, it was pointed out in [12,13,17,31] that this is equivalent to the following sparse L^{1}regularization formulation:
The L^{1}regularization term encourages the sparsity [32] of RKHSnorm terms, and thus indicates the relative importance of data sources. It was shown in [13] that the standard L^{2}regularization is equivalent to the use of uniformly weighted kernel weights λ, i.e. for any ℓ ∈ ℕ_{m}. Recently, Ye et al. [18] proposed an appealing kernel learning approach based on regularized kernel discriminant analysis. This can similarly be shown to be equivalent to a sparse L^{1}regularization formulation with a least square loss, see Appendix 1 for more details.
Informationtheoretic Data Integration
In this paper we adopt a novel informationtheoretic approach to learn the kernel combinatorial weights. The main idea is to quantify the similarity between K_{λ }and K_{y }through a KullbackLeibler (KL) divergence or relative entropy term [2428]. This approach is based on noting that these kernel matrices encode the similarity of data objects within their respective input and label data spaces. Furthermore, there is a simple bijection between the set of distance measures in these data spaces and the set of zeromean multivariate Gaussian distributions [25]. Using this bijection, the difference between two distance measures, parameterized by K_{λ }and K_{y}, can be quantified by the relative entropy or KullbackLeibler (KL) divergence between the corresponding multivariate Gaussians. Matching kernel matrices K_{λ }and K_{y }can therefore be realized by minimizing a KL divergence between these distributions and we will exploit this approach below in the context of multiple kernel learning.
Kernel matrices are generally positive semidefinite and thus can be regarded as the covariance matrices of Gaussian distributions. As described in [24], the KullbackLeibler (KL) divergence (relative entropy) between a Gaussian distribution (0, K_{y}) with the output covariance matrix K_{y }and a Gaussian distribution (0, K_{x}) with the input kernel covariance matrix K_{x }is:
where, for any square matrix B, the notation Tr(B) denotes its trace. The a priori choice of the output matrix K_{y }will be discussed later. Though KL ((0, K_{y})(0, K_{x})) is nonconvex w.r.t. K_{x}, it has a unique minimum at K_{x }= K_{y }if K_{y }is positive definite, suggesting that minimizing the above KLdivergence encourages K_{x }to approach K_{y}. If the input kernel matrix K_{x }is represented by a linear combination of m candidate kernel matrices, i.e. , the above KLdivergence based kernel learning is reduced to the following formulation:
where I_{n }denotes the n × n identity matrix and σ > 0 is a supplemented small parameter to avoid the singularity of K_{λ}.
Since the KLdivergence is not symmetric with respect to K_{y }and K_{λ}, another alternative approach to matching kernel matrices is given by
where parameter σ > 0 is added to avoid the singularity of K_{y}. If there is no positive semidefiniteness restriction over K_{ℓ}, this formulation is a wellknown convex maximumdeterminant problem [33] which is a more general formulation than semidefinite programming (SDP), its implementation is computationally intensive, and thus cannot be extended to largescale problems according to [33]. However, formulation (5) has a special structure here: λ_{ℓ }is nonnegative and all candidate kernel matrices are positive semidefinite. Hence, we can solve this problem by a simple projected gradient descent method, see below for more details.
The KLdivergence criterion for kernel integration was also successfully used in [27,28] which formulated the problem of supervised network inference as a kernel matrix completion problem. In terms of information geometry [34], formulation (4) corresponds to finding the mprojection of K_{y }over an eflat submanifold. The convex problem (5) can be regarded as finding the eprojection of K_{y }over a mflat submanifold. In [26], formulation (4) was developed for learning an optimal linear combination of diffusion kernels for biological networks. A gradientbased method was employed in [26] to learn a proper linear combination of diffusion kernels. This optimization method largely relies on the special property of all candidate diffusion kernel matrices enjoying the same eigenvectors and the gradientbased learning method could be a problem if we deal with general kernel matrices. In the next section, we propose to solve the general kernel learning formulation (4) using a difference of convex optimization method.
The formulation (4) also has a close relation with Gaussian Process regression [35]. A Gaussian process f can be fully specified by giving the covariance matrix for any finite set of zeromean random variables f = {f(x_{i}): i ∈ ℕ_{m}}. The relation between the inputs x = {x_{i }: i ∈ ℕ_{n}} and outputs y = {y_{i }: i ∈ ℕ_{m}} is realized by the latent variable f as follows:
where I_{n }denotes the n × n identity matrix and the latent random variable f = (f (x_{1}, ..., f (x_{n}))) is distributed as a Gaussian process prior. The Gaussian process prior can be fully specified by a kernel K with a random covariance matrix associated with random samples x = {x_{i}: i ∈ ℕ_{n}}. Specifically, it can be written as fx ~(f 0, K_{λ}). We assume a uniform distribution over λ, i.e. a Dirichlet prior distribution with α_{0 }= 1. If we let K_{y }= yy^{⊤ }in the objective function of formulation (4), then one can easily check that, up to a constant term, the objective function in formulation (4) is the negative of the log likelihood of Gaussian process regression, and maximizing the log likelihood is equivalent to the minimization problem (4).
Optimization Formulation
We now turn our attention to optimization approaches for the KLdivergence based kernel learning formulations (4) and (5). In particular, we show that formulation (5) can be approached by a projected gradient descent method and (4) can be solved by a difference of convex algorithm (DCA) [36] which, for linear constraint conditions, reduces to the special case of a concave convex procedure (CCCP) [37]. To this end, let
and
Theorem 1 Let functions g and f be defined by (6) and (7). Then, both f and g are convex with respect to λ ∈ Δ. Moreover, problem (5) is convex and problem (4) is a difference of convex problem, i.e.
Proof It suffices to prove the convexity of f and g. To this end, from [38] we observe that functions – log C and Tr(K_{y}C^{1}) are convex with respect to positive semidefinite matrices C. Hence, f and g are convex with respect to λ ∈ Δ. This completes the proof of the theorem.
For simplicity we refer to the KLdivergence kernel learning formulation (4) as MKLdivdc since it is a difference of convex problem and refer to formulation (5) as MKLdivconv since it is a convex problem.
Projected Gradient Descent Method for MKLdivconv
We propose a projected gradient descent (PGD) method to solve problem (5). The idea of this method is to alternately implement a gradient descent and then a projection to the feasible domain, see e.g. [39]. Recall the derivative of the log determinant,(see e.g. the matrix cookbook http://matrixcookbook.com/ webcite
With a little abuse of notation, we also denote by
the objective function of problem (5). Consequently, its gradient is given by
Then, at iteration step t the gradient descent step is realized by
where η > 0 is a prescribed step size. The projection of β to the feasible domain Δ can be written as the following quadratic programming problem
The theoretical convergence rate of the projected gradient descent method is generally of complexity where t is the iteration number and L is the Lipschitz constant of the gradient function defined by (10), see e.g. [39]. Here, the Lipschitz constant L is bounded by the largest eigenvalue of the Hessian of the objective function defined, for any i, j ∈ ℕ_{m}, by
Since ℒ is convex, the Hessian ℋ(ℒ) is positive semidefinite and thus
where ·_{Fro }denotes the Frobenious norm of a matrix. Hence, the projected gradient descent algorithm could take longer time to become convergent if the value of σ is very small.
Difference of Convex Algorithm for MKLdivdc
By Theorem 1, problem (4) is a difference of convex problem. We propose to solve this problem by a concave convex procedure (CCCP) [36,37]. This procedure iteratively solves the following convex problem:
where, for any j ∈ ℕ_{m}, the derivative of the log determinant is given by equation (9). Before we continue the main discussion, let us first note an interesting property of CCCP. By the definition of λ^{(t+1)}, we know that
Since g is convex, we have that
Consequently,
which means that the objective value ℒ(λ^{(t)}) monotonically decreases with each iteration. Consequently, we can use the relative change of the objective function as a stopping criterion. Local convergence of the DCA algorithm is proven in [36] (Lemma 3.6, Theorem 3.7). Tao and An [36] state that the DCA often converges to the global solution. Overall, the DC programming approach to MKLdivdc can be summarized as follows.
• Given a stopping threshold ε
• Initialize λ^{(1)}, e.g. for any ℓ ∈ ℕ_{m}
• Given the solution λ^{(t) }at step t, for step t + 1, first compute ▽g(λ^{(t)}) by equation (9). Then, compute solution λ^{(t+1) }of convex subproblem (13).
• Stop until the relative change where ε is a stopping threshold
SILP Formulation for the Convex Subproblem (13)
We now turn to the solution of the convex subproblem (13). To see this, first decompose the output matrix K_{y }into the form K_{y }= AA^{⊤}, e.g. by eigendecomposition. Here, A is an n × r matrix with r = rank(A) which always exists since K_{y }is positive semidefinite. Hence, by introducing an auxiliary matrix α ∈ ℝ^{n × r}, we observe, for any positive definite matrix C, that
Applying the above equality with , up to a constant, equation (13) is equivalent to the augmented problem:
Equivalently, by the minmax theorem (see e.g. [38])
To solve the subproblem (15), we can formulate it as a quadratically constrained quadratic programming (QCQP) problem as in [11]. Here we formulate the problem in (15) as a semiinfinite linear programming (SILP) problem [14,40] since SILP usually has better scalability compared to QCQP. To see this, let , and . Then, letting , we can rewrite (15) as a SILP problem:
In (16), there are an infinite number of constraints (indexed by a), indicative of a semiinfinite linear programming (SILP) problem. The SILP task can be solved by an iterative column generation algorithm (or exchange method) which is guaranteed to converge to a global optimum. A brief description of the column generation method is illustrated in Appedix 2.
Alternatively we could apply the projected gradient descent (PGD) method in the above subsection directly to the convex subproblem (13). However, the gradient function of its objective function involves the matrix . In analogy to the argument of inequality (12), the Lipschitz constant of the gradient of the objective function in (13) is very large when the value of σ is very small, and thus the projected gradient descent algorithm could take longer to become convergent. Hence, this could make the overall DC programming unacceptable slow. In contrast, in the SILP formulation (16) we introduce the auxiliary variables α to avoid the matrix . In addition, the gradient descent algorithm generally needs to determine the step size η according to the value of σ, see also discussion in the experimental section.
Prior Choice of the Output Kernel Matrix
The choice of the output kernel matrix K_{y }will depend on the problem considered. We first consider a multiclass classification for the specific task of protein fold recognition. In this case, we preprocess the output labels using a oneagainstall strategy. In particular, for a Cclass classification we recast the outputs y = {y_{i }: i ∈ ℕ_{n}} as (y_{i1}, ..., y_{iC}) such that y_{ip }= 1 if x_{i }is in class p and otherwise 1. Hence the outputs are represented by an n × C indicator matrix Y = (y_{ip})_{i, p }whose pth column vector is denoted by y_{p}. Then, taking K_{y }= YY^{⊤}, formulation (4) can be extended to the joint optimization problem
and formulation (5) can be written as
For the protein fold recognition and yeast protein function prediction projects discussed below, we choose K_{y }= YY^{⊤ }as stated.
In general, though, K_{y }might encode a known structural relationship between labels. For example, in supervised gene or protein network inference (see e.g. [41,42]) the output information corresponds to an adjacency (square) matrix A where A_{ij }= 1 means there is an interaction between gene or protein pair (e_{i}, e_{j}) of an organism, otherwise A_{ij }= 0. In this case, the output kernel matrix K_{y }can potentially be chosen as the graph Laplacian defined as L = diag(A1)  A, where 1 is the vector of all ones. It can also be formulated as a diffusion kernel [43] defined by , where hyperparameter β > 0. Other potential choices of K_{y }can be found in [19,20] for multilabeled datasets.
Results and Discussion
We mainly evaluate MKLdiv methods (MKLdivdc and MKLdivconv) on protein fold recognition, and then consider an extension to the problem of yeast protein function prediction. In these tasks we first compute the kernel weights by MKLdiv and then feed these into a oneagainstall multiclass SVM to make predictions. The tradeoff parameter in the multiclass SVM is adjusted by 3fold cross validation over the training dataset. For all experiments with MKLdivdc, we choose σ = 10^{5 }and for MKLdivconv, we tune σ = {10^{5}, ..., 10^{1}} using cross validation. In both methods, we use a stopping criterion of ε = 10^{5 }and initialize the kernel weight λ by setting for any ℓ ∈ ℕ_{m }where m is the number of candidate kernel matrices.
Synthetic Data
We first validated the proposed MKLdiv algorithms on a simple threeclass dataset illustrated in subfigure (a) of Figure 1. As in [11], we use a Gaussian kernel with unit variance, a polynomial kernel of order two and a linear kernel. In this case we demonstrate the effect of our approaches on combining kernel matrices derived from a single data source. Subfigures (e) and (f) of Figure 1 illustrate the kernel weights learned by MKLdivdc and MKLdivconv. In particular, MKLdivdc successfully picked up the Gaussian kernel as the most dominant kernel, which is more reasonable than MKLdivconv. Subfigures (b) and (c) of Figure 1 show the relative change of objective function values versus iteration, i.e. (ℒ(λ^{(t1)})  ℒ(λ^{(t)}))/ℒ(λ^{(t)}), of MKLdivdc and MKLdivconv. We can see that the DC algorithm of MKLdivdc converges quickly to a local minimum while the projected gradient descent algorithm converges a little slower to a global minimum. However, MKLdivdc needs more time per iteration than MKLdivconv since MKLdivdc needs to solve the subproblem (13) at each iteration. As mentioned before, the subproblem (13) can be solved by either semiinfinite linear programming (SILP) or a projected gradient descent (PGD) method. To see their convergence, in subfigure (d) of Figure 1 we plot the relative changes of the objective function in subproblem (13) when for ℓ ∈ ℕ_{m}. We can see from subfigure (d) that the PGD approach converges faster in the beginning but stalls at a higher precision and the SILP method converges faster at higher precision.
Figure 1. Synthetic data verification. Synthetic data verification of MKLdiv: (a) depiction of the threecircle dataset; (b) relative change of objective values of MKLdivdc versus iteration number of CCCP; (c) relative change of objective values of MKLdivconv versus iteration number of projected gradient descent (PGD) method; (d) relative change of objective values of subproblem (13) by SILP (dishline) and PGD methods; (e) kernel weights learned by MKLdivdc; (f) kernel weights learned by MKLdivconv.
Protein Fold Recognition
Next we evaluated MKLdiv on a wellknown protein fold prediction dataset [3]. This benchmark dataset (based on SCOP PDB40D) has 27 SCOP fold classes with 311 proteins for training and 383 for testing. This dataset was originally proposed by Ding and Dubchak [3] and it has 313 samples for training and 385 for testing. There is less than 35% sequence identity between any two proteins in the training and test set. We follow Shen and Chou [4] who proposed to exclude two proteins from the training and test datasets due to a lack of sequence information. We compare our MKLdiv methods with kernel learning based on oneagainstall multiclass SVM using the SimpleMKL software package [44], kernel learning for regularized discriminant analysis (MKLRKDA) [18]http://www.public.asu.edu/~jye02/Software/DKL/ webcite and a probabilistic Bayesian model for kernel learning (VBKC) [21]. The tradeoff parameters in SimpleMKL and MKLRKDA were also adjusted by 3fold cross validation on the training set.
Description of the Fold Discriminatory Data Sources
As listed in Table 1, there are a total of 12 different fold discriminatory data sources available: Amino Acid Composition (C), Predicted Secondary Structure (S), Hydrophobicity (H), Polarity (P), van der Waals volume (V), Polarizability (Z), PseAA λ = 1 (L1), PseAA λ = 4 (L4), PseAA λ = 14 (L14), PseAA λ = 30 (L30), SW with BLOSUM62 (SW1) and SW with PAM50 (SW2). The first six data sources were originally from [3]. Four data sources using different dimensions of pseudoamino acid composition (PseAA) were introduced in [4] to replace the aminoacid composition. The last two data sources used in [21] are derived from a pairwise kernel [45] for local sequence alignment based on SmithWaterman scores.
Table 1. Performance with individual and all data sources
As in [21], we employ linear kernels (SmithWaterman scores) for SW1 and SW2 and second order polynomial kernels for the other data sources. Ding and Dubchak [3] conducted an extensive study on the use of various multiclass variants of standard SVMs and neural network classifiers. For these authors the best test set accuracy (TSA) was 56%, and the most informative among their six data sources (CSHPVZ) were aminoacid composition (C), the predicted secondary structure (S) and hydrophobicity (H). Shen and Chou [4] introduced four additional PSeAA data sources to replace the amino acid composition (C) and raised test performance to 62.1%. The latter authors used an ad hoc ensemble learning approach involving a combination of multiclass k nearest neighbor classifiers individually trained on each data source. Recently, test performance was greatly improved by Damoulas and Girolami [21] using a Bayesian multiclass multikernel algorithm. They reported a best test accuracy of 70% on a single run.
Performance with Individual and All Data Sources
We ran MKLdivdc, MKLdivconv, SimpleMKL and MKLRKDA on the overall set of 12 data sources, also evaluating performance on a uniformly weighted (averaged) composite kernel in addition to individual performance on each separate data source. In Table 1 we report the test set accuracy on each individual data source. The performance of MKLdivdc and MKLdivconv inclusive of all data sources achieves a test set accuracy of 73.36% and 71.01% respectively, consistently outperforming all individual performances and the uniformly weighted composite kernel (68.40%). Moreover, individual performance for MKLdivdc, SimpleMKL and MKLRKDA indicates that the most informative data sources are local sequence alignments (SW1 and SW2) and the amino acid composition (C). The performance with individual data sources for MKLdivdc, MKLdivconv, and SimpleMKL are almost the same since, for a fixed kernel, they use the same oneagainstall multiclass SVM.
From Table 1, performances of MKLdivdc and MKLdivconv with all the available data sources achieve test set accuracies of 73.36% and 71.01%, both of which outperform the stateofart performance 70% on a single run reported in [21] and other kernel learning methods including SimpleMKL (66.57%) and MKLRKDA (68.40%). The performance of the uniformly weighted kernel is 68.14% which is better than the performance 66.57% of SimpleMKL. This indicates that sparse L^{1}regularization does not necessarily yield better performance. The kernel weights λ of MKLdivdc, SimpleMKL, and MKLRKDA are shown in subfigures (b), (e) and (g) of Figure 2 which indicates that Amino Acid Composition (C), predicted secondary structure (S), Hypdrophobicity (H), and the last two data sources SW1 and SW2 are the most informative data sources, and the remaining data sources of H, P, V, and PseAA are less informative. As depicted in the subfigure (b) of Figure 2, MKLdivdc and MKLdivconv include some less informative data sources such as P, Z, L1, L4, L14, L30 etc., with small (but not zero) kernel weights. In contrast, as shown in (e) and (g) of Figure 2, SimpleMKL and MKLRKDA completely discard these less informative data sources. However, as shown in (d) and (f) of Figure 2, SimpleMKL and MKLRKDA achieve poorer performance, less than 70%, while MKLdivdc achieves 73.36% and MKLdivconv achieves 71.01%. This suggests that MKLdivdc provides a more reasonable balance over the entire set of data sources. This observation also suggests that achieving a sparsity among kernel weights does not necessarily guarantee good generalization performance since some available data sources may be weakly informative but may still carry some useful additional information.
Figure 2. Performance with all data sources on protein fold recognition. Test set accuracy of individual (bars) and all data sources (horizontal lines) on the protein fold recognition dataset: (a) MKLdivdc and MKLdivconv, where the solid line is the performance of MKLdivdc and the stardashed line is the performance of MKLdivconv; (d) SimpleMKL; (f) MKLRKDA. Kernel weights: (b) MKLdivdc, (c) MKLdivconv, (e) SimpleMKL and (g) MKLRKDA.
Performance with Sequential Addition of Data Sources
As mentioned above, the kernel weights learned by MKLdiv on all the data sources can provide useful insights into the significance of informative data sources. Hence, we further investigated the effect of sequentially adding data sources based on information from learned kernel weights in Tables 2 and 3. Without loss of generality, we take the kernel weights learned by MKLdivdc as an example.
Table 2. Effects of sequentially adding data sources
Table 3. Effects of sequentially adding data sources (continued)
We first report in Table 2 the effect of sequentially adding the sources in the order which was used in [3] and [21] and MKLdivdc and MKLdivconv consistently outperform the competitive kernel learning methods VBKC, SimpleMKL, MKLRKDA and the best performing SVM combination methodology stated in [3]. As suggested by the kernel weights of MKLdivdc in the subfigure (b) of Figure 2, the sequence alignment based data source SW1 is most informative, then S, then SW2 and so on. Hence, in Table 3 we further report the effect of sequentially adding data sources in this rank order. As shown in Table 3, there is a significant improvement over SW1SW2 in MKLdivdc when we sequentially add the data sources of amino acid composition (C) and predicted secondary structure (S). The performance of MKLdivdc keeps increasing until we include CSHPZ, giving the best performance of 75.19%. Although according to [4], the PseAA data sources are believed to contain more information than the conventional amino acid composition. The same behaviour appears for MKLdivconv. However, the MKLdivdc performance degenerates if we continue to add PseAA composition data sources and the same behaviour appears for MKLdivconv. Similar observations were made by [21] which suggests that PseAA measurements may carry noncomplementary information with the conventional amino acid compositions.
With regard to the best performance of MKLdivdc with the feature set SW1SW2CSHPZ, we display the corresponding kernel weights in Figure 3. We can see in Figure 3 that SimpleMKL and MKLRKDA almost eliminate the informative feature set HPZ while MKLdivdc and MKLdivconv include them into the composite kernel. The sparse L^{1}regularization of SimpleMKL and MKLRKDA accounts for the sparse weights of SimpleMKL and MKLRKDA.
Figure 3. Kernel weights on dominant data sources for protein fold recognition. Kernel weights on the dominant data sources SW1SW2CSHPZ which yields the best prediction on the protein fold recognition dataset: (a) MKLdivdc, (b) MKLdivconv, (c) SimpleMKL and (d) MKLRKDA.
Comparison of Running Time
To investigate the runtime efficiency of MKLdiv on protein fold recognition dataset, we list their CPU time in Tables 2 and 3. The running time (in seconds) is the term inside the parenthesis. The SILP approach for MKLRKDA is very efficient while SimpleMKL takes a bit longer. The reason could be that MKLRKDA essentially used the leastsquare loss for multiclass classification in contrast to the oneagainstall SVM used in SimpleMKL. Generally, more time is required to run the interior method for oneagainstall SVM than directly computing the solution of the leastsquare regression. The projected gradient descent method for MKLdivconv is also slower than MKLRKDA. It is to be expected that MKLdivconv converges faster than MKLdivdc since the DC algorithm for MKLdivdc is nonconvex and it needs to solve the subproblem (13) in each iteration of CCCP. Nevertheless, the price we paid in running time for MKLdivdc is worthwhile given its significantly better performance on the protein fold prediction problem.
Sensitivity against Parameter σ
The initial purpose of introducing σ is to avoid the singularity of the input kernel matrix or the output kernel matrix. However, in practice we found that, in the convex formulation MKLdivconv, values of σ have a great influence on performance for protein fold recognition. Hence, when we ran MKLdivconv, we always did cross validation over the training set to select the parameter σ. To see how sensitive the test set accuracy is with respect to σ, in Figure 4 we depicted the test set accuracy versus values of σ. In Figure 4 we can observe that the test set accuracy of MKLdivdc is relatively stable for small values of σ's. However, this is not the case for MKLdivconv and generally suggests that the parameter σ has a great impact on performance of MKLdivconv. This could be because the output kernel matrix K_{y }= YY^{⊤ }is of low rank (rank one in binary classification) and thus adding a small matrix σI_{n }in the formulation MKLdivconv could dramatically change the information of the output kernel matrix. In contrast, we can reasonably assume the input kernel matrices are nonsingular or not of low rank and the effect of adding a small matrix σI_{n }in the formulation MKLdivdc can be ignored.
Figure 4. Sensitivity against parameter σ for protein fold recognition. Test set accuracy versus different values of σ on the protein fold recognition dataset: (a) MKLdivdc and (b) MKLdivconv.
Extension of Investigation to Yeast Protein Classification
We next extend our investigation of MKLdivdc and MKLdivconv on a yeast membrane protein classification problem [23]. This binary classification task has 2316 examples derived from the MIPS comprehensive Yeast Genome Database (CYGD) (see [46]). There are eight kernel matrices http://noble.gs.washington.edu/proj/sdpsvm/ webcite. The first three kernels (K_{SW}, K_{B}, and K_{Pfam}) are respectively designed to measure the similarity of protein sequences using BLAST, SmithWaterman pairwise sequence comparison algorithms and a generalization of pairwise comparison method derived from hidden Markov models. The fourth sequencebased kernel matrix (K_{FFT}) incorporates information about hydrophobicity which is known to be useful in identifying membrane proteins, computed by Fast Fourier Transform. The fifth and sixth kernel matrices (K_{LI}, K_{D}) are respectively derived from linear and diffusion kernels based on proteinprotein interaction information. The seventh kernel matrix (K_{E}) is a Gaussian kernel encoding gene expression data. Finally, we added a noise kernel matrix K_{Ran }generated by first generating random numbers and then using a linear kernel.
The performance of MKLdivdc and MKLdivconv is evaluated by 10 random partitions of the data into a training and test set in a proportion of 4: 1. We report the receiver operating characteristic (ROC) score, which measures the overall quality of the ranking induced by the classifier, rather than the quality of a single point in that ranking. The first subfigure of Figure 5 shows the performance with individual kernels and the performance of MKLdivdc (the third to last bar), MKLdivconv (the next to last bar), and the uniformly weighted kernel (last bar). Specifically, MKLdivdc yields a ROC score of 0.9189 ± 0.0171 which is competitive with the result in [23]. MKLdivconv, however, achieved a ROC score of 0.9016 ± 0.0161 which is worse than MKLdivdc. The performance of MKLdivdc is also slightly better than the performance of the uniformly weighted kernel 0.9084 ± 0.0177 excluding the noise kernel and 0.8979 ± 0.0120 including the noise kernel. We also plot the kernel weights on (b) and (c) of Figure 5. As expected, in MKLdivdc the BLAST kernel (K_{B}) derived from the protein sequence similarity comparison is very informative which is consistent with [23]. The derived kernel weights also show that the interactionbased diffusion kernel is more informative than the expression kernel, which is consistent with [23]. Also, it is interesting to note that MKLdivdc shows that the noise kernel (K_{Ran}) is least informative. This is indicated by its individual ROC score: a ROC score around 0.5 corresponds to random ranking. The kernel weights of MKLdivconv indicate that the diffusion kernel (D) is the most important data source, and also suggest that Pfam and FFT are almost noninformative regardless of their good individual performances. For the kernel weights, MKLdivdc is more reasonable than MKLdivconv since MKLdivdc is more consistent with the individual data source's performance and MKLdivdc outperforms MKLdivconv using all data sources.
Figure 5. Performance of MKLdiv on yeast membrane protein. Performance on the yeast membrane protein function dataset: (a) average ROC score for individual data sources, using MKLdivdc and MKLdivconv, where the third bar to last (Alldc) is MKLdivdc, the second bar to last (Allconv) is MKLdivconv and the last bar (Averg) is the performance using uniformly weighted kernels. Kernel weights: (b) MKLdivdc and (c) MKLdivconv.
Conclusion
In this paper we developed a novel informationtheoretic approach to learning a linear combination of kernel matrices based on the KLdivergence [2428], especially focused on the protein fold recognition problem. Based on the different position of the input kernel matrix and the output kernel matrix in the KLdivergence objective, there are two formulations. The first one is a difference of convex (DC) problem termed MKLdivdc and the second formulation is a convex formulation called MKLdivconv. The sparse formulation for kernel learning based on discriminant analysis [18] was also established. Our proposed methods are able to achieve stateoftheart performance on the SCOP PDB40D benchmark dataset for protein fold recognition problem. In particular, MKLdivdc further improves the fold discrimination accuracy to 75.19% which is a more than 5% improvement over a competitive Bayesian probabilistic approach [21], SVM marginbased kernel learning methods [11], and the kernel learning based on discriminant analysis [18]. We further extended the investigation to the problem of yeast protein function prediction.
Generally, it is difficult to determine which criterion is better for multiple kernel combination since this problem is highly datadependent. For the informationtheoretic approaches MKLdivdc and MKLdivconv, although MKLdivdc is not convex and its DC algorithm tends to find a local minima, in practice we would recommend MKLdivdc for the following reasons. Firstly, as mentioned above MKLdivdc has a close relation with the kernel matrix completion problem using information geometry [27,28] and the maximization of the log likelihood of Gaussian Process regression [35], which partly explains the success of MKLdivdc. Secondly, we empirically observed that MKLdivdc outperforms MKLdivconv in protein fold recognition and yeast protein function prediction. Finally, as we showed in Figure 4, the performance of MKLdivconv is quite sensitive to the parameter σ and the choice of σ remains a challenging problem. MKLdivdc is relatively stable with respect to small values of σ and we can fix σ to be a very small number e.g. σ = 10^{5}. In future, we are planning to empirically compare performance with other existing kernel integration formulations on various datasets, and discuss convergence properties of the DC algorithm for MKLdivdc based on the theoretical results of [36].
Authors' contributions
YY and CC conceived the project. YY proposed and implemented the method, drafted the manuscript. KH joined the project and participated in the design of the study. All authors read and improved the manuscript.
Appendix
Appendix 1 – Column generation method for SILP
Here we briefly describe the column generation method (see e.g. [40]) for SILP (16) to solve the subproblem (15), i.e.
where , and S_{0}(α) = 2Tr(α^{⊤ }α). The basic idea is to compute the optimum (λ, γ) by linear programming for a restricted subset of constraints, and update the constraint subset based on the obtained suboptimal (λ, γ). More precisely, given restricted constraints {α_{p }: p = 1, ..., P}, first we find the intermediate solution (λ, γ) by the following linear programming optimization with P linear constraints
This problem is often called the restricted master problem. Then, we find the next constraint with the maximum violation for the given intermediate solution (λ, γ), i.e.
If the optimizer α * of the above equation satisfies then the current intermediate solution (λ, γ) is optimal for the optimization (19). Otherwise α* should be added to the restriction set. We repeat the above iteration until convergence which is guaranteed to be globally optimal, see e.g. [14,40]. In a similar fashion to the convergence criterion in [14], the algorithm stops when
For instance, the threshold ε is usually chosen to be 5 × 10^{4}.
Appendix 2 – Sparse formulation of kernel learning based on discriminant analysis
In this appendix we show that kernel learning for regularized discriminant analysis [18] is closely related to sparse regularization. To see this, consider the following algorithm
Using the fact [31] that min , the above equation is identical to
The equivalence between the above algorithm and RKDA kernel learning becomes clear if we formulate its dual problem as follows:
Theorem 2 Let , I_{n }be the identity matrix and 1_{n }be an ndimensional column vector of all ones. Define , and for any i ∈ ℕ_{n}. Then, the dual problem of algorithm (22) can be written as
where .
Proof: Taking the minimization of b first, algorithm (22) yields . Then, algorithm (22) can be further rewritten as
Here, for any ℓ and i, which can be further represented by . Then, letting for any i and solving the standard Lagrangian formulation of (23) with Lagrangian variables α yields
Now, replacing α_{i }by μα_{i }and letting completes the argument. □
Let n_{ }and n_{+ }denote the number of samples in class +1 and 1. If we redefine the class indicator output y, for any i ∈ ℕ_{n }by y_{i }= if x_{i }is in class +1 otherwise , then the class indicator output reduces to the vector a defined in [18] for binary classification, i.e.
Now we turn our attention to multiclass classification. To this end, consider
Using the above argument for binary classification it is easy to check its dual problem is as follows
where . Let n_{c }denote the number of samples in class c. If we redefine the class indicator matrix Y, for any i ∈ ℕ_{n }and c ∈ ℕ_{C }by if y_{i }= c, otherwise , then the class indicator matrix reduces to the matrix H defined in [18] for multiclass classification, i.e.
Now we can see that the dual problem of algorithm (24) is exactly the same as the formulation (see equation (36) in [18]) of RKDA kernel learning.
Acknowledgements
We would like to thank the referees for their constructive comments and suggestions which greatly improve the paper. We also thank Prof. Mark Girolami, Dr. Theodoros Damoulas and Dr. Simon Rogers for stimulating discussions. This work is supported by EPSRC grant EP/E027296/1.
References

Baker D, Sali A: Protein structure prediction and structural genomics.
Science 294:9396. PubMed Abstract  Publisher Full Text

Jones DT, et al.: A new approach to protein fold recognition.
Nature 1992, 358:8689. PubMed Abstract  Publisher Full Text

Ding C, Dubchak I: Multiclass protein fold recognition using support vector machines and neural networks.
Bioinformatics 2001, 17:349358. PubMed Abstract  Publisher Full Text

Shen HB, Chou KC: Ensemble classifier for protein fold pattern recognition.
Bioinformatics 2006, 22:17171722. PubMed Abstract  Publisher Full Text

Andreeva A, et al.: SCOP database in 2004: refinements integrate strucuture and sequence family data.
Nucleic Acids Res 2004, 32:226229. Publisher Full Text

Lo Conte L, et al.: SCOP: a structural classification of protein database.
Nucleic Acids Res 2000, 28:257259. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Chou K, Zhang C: Prediction of protein structural classes.
Crit Revi Biochem Mol Biol 1995, 30:275349. Publisher Full Text

Dubchak I, et al.: Prediction of protein folding class using global description of amino acid sequence.
Proc Natl Acad Sci 1995, 92:87008704. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Schölkopf B, Smola AJ: Learning with Kernels. The MIT Press, Cambridge, MA, USA; 2002.

ShaweTaylor J, Cristianini N: Kernel Methods for Pattern Analysis. Cambridge university press; 2004.

Lanckriet GRG, Cristianini N, Bartlett PL, Ghaoui LE, Jordan ME: Learning the kernel matrix with semidefinite programming.

Bach F, Lanckriet GRG, Jordan MI: Multiple kernel learning, conic duality and the SMO algorithm.
Proceedings of the Twentyfirst International Conference on Machine Learning (ICML) 2004.

Rakotomamonjy A, Bach F, Canu S, Grandvalet Y: More efficiency in multiple kernel learning.
Proceedings of the 24th International Con ference on Machine Learning (ICML) 2007.

Sonnenburg S, Rätsch G, Schäfer C, Schölkopf B: Large scale multiple kernel learning.

Bach F: Consistency of the group Lasso and multiple kernel learning.

Ying Y, Zhou DX: Learnability of Gaussians with flexible variances.

Lin Y, Zhang H: Component selection and smoothing in multivariate nonparametric regression.
Annals of Statistics 2006, 34:22722297. Publisher Full Text

Ye J, Ji S, Chen J: Multiclass discriminant kernel learning via convex programming.

Ye J, et al.: Heterogeneous data fusion for Alzheimer's disease study.
The Fourteenth ACM SIGKDD International Conference On Knowledge Discovery and Data Mining (SIGKDD) 2008.

Ji S, Sun L, Jin R, Ye J: Multilabel multiple kernel learning.
The TwentySecond Annual Conference on Neural Information Processing Systems (NIPS) 2008.

Damoulas T, Girolami M: Probabilistic multiclass multikernel learning: On protein fold recognition and remote homology detection.
Bioinformatics 2008, 24(10):12641270. PubMed Abstract  Publisher Full Text

Girolami M, Zhong M: Data integration for classification problems employing gaussian process priors. In Advances in Neural Information Processing Systems. Cambridge, MA: MIT Press; 2007.

Lanckriet GRG, De Bie T, Cristianini N, Jordan MI, Noble WS: A statistical framework for genomic data fusion.
Bioinformatics 2004, 20(16):26262635. PubMed Abstract  Publisher Full Text

Lawrence ND, Sanguinetti G: Matching kernels through KullbackLeibler divergence minimization. In Technical Report CS04–12. Department of Computer Science, University of Sheffield; 2004.

Davis JV, Kulis B, Jain P, Sra S, Dhillon IS: Informationtheoretic metric learning.
Proceedings of the 24th International Conference on Machine Learning 2007.

Sun L, Ji S, Ye J: Adaptive diffusion kernel learning from biological networks for protein function prediction.
BMC Bioinformatics 2008, 9:162. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Tsuda K, Akaho S, Asai K: The em algorithm for kernel matrix completion with auxiliary data.
Journal of Machine Learning Research 2003, 4:6781. Publisher Full Text

Kato T, Tsuda K, Asai K: Selective integration of multiple biological data for supervised network inference.
Bioinformatics 2005, 21:24882495. PubMed Abstract  Publisher Full Text

Aronszajn N: Theory of reproducing kernels.
Trans Amer Math Soc 1950, 68:337404. Publisher Full Text

Cristianini N, ShaweTaylor J, Elisseeff A: On kerneltarget alignment. In Advances in Neural Information Processing Systems 14. Cambridge, MA: MIT Press; 2002.

Micchelli CA, Pontil M: Learning the kernel function via regularization.

Hastie T, Tibshirani R, Friedman J: Elements of Statistical Learning: Data Mining, Inference, and Prediction. New York, Springer; 2001.

Vandenberghe L, Boyd S, Wu S: Determinant maximization with linear matrix inequality constraints.
SIAM J Matrix Analysis and Applications 1998, 19:499533. Publisher Full Text

Amari S: Information geometry of the EM and em algorithms for neural networks.
Neural Networks 1995, 8:13791408. Publisher Full Text

Rasmussen CE, Williams C: Gaussian Processes for Machine Learning. the MIT Press; 2006.

Tao PD, An LTH: A D.C. optimization algorithm for solving the trust region subproblem.
SIAM J Optim 1998, 8:476505. Publisher Full Text

Yuille AL, Rangarajan A: The concave convex procedure.
Neural Computation 2003, 15:915936. PubMed Abstract  Publisher Full Text

Borwein JM, Lewis AS: Convex Analysis and Nonlinear Optimization: Theory and Examples. CMS Books in Mathematics, SpringerVerlag, New York; 2000.

Nesterov Y: Introductory Lectures on Convex Optimization: A Basic Course. Springer; 2003.

Hettich R, Kortanek KO: Semiinfinite programming: theory, methods, and applications.
SIAM Review 1993, 3:380429. Publisher Full Text

Bleakley K, Biau G, Vert JP: Supervised reconstruction of biological networks with local models.
Bioinformatics 2007, 23:5765. PubMed Abstract  Publisher Full Text

Yamanishi Y, Vert JP, Kanehisa M: Protein network inference from multiple genomic data: a supervised approach.
Bioinformatics 2004, 20:363370. Publisher Full Text

Kondor RI, Lafferty JD: Diffusion kernels on graphs and other discrete structures.
Proceedings of the Nineteenth International Conference on Machine Learning 2002.

The SimpleMKL Toolbox [http://asi.insarouen.fr/enseignants/~arakotom/code/mklindex.html] webcite

Liao L, Noble WS: Combining pairwise sequence similarity and support vector machine for detecting remote protein evolutionary and structural relationships.
J Comput Biol 2003, 6:857868. Publisher Full Text

Mewes HW, et al.: MIPS: a database for genomes and protein sequences.
Nucleic Acids Res 2000, 28:3740. PubMed Abstract  Publisher Full Text  PubMed Central Full Text