1 Introduction↩︎

Learning from biological sequences is important for a variety of scientific fields such as evolution [1] or human health [2]. In order to use classical statistical models, a first step is often to map sequences to vectors of fixed size, while retaining relevant features for the considered learning task. For a long time, such features have been extracted from sequence alignment, either against a reference or between each others [3]. The resulting features are appropriate for sequences that are similar enough, but they become ill-defined when sequences are not suited to alignment. This includes important cases such as microbial genomes, distant species, or human diseases, and calls for alternative representations [4].

String kernels provide generic representations for biological sequences, most of which do not require global alignment  [5]. In particular, a classical approach maps sequences to a huge-dimensional feature space by enumerating statistics about all occuring subsequences. These subsequences may be simple classical \(k\)-mers leading to the spectrum kernel [6], \(k\)-mers up to mismatches [7], or gap-allowing subsequences [8]. Other approaches involve kernels based on a generative model [9], [10], or based on local alignments between sequences [11] inspired by convolution kernels [12], [13].

The goal of kernel design is then to encode prior knowledge in the learning process. For instance, modeling gaps in biological sequences is important since it allows taking into account short insertion and deletion events, a common source of genetic variation. However, even though kernel methods are good at encoding prior knowledge, they provide fixed task-independent representations. When large amounts of data are available, approaches that optimize the data representation for the prediction task are now often preferred. For instance, convolutional neural networks [14] are commonly used for DNA sequence modeling [15]–[17], and have been successful for natural language processing [18]. While convolution filters learned over images are interpreted as image patches, those learned over sequences are viewed as sequence motifs. RNNs such as long short-term memory networks (LSTMs) [19] are also commonly used in both biological [20] and natural language processing contexts [21], [22].

Motivated by the regularization mechanisms of kernel methods, which are useful when the amount of data is small and are yet imperfect in neural networks, hybrid approaches have been developed between the kernel and neural networks paradigms [23]–[25]. Closely related to our work, the convolutional kernel network (CKN) model originally developed for images [26] was successfully adapted to biological sequences in [27]. CKNs for sequences consist in a continuous relaxation of the mismatch kernel: while the latter represents a sequence by its content in \(k\)-mers up to a few discrete errors, the former considers a continuous relaxation, leading to an infinite-dimensional sequence representation. Finally, a kernel approximation relying on the Nyström method [28] projects the mapped sequences to a linear subspace of the RKHS, spanned by a finite number of motifs. When these motifs are learned end-to-end with backpropagation, learning with CKNs can also be thought of as performing feature selection in the—infinite dimensional—RKHS.

In this paper, we generalize CKNs for sequences by allowing gaps in motifs, motivated by genomics applications. The kernel map retains the convolutional structure of CKNs but the kernel approximation that we introduce can be computed using a recurrent network, which we call recurrent kernel network (RKN). This RNN arises from the dynamic programming structure used to compute efficiently the substring kernel of [8], a link already exploited by [29] to derive their sequence neural network, which was a source of inspiration for our work. Both our kernels rely on a RNN to build a representation of an input sequence by computing a string kernel between this sequence and a set of learnable filters. Yet, our model exhibits several differences with [29], who use the regular substring kernel of [8] and compose this representation with another non-linear map—by applying an activation function to the output of the RNN. By contrast, we obtain a different RKHS directly by relaxing the substring kernel to allow for inexact matching at the compared positions, and embed the Nyström approximation within the RNN. The resulting feature space can be interpreted as a continuous neighborhood around all substrings (with gaps) of the described sequence. Furthermore, our RNN provides a finite-dimensional approximation of the relaxed kernel, relying on the Nyström approximation method [28]. As a consequence, RKNs may be learned in an unsupervised manner (in such a case, the goal is to approximate the kernel map), and with supervision with backpropagation, which may be interpreted as performing feature selection in the RKHS.

1.0.0.1 Contributions.

In this paper, we make the following contributions:
 \(\bullet\) We generalize convolutional kernel networks for sequences [27] to allow gaps, an important option for biological data. As in [27], we observe that the kernel formulation brings practical benefits over traditional CNNs or RNNs [20] when the amount of labeled data is small or moderate.
 \(\bullet\) We provide a kernel point of view on recurrent neural networks with new unsupervised and supervised learning algorithms. The resulting feature map can be interpreted in terms of gappy motifs, and end-to-end learning amounts to performing feature selection.
 \(\bullet\) Based on [30], we propose a new way to simulate max pooling in RKHSs, thus solving a classical discrepancy between theory and practice in the literature of string kernels, where sums are often replaced by a maximum operator that does not ensure positive definiteness [11].

2 Background on Kernel Methods and String Kernels↩︎

