September 4, 2009

12:3

WSPC - Proceedings Trim Size: 9.75in x 6.5in

23-Sato1

GRADIENT-BASED OPTIMIZATION OF HYPERPARAMETERS
FOR BASE-PAIRING PROFILE LOCAL ALIGNMENT KERNELSKENGO SATO1,2,3
sato-kengo@aist.go.jp

YUTAKA SAITO3
saito@dna.bio.keio.ac.jp

YASUBUMI SAKAKIBARA3
yasu@bio.keio.ac.jp
1 Japan Biological Informatics Consortium (JBIC), 2–45 Aomi, Koto-ku, Tokyo 135–

8073, Japan
2 Computational Biology Research Center (CBRC), National Institute of Advanced

Industrial Science and Technology (AIST), 2–42 Aomi, Koto-ku, Tokyo 135–0064,
Japan
3 Department of Biosciences and Informatics, Keio University, 3–14–1 Hiyoshi,
Kohoku-ku, Yokohama, Kanagawa 223–8522, Japan
We have recently proposed novel kernel functions, called base-pairing proﬁle local alignment (BPLA) kernels for discrimination and detection of functional RNA sequencesusing SVMs. We employ STRAL’s scoring function which takes into account sequencesimilarities as well as upstream and downstream base-pairing probabilities, which enables us to model secondary structures of RNA sequences. In this paper, we develop amethod for optimizing hyperparameters of BPLA kernels with respect to discrimination
accuracy using a gradient-based optimization technique. Our experiments show that the
proposed method can ﬁnd a nearly optimal set of parameters much faster than the gridsearch on all parameter combinations.
Keywords:
RNAs.

support vector machines; kernel functions; gradient decent; non-coding

1. Introduction
Recent research has discovered crucial roles of non-coding RNAs (ncRNAs) in cells,
including post-transcriptional gene regulation, maturation of mRNAs, rRNAs andtRNAs. Most functional ncRNAs form secondary structures with hydrogen bonds
such as Watson-Crick base-pairs (A-U and G-C) and Wobble base-pairs (G-U). It is
well-known that functions of ncRNAs are deeply related to their secondary structures rather than primary sequences. For example, transfer RNAs registered in the
Rfam database [?] have only 44 % average sequence identity, whereas they arefolded into conserved clover leaf secondary structures. Therefore, ncRNAs are modeled by means of their secondary structures for computational ncRNA analysis. Oneof such models is stochastic context-free grammars (SCFGs), which can represent
RNA secondary structures without pseudoknots [?, ?]. These grammatical methods
have succeeded in modeling typical secondary structures of RNAs, and are com-

September 4, 2009

2

12:3

WSPC - Proceedings Trim Size: 9.75in x 6.5in

23-SatoK. Sato, Y. Saito, & Y. Sakakibara

monly used for structural alignment of RNA sequences. However, such stochastic
models are not capable of discriminating the member sequences of RNA families
from nonmembers with suﬃciently high accuracy to detect non-coding RNA regionsin genome sequences.
On the other hand, a compelling attraction of support vector machines (SVMs)
and other kernel methods has been increasing in many research ﬁelds, including
bioinformatics [?], because of their robustness. For example, Saigo et al. have proposed local alignment kernels for amino acid sequences [?, ?], and their kernels
with SVMs have outperformed other state-of-the-art methods in benchmarks forremote homology detection of amino acid sequences. Therefore, we have been considering using kernel methods, including SVMs, for the analysis of functional ncRNAs [?, ?, ?].
We have recently proposed novel kernel functions, called base-pairing proﬁle localalignment (BPLA) kernels for discrimination and detection of functional RNA sequences using SVMs [?]. We extend the concept of local alignment kernels in such a
way that it can handle RNA sequences using the scoring function in STRAL, which
is a practical aligner for non-coding RNAs [?]. The local alignment kernels measure
the similarity between two sequences by summing the scores over all possible localalignments with gaps, whereas STRAL’s scoring function takes into account sequence similarities as well as upstream and downstream base-pairing probabilities,
which enables us to model secondary structures of RNA sequences. Our experiments
have proved powerful discrimination ability of BPLA kernels for functional RNA
families. Furthermore, we have demonstrated the applicability of our kernels tothe problem of genome-wide search of snoRNA families in the C. elegans genome,
and have conﬁrmed that the expression is valid in 14 out of 48 of our predicted
candidates by using qRT-PCR (p-value = 2.8 × 10−32 ).
Our previous experiments showed that parameters for the BPLA kernels arequite diﬀerent depending on what kind of data (ﬂanking regions or not) and whatfamily of sequences should be trained. Our previous study employed a brute force
approach in which optimal parameters for each case were calibrated by the gridsearch on all parameter combinations from selected search space for each parameter.
It is easily predicted that this approach should require exponential time in thenumber of parameters to be optimized.
Therefore, we propose a gradient-based parameter optimization for BPLA kernels. This work is based on Keerthi et al. [?], in which gradients of accuracy measures(e.g. AUC) with respect to hyperparameters of kernel functions are calculated, thenthese hyperparameters are optimized with a numerical optimization algorithm suchas the quasi-Newton’s method. We develop a dynamic programming algorithm forcalculating the gradients of BPLA kernels with respect to their parameters to be optimized by the Keerthi’s method. Our experiments show that the proposed method
can ﬁnd a nearly optimal set of parameters much faster than the grid search on allparameter combinations.
This paper is organized as follows: in section 2, we ﬁrst introduce BPLA kernels

September 4, 2009

12:3

WSPC - Proceedings Trim Size: 9.75in x 6.5in

23-Sato

BPLA Kernels

3

for functional RNA analysis and dynamic programming algorithms for calculating
BPLA kernels and their gradient with respect to their hyperparameters. Then, we
brieﬂy describe the Keerthi’s gradient-based optimization method. In section 3,
computational experiments are shown to conﬁrm the eﬃciency of our method, andwe discuss our framework in section 4. Finally, we conclude this paper in section 5.
2. Methods
2.1. Base-pairing proﬁle local alignment kernels
Local alignment kernels can calculate the similarity between a pair of sequences
by taking into account only sequence homology [?]. However, since it is well-known
that ncRNAs form secondary structures which are related to their functions, it isnecessary to consider the secondary structures of ncRNAs in order to compare two
RNA sequences. To do this, we have proposed stem kernels [?] and their approximation using directed acyclic graphs (DAGs) of ncRNAs [?]. The stem kernels canconsider RNA secondary structures strictly, but require huge computational time forpractical problems. Diﬀering from them, BPLA kernels employ a kind of summary
information of secondary structures, called base-pairing proﬁles, as well as STRAL,
which is a practical aligner for non-coding RNAs [?].
For each sequence, we ﬁrst calculate a base-pairing probability matrix using
the McCaskill algorithm [?] or the inside-outside algorithm [?]. The base-pairing
probability matrix for a sequence x consists of the base-pairing probabilities Pij
that the i-th and the j-th nucleotides of x form a base pair, which is deﬁned as:∑p(y|x)Iij (y),(1)
Pij = E[Iij |x] =y∈Y(x)

where Y(x) is an ensemble of all possible secondary structures of x, p(y|x) is theposterior probability of a secondary structure y given x, and Iij (y) is an indicatorfunction, which equals 1 if the i-th and the j-th nucleotides form a base-pair in y and0 otherwise. In this study, we employ the Vienna RNA package [?] for computing
the expected counts (1) using the McCaskill algorithm.
Subsequently, for each position i of x, the base-pairing probabilities are sum∑
marized into tree kinds of sums: the probability Pilef t = j>i Pij that a pair is∑formed with one of the downstream nucleotides, the probability Piright = j<i Pji
that a pair is formed with one of the upstream nucleotides, and the probability
Piunpair = 1 − (Pilef t + Piright ) that the nucleotide is unpaired. A probability distribution consisting of these three probabilities is called a base-pairing proﬁle [?].
In accordance with STRAL [?], we deﬁne the match score between two nucleotides xi and yj using base-pairing proﬁles as follows:
S(xi , yj ) =α(Sstruct ) + Sseq(√
) √√lef t lef tright right=α
Pxi Pyj + Pxi Pyj+ Pxunpair
Pyunpair
S 0 (xi , yj ), (2)ijSeptember 4, 2009