Kernel methods consist in mapping data points living in a set \({\mathcal{X}}\) to a possibly infinite-dimensional Hilbert space \({\mathcal{H}}\), through a mapping function \(\Phi: {\mathcal{X}}\to {\mathcal{H}}\), before learning a simple predictive model in \({\mathcal{H}}\) [31]. The so-called kernel trick allows to perform learning without explicitly computing this mapping, as long as the inner-product \(K({\mathbf{x}},{\mathbf{x}}')=\langle \Phi({\mathbf{x}}), \Phi({\mathbf{x}}')\rangle_{\mathcal{H}}\) between two points \({\mathbf{x}}, {\mathbf{x}}'\) can be efficiently computed. Whereas kernel methods traditionally lack scalability since they require computing an \(n\times n\) Gram matrix, where \(n\) is the amount of training data, recent approaches based on approximations have managed to make kernel methods work at large scale in many cases [28], [32].

For sequences in \({\mathcal{X}}={\mathcal{A}}^*\), which is the set of sequences of any possible length over an alphabet \({\mathcal{A}}\), the mapping \(\Phi\) often enumerates subsequence content. For instance, the spectrum kernel maps sequences to a fixed-length vector \(\Phi({\mathbf{x}}) = \left(\phi_u({\mathbf{x}})\right)_{u\in{\mathcal{A}}^k}\), where \({\mathcal{A}}^k\) is the set of \(k\)-mers—length-\(k\) sequence of characters in \({\mathcal{A}}\) for some \(k\) in \(\mathbb{N}\), and \(\phi_u({\mathbf{x}})\) counts the number of occurrences of \(u\) in \({\mathbf{x}}\) [6]. The mismatch kernel [7] operates similarly, but \(\phi_u({\mathbf{x}})\) counts the occurrences of \(u\) up to a few mismatched letters, which is useful when \(k\) is large and exact occurrences are rare.

2.1 Substring kernels↩︎

As [29], we consider the substring kernel introduced in [8], which allows to model the presence of gaps when trying to match a substring \(u\) to a sequence \({\mathbf{x}}\). Modeling gaps requires introducing the following notation: \({\mathcal{I}}_{{\mathbf{x}},k}\) denotes the set of indices of sequence \({\mathbf{x}}\) with \(k\) elements \((i_1,\dots,i_k)\) satisfying \(1\le i_1<\cdots< i_k\le |{\mathbf{x}}|\), where \(|{\mathbf{x}}|\) is the length of \({\mathbf{x}}\). For an index set \({\mathbf{i}}\) in \({\mathcal{I}}_{{\mathbf{x}},k}\), we may now consider the subsequence \({\mathbf{x}}_{{\mathbf{i}}} = ({\mathbf{x}}_{i_1},\dots,{\mathbf{x}}_{i_k})\) of \({\mathbf{x}}\) indexed by \({\mathbf{i}}\). Then, the substring kernel takes the same form as the mismatch and spectrum kernels, but \(\phi_u({\mathbf{x}})\) counts all—consecutive or not—subsequences of \({\mathbf{x}}\) equal to \(u\), and weights them by the number of gaps. Formally, we consider a parameter \(\lambda\) in \([0,1]\), and \(\phi_u({\mathbf{x}}) = \sum_{{\mathbf{i}}\in{\mathcal{I}}_{{\mathbf{x}},k}} \lambda^{{\text{gaps}}({\mathbf{i}})}\delta(u, {\mathbf{x}}_{{\mathbf{i}}})\), where \(\delta(u, v)=1\) if and only if \(u=v\), and \(0\) otherwise, and \({\text{gaps}}({\mathbf{i}}) :=i_k - i_1 - k + 1\) is the number of gaps in the index set \({\mathbf{i}}\). When \(\lambda\) is small, gaps are heavily penalized, whereas a value close to \(1\) gives similar weights to all occurrences. Ultimately, the resulting kernel between two sequences \({\mathbf{x}}\) and \({\mathbf{x}}'\) is \[\label{eq:substring95kernel} {\mathcal{K}}^s({\mathbf{x}},{\mathbf{x}}'):=\sum_{{\mathbf{i}}\in{\mathcal{I}}_{{\mathbf{x}},k}} \sum_{{\mathbf{j}}\in{\mathcal{I}}_{{\mathbf{x}}',k}} \lambda^{{\text{gaps}}({\mathbf{i}})}\lambda^{{\text{gaps}}({\mathbf{j}})} \delta\left({\mathbf{x}}_{{\mathbf{i}}}, {\mathbf{x}}'_{{\mathbf{j}}}\right).\tag{1}\] As we will see in Section 3, our RKN model relies on (1 ), but unlike [29], we replace the quantity \(\delta({\mathbf{x}}_{{\mathbf{i}}}, {\mathbf{x}}'_{{\mathbf{j}}})\) that matches exact occurrences by a relaxation, allowing more subtle comparisons. Then, we will show that the model can be interpreted as a gap-allowed extension of CKNs for sequences. We also note that even though \({\mathcal{K}}^s\) seems computationally expensive at first sight, it was shown in [8] that (1 ) admits a dynamic programming structure leading to efficient computations.

2.2 The Nyström method↩︎

When computing the Gram matrix is infeasible, it is typical to use kernel approximations  [28], [32], consisting in finding a \(q\)-dimensional mapping \(\psi: {\mathcal{X}}\to {\mathbb{R}}^q\) such that the kernel \(K({\mathbf{x}},{\mathbf{x}}')\) can be approximated by a Euclidean inner-product \(\langle \psi({\mathbf{x}}), \psi({\mathbf{x}}')\rangle_{{\mathbb{R}}^q}\). Then, kernel methods can be simulated by a linear model operating on \(\psi({\mathbf{x}})\), which does not raise scalability issues if \(q\) is reasonably small. Among kernel approximations, the Nyström method consists in projecting points of the RKHS onto a \(q\)-dimensional subspace, allowing to represent points into a \(q\)-dimensional coordinate system.

Specifically, consider a collection of \(Z=\{{\mathbf{z}}_1,\ldots,{\mathbf{z}}_q\}\) points in \({\mathcal{X}}\) and consider the subspace \[{\mathcal{E}}= \text{Span}(\Phi({\mathbf{z}}_1), \ldots, \Phi({\mathbf{z}}_q))~~~\text{and define}~~~\psi({\mathbf{x}}) = K_{ZZ}^{-\frac{1}{2}}K_Z({\mathbf{x}}),\] where \(K_{ZZ}\) is the \(q \times q\) Gram matrix of \(K\) restricted to the samples \({\mathbf{z}}_1, \ldots, {\mathbf{z}}_q\) and \(K_Z({\mathbf{x}})\) in \({\mathbb{R}}^q\) carries the kernel values \(K({\mathbf{x}}, {\mathbf{z}}_j), j=1, \ldots, q\). This approximation only requires \(q\) kernel evaluations and often retains good performance for learning. Interestingly as noted in [26], \(\langle \psi({\mathbf{x}}), \psi({\mathbf{x}}')\rangle_{{\mathbb{R}}^q}\) is exactly the inner-product in \({\mathcal{H}}\) between the projections of \(\Phi({\mathbf{x}})\) and \(\Phi({\mathbf{x}}')\) onto \({\mathcal{E}}\), which remain in \({\mathcal{H}}\).

When \({\mathcal{X}}\) is a Euclidean space—this can be the case for sequences when using a one-hot encoding representation, as discussed later— a good set of anchor points \({\mathbf{z}}_j\) can be obtained by simply clustering the data and choosing the centroids as anchor points [33]. The goal is then to obtain a subspace \({\mathcal{E}}\) that spans data as best as possible. Otherwise, previous works on kernel networks [26], [27] have also developed procedures to learn the set of anchor points end-to-end by optimizing over the learning objective. This approach can then be seen as performing feature selection in the RKHS.

3 Recurrent Kernel Networks↩︎

With the previous tools in hand, we now introduce RKNs. We show that it admits variants of CKNs, substring and local alignment kernels as special cases, and we discuss its relation with RNNs.

3.1 A continuous relaxation of the substring kernel allowing mismatches↩︎

From now on, and with an abuse of notation, we represent characters in \({\mathcal{A}}\) as vectors in \({\mathbb{R}}^d\). For instance, when using one-hot encoding, a DNA sequence \({\mathbf{x}}=({\mathbf{x}}_1,\ldots,{\mathbf{x}}_m)\) of length \(m\) can be seen as a 4-dimensional sequence where each \({\mathbf{x}}_j\) in \(\{0, 1\}^4\) has a unique non-zero entry indicating which of \(\{A,C,G,T\}\) is present at the \(j\)-th position, and we denote by \({\mathcal{X}}\) the set of such sequences. We now define the single-layer RKN as a generalized substring kernel 1 in which the indicator function \(\delta({\mathbf{x}}_{\mathbf{i}},{\mathbf{x}}'_{\mathbf{j}})\) is replaced by a kernel for \(k\)-mers: \[{\mathcal{K}}_k({\mathbf{x}},{\mathbf{x}}'):=\sum_{{\mathbf{i}}\in{\mathcal{I}}_{{\mathbf{x}},k}} \sum_{{\mathbf{j}}\in{\mathcal{I}}_{{\mathbf{x}}',k}} \lambda_{{\mathbf{x}},{\mathbf{i}}}\lambda_{{\mathbf{x}},{\mathbf{j}}} e^{-\frac{\alpha}{2}\|{\mathbf{x}}_{\mathbf{i}}- {\mathbf{x}}'_{\mathbf{j}}\|^2},\label{eq:kRKN}\tag{2}\] where we assume that the vectors representing characters have unit \(\ell_2\)-norm, such that \(e^{-\frac{\alpha}{2}\|{\mathbf{x}}_{\mathbf{i}}- {\mathbf{x}}'_{\mathbf{j}}\|^2}=e^{\alpha( \langle {\mathbf{x}}_{\mathbf{i}}, {\mathbf{x}}'_{\mathbf{j}}\rangle - k)}=\prod_{t=1}^k e^{\alpha\left(\langle {\mathbf{x}}_{i_t},{\mathbf{x}}'_{j_t}\rangle -1\right)}\) is a dot-product kernel, and \(\lambda_{{\mathbf{x}},{\mathbf{i}}}=\lambda^{{\text{gaps}}({\mathbf{i}})}\) if we follow (1 ).

For \(\lambda=0\) and using the convention \(0^0=1\), all the terms in these sums are zero except those for \(k\)-mers with no gap, and we recover the kernel of the CKN model of [27] with a convolutional structure—up to the normalization, which is done \(k\)-mer-wise in CKN instead of position-wise.

Compared to 1 , the relaxed version 2 accommodates inexact \(k\)-mer matching. This is important for protein sequences, where it is common to consider different similarities between amino acids in terms of substitution frequency along evolution [34]. This is also reflected in the underlying sequence representation in the RKHS illustrated in Figure 1: by considering \(\varphi(.)\) the kernel mapping and RKHS \({\mathcal{H}}\) such that \(K({\mathbf{x}}_{\mathbf{i}}, {\mathbf{x}}'_{\mathbf{j}}) = e^{-\frac{\alpha}{2}\|{\mathbf{x}}_{\mathbf{i}}- {\mathbf{x}}'_{\mathbf{j}}\|^2} = \langle \varphi({\mathbf{x}}_{\mathbf{i}}), \varphi({\mathbf{x}}'_{\mathbf{j}}) \rangle_{{\mathcal{H}}}\), we have \[{\mathcal{K}}_k({\mathbf{x}},{\mathbf{x}}')= \left\langle \sum_{{\mathbf{i}}\in{\mathcal{I}}_{{\mathbf{x}},k}} \lambda_{{\mathbf{x}},{\mathbf{i}}} \varphi({\mathbf{x}}_{\mathbf{i}}) , \sum_{{\mathbf{j}}\in{\mathcal{I}}_{{\mathbf{x}}',k}} \lambda_{{\mathbf{x}},{\mathbf{j}}} \varphi({\mathbf{x}}'_{\mathbf{j}})\right\rangle_{{\mathcal{H}}}. \label{eq:kRKNconv}\tag{3}\] A natural feature map for a sequence \({\mathbf{x}}\) is therefore \(\Phi_k({\mathbf{x}}) = \sum_{{\mathbf{i}}\in{\mathcal{I}}_{{\mathbf{x}},k}} \lambda_{{\mathbf{x}},{\mathbf{i}}} \varphi({\mathbf{x}}_{\mathbf{i}})\): using the RKN amounts to representing \({\mathbf{x}}\) by a mixture of continuous neighborhoods \(\varphi({\mathbf{x}}_{\mathbf{i}}): {\mathbf{z}}\mapsto e^{-\frac{\alpha}{2}\|{\mathbf{x}}_{\mathbf{i}}- {\mathbf{z}}\|^2}\) centered on all its \(k\)-subsequences \({\mathbf{x}}_{\mathbf{i}}\) , each weighted by the corresponding \(\lambda_{{\mathbf{x}},{\mathbf{i}}}\) (e.g., \(\lambda_{{\mathbf{x}},{\mathbf{i}}} = \lambda^{{\text{gaps}}({\mathbf{i}})}\)). As a particular case, a feature map of CKN [27] is the sum of the kernel mapping of all the \(k\)-mers without gap.

3.2 Extension to all \(k\)-mers and relation to the local alignment kernel↩︎

Dependency in the hyperparameter \(k\) can be removed by summing \({\mathcal{K}}_k\) over all possible values: \[{\mathcal{K}}_{\textrm{sum}}({\mathbf{x}},{\mathbf{x}}'):= \sum_{k=1}^\infty{\mathcal{K}}_k({\mathbf{x}},{\mathbf{x}}') = \sum_{k=1}^{\max(|{\mathbf{x}}|, |{\mathbf{x}}'|)}{\mathcal{K}}_k({\mathbf{x}},{\mathbf{x}}').\] Interestingly, we note that \({\mathcal{K}}_\text{sum}\) admits the local alignment kernel of [11] as a special case. More precisely, local alignments are defined via the tensor product set \({\mathcal{A}}_k({\mathbf{x}},{\mathbf{x}}'):={\mathcal{I}}_{{\mathbf{x}},k}\times{\mathcal{I}}_{{\mathbf{x}}',k}\), which contains all possible alignments of \(k\) positions between a pair of sequences \(({\mathbf{x}}, {\mathbf{x}}')\). The local alignment score of each such alignment \(\pi=({\mathbf{i}},{\mathbf{j}})\) in \({\mathcal{A}}_k({\mathbf{x}},{\mathbf{x}}')\) is defined, by [11], as \(S({\mathbf{x}},{\mathbf{x}}', \pi):=\sum_{t=1}^{k} s({\mathbf{x}}_{i_t}, {\mathbf{x}}'_{j_t})-\sum_{t=1}^{k-1}[g(i_{t+1} - i_t - 1)+g(j_{t+1} - j_t - 1)]\), where \(s\) is a symmetric substitution function and \(g\) is a gap penalty function. The local alignment kernel in [11] can then be expressed in terms of the above local alignment scores (Thrm. 1.7 in [11]): \[\label{eq:LA} K_{LA}({\mathbf{x}},{\mathbf{x}}')=\sum_{k=1}^{\infty}K_{LA}^{k}({\mathbf{x}},{\mathbf{x}}'):=\sum_{k=1}^{\infty}\sum_{\pi\in{\mathcal{A}}_k({\mathbf{x}},{\mathbf{x}}')}\exp(\beta {\mathcal{S}}({\mathbf{x}},{\mathbf{x}}',\pi))~~~\text{for some}~\beta > 0.\tag{4}\] When the gap penalty function is linear—that is, \(g(x)=cx\) with \(c>0\), \(K_{LA}^{k}\) becomes \(K_{LA}^{k}({\mathbf{x}},{\mathbf{x}}')=\sum_{\pi\in{\mathcal{A}}_k({\mathbf{x}},{\mathbf{x}}')}\exp(\beta {\mathcal{S}}({\mathbf{x}},{\mathbf{x}}',\pi)) =\sum_{({\mathbf{i}},{\mathbf{j}})\in{\mathcal{A}}_k({\mathbf{x}},{\mathbf{x}}')} e^{-c\beta {\text{gaps}}({\mathbf{i}})} e^{-c\beta {\text{gaps}}({\mathbf{j}})} \prod_{t=1}^k e^{\beta s ({\mathbf{x}}_{i_t},{\mathbf{x}}'_{j_t})}\). When \(s({\mathbf{x}}_{i_t},{\mathbf{x}}'_{j_t})\) can be written as an inner-product \(\langle \psi_s({\mathbf{x}}_{i_t}),\psi_s({\mathbf{x}}'_{j_t})\rangle\) between normalized vectors, we see that \(K_{LA}\) becomes a special case of 2 —up to a constant factor—with \(\lambda_{{\mathbf{x}},{\mathbf{i}}}=e^{-c\beta{\text{gaps}}({\mathbf{i}})}\), \(\alpha={\beta }\).

This observation sheds new lights on the relation between the substring and local alignment kernels, which will inspire new algorithms in the sequel. To the best of our knowledge, the link we will provide between RNNs and local alignment kernels is also new.

3.3 Nyström approximation and recurrent neural networks↩︎

As in CKNs, we now use the Nyström approximation method as a building block to make the above kernels tractable. According to (3 ), we may first use the Nyström method described in Section 2.2 to find an approximate embedding for the quantities \(\varphi({\mathbf{x}}_{\mathbf{i}})\), where \({\mathbf{x}}_{\mathbf{i}}\) is one of the \(k\)-mers represented as a matrix in \({\mathbb{R}}^{k \times d}\). This is achieved by choosing a set \(Z=\{{\mathbf{z}}_1,\ldots,{\mathbf{z}}_q\}\) of anchor points in \({\mathbb{R}}^{k \times d}\), and by encoding \(\varphi({\mathbf{x}}_{\mathbf{i}})\) as \(K_{ZZ}^{-{1}/{2}}K_Z({\mathbf{x}}_{\mathbf{i}})\)—where \(K\) is the kernel of \({\mathcal{H}}\). Such an approximation for \(k\)-mers yields the \(q\)-dimensional embedding for the sequence \({\mathbf{x}}\): \[\psi_k({\mathbf{x}}) = \sum_{{\mathbf{i}}\in{\mathcal{I}}_{{\mathbf{x}},k}} \lambda_{{\mathbf{x}},{\mathbf{i}}} K_{ZZ}^{-\frac{1}{2}} K_Z({\mathbf{x}}_{\mathbf{i}}) = K_{ZZ}^{-\frac{1}{2}} \sum_{{\mathbf{i}}\in{\mathcal{I}}_{{\mathbf{x}},k}} \lambda_{{\mathbf{x}},{\mathbf{i}}}K_Z({\mathbf{x}}_{\mathbf{i}}).\label{eq:psik}\tag{5}\] Then, an approximate feature map \(\psi_{\text{sum}}({\mathbf{x}})\) for the kernel \({\mathcal{K}}_\text{sum}\) can be obtained by concatenating the embeddings \(\psi_1({\mathbf{x}}),\ldots,\psi_k({\mathbf{x}})\) for \(k\) large enough.

3.3.0.1 The anchor points as motifs.

The continuous relaxation of the substring kernel presented in (2 ) allows us to learn anchor points that can be interpreted as sequence motifs, where each position can encode a mixture of letters. This can lead to more relevant representations than \(k\)-mers for learning on biological sequences. For example, the fact that a DNA sequence is bound by a particular transcription factor can be associated with the presence of a T followed by either a G or an A, followed by another T, would require two \(k\)-mers but a single motif [27]. Our kernel is able to perform such a comparison.

3.3.0.2 Efficient computations of \({\mathcal{K}}_k\) and \({\mathcal{K}}_\text{sum}\) approximation via RNNs.

A naive computation of \(\psi_k({\mathbf{x}})\) would require enumerating all substrings present in the sequence, which may be exponentially large when allowing gaps. For this reason, we use the classical dynamic programming approach of substring kernels [8], [29]. Consider then the computation of \(\psi_j({\mathbf{x}})\) defined in (5 ) for \(j=1,\ldots,k\) as well as a set of anchor points \(Z_k=\{{\mathbf{z}}_1,\ldots,{\mathbf{z}}_q\}\) with the \({\mathbf{z}}_i\)’s in \({\mathbb{R}}^{d \times k}\). We also denote by \(Z_j\) the set obtained when keeping only \(j\)-th first positions (columns) of the \({\mathbf{z}}_j\)’s, leading to \(Z_j=\{[{\mathbf{z}}_1]_{1:j},\ldots,[{\mathbf{z}}_q]_{1:j}\}\), which will serve as anchor points for the kernel \({\mathcal{K}}_j\) to compute \(\psi_j({\mathbf{x}})\). Finally, we denote by \({\mathbf{z}}_i^j\) in \({\mathbb{R}}^d\) the \(j\)-th column of \({\mathbf{z}}_i\) such that \({\mathbf{z}}_i=[{\mathbf{z}}_i^1,\ldots,{\mathbf{z}}_i^k]\). Then, the embeddings \(\psi_1({\mathbf{x}}),\ldots,\psi_k({\mathbf{x}})\) can be computed recursively by using the following theorem:

A proof is provided in Appendix 5 and is based on classical recursions for computing the substring kernel, which were interpreted as RNNs by [29]. The main difference in the RNN structure we obtain is that their non-linearity is applied over the outcome of the network, leading to a feature map formed by composing the feature map of the substring kernel of [8] and another one from a RKHS that contains their non-linearity. By contrast, our non-linearities are built explicitly in the substring kernel, by relaxing the indicator function used to compare characters. The resulting feature map is a continuous neighborhood around all substrings of the described sequence. In addition, the Nyström method yields an orthogonalization factor \(K_{ZZ}^{-{1}/{2}}\) to the output \(K_Z({\mathbf{x}})\) of the network to compute our approximation, which is perhaps the only non-standard component of our RNN. This factor provides an interpretation of \(\psi({\mathbf{x}})\) as a kernel approximation. As discussed next, it makes it possible to learn the anchor points by \(k\)-means, see [27], which also makes the initialization of the supervised learning procedure simple without having to deal with the scaling of the initial motifs/filters \({\mathbf{z}}_j\).

3.3.0.3 Learning the anchor points \(Z\).

We now turn to the application of RKNs to supervised learning. Given \(n\) sequences \({\mathbf{x}}^1,\dots,{\mathbf{x}}^n\) in \({\mathcal{X}}\) and their associated labels \(y^1,\dots,y^n\) in \({\mathcal{Y}}\), e.g., \({\mathcal{Y}}=\{-1,1\}\) for binary classification or \({\mathcal{Y}}={\mathbb{R}}\) for regression, our objective is to learn a function in the RKHS \({\mathcal{H}}\) of \({\mathcal{K}}_k\) by minimizing \[\min_{f\in{\mathcal{H}}} \frac{1}{n}\sum_{i=1}^n L(f({\mathbf{x}}^i),y^i)+\frac{\mu}{2} \|f\|^2_{{\mathcal{H}}},\] where \(L:{\mathbb{R}}\times{\mathbb{R}}\to{\mathbb{R}}\) is a convex loss function that measures the fitness of a prediction \(f({\mathbf{x}}^i)\) to the true label \(y^i\) and \(\mu\) controls the smoothness of the predictive function. After injecting our kernel approximation \({\mathcal{K}}_k({\mathbf{x}},{\mathbf{x}}')\simeq\langle \psi_k({\mathbf{x}}), \psi_k({\mathbf{x}}') \rangle_{{\mathbb{R}}^q}\), the problem becomes \[\label{eq:sup} \min_{{\mathbf{w}}\in{\mathbb{R}}^q} \frac{1}{n}\sum_{i=1}^n L\left(\langle \psi_k({\mathbf{x}}^i),{\mathbf{w}}\rangle ,y^i\right) + \frac{\mu}{2}\|{\mathbf{w}}\|^2.\tag{6}\]

Following [26], [27], we can learn the anchor points \(Z\) without exploiting training labels, by applying a \(k\)-means algorithm to all (or a subset of) the \(k\)-mers extracted from the database and using the obtained centroids as anchor points. Importantly, once \(Z\) has been obtained, the linear function parametrized by \({\mathbf{w}}\) is still optimized with respect to the supervised objective 6 . This procedure can be thought of as learning a general representation of the sequences disregarding the supervised task, which can lead to a relevant description while limiting overfitting.

Another strategy consists in optimizing (6 ) jointly over \((Z,{\mathbf{w}})\), after observing that \(\psi_k({\mathbf{x}}) = K_{ZZ}^{-{1}/{2}} \sum_{{\mathbf{i}}\in{\mathcal{I}}_{{\mathbf{x}},k}} \lambda_{{\mathbf{x}},{\mathbf{i}}}K_Z({\mathbf{x}}_{\mathbf{i}})\) is a smooth function of \(Z\). Learning can be achieved by using backpropagation over \((Z,{\mathbf{w}})\), or by using an alternating minimization strategy between \(Z\) and \({\mathbf{w}}\). It leads to an end-to-end scheme where both the representation and the function defined over this representation are learned with respect to the supervised objective 6 . Backpropagation rules for most operations are classical, except for the matrix inverse square root function, which is detailed in Appendix 6. Initialization is also parameter-free since the unsupervised learning approach may be used for that.

3.4 Extensions↩︎

3.4.0.1 Multilayer construction.

In order to account for long-range dependencies, it is possible to construct a multilayer model based on kernel compositions similar to [29]. Assume that \({\mathcal{K}}^{(n)}_k\) is the \(n\)-th layer kernel and \(\Phi_k^{(n)}\) its mapping function. The corresponding \((n+1)\)-th layer kernel is defined as \[{\mathcal{K}}_k^{(n+1)} ({\mathbf{x}},{\mathbf{x}}')=\sum_{{\mathbf{i}}\in{\mathcal{I}}_{{\mathbf{x}},k},{\mathbf{j}}\in{\mathcal{I}}_{{\mathbf{x}}',k}}\lambda^{(n+1)}_{{\mathbf{x}},{\mathbf{i}}}\lambda^{(n+1)}_{{\mathbf{x}}',{\mathbf{j}}} \prod_{t=1}^k K_{n+1}(\Phi_k^{(n)}({\mathbf{x}}_{1:i_t}),\Phi_k^{(n)}({\mathbf{x}}'_{1:j_t})), \label{eq:multilayer}\tag{7}\] where \(K_{n+1}\) will be defined in the sequel and the choice of weights \(\lambda^{(n)}_{{\mathbf{x}},{\mathbf{i}}}\) slightly differs from the single-layer model. We choose indeed \(\lambda^{(N)}_{{\mathbf{x}},{\mathbf{i}}} = \lambda^{{\text{gaps}}({\mathbf{i}})}\) only for the last layer \(N\) of the kernel, which depends on the number of gaps in the index set \({\mathbf{i}}\) but not on the index positions. Since (7 ) involves a kernel \(K_{n+1}\) operating on the representation of prefix sequences \(\Phi_k^{(n)}({\mathbf{x}}_{1:t})\) from layer \(n\), the representation makes sense only if \(\Phi_k^{(n)}({\mathbf{x}}_{1:t})\) carries mostly local information close to position \(t\). Otherwise, information from the beginning of the sequence would be overrepresented. Ideally, we would like the range-dependency of \(\Phi_k^{(n)}({\mathbf{x}}_{1:t})\) (the size of the window of indices before \(t\) that influences the representation, akin to receptive fields in CNNs) to grow with the number of layers in a controllable manner. This can be achieved by choosing \(\lambda^{(n)}_{{\mathbf{x}},{\mathbf{i}}}=\lambda^{|{\mathbf{x}}|-i_1-k+1}\) for \(n < N\), which assigns exponentially more weights to the \(k\)-mers close to the end of the sequence.

For the first layer, we recover the single-layer network \({\mathcal{K}}_k\) defined in 2 by defining \(\Phi_k^{(0)}({\mathbf{x}}_{1:i_k}) \!=\! {\mathbf{x}}_{i_k}\) and \(K_1({\mathbf{x}}_{i_k},{\mathbf{x}}'_{j_k}) = e^{\alpha(\langle {\mathbf{x}}_{i_k},{\mathbf{x}}'_{j_k}\rangle -1)}\). For \(n>1\), it remains to define \(K_{n+1}\) to be a homogeneous dot-product kernel, as used for instance in CKNs [26]: \[K_{n+1}({\mathbf{u}}, {\mathbf{u}}') = \| {\mathbf{u}}\|_{{\mathcal{H}}_n} \| {\mathbf{u}}\|_{{\mathcal{H}}_n} \kappa_n\left (\left\langle \frac{{\mathbf{u}}}{\|{\mathbf{u}}\|_{{\mathcal{H}}_n}}, \frac{{\mathbf{u}}'}{\|{\mathbf{u}}'\|_{{\mathcal{H}}_n}} \right\rangle_{{\mathcal{H}}_n}\right)~~~~\text{with}~~~\kappa_n(t) = e^{\alpha_n (t-1)}.\label{eq:dotprod}\tag{8}\] Note that the Gaussian kernel \(K_1\) used for 1st layer may also be written as (8 ) since characters are normalized. As for CKNs, the goal of homogenization is to prevent norms to grow/vanish exponentially fast with \(n\), while dot-product kernels lend themselves well to neural network interpretations.

As detailed in Appendix 7, extending the Nyström approximation scheme for the multilayer construction may be achieved in the same manner as with CKNs—that is, we learn one approximate embedding \(\psi_k^{(n)}\) at each layer, allowing to replace the inner-products \(\langle \Phi_k^{(n)}({\mathbf{x}}_{1:i_t}),\Phi_k^{(n)}({\mathbf{x}}'_{1:j_t})\rangle\) by their approximations \(\langle \psi_k^{(n)}({\mathbf{x}}_{1:i_t}),\psi_k^{(n)}({\mathbf{x}}'_{1:j_t})\rangle\), and it is easy to show that the interpretation in terms of RNNs is still valid since \({\mathcal{K}}^{(n)}_k\) has the same sum structure as (2 ).

3.4.0.2 Max pooling in RKHS.

Alignment scores (e.g. Smith-Waterman) in molecular biology rely on a max operation—over the scores of all possible alignments—to compute similarities between sequences. However, using max in a string kernel usually breaks positive definiteness, even though it seems to perform well in practice. To solve such an issue, sum-exponential is used as a proxy in [35], but it leads to diagonal dominance issue and makes SVM solvers unable to learn. For RKN, the sum in 3 can also be replaced by a max \[{\mathcal{K}}^{\text{max}}_k({\mathbf{x}},{\mathbf{x}}')= \left\langle \max_{{\mathbf{i}}\in{\mathcal{I}}_{{\mathbf{x}},k}} \lambda_{{\mathbf{x}},{\mathbf{i}}} \psi_k({\mathbf{x}}_{\mathbf{i}}) , \max_{{\mathbf{j}}\in{\mathcal{I}}_{{\mathbf{x}}',k}} \lambda_{{\mathbf{x}},{\mathbf{j}}} \psi_k({\mathbf{x}}'_{\mathbf{j}})\right\rangle, \label{eq:kRKNmax}\tag{9}\] which empirically seems to perform well, but breaks the kernel interpretation, as in [35]. The corresponding recursion amounts to replacing all the sum in ?? by a max.

An alternative way to aggregate local features is the generalized max pooling (GMP) introduced in [30], which can be adapted to the context of RKHSs. Assuming that before pooling \({\mathbf{x}}\) is embedded to a set of \(N\) local features \((\varphi_1, \dots, \varphi_N)\in{\mathcal{H}}^N\), GMP builds a representation \(\varphi^{{\text{gmp}}}\) whose inner-product with all the local features \(\varphi_i\) is one: \(\langle \varphi_i, \varphi^{{\text{gmp}}} \rangle_{{\mathcal{H}}}=1,~\text{for}~i=1,\dots,N\). \(\varphi^{{\text{gmp}}}\) coincides with the regular max when each \(\varphi\) is an element of the canonical basis of a finite representation—i.e., assuming that at each position, a single feature has value 1 and all others are \(0\).

Since GMP is defined by a set of inner-products constraints, it can be applied to our approximate kernel embeddings by solving a linear system. This is compatible with CKN but becomes intractable for RKN which pools across \(|{\mathcal{I}}_{{\mathbf{x}},k}|\) positions. Instead, we heuristically apply GMP over the set \(\psi_k({\mathbf{x}}_{1:t})\) for all \(t\) with \(\lambda_{{\mathbf{x}},{\mathbf{i}}}=\lambda^{|{\mathbf{x}}|-i_1-k+1}\), which can be obtained from the RNN described in Theorem 1. This amounts to composing GMP with mean poolings obtained over each prefix of \({\mathbf{x}}\). We observe that it performs well in our experiments. More details are provided in Appendix 8.

4 Experiments↩︎

We evaluate RKN and compare it to typical string kernels and RNN for protein fold recognition. Pytorch code is provided with the submission and additional details given in Appendix 9.

4.1 Protein fold recognition on SCOP 1.67↩︎

Sequencing technologies provide access to gene and, indirectly, protein sequences for yet poorly studied species. In order to predict the 3D structure and function from the linear sequence of these proteins, it is common to search for evolutionary related ones, a problem known as homology detection. When no evolutionary related protein with known structure is available, a—more difficult—alternative is to resort to protein fold recognition. We evaluate our RKN on such a task, where the objective is to predict which proteins share a 3D structure with the query [36].

Here we consider the Structural Classification Of Proteins (SCOP) version 1.67 [37]. We follow the preprocessing procedures of [38] and remove the sequences that are more than 95% similar, yielding 85 fold recognition tasks. Each positive training set is then extended with Uniref50 to make the dataset more balanced, as proposed in [20]. The resulting dataset can be downloaded from http://www.bioinf.jku.at/software/LSTM_protein. The number of training samples for each task is typically around 9,000 proteins, whose length varies from tens to thousands of amino-acids. In all our experiments we use logistic loss. We measure classification performances using auROC and auROC50 scores (area under the ROC curve and up to 50% false positives).

For CKN and RKN, we evaluate both one-hot encoding of amino-acids by 20-dimensional binary vectors and an alternative representation relying on the BLOSUM62 substitution matrix [34]. Specifically in the latter case, we represent each amino-acid by the centered and normalized vector of its corresponding substitution probabilities with other amino-acids. The local alignment kernel 4 , which we include in our comparison, natively uses BLOSUM62.

4.1.0.1 Hyperparameters.

We follow the training procedure of CKN presented in [27]. Specifically, for each of the \(85\) tasks, we hold out one quarter of the training samples as a validation set, use it to tune \(\alpha\), gap penalty \(\lambda\) and the regularization parameter \(\mu\) in the prediction layer. These parameters are then fixed across datasets. RKN training also relies on the alternating strategy used for CKN: we use an Adam algorithm to update anchor points, and the L-BFGS algorithm to optimize the prediction layer. We train 100 epochs for each dataset: the initial learning rate for Adam is fixed to 0.05 and is halved as long as there is no decrease of the validation loss for 5 successive epochs. We fix \(k\) to 10, the number of anchor points \(q\) to 128 and use single layer CKN and RKN throughout the experiments.

4.1.0.2 Implementation details for unsupervised models.

The anchor points for CKN and RKN are learned by k-means on 30,000 extracted \(k\)-mers from each dataset. The resulting sequence representations are standardized by removing mean and dividing by standard deviation and are used within a logistic regression classifier. \(\alpha\) in Gaussian kernel and the parameter \(\lambda\) are chosen based on validation loss and are fixed across the datasets. \(\mu\) for regularization is chosen by a 5-fold cross validation on each dataset. As before, we fix \(k\) to 10 and the number of anchor points \(q\) to 1024. Note that the performance could be improved with larger \(q\) as observed in [27], at a higher computational cost.

4.1.0.3 Comparisons and results.

The results are shown in Table 1. The blosum62 version of CKN and RKN outperform all other methods. Improvement against the mismatch and LA kernels is likely caused by end-to-end trained kernel networks learning a task-specific representation in the form of a sparse set of motifs, whereas data-independent kernels lead to learning a dense function over the set of descriptors. This difference can have a regularizing effect akin to the \(\ell_1\)-norm in the parametric world, by reducing the dimension of the learned linear function \(w\) while retaining relevant features for the prediction task. GPkernel also learns motifs, but relies on the exact presence of discrete motifs. Finally, both LSTM and [29] are based on RNNs but are outperformed by kernel networks. The latter was designed and optimized for NLP tasks and yields a \(0.4\) auROC50 on this task.

RKNs outperform CKNs, albeit not by a large margin. Interestingly, as the two kernels only differ by their allowing gaps when comparing sequences, this results suggests that this aspect is not the most important for identifying common foldings in a one versus all setting: as the learned function discriminates on fold from all others, it may rely on coarser features and not exploit more subtle ones such as gappy motifs. In particular, the advantage of the LA-kernel against its mismatch counterpart is more likely caused by other differences than gap modelling, namely using a max rather than a mean pooling of \(k\)-mer similarities across the sequence, and a general substitution matrix rather than a Dirac function to quantify mismatches. Consistently, within kernel networks GMP systematically outperforms mean pooling, while being slightly behind max pooling.

Additional details and results, scatter plots, and pairwise tests between methods to assess the statistical significance of our conclusions are provided in Appendix 9. Note that when \(k=14\), the auROC and auROC50 further increase to 0.877 and 0.636 respectively.

4.2 Protein fold classification on SCOP 2.06↩︎

We further benchmark RKN in a fold classification task, following the protocols used in [40]. Specifically, the training and validation datasets are composed of 14699 and 2013 sequences from SCOP 1.75, belonging to 1195 different folds. The test set consists of 2533 sequences from SCOP 2.06, after removing the sequences with similarity greater than 40% with SCOP 1.75. The input sequence feature is represented by a vector of 45 dimensions, consisting of a 20-dimensional one-hot encoding of the sequence, a 20-dimensional position-specific scoring matrix (PSSM) representing the profile of amino acids, a 3-class secondary structure represented by a one-hot vector and a 2-class solvent accessibility. We further normalize each type of the feature vectors to have unit \(\ell_2\)-norm, which is done for each sequence position. More dataset details can be found in [40]. We use mean pooling for both CKN and RKN models, as it is more stable during training for multi-class classification. The other hyperparameters are chosen in the same way as previously. More details about hyperparameter search grid can be found in Appendix E.

The accuracy results are obtained by averaging 10 different runs and are shown in Table 2, stratified by prediction difficulty (family/superfamily/fold, more details can be found in [40]). By contrast to what we observed on SCOP 1.67, RKN sometimes yields a large improvement on CKN for fold classification, especially for detecting distant homologies. This suggests that accounting for gaps does help in some fold prediction tasks, at least in a multi-class context where a single function is learned for each fold.

Acknowledgments↩︎

We thank the anonymous reviewers for their insightful comments and suggestions. This work has been supported by the grants from ANR (FAST-BIG project ANR-17-CE23-0011-01), by the ERC grant number 714381 (SOLARIS), and ANR 3IA MIAI@Grenoble Alpes.

5 Nyström Approximation for Single-Layer RKN↩︎

We detail here the Nytröm approximation presented in Section 3.3, which we recall here for a sequence \({\mathbf{x}}\): \[\psi_k({\mathbf{x}}) = K_{ZZ}^{-1/2} \sum_{{\mathbf{i}}\in{\mathcal{I}}_{{\mathbf{x}},k}} \lambda_{{\mathbf{x}},{\mathbf{i}}}K_Z({\mathbf{x}}_{\mathbf{i}}).\label{eq:psik95app}\tag{10}\] Consider then the computation of \(\psi_j({\mathbf{x}})\) defined in (10 ) for \(j=1,\ldots,k\) given a set of anchor points \(Z_k=\{{\mathbf{z}}_1,\ldots,{\mathbf{z}}_q\}\) with the \({\mathbf{z}}_i\)’s in \({\mathbb{R}}^{d \times k}\). Given the notations introduced in Section 3.3, we are now in shape to prove Theorem 1.

6 Back-propagation for Matrix Inverse Square Root↩︎

In Section 3.3, we have described an end-to-end scheme to jointly optimize \(Z\) and \({\mathbf{w}}\). The back-propagation of \(Z\) requires computing that of the matrix inverse square root operation as it is involved in the approximate feature map of \({\mathbf{x}}\) as shown in 10 . The back-propagation formula is given by the following proposition, which is based on an errata of [26] and we include it here for completeness.

When doing back-propagation, one is usually interested in computing a quantity \(\bar{\mathbf{A}}\) such that given \(\bar{\mathbf{B}}\) (with appropriate dimensions), we have \[\langle \bar{\mathbf{B}}, d(\mathbf{A}^{-{\frac{1}{2}}}) \rangle_F = \langle \bar{\mathbf{A}}, d\mathbf{A}\rangle_F,\] see [41], for instance. Here, \(\langle, \rangle_F\) denotes the Frobenius inner product. Then, it is easy to show that \[\bar{\mathbf{A}} = - \mathbf{U}(\mathbf{F}\circ ( \mathbf{U}^\top \bar{\mathbf{B}} \mathbf{U})) \mathbf{U}^\top.\]

7 Multilayer Construction of RKN↩︎

For multilayer RKN, assume that we have defined \({\mathcal{K}}^{(n)}\) the \(n\)-th layer kernel. To simplify the notation below, we consider that an input sequence \({\mathbf{x}}\) is encoded at layer \(n\) as \({\mathbf{x}}^{(n)}:=(\Phi_k^{(n)}({\mathbf{x}}_1), \Phi_k^{(n)}({\mathbf{x}}_{1:2}),\dots, \Phi_k^{(n)}({\mathbf{x}}) )\) where the feature map at position \(t\) is \({\mathbf{x}}^{(n)}_t=\Phi_k^{(n)}({\mathbf{x}}_{1:t})\). The \((n+1)\)-layer kernel is defined by induction by \[\label{eq:rkn95multi95app} {\mathcal{K}}_k^{(n+1)} ({\mathbf{x}},{\mathbf{x}}')=\sum_{{\mathbf{i}}\in{\mathcal{I}}_{{\mathbf{x}},k},{\mathbf{j}}\in{\mathcal{I}}_{{\mathbf{x}}',k}}\lambda^{(n)}_{{\mathbf{x}},{\mathbf{i}}}\lambda^{(n)}_{{\mathbf{x}}',{\mathbf{j}}} \prod_{t=1}^k K_{n+1}({\mathbf{x}}^{(n)}_{i_t},{\mathbf{x}}'^{(n)}_{j_t}),\tag{11}\] where \(K_{n+1}\) is defined in (8 . With the choice of weights described in Section 3.4, the construction scheme for an \(n\)-layer RKN is illustrated in Figure 2

The Nyström approximation scheme for multilayer RKN is straightforward by inductively applying the Nytröm method to the kernels \({\mathcal{K}}^{(1)},\dots,{\mathcal{K}}^{(n)}\) from bottom to top layers. Specifically, assume that \({\mathcal{K}}^{(n)}({\mathbf{x}},{\mathbf{x}}')\) is approximated by \(\langle \psi^{(n)}_k({\mathbf{x}}), \psi^{(n)}_k({\mathbf{x}}')\rangle_{{\mathbb{R}}^{q_n}}\) such that the approximate feature map of \({\mathbf{x}}^{(n)}\) at position \(t\) is \(\psi^{(n)}_k({\mathbf{x}}_{1:t})\). Now Consider a set of anchor points \(Z_k=\{{\mathbf{z}}_1,\ldots,{\mathbf{z}}_{q_{n+1}}\}\) with the \({\mathbf{z}}_i\)’s in \({\mathbb{R}}^{q_n \times k}\) which have unit norm at each column. We use the same notations as in single-layer construction. Then very similar to the single-layer RKN, the embeddings \((\psi^{(n+1)}_j({\mathbf{x}}_{1:t}^{(n)}))_{1=1,\dots,k, t=1,\dots,|{\mathbf{x}}^{(n)}|}\) are given by the following recursion

8 Generalized Max Pooling for RKN↩︎

Assume that a sequence \({\mathbf{x}}\) is embedded to \((\varphi_1, \dots,\varphi_n)\in{\mathcal{H}}^n\) local features, as in Section 3.4. Generalized max pooling (GMP) looks for a representation \(\varphi^{{\text{gmp}}}\) such that the inner product between this vector and all the local representations is one: \(\langle \varphi_i, \varphi^{{\text{gmp}}} \rangle_{{\mathcal{H}}}=1,~\text{for }i=1,\dots,n.\) Assuming that each \(\varphi_i\) is now represented by a vector \(\psi_i\) in \({\mathbb{R}}^q\), the above problem can be approximately solved by search an embedding vector \(\psi^{{\text{gmp}}}\) in \({\mathbb{R}}^q\) such that \(\langle \psi_i, \psi^{{\text{gmp}}} \rangle =1\) for \(i=1,\dots,n\). In practice, and to prevent ill-conditioned problems, as shown in [30], it is possible to solve a ridge regression problem: \[\psi^{{\text{gmp}}}=\mathop{\mathrm{arg\,min}}_{\psi\in{\mathbb{R}}^q} \| \Psi ^{\top} \psi - \mathbf{1} \|^2 + \gamma \| \psi\|^2,\] where \(\Psi=[\psi_1,\dots,\psi_n]\in{\mathbb{R}}^{q\times n}\) and \(\mathbf{1}\) denotes the \(n\)-dimensional vectors with only 1 as entries. The solution is simply given by \(\psi^{{\text{gmp}}}=(\Psi \Psi^{\top}+\gamma I )^{-1}\Psi \mathbf{1}\). It requires inverting a \(q\times q\) matrix which is usually tractable when the number of anchor points is small. In particular, when \(\psi_i=K_{ZZ}^{-1/2} K_Z({\mathbf{x}}_i)\) the Nyström approximation of a local feature map, we have \(\Psi=K_{ZZ}^{-1/2} K_{ZX}\) with \([K_{ZX}]_{ji}=K({\mathbf{z}}_j,{\mathbf{x}}_i)\) and thus \[\psi^{{\text{gmp}}}=K_{ZZ}^{\frac{1}{2}} (K_{ZX}K_{ZX}^{\top}+\gamma K_{ZZ} )^{-1} K_{ZX}\mathbf{1}.\]

9 Additional Experimental Material↩︎

In this section, we provide additional details about experiments and scatter plots with pairwise statistical tests.

9.1 Protein fold recognition on SCOP 1.67↩︎

9.1.0.1 Hyperparameter search grids.

Here, we first provide the grids used for hyperparameter search. In our experiments, we use \(\sigma\) instead of \(\alpha\) such that \(\alpha=1/k \sigma^2\). The search range is specified in Table 3.

9.1.0.2 Comparison of unsupervised CKNs and RKNs.

Then, we provide an additional table of results to compare the unsupervised models of CKN and RKN. In this unsupervised regime, mean pooling perform better than max pooling, which is different than what we have observed in the supervised case. RKN tend to work better than CKN, while RKN-sum—that is, using the kernel \({\mathcal{K}}_{\text{sum}}\) instead of \({\mathcal{K}}_k\), works better than RKN.

9.1.0.3 Study of filter number \(q\) and size \(k\).

Here we use max pooling and fix \(\sigma\) to 0.4 and \(\lambda\) to 0.1. When \(q\) varies \(k\) is fixed to 10 and \(q\) is fixed to 128 when \(k\) varies. We show here the performance of RKN with different choices of \(q\) and \(k\). The gap hyperparameter \(\lambda\) is chosen optimally for each \(q\) and \(k\). The results are shown in Figure 3.

9.1.0.4 Discussion about complexity.

Performing backpropgation with our RKN model has the same complexity has a performing a similar step within a recurrent neural network, up to the computation of the inverse square root matrix \(K_{ZZ}^{-1/2}\), which has complexity \(O(q^3)\). When \(q\) is reasonably small \(q=128\) in our experiments, such a complexity is negligible. For instance, one forward pass with a minibatch of \(b=128\) sequences of length \(m\) yields a complexity \(O(k^2 m b q)\), which can typically be much greater than \(q^3\).

9.1.0.5 Computing infrastructures.

Experiments were conduced by using a shared GPU cluster, in large parts build with Nvidia gamer cards (Titan X, GTX1080TI). About 10 of these GPUs were used simultaneously to perform the experiments of this paper.

9.1.0.6 Scatter plots and statistical testing.

Even though each method was run only one time for each task, the 85 tasks allow us to conduct statistical testing when comparing two methods. In Figures 4 and 5, we provide pairwise comparisons allowing us to assess the statistical significance of various conclusions drawn in the paper. We use a Wilcoxon signed-rank test to provide p-values.

9.2 Protein fold classification on SCOP 2.06↩︎

9.2.0.1 Hyperparameter search grids.

We provide the grids used for hyperparameter search, shown in Table 5.

9.2.0.2 Complete results with error bars.

The classification accuracy for CKNs and RKNs on protein fold classification on SCOP 2.06 are obtained by averaging on 10 runs with different seeds. The results are shown in Table ¿tbl:fig:scop20695errbar? with error bars.

References↩︎