4

12:3

WSPC - Proceedings Trim Size: 9.75in x 6.5in

23-SatoK. Sato, Y. Saito, & Y. Sakakibara

Initialization:for i ∈ {0, . . . , |x|} and j ∈ {0, . . . , |y|} do
M (i, 0) = IX (i, 0) = IY (i, 0) = RX (i, 0) = RY (i, 0) = 0
M (0, j) = IX (0, j) = IY (0, j) = RX (0, j) = RY (0, j) = 0
end forIteration:for i ∈ {1, . . . , |x|} and j ∈ {1, . . . , |y|} do
M (i, j) = eβS(xi ,yj ) (1 + IX (i − 1, j − 1) + IY (i − 1, j − 1) + M (i − 1, j − 1))
IX (i, j) = eβg M (i − 1, j) + eβe IX (i − 1, j)
IY (i, j) = eβg (M (i, j − 1) + IX (i, j − 1)) + eβe IY (i, j − 1)
RX (i, j) = M (i − 1, j) + RX (i − 1, j)
RY (i, j) = M (i, j − 1) + RX (i, j − 1) + RY (i, j − 1)
end forTermination:
K(x, y) = 1 + RX (|x|, |y|) + RY (|x|, |y|) + M (|x|, |y|)Fig. 1.

The dynamic programming of the BPLA kernel.

where α is a weight parameter for structural information, and S 0 is a substitution
scoring function between two nucleotides. In this study, we use a modiﬁed RIBOSUM 85-60 [?] as the substitution scoring matrix, in which its smallest eigenvalue issubtracted from each of its diagonal elements in order to make it positive semidefinite. Then, we employ equation (2) as the match scores for the local alignment
kernels between two sequences x and y, which is deﬁned as follows:∑
K(x, y) =exp βs(x, y, π),(3)π∈Π(x,y)where β ≥ 0 is a concentration parameter, Π(x, y) is a set of all possible localalignments of x and y, and s(x, y, π) is a score of alignment π deﬁned by:∑na,b (x, y, π) · S(a, b) − ng (x, y, π) · g − ne (x, y, π) · e,(4)s(x, y, π) =a,bwhere na,b (x, y, π) is the number of times that nucleotides a is aligned with nucleotides b in π, ng (x, y, π) and ne (x, y, π) are the numbers of gap opens and gap
extensions, respectively, and g and e are penalties for gap opens and gap extensions, respectively. We call equation (3) with the match scoring function (2) the
base-pairing proﬁle local alignment (BPLA) kernel. Equation (3) can be computedby the dynamic programming technique shown as Figure 1.
The hyperparameter optimization algorithm, described in Section 2.2, employs
the gradient of K(x, y) with respect to parameters α, β, g and e, which can also becomputed by the dynamic programming. Figure 2 shows the dynamic programmingof calculating the derivative with respect to α. The derivatives with respect to theother parameters can be computed as well as α.
In order to avoid length bias, we usually normalize kernel values. Therefore, wewould like to optimize hyperparameters with respect to the normalized kernels Kn
by using its derivative Kn0 :September 4, 2009

12:3

WSPC - Proceedings Trim Size: 9.75in x 6.5in

23-Sato

BPLA Kernels

5

Initialization:for i ∈ {0, . . . , |x|} and j ∈ {0, . . . , |y|} do
∂
∂
∂
∂
∂
∂α M (i, 0) = ∂α IX (i, 0) = ∂α IY (i, 0) = ∂α RX (i, 0) = ∂α RY (i, 0) = 0
∂
∂
∂
∂
∂
∂α M (0, j) = ∂α IX (0, j) = ∂α IY (0, j) = ∂α RX (0, j) = ∂α RY (0, j) = 0
end forIteration:for i ∈ {1, . . . , |x|} and j„∈ {1, . . . , |y|} do«qq
lef t lef tright right
∂
P
P
M
(i,
j)
=
βM
(i,
j)
P
P
+
xi
xi
yj
yj
∂α
∂
∂
∂
+eβS(xi ,yj ) ( ∂αIX (i − 1, j − 1) + ∂α
IY (i − 1, j − 1) + ∂αM (i − 1, j − 1))
βg ∂
βe ∂
∂
I
(i,
j)
=
e
M
(i
−
1,
j)
+
e
I
(i
−
1,
j)
∂α X
∂α
∂α X
βg ∂
βe ∂
∂
∂
∂α IY (i, j) = e ( ∂α M (i, j − 1) ∂α IX (i, j − 1)) + e ∂α IY (i, j − 1)∂∂∂
∂α RX (i, j) = ∂α M (i − 1, j) + ∂α RX (i − 1, j)∂∂∂∂∂α RY (i, j) = ∂α M (i, j − 1) + ∂α RX (i, j − 1) + ∂α RY (i, j − 1)

end forTermination:∂∂∂∂∂α K(x, y) = ∂α RX (|x|, |y|) + ∂α RY (|x|, |y|) + ∂α M (|x|, |y|)Fig. 2.

The dynamic programming of the derivative of the BPLA kernel with respect to α.

K(x, y)Kn (x, y) = √K(x, x)K(y, y)()
Kn (x, y) K 0 (x, x) K 0 (y, y)
K 0 (x, y)−
Kn0 (x, y) = √+.
2
K(x, x)K(y, y)K(x, x)K(y, y)2.2. Gradient-based optimization for SVMs
Tuning hyperparameters of kernel functions is known as an important task toachieve enough accuracy in practical problems. Although brute force approaches
like the grid search are usually employed for such tasks, huge computational time
will obviously be required if the search space of hyperparameters is large. Keerthi
et al. have proposed a more eﬃcient method which optimizes hyperparameters via
the gradient-based optimization [?].
We assume that n training data {(xi , yi )}ni=0 (yi ∈ {−1, +1}) are given. In SVM
training, we optimize the parameters w and b which satisfy the following condition:∑1
min ||w||2 + Cl(oi , yi ),w,b 2iwhere l is a loss function such as the hinge-loss l(oi , yi ) = max{0, 1 − yi oi }, C isa weight for the loss function, and oi is the output from the SVM for the trainingdata:∑oi = w · φ(xi ) − b =αj yj K(xi , xj ) − b.
j

September 4, 2009

6

12:3

WSPC - Proceedings Trim Size: 9.75in x 6.5in

23-SatoK. Sato, Y. Saito, & Y. Sakakibara

Here, K(xi , xj ) = φ(xi ) · φ(xj ) is called a kernel function which can implicitly deﬁne
a high-dimensional mapping φ. The terms αj s constitute a dual representation forw. We can obtain a linear system from the Karush-Kuhn-Tucker (KKT) condition:
( )
α
P
= q.
(5)
b
The training data and their coeﬃcients α are divided into three parts: a part I0 for
αi = 0, a part Iu for 0 < αi < C, and a part Ic for αi = C. The linear system (5)can be rewritten as follows:)( ) ()(
eu − Ωuc αcαuΩuu −yu,=α0 = 0, αc = Cec ,−yuT 0bycT αcwhere Ωij = yi yj Kij and eu is a unit vector. These α and b can be computed bysolving quadratic programming.
For given validation data {(x̄l , ȳl )}, we can calculate:∑ōl =αi yi K(x̄l , xi ) − b,iand can rewrite it with a matrix form:ōl = ψlT β,
where β is a vector containing α and b, and ψl is a vector containing yi Kli , i =1, . . . , n and −1 as the last element (corresponding to b).
Our objective is to maximize f (which is a kind of accuracy measures such as
speciﬁcity, sensitivity, F-value, or AUC) with respect to hyperparameters θ of agiven kernel function K. Here, f˙ denotes a derivative of f with respect to θ, thatis, f˙ = ∂f∂θ , which can be computed as:
)T
(
(
)
)
(
)
∑
∑ (
∑
T
δl ψ̇lβ,δl ψ T P −1 q̇ − Ṗ β + ψ̇l β = dT q̇ − Ṗ β +δl ō˙ l =
f˙ =
l

l

l

l

(6)
where d is the solution of
PTd =

(
∑

)
δl ψl

,

(7)

l

and
δl =

∂f
.
∂ ōl

(8)

The steps for calculating the derivative f˙ is as follows. First, compute δl from
(8). Then, solve (7) for d. Equation (7) can eﬃciently be solved by the conjugate
gradient method. Then, compute f˙ for each θ from (6). After f and its gradientsf˙ are computed, we can apply the gradient-based optimization by the L-BFGS [?](or the L-BFGS-B [?] with box-constraints).

September 4, 2009

12:3

WSPC - Proceedings Trim Size: 9.75in x 6.5in

23-Sato

BPLA Kernels

7

We optimize the hyperparameters of BPLA kernels by means of AUC, that is,
f = auc deﬁned as:∑∑1
sign(ōi − ōj ),(9)
auc = + −n n i:ȳ =+1 j:ȳ =−1
i

+

j

−

where n and n are the numbers of the positive and negative examples in thevalidation dataset, respectively. The step function sign(x) is approximated by using1
the sigmoidal function s(x) = 1+exp(−σx)
to make (9) diﬀerentiable.
3. Results
In order to illustrate the eﬃciency of our gradient-based optimization, we compared
the proposed method with the brute-force grid search. In each experiment, we ﬁrst
calculated AUC for all the parameter combinations in the predeﬁned search space,
and then performed our optimization procedure with the corresponding box constraints to evaluate whether it could approach the AUC value at the optimal grid
point. We used the same data set as in our previous paper [?], which includes ﬁve
ncRNA families taken from the Rfam database [?]. For each family, 100 sequenceswere chosen as positive samples in such a way that the pairwise identity is not above
80% for any pair of sequences, and as negative samples, 100 randomly shuﬄed sequences with the same dinucleotide composition as the positives were generated.
AUC was calculated using four-fold cross-validation, and its gradient was computedby summing over the four data partitions.
In our ﬁrst experiment, we attempted to simultaneously optimize all the fourparameters α, β, g, and e in BPLA kernels, as well as C in the loss functionsof SVMs. To do this, we deﬁned the search space with the ﬁve-dimensional box
constraints: 0 < α ≤ 20, 0 < β ≤ 0.2, 0 < g ≤ 30, 0 < e ≤ 1, and 0 < C ≤40. Subsequently, we carried out the grid search in which we tried six values foreach parameter at even intervals, i.e. α ∈ {1 × 10−9 , 4, 8, 12, 16, 20}, β ∈ {1 ×
10−9 , 0.04, 0.08, 0.12, 0.16, 0.20}, and similarly for the others. Then, we performed
the gradient-based optimization. Since the proposed method might fall into one of
local maxima, resulting AUC depends on initial values of the parameters. Therefore,
we executed the optimization procedure from various start points. These start points
were chosen by combining the ﬁrst or third quantile of each parameter, i.e. α ∈{5, 15}, β ∈ {0.05, 0.15}, and similarly for the others, and we evaluated the mean
AUC value of the 25 trials.
Table 1 presents the experimental results. Although the grid search attained
the best discrimination for all the ncRNA families tested, it took the huge computational time because of brute-force calculation for 7776 (= 65 ) combinationsof the parameters. Note that even the relatively coarse grid with only six values
for one dimension required such massive computation. This clearly shows that thegrid search is not a practical way to adapt the parameters of BPLA kernels against
various kinds of data set. On the other hand, our gradient-based optimization successfully found the nearly optimal set of the parameters much faster than the grid

September 4, 2009

8

12:3

WSPC - Proceedings Trim Size: 9.75in x 6.5in

23-SatoK. Sato, Y. Saito, & Y. Sakakibara
Table 1.
search.

Comparison between the gradient-based optimization and the grid

Family5S rRNAtRNA
C/D snoRNA
H/ACA snoRNA
miRNA

gradient-based optimization
N
AUC
time (h)240.997 (0.956)
1.8
15
0.993 (0.957)
0.6
25
0.843 (0.773)
1.8
23
0.879 (0.832)
3.5
20
0.995 (0.971)
1.1

N
7776
7776
7776
7776
7776

grid search
AUC
time (h)
0.997
149.0
0.995
67.5
0.887
124.8
0.906
248.0
0.995
89.0

Note: Family: name of the target ncRNA family. N: number of parameter sets tried.
AUC: area under the ROC curve, the value in parentheses is AUC calculated at thestart point of optimization. time: computational time on a 2.8 GHz AMD Opteronprocessor. In the columns of the gradient-based optimization, mean values of 25
trials from the diﬀerent start points are shown.

search. The AUC was substantially improved compared to the start point on average, and convergence was obtained within a few dozen of iterations. In particular,
we achieved the comparable discrimination ability to the grid search for 5S rRNAs,
tRNAs and miRNAs showing the eﬀectiveness of the proposed method.
The second experiment was organized to visually demonstrate how our method
works in the parameter space. The idea was that we would show a trajectory of
optimization on the heat map of AUC. For this purpose, we focused only on α andβ, and ﬁxed the other parameters. Thus, we deﬁned the same box constraints for
α and β as in the ﬁrst experiment, whereas we consistently used g = 27, e = 1,
and C = 1 for the others. We calculated AUC over the ﬁner grid with 100 values
for each of α and β, and initiated the gradient-based optimization. Figure 3 shows
the corresponding heat map and the trajectory plot. We present the typical results
for the C/D snoRNA data set. The parameters α and β were moved sharply tothe direction of higher AUC, and well calibrated along with the iterations. Indeed,
the optimization in this case reached the region of highest AUC, and successfully
converged at the global optima in the parameter space.
4. Discussion
Our current implementation of the hyperparameter optimizer for BPLA kernels
does not optimize the substitution scoring function for RNAs (S 0 in (2)) althoughderivatives of substitution scores can easily be calculated. In fact, Saigo et al. haveproposed a method for optimizing substitution scores for local alignment kernels
using the gradient decent with respect to a conﬁdence score for remote homology
search [?] (not with respect to the SVM classiﬁcation). However, the optimization of
the substitution scoring matrix using the L-BFGS-B, which our method employed,cannot hold the positive semideﬁniteness of S 0 and BPLA kernels, which is theessential basis of kernel methods. This kind of optimization problems can be solved
by a method for semideﬁnite programming (SDP) [?]. Because of small eﬀect ofchanging substitution scoring matrices in our method (data not shown), we ﬁxed
the substitution scoring matrix to a modiﬁed RIBOSUM 85-60 [?], in which its

September 4, 2009

12:3

WSPC - Proceedings Trim Size: 9.75in x 6.5in

23-Sato

BPLA Kernels

9

AUC
0.9

0.8
0.15
0.7



0.10

0.6

0.5
0.05
0.4

0.00

0.3
0

5

10

15

Fig. 3. Trajectory of parameter optimization on the heat map of AUC. AUC values evaluated for
the C/D snoRNA data set are represented by gradation over the parameter space. The optimization
was initiated from (α, β) = (6.50, 0.05), and terminated at (α, β) = (1.01, 0.09) after 35 iterations.

smallest eigenvalue is subtracted from each of its diagonal elements in order tomake it positive semideﬁnite.
5. Conclusion
We have developed a method for optimizing hyperparameters of BPLA kernels withrespect to discrimination accuracy using a gradient-based optimization technique.
Our experiments showed that the proposed method can ﬁnd a nearly optimal setof parameters much faster than the grid search on all parameter combinations.
We are planing to apply BPLA kernels to comparative analysis of functionalRNAs, which can consider the sequence conservation between species as well as the
structural conservation. For this purpose, the proﬁle-proﬁle comparison should be
implemented.
Acknowledgment
This work was supported in part by a grant from “Functional RNA Project” fundedby the New Energy and Industrial Technology Development Organization (NEDO)
of Japan, and was also supported in part by Grant-in-Aid for Scientiﬁc Research on
Priority Area “Comparative Genomics” No. 17018029 from the Ministry of Education, Culture, Sports, Science and Technology of Japan. We thank our colleagues

September 4, 2009

10

12:3

WSPC - Proceedings Trim Size: 9.75in x 6.5in

23-SatoK. Sato, Y. Saito, & Y. Sakakibara

from the RNA Informatics Team at the Computational Biology Research Center(CBRC) for fruitful discussions.