Tanimoto Random Features for Scalable Molecular Machine Learning
Authors: Austin Tripp
Sergio Bacallado
Sukriti Singh
José Miguel Hernández-Lobato
Venue: NeurIPS 2023 poster (2023)
Keywords: Tanimoto
Kernel
MinMax
Gaussian process
molecule
chemistry
random features
TL;DR: We propose two kinds of random features for the Tanimoto/Jaccard kernel and apply it to molecule data.

Abstract: The Tanimoto coefficient is commonly used to measure the similarity between molecules represented as discrete fingerprints,
either as a distance metric or a positive definite kernel. While many kernel methods can be accelerated using random feature approximations, at present there is a lack of such approximations for the Tanimoto kernel. In this paper we propose two kinds of novel random features to allow this kernel to scale to large datasets, and in the process discover a novel extension of the kernel to real-valued vectors. We theoretically characterize these random features, and provide error bounds on the spectral norm of the Gram matrix. Experimentally, we show that these random features are effective at approximating the Tanimoto coefficient of real-world datasets
and are useful for molecular property prediction and optimization tasks. Future updates to this work will be available at http://arxiv.org/abs/2306.14809.

Tanimoto Random Features for Scalable Molecular Machine Learning
Austin Tripp
Sergio Bacallado Unversity of Cambridge Sb21160cam. ac
Unversity of Cambridge ajt2120ccmm ac uk
Sukriti Singh University of Cambridge 8s29710cam ac.uk uk
José Miguel Hernández-Lobato
University of Cambridge jmh2338cam. ac uk
Abstract
The Tanimoto coefficient is commonly used to measure the similarity between molecules represented as discrete fingerprints, either as a distance metric or a positive definite kernel. While many kernel methods can be accelerated using random feature approximations, at present there is a lack of such approximations for the Tanimoto kernel. In this paper propose two kinds of novel random features to allow this kernel to scale to large datasets, and in the process discover a novel extension of the kernel to real-valued vectors. We theoretically characterize these random features, and provide error bounds on the spectral norm of the Gram matrix. Experimentally, approximating the Tanimoto coefficient of real-world datasets and are useful for molecular property prediction and optimization tasks. Future updates to this work will be available at http: // arxiv. org/ ab///3////22./////////// 14809.
show that these random features effective
Introduction
In recent years there have been notable advances in the use of machine learning (ML) for drug discovery, including molecule generation and property prediction (Dara et al., 2022). Despite ceaseless progress in deep learning, conventional methods such as support vector machines or random forest trained on molecular fingerprints are still competitive in the low-data regime (Walters and Barzilay, 2020; Stanley et al., 2021). These fingerprints essentially encode fragments from molecule into a sparse vector, thereby compactly representing a large number of molecular substructures (David et al., 2020). They are extensively used in virtual screening for substructure and similarity searches as well as an input for MI. models (Cereto-Masssauué et al., 2015; Granda et al., 2018).
The Tanimoto coefficient (also known as the Jaccard index) stands out as a natural way to compare function on sets T3s or
such fingerprints. This coefficient is most commonly expressed function of non-negative vectors TMa (Jaccard, 1912: Tanimoto, 1958; Ralaivola et al., 2005; Costa, 2021; Tan al., 2016):
X X'(XXX') XUX
E. min((,,))) x() E. max(as,
XXX'
TasM(2,2))
I, r' E R (1)
=
and are binary indicator vectors representing sets and B respectively, then Ts B) T Tara (a, b). Therefore TMM be viewed generalization of Ts; for this reason it is sometimes called the "weighted Jaccard coefficient" or min-max coefficient.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

The Tanimoto coefficient is widely used in machine learning and cheminformaties to compute similarities between molecular fingerprints (Bajusz et al., 2015; O' Boylle and Sayle, 2016; Miranda- Quintana et al., 2021), chiefly because of the following properties:
1. Clear Interpretation: The value of I TMM((a (z, x') represents the degree of overlap between x and x' and is always between 0 and [xx x²) = only when x x'.
2. Kernel: T' (ss)) is positive definite (Gower, 1971; Ralaivola et al., 2005), meaning can be used as the kernel for algorithms like support vector machines or Gaussian processes.
3. Metric: 1 a (z, x²) is valid distance metric (typically called Jaccard/Soergel distance) and can therefore be used in nearest- neighbhu and clustering algorithms (Marczewski and Steinhaus, 1958; Levandowsky and Winter, 1971).
In this paper, present and characterize two efficient low-rank approximations for large matrices of Tanimoto coefficients. The first method, presented in section 3, uses a random hash function to index random tensor and enjoys exceptionally low variance. The second method, presented in section 4, uses a power series expansion of the Tanimoto similarity for binary vectors. This line of research also unexpectedly led to the discovery of new generalization of the Tanimoto coefficient to arbitrary vectors in RR TD,, which is also kernel and can be used to form a distance metric. In section we demonstrate experimentally that our random features are effective at approximating Tanimoto matrices of real-world fingerprint data and demonstrate its application to molecular property prediction and optimization problems.
Background: kernel methods and random features
Kernel methods are a broad class of machine learning algorithms which make predictions using positive definite kernel function k X x X' + I (Schölkopf et al., 2002). Common methods in this class are support vector machines (Cortes and Vapnik, 1995) and Gaussian processes (Williams and Rasmussen, 2006). Given dataset of data points, training most kernel methods requires computing the n X kernel matrix! K,j k(x(()) 206) (with O(n²) time complexity) and possibly inverting (with (n²²) time complexity). Because of this, applying kernel methods to large datasets generally requires approximations.
Given a kernel k, random features map is random function f: X' r> RMM with the property that f()) - f(x") approximates k(x, x') for every pair x, x' E X. The approximation is often exact in expectation:
By [f(x) f(x)] k(x,x') for all x. x, E &.
Random features allow the kernel matrix to be approximated as R, = [(x(1)). [(x Because this matrix has rank most M, this approximation generally reduces the cost of (nn))/O(m²) computations to (OMM) i.e. at most linear in n.
The seminal work of Rahimi and Recht (2007), which coined the term random features, gave general method based on Fourier analysis to construct random features for any Bochner or stationary kernel, for which k(xxxx²) xx) is function of x x'. This class includes many common kernels including the RBF and Matérn kernels, but excludes T- Subsequent works have proposed random features for other kernels including the polynomial kernel and the arc-cosine kernel (Liu al., 2021). However, there is no general formula to define random features for non-stationary kernels, such as TMM-
A random features map is sometimes called data-oblivious sketch, to distinguish it from other data-dependent low-rank approximation methods which depend on given dataset ),.)),.,, Examples data-dependent low rank sketches are the Nyström approximation and leverage-score sampling (Drineas et al., 2005, 2012). Although data-dependent methods may result in lower approximation errors for a given dataset, data-oblivious sketches are naturally parallelizable and useful in cases where the dataset changes over time (e.g. streaming optimization) or for ultra-large datasets which may not fit in memory.
We write 2; to refer the äth element of vector and (0) to denote the ith vector list.

---

3 Low-variance random features for Tanimoto and MinMax kernels
Outside of chemistry, the Tanimoto coefficient has been widely used measure the similarity between text documents and rank results in search engines. To quickly find documents with high similarity to user's s query, many prior works have studied random hashes for the Tanimoto coefficient, i.e. family of random functions X + {l...., , such that Pa (h(x) = h(z')) = TMM((Z,x')
(3) Although initially these hashes were only applicable to binary inputs (Broder, 1997; Broder et al., 1998; Charikar, 2002), more recent work has produced efficient random hashes for arbitrary non- negative vectors (Manasse et al., 2010; Ioffe, 2010; Shrivastava, 2016). this section we propose a novel family of low-variance random features for Taras (and by extension Ts) which is based on random hashes.
It is important to clarify that although the definition of random hashes in equation 3 resembles the definition of random features in equation they are actually distinct. Random hash functions output discrete objects (typically an integer or tuple of integers) whose probability of equality is TM while random features must output vectors in RM whose expected inner product is T TMM- M- If random hash naive approach may be to use a K -K-dimensional indicator vector as a random feature. Because hash equality is a binary outcome, the variance of such random features would be as (1 - as).). Realistic hash functions like that of loffe (2010) use K 2 103, implying a "variance per feature" of = 10?, which is undesirably high.
maps to {I....,'' K}.
Our main insight is that low-variance scalar random features can be created by using random hash to index suitably distributed random vector. In the following theorem, we show that a vector of i..dd samples from any distribution with the correct first and second moments can be combined with random hashes to produce random features for T.M-
Theorem 3.1. Let h X + y be a random hash for TarM satisfying equation 3, with | = K. Furthermore, let & be a random variable such that [E[{] O and TB[{2] 1, and let = -- Ex] be a vector of independent copies of &. Then the ID random features
$ (2) = En(x)
(4)
estimate Tars without bias: Ex, (Z) & (x')) = TaSM (2r, x²), and with variance
VE, (0z,((z) (2,((z)) = 1+TMM 2') (EK-) - 1- TTMM 2')) 2 (I, x')² (5) Furthermore, the lower bound is tight and achieved when & is Rademacher distributed (i.e. uniform in 13).
The proof is given in Appendix D.I. This theorem shows that Rademacher & yields the smallest possible variance in the class of random features defined in eq. (4).
These random features have many desirable properties. First, unlike random features for many other kernels such as the Gaussian kernel (Liu et al., 2021), the variance does not depend on the dimension of the input data or norms of the input vectors. Second, because these random features are 1 I-dimensional scalars, M independent random feature functions can be concatenated to produce M-dimensional random feature vectors with variance at most 1/M. This suggests that as few as s 10s random features could be used in practical problems. Third, although each instance of = requires storing dimensional random vector, if Ş is chosen to be Rademacher distributed, then each entry can be stored with single bit, requiring just R 100 kB of memory when K = 10
One disadvantage of these random features is that they are not continuous or differentiable with respect to their inputs. For applications such as Bayesian optimization which require optimizing over model inputs this would create difficulties as gradient-based optimization could no longer be done. It was this disadvantage which motivated us to search for other random features, leading to the discoveries in the following section.
4 Tanimoto dot product kernel and its random features
Ralaivola et al. (2005) gave a definition for the Tanimoto coefficient involving dot products:
T x
Tpr(o,x))
(6)
|
+ | - x
3

---

with TDP(2,2²) x²) = 1 when x' = 0. It is easy to check that TDP(2, x²) = TMM (z, x²) on binary vectors, which was used by Ralaivola et al. (2005) to prove that Dop is a kernel on the space {0, 1] referencing prior work by Gower (1971). However, TD not identical TasM for general inputs x, 6 E 1850 Here, we give the first proof that TDP is positive definite function in g and thus, also a valid kernel in this space.
Theorem 4.1. For x, # 0 in Rd, we have
Top(a, =
Σ (x Co-at x') (1<=112 + || rml
where the series is absolutely convergent. The function Dp is positive definite kernel in R
It has been noticed previously that, unlike 1 - TAS Mr the function 1 Tp is not a distance metric on non-binary inputs (Kosub, 2019). Indeed, when d 1, the inputs {1, 2, violate the triangle inequality. However, we can easily derive a distance metric from TDr.
Corollary 4.2. dDP(z,x')= = WV--TDP(a x²) corresponds to the RKHS norm of the function [TDp( - (2'',, and is therefore valid distance metric on R4
Proofs are given in Appendix D.2. These results imply that Tor, like TM is an extension of the set-valued Tanimoto coefficient (equation 1) to real vectors and can be used as a substitute for TMa in machine learning algorithms that require kernll distance metrie. Unlike TM the kernel Top is differentiable everywhere with respect to its inputs. It can also be computed in batches using matrix-matii matrix multiplication, allowing for efficient vectorized computation.
We now consider producing a random features approximation Tpp for large-scale applications. Motivated by the close relationship between TD and TMM, may be tempted to find a random hash for TDp and apply the techniques developed in section 3. Unfortunately, we are able to prove that this is not possible.
Proposition 4.3. There exists no random hash function for TDp over non-binary vectors.
Proof. Charikar (2002) proved that if s(x, is a similarity function for which there exists random hash, then s(x, x') must satisfy the triangle inequality (see their Lemma 1). Because 1 - TDP(x, ar, x²) does not satisfy the triangle inequality it follows by contradiction that there does not exist a random hash for TDp.
Therefore producing random features for Tor will require another approach. In the remainder of this section we present a framework to produce random features for Tor by directly approximating its power series (equation 7). We first describe method to produce random features for [|æ=|| ||2/+ ) (4.1). Then we describe how these features can be combined with existing random features for the polynomial kernel approximate Tpp's 's truncated power series (4.2-4.3). Lastly, we present an error bound for the kernel matrix of a dataset in the spectral norm, showing that the required dimension for the sketch scales optimally with the stable rank of the kernel matrix (4.4).
4.1 Random features for the "prefactor" (1|==12 + ||
In this section we present random feature map for the positive definite kernel (a. x²) + (||=||2- || [|2--r, which we will refer to as the the prefactor. We defer all proofs to Appendix D.3. We begin with the following lemma, which defines scalar random features for the prefactor:
Lemma 4.4. If 7 Gamma (s, c) (where e a rate parameter), then
(r,z(x) (11/2-11=12) ((r-))/2
/ (===
(8)
an unbiased scalar random feature for the prefactor (||x||2 + 1|=' for all s,c 0.
Although independent copies of Z could be combined to form an M-dimensional sketch, we instead
propose to use a dependent point set Zi...., Zm where each element Z, has a Gamma(s, c) distribu- tion whilst maximally covering the real line. This is well-established Quasi-Monte Carlo (QMC) technique which generally attains lower variance. We define our M -ddimensional QMC features in the following lemma:

---

Lemma 4.5. Let Ys.c be the inverse cumulative distribution function of Gamma(s, (s, c) random variable. Fix M, N, u E (0, 1), c, > O and let u, = 2/M 1u+ i/M] for ...... ............................... M. Define (x)) = (du.r, (z). ..... Pu.r, (x)), where:
(ur,(()) 1 Tss) -1/2)a/))))) (U:))( T(r)
VM
If u ~ U(0, 1) then our (x) forms unbiased random features of the prefactor (1|=|1 + ||-"
Although the random features are unbiased for all c > O, the value of these parameters will impact the error. We show that if 8, e are suitably tuned, then the relative error can be bounded:
Lemma 4.6. Let x ,!)),..., (m) E Rd with min, 112 enror
1| 2 C. and fix u E [o, 1]. Define the relative
E )u)( (x))) (1 12 + |,)) (|1=(0)]22 + (2)
(10)
c = 2 -r,, then for some constant C independent of this error satisfies,
2 I(r()sTS (1/e) C((M) M T(r)
max 15tJ5h
|
(11)
Together, these lemmas suggest random features for the prefactor can be created by first estimating (the minimum ratio of norms of input vectors), then using the random features from Lemma 4.5 with the values of 8, specified in Lemma 4.6.
4.2 A framework to produce random features for TDp
There are straightforward rules for producing random features for sums and products of kernels whose individual random features are known (Duvenaud, 2014, sec. 2.6.2). Random features for kernels k,, k2 can be concatenated (denoted @) form random features for the sum kernel h + l2- while their tensor product² (denoted ()) forms random features for the product kernel K1 k2. Our strategy produce features for T' is to combine random features for the "prefactor" (presented in section 4.1) with random features for the polynomial kernel to produce random features for TDP's 's power series (equation 7) truncated at R terms.
......,, R, let o, be a middiiinsi random features map for the prefactor [||= 12 + ||= and let vé, be a m mim-dimmnsion random features map for (x =')". The function:
Fix R E N, and for = 1.....
$ (x(z) = e (x) @ (x))
(12)
therefore random feature estimate for TDp power series, truncated at R terms. Unfortunately. these random features have dimension M = D m, which depends on the product of the random features dimension of é, and ti,,. Furthermore, the dimension m.. of the random features vi required to approximate the polynomial kernel (xr x²)" with good accuracy can scale poorly with r. For even modest values of m,, mi, the resulting value of M will likely be prohibitively large.
To remedy this, we turn to recent works which propose powerful linear maps to approximate tensor products with a lower-dimensional vector. Assuming (1)), (2) e ıdd and ,((2, Rd, these maps are effectively random matrices II E (d4dd)) which exhibit subspace embedding property whereby (1) @ (-) [I
& y2)) concentrates sharply around (z @ z(()). (g(2)) @ y())).
Critically, the product [1(x((1 @ (2)) can be computed without instantiating either matrix II or the tensor product (1) @ (2), Examples of such methods include TENSORSKETCH and TENSORSRHT (Pagh, 2013; Pham and Pagh, 2013; Ahle al., 2020), but for generality we will simply refer to these methods as SKETCH. Defining a series of sketches SKETCH, m x mRm + m., we can modify & at from equation into: TT((T) = @ SKETCH, (Ør (2). 14- (x))
(13)
For vectors E dd and y d2 detine the tensor product y veo(xy") E Rdid2
5

---

which has output dimension M = - mr, i.. without any pathological dependencies on the dimensions of (x),, (z). However, because these features approximate truncated power series, they will be biased downward due to the monotonicity of the power series (7). We propose two bias correction techniques to potentially improve empirical accuracy. One approach to normalize the random features such that ((x) ((x)) = 1 = (1r x)) for all E R i.e., such that the diagonal entries of the kernel matrix K are estimated exactly. A second approach is based on sketching the
residual of the power series. These approaches are presented in detail in Appendix E.
4.3 Implementing the random features
Instantiating the random features from the previous subsection (13) requires making concrete choices for R. SKETCH,,, my, my my Pr, the for all r, and choosing bias correction technique. There are many reasonable choices for SKETCH,, such as TENSORSKETCH (Pham and Pagh, 2013) and TENSORSRHT (Ahle et al., 2020). These sketches generally allow my; m/., My be chosen freely (although naturally error will increase as m, decreases). Many of these sketches can also be used as random features for the polynomial kernel ti,- (Wacker al., 2022), either directly or as part of more complex algorithms like TREESKETCH (Ahle et al., 2020) or complex- to-real sketches (Wacker al., 2023). The QMC random features from section 4.1 can be used for the prefactor Φr, with the parameters e chosen based on the anticipated norms of the input vectors.
The only remaining inputs are R (the number of power series terms to approximate) and m..,...., mn (how many random features to use for each term). Assuming a fixed dimension M for the final random features, this choice involves a bias variance trade-off, as a higher value of Rt will reduce bias but require each term in the power series to have fewer features, thereby increasing variance. Intuitivelyy because the terms of the power series decrease monotonically the variance of terms for small r is likely to dominate the overall variance, and therefore we surmise that a decreasing sequence for (m, 1 will be the best choice. Ultimately however we do not have theoretical results to dictate this choice in practice. We will evaluate these choices empirically in section 6.
4.4 Asymptotic error bound for Tor random features
Because many kernel methods use kernel matrices as linear operators, it is natural to examine the error of kernel approximations in the operator norm. Previous works have produced asymptotic error bounds of the random feature dimension m required to achieve a relative approximation error of in the operator norm. Defining (ir(K) [[(K)/[K]] Cohen et al. (2015) show that m (i ()))s(()) /s") essentially optimal for data-oblivious random features of linear kernels, even though it is possible to eliminate logarithmic factors. Our main theoretical result is that with the correct choices base random features, the random features for TDP presented in section 4.2 achieve similar scaling. We now state this as theorem.
Theorem 4.7. For any 1. let a!(1), ,(!,),..,,,,))) E Rd be set of inputs with min; (a max, | S. Let
K be the matrix with entries Ks,j : Rdd + 1. with m (((r(K)/=²)) such that
.o.)). For all 0, there exists an oblivious sketch
1
P& (1I Kllop =1|A[1oog llop)
(14)
poly(n)
where Riy = ((x(0)) (()())). Furthermore, the sketch can be computed in time Ö(r (K)ns (ee()))= - + n -3).
The random features in the theorem follow equation 13, with specific choices for M,R, and (SKETCH,, m, Φr, t ) given in the proof in Appendix D.4. Theorem 4.7 essentially suggests that, with the correct settings, the error of the random features proposed in section 4.2 scales as well as one could reasonably expect for kernel of this type. We would highlight that the computational cost of the sketch sub-quadratic n, and compares favourably with the cost of data-dependent low-rank approximation methods.

---

2,0
TMM Variance
Tane MSE
1.5
10
10
0.5
10
10
10
10
10"
0.2
0.6
1.0
Number of random features
T T (x, xx)
Rademacher, count
Gaussian, count
Rademacher, binury
Gaussian, binary
Figure I: Left: MSE of TM matrix reconstruction as a function of number of random features (median over trials, shaded regions first/third quartiles). Right: empirical variance of scalar Tar random feature estimates for M = 109, closely matching theoretical predictions (dashed lines).
Related work
Our work on random features fits into large body of literature random features for kernels (Liu et al., 2021). The majority of work in this area focuses on stationary kernels (i.e. k(xxx²) f(x x z')). because the Fourier transform can be applied to any stationary kernel produce random features in systematic way (Rahimi and Recht, 2007). There is however no analogous universal formula to produce random features for non-stationary kernels like Tим and TDP: therefore each kernel requires a bespoke approach. Although our random features are novel, they build upon ideas present prior works. Our random features for TMи critically rely on previously-pee propposed random hashes for TMм- Our approach to create random features for Tor via approximating its power series follows was inspired by the random features for the Gaussian kernel from Cotter al. (2011), which were subsequently improved upon by Ahle et al. (2020). Similar techniques have also been used to create random features for neural tangent kernel (Zandieh et al., 2021). However, to the best of our knowledge no prior works have proposed random features specifically for the Tanimoto kernel or its variants.
Other works have proposed other types of scalable approximations for Tanimoto coefficients which are based on random features. Haque and Pande (2010) propose SCISSORS, an optimization- based approach to estimate Tanimoto coefficients which is akin to a data-dependent sketch. A large number of works use hash-based techniques to find approximate nearest neighbours with the Tanimoto distance metric (Nasr et al., 2010; Kristensen et al., 2011; Tabei and Tsuda, 2011; Anastasiu and Karypis, 2017). Although these techniques are useful for information retrieval, unlike random features they cannot be used to directly scale kernel methods to larger datasets.
6 Experiments
this section we apply the techniques in this paper to realistic datasets molecular fingerprints. All experiments were performed in python using the numpy (Harris et al., 2020), pytorch (Paszke et al., 2019), gpytorch (Gardner et al., 2018), and rdkit (Landrum et al., 2023) packages. Molecules represented with both binary (B) and count (C) Morgan fingerprints (Rogers and Hahn, 2010) of dimension 1024 (additional details in Appendix F. F.1). These vectors indicate the presence (B) or count (C) of different subgraphs in molecule. Code to reproduce all experiments is available at: https:/ggithub.com wtttt//ooo///////taaaaaaaa T//annioto--ranneeoo www.oto--nnooooottoraoeeeeee nestturessnuurip2233
6.1 Errors of random features on real datasets
Here we study the error of approximating matrices of Tanimoto coefficients using our random features, with the general goal of verifying claims sections 3-4 on realistic dataset of molecules. We choose to study sample of 1000 small organic molecules from the GuacaMol dataset (Brown et al., 2019; Mendez al., 2019) which exemplify the types of molecules typically considered in drug discovery projects. We use both binary (B) and count (C) fingerprints of radius 2.

---

2-1 (count FP)
3 (count FF)
lc+02
10
count FP, ml
Median log 10
Prefactor MSE
count FP. r
dc
10"
binary FP.
-1.5
0.01
binary FP:
-3.0
(binary FP)
r=3 (binary FP)
10"
-.55
le102
6.0
0.01
-9.0
0.01
1er02
0.01
let02
10
10
10
10
s/s
s/s
Number of random features
Figure 2: Left: Contour plots of MSE for prefactor random features with M = 10 with varying 8, Right: MSE number of prefactor random features with s, e values from Lemma 4.6. As in Figure 1, lines are medians over 5 trials, shaded regions are first/third quartiles.
0.0020
Tor MSE (p=-1)
10
0.0015
TOP MSE
10
0.0010
10
10
10"
10"
10
Feature allocation p
Number of random fcatures
hias correction
normalize
sketch
Figure 3: Left: MSE for M 10a dimensional random features when allocating features by x yP, 1....,4. 4. Right: MSE of Trp random features using R = 4, p = -1 and various bias correction strategies. Both subplots use binary fingerprints (the equivalent plot for count fingerprints Figure F.2). As in Figure 1, lines medians over 5 trials, shaded regions first/third quartiles.
First, we investigate the random features for TMM proposed in section We instantiate these features using the random hash from loffe (2010) (explained further in Appendix F.2) with E both Gaussian and Rademacher distributed. The results are shown in Figure The left subplot shows the median mean squared error (MSB), i.e. Ei, (Thrr (x) x(5ر ((x(1) )(x()))² 21
function of the random feature dimension M. As expected for Monte Carlo estimator, the square error decreases with (1/M) (i.e. increasing the number of random features by 10 reduces the MSE by a factor of 10). As predicted by Theorem 3.1, the estimation error seems to depend only the distribution of E and not on the input vectors themselves: therefore the етrог curves for count and binary fingerprints overlap completely. The error is lowest when =\t is Rademacher distributed, although the empirical difference error seems small. The right subplot looks at the variance across across scalar random features, showing close matching with the predictions of Theorem 3.1. Overall these features behave exactly as expected.
Next, we investigate the random features for TDp from section 4. These features are more complex, so we start by studying the random features for the "prefactor" from section 4.1. Recall that these features had free parameters e > 0. Fixing the number of features M = 10*, Figure 2 (left) shows the MSE for the r and r terms for both the binary and count fingerprints (which have different norms) as function of s and c. The values which minimize the relative error bound from lemma 4.6, denoted 8", c², seem to lie broad plateau low error all settings, suggesting that these values of e are prudent choice. Using these values of Figure 2 (right) shows the MSE with respect the number of features M. expected for QMC method, the error dependence appears to be quadratic O(1/M²) (i.e. ten-fold increase in M reduces MSE by 100-fold).

---

Table 1: Average log probability of test set labels with various approximate GPs for 5 targets from DOCKSTRING dataset (García-Ortegón et al., 2022). + values are standard deviations over 5 trials.
KERNEL
METHOD
ESR2
F2
KIT
PARPI
PGR
Tur
RAND SUBSET GP
-1.084=0.004
-0.951=0.002
-1.09420.002
-0.999
-1.1834:0.00
SVGP
-0.908
-0.502-1.005
-0.846:200022
-0.606:000005
-1.0301:0.005
RFGP (E RAD.)
-0.954 0.0099
-0.658+0.011
-1.222400004
-0.968 +00037
-1.027+0.023
RFGP (E GAUSS.)
-0.956=0.010
-0.663=400155
-1.230==0.048
-0.967=0.036
-1.124+0.025
Top RAND SUBSET GP
-1.073=0.002
-0.940=0.001
-1.077=0.002
-0.988=0.000 0.00
-1.187=0.006
SVGP
-0.880=0.004
-0.459+0.00 002
-0.804+0.002
-0.568=0.000
-1.010+0.000
RFGP (PLAIN)
-0.902=0.003
-0.513 004
-0.979=0.015
-0.690
-1.029=0.002
REGP (NORM)
0..002000003
0..115
-0,98010.015
-0.69110.022
-10028
REGP (SKETCH)
-0.904 + 0.002
0.505
0.971
-0.600+000221
-1.030+0.002
Because it is implemented in Salkit-leary (Pedregosa et al., 2011), we use TENSORSKETCH (Pagh, 2013) both as the polynomial random feature map and to combine the polynomial and prefactor random features. We fix the number of random features for the prefactor to be 10 (recall it can be chosen freely without impacting the final random feature dimension). Figure F. shows the MSE of approximating both (x x²) and ((x x²)/ ([1=( + |=="|")))" In both cases, the MSE decreases approximately with ((/M)). Finally, we empirically examine how to allocate M random features across R terms. Using R=4, 4, Figure 3 (left) shows that allocating most of the features to the terms with small R results lower error. We therefore heuristically suggest allocating features according to m, r- Recall that truncating the power series biases the random features downward, and in section 4.2 two bias correction techniques were proposed. Figure 3 (right) studies the overall MSE for the plain features and both bias correction techniques. It appears that, in practice, neither technique is particularly helpful (normalization in fact appears harmful for large M). All techniques show an error dependence of approximately O1/MM).
6.2 Molecular property prediction and uncertainty quantification
To evaluate their efficacy in practice, we use our random features to approximate large-scale Gaussian processes (GPs) (Williams and Rasmussen, 2006) for molecular property prediction. Specifically, we study tasks from the DOCKSTRING benchmark which entail predicting protein binding affinity from a molecular graph structure (Garefa-Ortegoe et al., 2022). Each task contains 250k molecules, making exact GP regression infeasible. We represent molecules with count fingerprints of radius 1.
We use M = 5000 random features for all methods. We compare to two approximate GP baselines. The first is exact GP a random subset of size M. Since this approach ignores most of the dataset, one should expect reasonable approximate GP to perform better. The second is a sparse variational GP (SVGP) which approximates the dataset using M pseudo-data points Z (Titsias, 2009; Hensman et al., 2013). The locations of Z are typically chosen based on the input dataset (we use K-means clustering), making this method effectively a data-dependent sketch. Accordingly, one might expect the performance of this approximation to be better than data-oblivious random features. Details of Gaussian process training are given appendix F.4.
Table 1 shows the average log probability of test set labels for all types of GP with TArm and Tdp kernels. Several trends evident. First, for each kemeel random feature GPs (RFGPs) consistently outperform random subset GPs, but underperform SVGP. Second, for each kernel the difference between the RFGP varieties is small (generally less than the standard deviation). Third, most targets TDP seems to perform better than the TM kernel. The reason for this is unclear. Similar trends can be see in the R2 metric (Table F.I). This suggests that the RFGPs in this paper can be used in large-scale regression, although it seems in practice that data-dependent approximations are more accurate.
6.3 Bayesian optimization in molecule space via Thompson sampling
Bayesian optimization (BO) uses a probabilistic surrogate model to guide optimization and is generally considered one of the most promising techniques for sample-efficient optimization (Shahriari et al., 2015). Because wet-lab experiments are expensive and time-consuming. there considerable interest in using BO for experiment design. Chemistry experiments are often done in large batches, and

---

Average F2 score
10
Time (s)
7.8
10
-8.0
10'
10¹
10"
10"
Number data points
Number of data points
Tine (exact)
Tirw (RFs)
Tor (exaci)
Tor (RFs)
Figure 4: Run time and docking scores of BO using exact and approximate Thompson sampling. Solid lines are means over 5 trials, shaded regions are standard errors.
therefore algorithms which use functions sampled from probabilistic model are of particular interest (Hernández-Loobe et al., 2017). To sample from a normal distribution N N(u, K) one typically transforms i.i.d. samples Z ~ (0, 1) via μ + K1/27. This requires computing K1/2, and thereby causes exact GP sampling to scale cubically the number of evaluation points. By approximating K z q 4, our random features allow for approximate sampling in linear time. In this section we apply this to a real-world dataset.
As a demonstration, we consider a single round of selecting 100 molecules from a random sub-sample of n molecules using Thompson sampling, procedure for Bayesian optimization which selects molecules that maximize a function sampled from GP prior using the TaM and TDP kernels. Similar to the setup from section 6.2, we molecules and labels from the DOCKSTRING dataset. Molecules are represented as count fingerprints. M = 5000 random features used, with Rademacher E for Tим and no bias correction for Tpr. All other implementation details the same in the previous subsection. Figure (left) shows that, as expected, exact Thompson sampling scales worse than approximate Thompson sampling with random features. Figure 4 (right) shows that using approximate instead of exact Thompson sampling does not seem to change the average F2 docking scores of the molecules chosen. This suggests that approximate Thompson sampling could fruitfully be applied to large datasets of molecules in Bayesian optimization tasks.
7 Discussion and conclusion
In this paper we presented two kinds of random features to estimate Tanimoto kernel matrices: one based on random hashes and another based on a power series expansion. To our knowledge, this is the first investigation into random features for the Tanimoto kernel. We theoretically analyze their approximation quality and demonstrate that they can effectively approximate the Tanimoto kernel matrices on realistic molecular fingerprint data. In the process we discovered a new Tanimoto-like kernel over all of Rd which is a promising substitute for the more established TMM regression and optimization tasks.
Despite promising theoretical and experimental results, our random features do have some limitations. We found that it was difficult efficiently vectorize the computation of the random features for TarM, making them undesirably slow to compute. For TDP, we were able to exhibit an error bound on the spectral norm which depends on certain choices for the base sketch and sketch dimensions; however, it is unclear whether these choices are optimal in practice. Nonetheless, in Appendix D.5 we prove that exact low-rank factorizations of TArM and TDP kernel matrices are not possible; this means that follow-up works could reduce but never climinate the approximation error.
We are most optimistic about the potential of our random features be applied in Bayesian opti- mization, in particular by enabling scalable approximate sampling from GP posteriors (Wilson et al., 2020). Although we briefly explored this technique in section 6.3, in the future could allow for sample-efficient Bayesian algorithms for complex tasks like Pareto frontier exploration and diverse optimization using the Bayesian algorithm execution framework (Neiswanger et al., 2021). These tasks are highly relevant to real-world drug discovery and there are scant new methods poised to solve them in a sample-efficient way. We hope that the methods presented in this paper enable impactful, large-scale applications of the Tanimoto kernel and its two extensions in chemoinformatics.
10

---

Acknowledgments and Disclosure of Funding
We thank Isaac Reid and Zhen Ning David Liu for helpful discussions. Austin Tripp acknowl- edges funding via a C T Taylor Cambridge International Scholarship and the Canadian Centennial Scholarship Fund. Sukriti Singh acknowledges funding from the UK Engineering and Physical Sciences Research Council. José Miguel Heránndez-Lobad acknowledges support from a Turing Al Fellowship under grant EP/VO23756/1.
Author contributions: the initial idea of using random hashes to produce random features for the Tanimoto kernel came from Sergio. Austin came across Ioffe (2010) and realized his proposed hash could be used to extend the scheme to I MM for general non-negative vectors. Austin derived and proved the variance statement from Theorem 3.1. Sergio proposed and proved Theorem 4. and Corollary 4.2 for the TDr kernel. Austin proposed and proved Proposition 4.3 after reading Charikar (2002). Sergio developed the random features for the prefactor (section 4.1) and the overall schema for the random features (section 4.2). Sergio proposed and proved Theorem 4.7. All experiments were designed and performed by Austin. Sukriti helped with some experiments which ultimately did not appear in this version of the manuscript. Sergio and José Miguel provided advising throughout the project. Writing was done jointly, but mostly by Austin and Sergio.
References
Ahle, D., Kapralov, M., Knudsen, J. B., Pagh, R., Velingker, A., Woodruff, D. P., and Zandich, A. (2020). Oblivious sketching of high-degree polynomial kernels. In Proceedings of the Fourteenth Annual ACM-SIAM Symposium on Discrete Algorithms, pages 141-160. SIAM.
Anastasiu, D. C. and Karypis, G. (2017). Efficient identification of tanimoto nearest neighbors: All-pairs similarity search using the extended jaccard coefficient. International Journal of Data Science and Analytics, 4:153-172.
Ando, T., Horn, R. A., and Johnson, C. (1987). The singular values of a hadamard product: A basic inequality. Linear and Multilinear Algebra, 21(4):345-365
Bajusz, D., Ráez, A., and Héberger, (2015). Why tanimoto index an appropriate choice for fingerprint-based similarity calculations? Journal of cheminformatics, 7(1):1-13.
Bishop, C. M. and Nasrabadi, N. M. (2006). Pattern recognition and machine learning, volume 4. Springer.
Broder, A. Z. (1997). On the resemblance and containment of documents. In Proceedings. Compres- sion and Complexity of SEQUENCES 1997 (Ca. No. /77B/00/71), pages 21-29. IEBE.
Broder, A. Z., Charikar, M., Frieze, A. M., and Mitzenmacher, M. (1998). Min-wise independent permutations. In Proceedings of the thirtieth anmual ACM symposium on Theory of computing.
pages 327-336.
Brown, N., Fiscato, M., Segler, M. H., and Vaucher, A. C. (2019). Guacamol: benchmarking models for de novo molecular design. Journal of chemical information and modeling, 1993):1096-1108
Cereto-Massagué A., Ojeda, M. J., Valls, C., Mulero, M., Gareia-Valle S., and Pujadas, G. (2015). Molecular fingerprint similarity search in virtual screening. Methods, 71:58-63.
Charikar, M. (2002). Similarity estimation techniques from rounding algorithms. In Proceedings of the foury-fourth annual ACM symposium on Theory of computing, pages 380-388.
Cohen, M. B., Nelson, J., and Woodruff, P. (2015). Optimal approximate matrix product in terms of stable rank. arXiv preprint 1r0iv:1507.02268.
Cortes, C. and Vapnik, V. (1995). Support-vector networks. Machine learning, 20:273-297.
Costa, d. F. (2021). Further generalizations of the jaccard index. arXiv preprint arXiv:2110.0940
Cotter, A., Keshet, J., and Srebro, N. (2011). Explicit approximations the gaussian kernel. arXiv preprint

---

Dara, S., Dhamercherla, S., Jadav, S. S., Babu, C. M., and Ahsan, M. J. (2022). Machine learning in drug discovery: review. Artificial Intelligence Review, 55(3):1447-1999
David, L., Thakkar, A., Mercado, R., and Engkvist, (2020). Molecular representations in ai ai-driven drug discovery: review and practical guide. Journal of Cheminformatics, 12(1):1-22.
Drineas, P., Magdon-Ismail,, M., Mahoney, M. W., and Woodruff, D. P. (2012). Fast approximation of matrix coherence and statistical leverage. The Journal of Machine Learning Research, 13(1):3475- 3506.
Drincas, P., Mahoney, M. W., and Cristianini, N. (2005). On the nyström method for approximating gram matrix for improved sernel-based learning. journal of machine learning research, 6(12).
Duvenaud, D. (2014). Automatic model construction with gaussian processes.
Garcia-Ortegón, M., Simm, G. N., Tripp. A. J., Hernández-Lobato, J. M., Bender, A., and Bacallado, S. (2022). Dockstring: easy molecular docking yields better benchmarks for ligand design. Journal of chemical information and modeling, 62((5))3486-3335
Gardner, J., Pleiss, G., Weinberger, K. Q-. Bindel, D., and Wilson, A. G. (2018). Gpytorch: Blackbox matrix gaussian process inference with gpu acceleration. Advances in neural information processing systems, 31.
Gilmer, J., Schoenholz, S. S.. Riley, F., Vinyals, O., and Dahl, G. E. (2017). Neural message passing for quantum chemistry. In International conference on machine learning, pages 1263-1272. PMLR.
Gower, J. C. (1971). A general coefficient of similarity and some of its properties. Biometrics, pages 857-871.
Granda, M., Donina, L., Dragone, V., Long, D.-L., and Cronin, (2018). Controlling an organic synthesis robot with machine learning to search for new reactivity. Nature, 559(7714):377-388
Haque, I. S. and Pande, V. (2010). Scissors: linear-algebraical technique to rapidly approximate chemical similarities. Journal of chemical information and modeling, 1066):
Harris, C. R., Millman, K. J., Van Der Walt, S. J., Gommers, R., Virtanen, P., Cournapeau, D., Wieser, E.. Taylor, J., Berg, S., Smith, N. J., et al. (2020). Arra programming with numpy. Nature, 585(7825):577332.
Hensman, J., Fusi, N., and Lawrence, N. D. (2013). Gaussian processes for big data. In Proceedings of the Twenty- Conference on Uncertainty in Artificial Intelligence, pages 282-290.
HemándezzLoobto, J. M., Requeima, J.. Pyzer-Knapp, E. O., and ispuru-Guzik, A. (2017). Parallel and distributed thompson sampling for large-scale accelerated exploration of chemical space. In International conference on machine learning, pages 1470-1479, PMLR.
Ioffe, S. (2010). Improved consistent sampling, weighted minhash and 11 sketching. In 2010 IEEE international conference on data mining, pages 246-255. IEEE.
Jaccard, P. (1912). The distribution of the flora in the alpine zone. 1. New phytologist, 11(2):37-50.
Kosub, S. (2019). A note on the triangle inequality for the jaccard distance. Pattern Recognition Letters, 120:36-38.
Kristensen, G., Nielsen, and Pedersen, N. (2011). Using inverted indices for accelerating lingo calculations. Journal of chemical information and modeling. 51(3):597-600.
Landrum, G., Tosco, P., Kelley, B., Ric. Cosgrove, D.,. sriniker, gedeck, Vianello, R., NadineSchneider, Kawashima, E., N, D., Jones, G., Dalke, A., Cole, B., Swain, M., Turk, S., AlexanderSavelyev
Vaucher, A., Wójcikowski, M., Take, I., Probst, D., Ujihara, K., Scalfani, V. F., guillaume godin, Lehtivarjo, J., Pahl, A., Walker, R., Berenger, F., jasondbiggs, and strets123 (2023). rdkit/rdkit: 2023_09__0 (q3 2023) release beta.
Levandowsky, M. and Winter, D. (1971). Distance between sets. Nature, 334(5333)334-35.
12

---

Liu, F., Huang, X., Chen, Y., and Suykens, A. (2021). Random features for kernel approximation: survey on algorithms, theory, and beyond. IEEE Transactions on Pattern Analysis and Machine Intelligence, 44(10):7128-7144
Manasse, M., McSherry, F., and Talwar, K. (2010). Consistent weighted sampling. Unpublished lechnical report) http://reeearch. microsoft. com/en-su//peplehmmaasss 2.
Marczewski, E. and Steinhaus, H. (1958). On a certain distance of sets and the corresponding distance of functions. In Colloquium Mathematicum, volume 6, pages 319-327. Instytut Matematyczny Polskiej Akademii Nauk.
Mendez, D., Gaulton, A., Bento, A. P., Chambers, J., Veij, M., Félix, E., Magariños, M. P., Mosquera, J. F., Mutowo, P., Nowotka, M., et al. (2019). Chembl: towards direct deposition of bioassay data. Nucleic acids research, 47(D1):D930-D990.
Miranda- Quintana, R. A., Bajusz, D., Rácz, A., and Héberger, K. (2021). Differential consistency analysis: which similarity measures can be applied in drug discovery? Molecular Informatics, 40(7):2060017.
Nasr, R., Hirschberg, D. S., and Baldi, (2010). Hashing algorithms and data structures for rapid searches of fingerprint vectors. Journal of chemical information and modeling, 50(8):1358-1368
Neiswanger, W., Wang, K. A., and Ermon, S. (2021). Bayesian algorithm execution: Estimating com- putable properties of black-box functions using mutual information. In International Conference Machine Learning, pages 5005-8015 PMLR.
Nelson, J. and Nguyên, H. (2013). Osnap: Faster numerical linear algebra algorithms via sparser subspace embeddings. In 2013 ieee 54th annual symposium on foundations of computer science, pages 117-126. IEEE.
O'Boyle, N. M. and Sayle, R. A. (2016). Comparing structural fingerprints using a literature-based similarity benchmark. Journal of cheminformatics, 8(1):1-14.
Pagh, R. (2013). Compressed matrix multiplication. ACM Transactions on Computation Theory (TOCT), 5(3):1117.
Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., et al. (2019). Pytorch: An imperative style, high-performance deep learning library. Advances in neural information processing systems, 32.
Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cournapeau, D., Brucher, M., Perrot, M., and Duchesmay, E. (2011) 1). Scikit-learn: Machine learning in Python. Journal of Machine Learning Research, 12:2825-2833
Pham, N. and Pagh, R. (2013). Fast and scalable polynomial kemnlls via explicit feature maps. In Proceedings of the 19th ACM SIGKDD international conference on Knowledge discovery and data mining, pages 239-247.
Rahimi, A. and Recht, B. (2007). Random features for large-scale kernel machines. Advances in neural information processing systems, 20.
Ralaivola, I.., Swamidass, S. J., Saigo, H., and Baldi, P. (2005). Graph kernels for chemical informatics. Neural networks, 1993-111.
Rogers, D. and Hahn, M. (2010). Extended-coomectiitty fingerprints. Journal of chemical information and modeling, 50(5):742-754
Schölkopf, B., Smola, A. J., Bach, F., et al. (2002). Learning with kernels: support vector machines, regularization, optimization, and beyond. MIT press.
Shahriari, B., Swersky, K., Wang, Z.,, Adams, R. P., and De Freitas, N. (2015). Taking the human out of the loop: A review of bayesian optimization. Proceedings of the IEEE, 104(1):148-155
13

---

Shrivastava, A. (2016). Simple and efficient weighted minwise hashing. Advances in Neural Information Processing Systems, 29.
Stanley, M., Bronskill, J. F, Maziarz, K., Misztela, H., Lanini, J., Segler, M., Schneider, N., and Brockschmidt, M. (2021). Fs-mol: A few-shot learning dataset of molecules. In Thitty-fith Conference on Neural Information Processing Systems Datasets and Benchmarks Track (Round 2).
Tabei, Y. and Tsuda, K. (2011). Sketchsort: Fast all pairs similarity search for large databases of molecular fingerprints. Molecular informatics, 30(9):01-8877
Tan, P.-N., Steinbach, M., and Kumar, V. (2016). Introduction to data mining. Pearson Education India.
Tanimoto, T. T. (1958). Elementary mathematical theory of classification and prediction.
Titsia,, M. (2009). Variational learning of inducing variables in sparse gaussian processes. In Artificial intelligence and statistics, pages 567-574. PMLR.
Wacker, J., Kanagawa, M., and Filippone, M. (2022). Improved random features for dot product kernels. arXiv preprint arXiv:2201.08712.
Wacker, J., Ohana, R., and Filippone, M. (2023). Complex-to-rent sketches for tensor products with applications to the polynomial kernel. In Ruiz, F., Dy, J., and van de Meent, J.-W., editors, Proceedings The 26th International Conference on Artificial Intelligence and Statistics, volume 206 of Proceedings of Machine Learning Research, pages 5181-5212. PMLR.
Walters, W. P. and Barzilay, R. (2020). Applications of deep learning in molecule generation and molecular property prediction. Accounts of chemical research, 54(2):263-270.
Williams, C. K. and Rasmussen, E. (2006). Gaussian processes for machine learning, volume MIT press Cambridge, MA.
Wilson, J., Borovitskiy, V., Terenin, A., Mostowsky, P., and Deisenroth, M. (2020). Efficiently sampling functions from gaussian process posteriors. In International Conference on Machine Learning, pages 10292-10302. PMLR.
Xiong, Z., Wang, D., Liu, X., Zhong, F., Wan, X., Li, X., Li, Z., Luo, X., Chen, K., Jiang, H., et al. (2019). Pushing the boundaries of molecular representation for drug discovery with the graph attention mechanism. Journal of medicinal chemistry, 63(16):8749 8760.
Zandieh, A., Han, I., Avron, H., Shoham, N., Kim, C., and Shin, (2021). Scaling neural tangent kernels via sketching and random features. In Ranzato, M., Beygelzimer, A., Dauphin, Y., Liang, P., and Vaughan, J. W., editors, Advances in Neural Information Processing Systems, volume 34, pages 1062-1073. Curran Associates, Inc.
14

---

Outline of appendices
Appendix A comments the items in the NeurIPS Paper Checklist.
Appendix B summarizes the notation used in this paper.
Appendix C explains TREESKETCH from Ahle et al. (2020), which is used in the proof of
Theorem 4.7.
Appendix D provides proofs for all theoretical results in the paper.
Appendix E explains possible techniques to correct the bias of the DD kernel which are mentioned briefly in section 4.
Appendix F gives additional details about the experiments from section and presents some supplementary results.
15

---

A Paper Checklist
Here explicitly comment on all areas of the NeurIPS paper checklist.
Claims The key claims in this paper are the creation of random features for the Tanimoto kernel and the creation of new continuous kernel. The claims are mainly supported theoretically with proofs. but we also test our random features experimentally and show that they are usable with real-world data.
Code of Ethics + broader impacts Read and acknowledged. Our work is fairly abstract and theoretical and we do not work with human-related data, so we do not foresee significant direct ethical impacts (positive or negative).
Limitations Our work is primarily theoretical, so the main strength (random features with low approximation error) is also the main weakness (the approximation error could plausibly be lower). Random features are well-studied as general approach so the strengths and weaknesses are generally undl-und (Liu et al., 2021). For this reason we did not include an explicit limitations section. However, we did mention some other limitations in section
Theory Our key theorems have minimal assumptions which are clearly stated. All theorems are proved in Appendix
Experiments Code to reproduce the experiments is available at: hwwp://iiiiTaaann--n-am----eeeee
The experiments in this paper have minimal details and the important de- tails are specified in Appendix All unspecified details should be easy find in the code hwtp:///wwwsssiaaaaaaaoooraaanna
Training details
Error bars
Our tables and figures include error bars.
Compute
The compute costs of the experiments in this paper were quite modest and were run single machine with no GPU usage. The approximate computational times were:
script produce all the plots took approximately 24 to run. The long runtime was mostly due to computing errors for many settings R and p.
Section 6.1:
Section 6.2: for TMMM SVGP and each RFGP inference took approximately h (total 3 h per target). For ToP, SVGP took 0.5 h while RFGP took 0.3 h (total 1.5 h per target). Multiplying this by targets and trials gives total compute time of 1112h.
Section 6.3: this experiment quite fast: total runtime perhaps h
Reproducibility Our contribution reproducible via the statements of our random features and our code. We aim for high standard of reproducibility: the exact commands reproduce the results stated explicitly in our code's README file, and the raw data behind all plots in this paper included alongside the code.
Safeguards
We believe there is nothing high-risk which we need to safeguard.
Licenses
We cite the assets which use in the paper.
Assets
We are not releasing assets.
Humans Subjects/ IRB Approvals Not applicable to our paper.
16

---

B Notation
We use | denote the Euclidean norm of vector æ. ||Al|lu for the spectral norm and ||||| for the Frobenius norm of a matrix A. I; denotes the ith element of a vector æ, while xx() denotes the ith vector in a list of vectors. We denote nnz(-) the number of non-zero entries in vector or matrix. For two functions f. 9 we say f(z) = (zg(z) if there is a constant C, such that 0 I f(z) S Cg(z) for large enough. Similarly, f() = ([g(2)) if there is a constant C, such that f(z) Cg(z) for large enough. Ô and f2 omit poly-logarithmic factors.
Throughout the paper, we use d for the dimension of input vectors, n for number samples in the dataset, and M for the dimension of the sketch or random features.
17

---

C Definition of TREESKETCH
For simplicity, we shall define TREESKETCH (Ahle al., 2020) for inputs which are r-fold tensor products (11) .... @ (r) where power of two, and each x:1) Rd is of the dimension d.
The main ingredients will two base sketches, one the leaves of tree, and one for internal nodes. The leaf sketch, Tase is random matrix in R The internal sketch, Stases is random matrix in _mxm", which can be rapidly applied to the tensor product of two vectors in _ It is possible to instantiate TreeSketch using simple sketches such as COUNTSKETCHI, for the leaves, and TENSORSKETCR for internal nodes. However, our theory assumes that Tase is OSNAP (Nelson and Nguyên, 2013), and Stuse is TENSORSHRT, both of which enjoy a useful spectral property.
Having picked the base sketches, we can define for any power of two 2 2, random map Q" : R + R" whose action on a tensor product v()) @.... @ vl1) with v(t) E R is given by the following recursion:
Q = QQ/2 (S9 @ vl) @ (1 (
Sa 2 (2
@ vl
Here, the matrices (St) are independent copies of Stwe- The action of TREESKETCH II" on x:)) Ø X (r) is then defined by
II' (x 8.. Q" (T) @ .... "...:"""))...."""""""" TT..I)(
where the matrices (T:) independent copies of Tasee Figure C. gives schematic view of the computational tree for tensor product with = 4.
Whilst the action of 1I was only defined for tensor products, this can be extended to arbitrary vectors in Rd Ahle et al. (2020) prove that the resulting linear sketch has many desirable subspace embedding properties, some of which are used in the proof of Theorem 4.7.
S7
St
S1
Ti
T2
T3
T
((1)
(x22)
(3
xx4)
Schematic view of the computational for TREESKETCH for a tensor product of four vectors (1) ... & (4).
Figure C.1:
18

---

D Proofs
D.1 Proof of Theorem 3.1
For simplicity, we drop the subscripts h, from probabilities and expectations. To show that the random features are unbiased, we re-write the expectation:
E [0),((>) (=,n(=')]
PP(((z) = h(2))) (Mo)))))))))))== 1(xx)
+1((h(2)) + ((x))()))))))))()))) )) (ncs)laa)) (/())(()) & h(z²))
-P(h(x) h(x))) h(x))
(E(()) = EM(rr)
+ +P(h(z) * h(2')) Em(a) E [M(O))
(EM(s): Ba(r)) independent)
-1-P(h(z) h P(((z) + h(z))
(assumed moments of {)
TTMM (1,=²))
(h(x) is an unbiased hash for TM
Because in general V[x] = = B[X]2 and we have (x) (=,((=)) = TMr (z, x'), we only need to compute EE (d=,((x) (x))(" This can be done using a similar decomposition:
PP(((>) = (( (NE, (x')) (M(x)) m(x')] E (VE,A())) (=,((')))
+P(h(x) # k(x'))
(dE,(()) 11(x)) & h(z'))
-TMM(I, )
((r)
+ +(1 - (1-Tarm (I, x
22 E
(
h(x)
= -TMM(E, x') ')E
+ (1 (MMM(E,
h(())
Here, the last step follows from the moments of E which were fixed by assumption. Subtracting Tar (a, x²)² from the above yields the expression for the variance.
By Jensen's 's inequality, F[(*)= 2 E[ and B[(")= = 1 by assumption, so R[c") > 1. Substituting E[S*| into the equation for (( (z") gives the lower bound and completes the proof.
D.2 Proof of Theorem 4.1 and Corollary 4.2
Take any pair of inputs x, y * 0 in Rd. By the Cauchy-Schwarz inequality, we have
1 2
||00|//1lvv
(15)
1
|
11y1
| |
Now, we can write the function TDP as
2.9 + Mll
x. y
TDp y) = || + |12112 x
(16)
2-G [[=]]] w
y/
(17)
||=:| + ||8/1
J.
|| 10
I
(18)
|0 12 + 1
2 y
(19)
|
+ |
where in the identity (18) we use the bound (15) to assert that the series bounded by 2- and thus absolutely convergent.
19

---

Furthermore, because (x. + y is positive definite kernel, and so is (x, v) +
e ([al"++loll>
| 2 101
each summand the power series is positive definite because it is a product of positive definite kernels. Hence, sums and limits of positive definite kernels positive definite, and the power series is convergent, TDP is positive definite the space I {0}. The extension to include the vector E Rd is straightforward. as any Gram matrix including this vector block diagonal. We conclude that Tpp is positive definite kernel in R4
To prove Corollary 4.2, note that by the Moore-Aronszan theorem, there exists RKHS of functions (H, ()x) on Red with reproducing kernel TDP, such that for any x, E Rd (TDp(z, (Dp +)H TDp(z.)) y). Then, observe that
[[Tpp(z.)) Tnp -) (Tnr(z,,) TTppy. Tnp(a,) (n.
(TDr(I,))TDD((( TDP(D,, )u (Torly, )), Tor
2TTrr((, (DD(U))))
= 2 (2.
where in the final identity, we use the fact that (D)(E, for all Rd. Dividing by 2 and taking square roots on both sides, we obtain
| Tor((, 2 T
V TDr(I,y),
where the RKHS norm the left is clearly distance metric.
D.3 Proofs and derivations for prefactor random features from section 4.1
D.3.1 Scalar random features (Lemma 4.4)
First, we present a derivation for the scalar random features for the prefactor (||=]|/2 + ||~"|1²)~~ from Lemma 4.4. As the kernel is function of ||=:|]2 and || with the one-dimensional case. We begin by observing that for any a, b > 0,
deal, without loss of generality,
1)² a b
bie t' -
01/2- alt (1/2-b1e (010-a)e11/222-10100
I(r)
al/2-alt (1/2 b)t c e ()-1)) e² c'tr-le l
dt
[(r)
0
(0 c² c-
c
(1/2-000
dt.
(20)
r(r)
r(s)
Defining the function
(a) = (1/2-a) 2/1-))/2 (c-
(
1/2
T(r)
and letting Z Gamma(s, c), recognise the right hand side of (20) the expectation of 4p (2(()) (b). Hence,
E(pz(a)sz(6)).
a +
This proves Lemma 4.4.
D.3.2
QMC random features (Lemma 4.5)
Next motivate and derive QMC random features for the prefactor, proving Lemma 4.5. It possible to sketch a with a vector (p(a) - (Pz, (q).... PZ (a)) where Zi,..., Zm are independent copies of Z. This makes the kernel approximation pl()) (p(b) a Monte Carlo estimator
20

---

of the expectation, with error decreasing with the dimension M of the sketch at the standard rate (MM-1/2)
However, shall instead consider Quasi-Monte Monte Carlo (QMC) estimator with error decreasing at the faster rate (M ((M Let Ta.o be the inverse cumulative distribution function Gamma c) random variable. With change of variables, we write the integral (20) as
1 a +
*T(s) T(r)
/ (--aatb-ehhoool (T (T)) du.
(21)
Fix 0, u E [0, 1]. and let u, (u + i/M) - Lu + i/M for I,..., M. We can now define features (u,r(a) E R as in Section 4.2:
c ((()
(a)
(a- (a-c/2)m..((u)) (Taro(:)))
VM
T(r)
A QMC estimator for (a + b) is given by the dot product ou (b). This estimator is unbiased: Lemma. If u is random and (0, 1), then for all a, > 0,
B(dd,r(a) our
+
Proof. By linearity of expectation,
(T() [(r)
B(dur(() ou. (b))
E
-atd -att--)rn (vs,o(us)))
Observing that each u, ~ U(0, 1), the result follows by Eq. (21).
This lemma equivalent to Lemma 4.5, thereby proving it.
D.3.3 Relative error bound (Lemma 4.6)
In the following lemma, establish basic Quasi Munsi-Monte Carlo error bound. Lemma D.1. For any [0, 1), e (a +c+b), <ss
1 b)r Qu,r(a) (b) S M r(r) eela+b-o
2 T(s) r-s
r-
Proof. Consider the QMC approximation of an integral fo' f(r)dx for function f of bounded variation, using regular net of n points. It is well known that the enror of this approximation bounded above by V(f)/n where V is the total variation norm. Hence,
1
1 TT))(s V())
(a b)r (a) (u) (b)
T(r)
0 < < 1. As 7s. is monotone increasing, f is unimodal tending to as - 0 and u - 1. Thus V( = maxogus1 ('u). The maximum of is attained where (u)) (r s)/(a (a + c). hence
for fu) e (a(b-- ==00ab--cy..000
r 8-8
v(f) = claa + b
We are interested in sketching (11,(() + 11)) for every pair of inputs x(1), x((( in dataset (,()),...,,,.1)) The previous lemma can be used to choose values and a S which minimise the relative error
Es, (z.((((0)) (x(()) (11 + 12))-))))-
(112 + ||
21

---

The kernel k(z, y) is invariant to scaling z and i.e., k(x, v) = k(Ez, Ey) for all E >> 0. Thus, if we are interested in approximating the Gram matrix for (x())....,Ix)) ,("), we may assume without loss of generality that max; ||xx(")|| = 1 and let < = min; || Here, S will act as condition number for the dataset, which will determine the error of the random feature approximation.
By Lemma D.I, we have
|Bal s (1a000122 113 2 TTsss
r r-8 8
|| c) M T(r) c(|la100112
2 (r()(r 8)"=
1-5
=11=(0)|2 + +1
M ((r)e(r=s)
||x10|12 - ||
2 (())(r- -) (1 | + ||l
-8
M ((r)e(r-a)
|10000|12 + |a:t5
The last two factors are log-convex in c/(1>/0111 ([llg)elll (1l] 0. Hence, for a fixed e, the product
(1 |2:00 ( |
||~00| ||
is maximised at either ||= + ||ا = 2 or ||ا0 + || 1|2 2Ç. Therefore,
2 T(s)(r - 5)-= | S M I(r)e(r-s)
(1- 2) 2 v (1 2
2C
Using again the log-convexity of y + (1 - y v we find that
11 (1-2) 2 2
8-r
is minimised as a function of c at c 28/r. And
I
2
2C
is minimised at 02 = 2çs/r. Recall that Lemma D.1 requires 0 S c S
||= + || for all j.j
or 0 S S 2C, which is satisfied by 02 when S r, but not necessarily by c1. Thus, is reasonable to set c = 02, which leads to
2 [()
- 2 ((s)(rr - -)"-""'"" [1- (r-(6))* (8))*² (-s v (r s
(6 () (1-
C8 (1-2)²
(2) +
|Bib
M ((r)e(r-8)
M ((r)e(r-a)
2 (ss))s
(==S+) x- 5-vv [
M ((r)
When (= (all inputs 20(9) have the same norm), setting s = r leads to an error bound of 4/n. However, small values of Ç can make the upper bound blow up unless we tune s suitably. For a given value of Ç, this upper bound could be maximised numerically over < s S r. However, setting = rC, we obtain
2 ((()(r//))> [(1=5) r(c-1)
|r,5 -
51
M ((r)
2 P(rC)s R
(1/e)T(S-1)
((==5)
22) c-1 5-5 v1
< M
I(T)
2 TGC)S-
(1/e)r(s-1)
M T(r)
22

---

where p =2 1.2 is defined by
(1-() 5-S 1-
11- 1---2)--1 5-1
P
max 0S431
This proves the first inequality Lemma 4.6.
To conclude the proof of Lemma 4.6, we study the behaviour of the upper bound as it diverges when
Ç 0. Using the asymptotic (xx) 1/2 + O, have
2 [(rO) (r/e) (MC) 2(r/e) [(r)
M
rl((r)
where 2(r/e) "p" -I(r) the proof.
is bounded for r 2 Therefore, the upper bound is O(M This completes
D.4 Proof of Theorem 4.7
We begin by defining the stable rank, a common notion of effective dimension for data matrix A _dx its associated Gram matrix A²"
Definition D.2. The stable rank sr(A) of data matrix A _dx= is given by
||
(A( = ||A||
For positive semi-definite matrix K E IR*x", define sir( Si(K)= Tr(K)/||
The stable rank of A is upper bounded by the rank, but it is insensitive to small singular values, so be much smaller than min(n, d)). Note that sr( concept of stable rank to general kernel matrix K.
(ir A). The function sir( K)) extends the
We now re-state theorem 4.7 for convenience, and provide proof:
,..., .(),,..,,(() Fdd be a set of inputs with min: || Theorem. For any 2 1, let x11),
max, !| 2 Ç. Let K
be the matrix with entries Kt,j (x(( x(6)). For all 0, there exists an oblivious sketch
Rd + n. with m ((5((K)))") such that
P& (1 Kllop 2 = (lop) S poly(n)
(22)
where K., = (x)) ((x())). Furthermore, the sketch can be computed in time Ö(sr(K))z ne-
nnz(X)e + n Ec3))
We begin by specifying the sketch 4. Define
((z) = 1 4'(x)
where (x:) (ou, & x 8r)), II" is TREESKETCH (see appendix C) and (x.)(x)) is the random feature expansion for the prefactor defined in Lemma 4.5. TREESKETCH (Ahle et al., 2020), must be instantiated with the base sketches OSNAP (Nelson and Nguyên, 2013) at the leaves (Tbsse), and TENSORSRHT (Ahle et al., 2020) internal nodes (Suase. The choices for R, and the dimensions of II" and purr will be made explicit in the proof.
For each r E [R], define matrices K", K", and K" in Raxn by
x11) (1)
KJ, = |000112 |1008111
KY, (o(Der (x(()))))) (x (z(r)r (x(1),
Ki = (x((11)) 2r (x))),
23

---

Then, by the triangle inequality, we have
n
||
K
lop
R
+
Σ |K"
K"'lop +
Σ R
K
R
K"
K" lop + IK" K" llop
Σ
Σ
K
lop
(23)
-=++1
We shall find high-probabillty bounds for each of the terms on the right in turn. For the last term, we can choose R C Coog(ss(( (()K)=))) with large enough, and apply the following upper bound
K lop v
K" ||
R+11
r=
R+1
|x 2000
1/2
Σ
max
iJEEm
(| + ||x(5) 122)5
R+1
Σ
n2 n/21/2
$ n2 (2-1)/2
ne
El K ol
(24)
3sT(K)
where the final inequality follows from the fact that sir(K) n2/ || op. as the kernel matrix K has ones on the diagonal.
The terms in the second sum of (23) may be written as
|| KM
dlop || K'
where @ denotes the Hadamard product and E is the error matrix defined in (10). Noting that E is symmetric, apply Corollary 11 in Ando et al. (1987) assert that
|| K - K" loop I max E |||l
By the triangle inequality, for each r 2 1, we have || S ||p Hence, with the choice M = Q -11) for the dimension of @u.,, and applying Lemma 4.6 we obtain
8||K
Σ RK KT
(25)
3
Finally, we turn to the first term in (23). Let A, be matrix with columns (ier & pur (x(()) for I.....n. n. We can rewrite the error question
R
K" op [I TI) A, A7 Ar
In Lemma D.3, we establish that ||A, < (1+ =)T((K) and || S (1 + +)+ lop. Let
8 1/poly(n) (n) denote the error tolerance. Choose the dimension of the sketch II++1 to be m, R² r (i(K) K) log² m(dvMr) and the sparsity parameter of OSNAP as s 2² log m(dvar)
Ω
R²,"
Then, Lemmata 32-34 Ahle et al. (2020) ensure that,
<||A||op = (()TT+*,,)"'T' A. AZ. Arllop 2
=|| 1lop
[RRT R'lpop >
3R
3R
poly(n)
Taking union bound, we obtain
R Σ
<l|K K||K<lop 1lop
P
RK
K" R'
(26)
poly(n)
Finally, combining the bounds (26), (25), and (24), we obtain
P
( K 1l 2 E||K|] loo)
< poly(n)
24

---

Dimension of the sketch e; The total dimension of the sketch is m = mi + .... + ma.
R²)*a(K) log (M(dvA)) with R C'log(r((K)=")) Hence, ignoring poly- where m, Ω o²
logarithmic factors, the sketch has dimension m (((r(()//²)
Runtime: We first estimate the runtime of applying "" Applying OSNAP with sparsity parameter 8 to a vector w E IRd takes ((snnz(w) (w)) time. Thus, applying r independent OSNAP sketches to each x(()) for E [n]] takes ((srnnz(X)) time. Applying OSNAP to (u.r(z(0) for i E [n] takes O(nsM) time. Then, applying one TENSORSRHT sketch takes Om,, logm)) time, and therefore applying all the necessary copies of Shase requires Orrm, log m.). In total, the sketch tioi may be computed in ( log m, + snnz(X) nsM). Adding these runtimes over E [R] gives us a total runtime which is Ósr m 2 + nnz(X)e + 3
Lemma D.3. For each r = 1,..., R, we have || - + [1+E)]|K) lop and |A,11 S (1+))T((K)
Proof. Note that AT A. K Hence, aim to show that || I (1 + e) ||lp and that Tr(K*) S (1 + (K). By the triangle inequality,
|| lop S 1lop + IK - op S (1 + K> S (1 + olpp
where the penultimate inequality was shown in the proof of Theorem 4.7, and the final inequality is due to the fact that K N K.
Then, by linearity of the trace,
Tr(k") (T(KK)+TT(E - K") = Tr(K*) K")
(27)
And applying once more Corollary 11 of Ando et al. (1987), we have
Tr( Σ G(EOKT S max
E
Σ (i((K) SE Σ (K())
iela
where a-(-) denotes the ith singular value. Using once more that K" 3 K and hence o, (K(K) S (K( for all E [n]. we obtain,
Tr(E O (O) <T((K).
Plugging this into (27) yields Tr(RT) < TT()).
D.5 Rank of Tanimoto Kernels
Lemma D.4. There does not exist an exact fi
map for either TMM or Tpp.
Proof. We will use proof by contradiction, first focusing TMr The setup for the contradiction as follows: suppose there existed an exact feature map f R + Mr such that TMM y) = f()) fy) Vx, y Rd This would imply that any kernel matrix between n points could be written an inner product
(M, = -fx))
Because f outputs M dimensional vectors this would imply the resulting matrix has rank at most M. Therefore, under this hypothesis it should not be possible to form kernel matrix of more than M inputs which is full-rank (invertible). Our contradiction will be to construct such matrix.
We now present a way to construct full-rank TMM kernel matrix with any number of points n. Con- sider the set of points (n Jiz1 where a = (2n)' E R. For any i, we have that T (M (a()) = and if then
Σ (a()) (a))) j
-00) (2n)
1
2
3
jzfi
meaning the kernel matrix strictly diagonally dominant and therefore non-singular. Since such a construction exists for any n, setting = M + contradicts the assumption of an A -iimensiional
25

---

feature map, and repeating this argument for all finite M proves that no lnite-dimensional feature map can exists for T.M
The proof for Top is almost identical. Consider the same sequence of points as in the preceding paragraph. For any i, we have that (DP(a , = and
Σ (Tr((a ,00))= = JF1
<1,
(2m) 1- (2n)
2n
jti
j+i
meaning the kernel matrix for Top is also strictly diagonally dominant, thereby also precluding the existence of an --dimensiionll feature map for finite M.
This result suggests that the best we can hope for is an approximate finite-dimensional feature map for both kernels, which is exactly what is provided this paper.
26

---

E
Methods to correct the bias of Tanimoto dot product random features
Here, we describe the bias correction strategies mentioned in Section 4.2 and tested experimentally in section 6. We begin by noting that when x' E R I x' 2 0, and therefore, the power series
TDp(r,z')= Σ (o y') (1|=]]² !|
monotone. Therefore, if we unbiased sketch for the truncated series
R
Σ (a + |1=-113)~
1=1
the one constructed in Section 4.2, the final sketch will be biased downward:
(()()) (xx)) < TDr(r,x")
for all vr x E IRd 20 Below we introduce two strategies to remedy this issue.
E.1 Bias correction strategy 1: normalize the features
Empirically, we observe that the highest bias occurs in the diagonal elements of the kernel matrix. As the kernel satisfies (2p x) = for all x E Rdd one possible correction is to normalize the sketch obtain
(xx) (x)]]
(x)
This remains an oblivious sketch, and ensures that the diagonal of the kernel matrix is estimated exactly, the expense of possibly introducing some bias in off-diagonal elements of K.
E.2
Bias correction strategy 2: sketch the residual
To simplify the algebra, lz,y
The power series for TDp can then re-written as:
++bl
TDp(2.
(
(fw.)) + (tow) R+1 + Σ (ta
r= R+2
R
Σ
(tw + (tz.w) s+1
Rt + (tw.)) R+1 Σ (tz.))
R
Σ (twyy" + (tx.y) R+1 (1 (tzz)
TDP(xx v))
(28)
trincation error
RR
Critically, the truncation error can be written in terms of the kernel value itself, without any infinite sums. This enables simple procedure to generate random features for both the truncated power series and the remainder:
1. Compute and store (z:) as random features for the truncated kernel R.R
2. Concatenate a single onto (xx) to produce features o+1 (xx) = (1, (<x)) which approx- imate the kernel + R using only single extra dimension. Treat this as approximate random features for the kernel +
3. Compute random features f,41(2) (a:) for t
4. Apply SKETCH to the tensor product o (xr (0 (x) to obtain A(x), which approximates random features of the truncation error R+1
(1 + TDP(2, v)).
2.
27

---

5. Concatenate the random features (x) and (x)), to obtain bias corrected random features (ac(x)
(z).
Overall, these random features are essentially a concatenation of the random features for the truncated power series with an additional random feature estimate of the truncation error, which is formed using both the random features for the first R terms and the random features for the (R + 1)th term. A nice property the procedure above that it re-uses the random features (z) to estimate the error.
28

---

Experimental details and Additional Results
F.1 Datasets and featurization
The 1000 molecules from GuacaMol are included in our code. The molecules from the DOCKSTRING dataset are available online
com/ com/dockstring datanst))ttttt..
Note that when using count fingerprints for Top we used the square root of the counts as the fingerprint. This was done for two reasons:
1. To reduce the norms of the vectors (and thereby increase ()-
2. To roughly make their interpretation consistent with TM i.c. the "weight" a fragment in which occurs times in the numeratoor/demmiinato of TDP is n times the weight of fragment which occurs once.
F.2 Details of aa random features
We use the random hash from Ioffe (2010) in our implementation. To hash vectors in Rd random
variables ri, G, [(2, 1), Be ~ (0, 1) are drawn i.i.d. for i I,..., D, and then the following are computed for all i:
fx = [(nx;))rr++ Bel,
(29)
(x) rx(f:(z) - B)),
(30)
a,(x) In - -B((z) Inr;
(31)
i' (xx) arg min a,(x),
(32)
I....
(r (i" (z),))(()))))))
(33)
Note that equation 29 uses the convention that In -0. These variables do not have a clear interpretation in isolation procedure produces random hash for TMM. The hash itself is formed of 2 integers: i* E and 44- E Σ. This unfortunately means that the number of possible outputs is potentially infinite, which would require us to potentially sample an arbitrarily large vector E. To avoid this, we first use python's 's built-in hash function map this pair of integers module 212 4096. This allows us to sample a small (finite) vector E. and although it introduces small amount of bias this does not appear to be an issue practice.
refer the reader to loffe (2010) for explanation of why this hashing
D}
single integer, then take the result
Elements of are always sampled i.i.d. from either a Rademacher or Gaussian distribution (the distribution should always be clear from context).
F.3 Additional results for Top random features
Figure F.I shows that the error for both polynomial random features from TENSORSKETCH and random features for individual terms in Tops power series decreases with O(1/M). Figure F.2 shows the overall error for count fingerprints.
F.4
Gaussian process training details from section 6.2
Our GP models use a constant mean and Gaussian noise. Specifically, the model of the observed labels y is:
f(X) gP (kxX,())))) v(X) N (f(
Therefore, the GP hyperparameters are three scalars:
The constant mean, P
NNote that the presentation of equations 29-33 be the exponential of equations 30/31: we wrote it this way because in practice working in log space avoids numerical stability issues. This is explicitly suggested in their paper.
33 differs slightly from loffe (2010) who defines yı and a; to
29

---

10a
Polynomial MSE
10
Top term MSE
10"
10
10
10
10
10
10
10
10
10
10'
10
10
10
Number of random features
Number of random features
binary FP.
binary FP, c
Figure F.1: MSEs for approximating (æ r')" (left) and ((x x / 2 + ||>")²²)) (right) using SKNSORSKETCH for various count and binary fingerprints function of the random feature dimension M.
10
0.0020
Tor MSE (p=
TDp MSE
0.0015
10
0.0010
10
0.0005
10
10"
10"
10
Feature allocation
Number of random features
biss correction
normalize
sketch orror
Figure F.2: Same as Figure but with count fingerprints.
The kernel amplitude/outputsaale
The observation noise 2²
GP performance will be greatly affected by the choice of kernel hyperparameters. To ensure that the difference in performance is not due differences in kernel hyperparameter settings fit them in consistent way for all methods. Specifically, for all methods, we start by fitting an exact GP to random subset of M 5000 data points by maximizing the marginal likelihood with L-BFGS. The different approximations are as follows:
gpytorch's 's exact GP inference is applied on random subset of M data points (a different subset then used fit the hyperparameters).
Random subset:
SVGP: acikit-learn' 's K-means clustering is run with M clusters to produce an initializa- tion of the inducing points. These inducing points are used initialize sparse variational GP (Hensman et al., 2013), implemented in gpytorch. The variational parameters are optimized via natural gradient descent with a learning rate of 10- and batch size of 2M = 10 000 for one pass through the dataset. Although the inducing point locations themselves could be further optimized with gradient descent, we chose not do so for this experiment.
RFGP: First, the training and test sets are converted into M dimensional random features. Then, the posterior equations for inference in Bayesian linear models are used to make predictions on the test set given the training set (Bishop and Nasrabadi, 2006, equations 3.49-3.51). The computation is done in specific order to avoid forming any n matrices (the Woodbury matrix identity is used extensively for this). This fairly clearly documented in the code.
30

---

Table E.l: Average R² score for approximate GPs on DOCKSTRING dataset. Attentive FP and MPNN results are taken from Carcfa-Ortegón et al. (2022). Other details are the same as Table 1.
KERNEL
METHOD
ESR2
KIT
PARP1
PGR
TMM
RAND SUBSET GP
0,514=0.00 0.002 0.578-0.001
0.810=1.002
0.695:0.0002
0.849:0.000
0.426 0.00
SVGP
0.861+0.000
0.749+00.00
0.889+00.000
0.542+0.000
RFGP (= RAD.)
0.518=0.002
0.838=0.001
0.703+0.002
0.864+:0.001
0.465=0.003
RFGP (2 GAUSS.)
0.517=0.002
0.837=0.000
0.70210.001
0.86410.001
0.467 0.004
TDP
RAND SUBSET GP
0,513=0.003
0.817=00001
0.696:000022
0.851+0.000
0.384=0.011
SVGP
0.581=0.001
0.865=00.000
0.753+000001
0.889+00000
0.543 0.002
RFGP (PLAIN)
0.546=0.001
0.852=0.001
0.716:10.002
0.876100000 0.000
0.512=0.002
RFGP (NORM)
0.546=0.001
0.852-10.001
0.715
0.876 0.00
0.513=0.002
REGP (SKETCI)
0.545=0.001
0.852=0.001
0.716:10.002
0.876+0.000
0.510=0.002
N/A
MPNN
0.506=0.001
0.798
0.755 005
0.81540.010
0.324
N/A
ÁTTENTIVE FP
0.627 0.010
0.800
00.008
0.910+0.002
0.678 0.008
Note that for Ta ar it was vital to implement the kernel as
| |
||
||
Tarm (z) x²)
(34)
|
||00||
||
+
instead of naive implementation which follows equation 1. This because such an implementation requires forming a tensor of shape N x M x d when calculating the kernel between N and M points in 1R (for example Tijk min(2;4 Vjk)) which can exceed memory limits for modest N, M, This identity allows us to use the relatively efficient torch. cdist function. Note that we did not invent this identity ourselves; we discovered it in loffe (2010).
Even with this implementation, M 5000 inducing points did not fit into GPU memory for T TM. so all experiments were run CPU.
E55 Metrics and additional results from section 6.2
Metrics Table reports average log probability, which calculated by first calculating the log probability of each test point individually (i.e. marginally, not jointly), then averaging these val- ues. We also report the coefficient of determination (R² score), calculated using the function sklearn. metmiis......... r2__r This measures only the error of the GP mean. A value of indicates perfect prediction, while a value of 0 can be achieved by predicting the sample mean for every data point.
Additional results Table F.1 reports the average R². Trends are similar to Table 1. We also include baselines for types of graph neural network: Attentive FP (Xiong et al., 2019) and MPNN (Gilmer et al., 2017). Although the performance of GP methods does not match that of Attentive FP, is often close, suggesting there is potential for approximate GPs to be competitive with graph neural networks for molecular property prediction.
31


Reviews:
[{"review_id": "NOWh0t23no", "reviewer": "Reviewer_stq5", "summary": "The paper derives two random feature maps for Tanimoto kernels of fingerprint data. The feature maps are applied in GP and BO settings to allow accelerated primal learning.", "questions": "- “In the early stage..”. Not sure I agree with this statement. Big datasets are used in all stages of drug development. It’s also not clear why deep learning methods struggle, and there is no citation. Can you explain?\n- \\cdot is undefined\n- What does random function mean? Does it change every time? Is it stochastic? Is it random? Random in what sense? What is a random feature map? The presentation lacks rigor.\n- What does eq 2 mean? What does E_f mean? What is the measure or probability of f? It seems that we take infinitely many feature maps, and require that they together average to a kernel’s feature map. This feels strange, where is this coming from? It does not seem to be RFFs, since in RFFs the expectation is over frequencies (ie taking m to infinity), not functions. Eq 2 doesn’t require that any particular f would approximate the correct kernel value, and thus the statement hatK = f.f is not true based on eq 2.\n- “There is no general formula to define random features for non-stationary kernels”. Not sure this is correct. Please check seminal papers Kom-Samo et al “Generalised spectral kernel”, and Shen et al “Harmonizable mixture kernels with variational Fourier features”; and followup works.\n- P_h is undefined\n- What is “discrete range”?\n- What is an “independent copy”? So we copy the number independently, resulting in just identical copies (surely all copies are independent)? Or is the idea that each is an iid sample of the random variable? Is \\xi scalar?\n- What does \\Xi_h(x) mean? I don’t understand where the h(x) goes. I think h(x) is a selector function that gives one choice out of K choices, and then we pick that \\xi from the list of \\xi’s? But that would be just \\xi_{h(x)}, so \\Xi_h(x) should then be something else. Maybe it’s a set of K \\xi’s that we somehow index? It would help a lot if you would describe the domain and size of all symbols and variables in the paper (eg. is \\phi(x) a scalar or vector?). It’s also annoying that \\cdot’s are used to what are seemingly scalar products (line 136, or maybe its not scalar? Ok, line 144 declares them as scalars.)\n- What is the E_\\Xi,h expectation over? Are there measures over \\Xi and h? What are they?\n- What is a “tight” lower bound?\n- Theorem 3.1. states that \\xi needs to have variance 1. Yet, the optimal \\xi is declared to be Unif(-1,1), which has variance of 1/3. This seems to nullify the entire Rademacher result. Can you explain?\n- Why is the variance at most 1/m?\n- I’m confused what the 1/m analysis actually says. Let’s assume we have million-dimensional input vectors. Now we seemingly only need 1000 features to accurately represent them. Surely this can’t be true. What if we have trillion input dimensions? I think this analysis has some caveats that are not made transparent. Can you explain? Furthermore, the eq 4 requires us to average both \\Xi and h’s. It seems that the text is forgetting the \\Xi part. I’m also now realising that I don’t understand what the “m” random features refer to. Is this m hashes, or m \\xi’s, or m \\Xis, or something else?\n- I don’t understand what the 1-bit remark means. What is an “entry” and how come it needs only one bit? Why do we get 100kB?\n- I don’t understand what is the contribution in sec 3. It seems that a new feature map is introduced. Ok, but how is it novel? What problem does it solve? What is its significance? How does it improve over earlier stuff? The text mentions low-variance, but there has not been any discussion of variance of earlier methods. How do we know the variance is then lower?\n- Sec 3 also feels vague: the hash has not been described. Surely the hash affects things and can’t be ignored? What if we have a really crappy hash that outputs constant value: surely then the line 136 would not be true. Please make all assumptions explicit.\n- Why eq 7 doesn’t support x=x’=0? The equation seems to work for it as well.\n- Are m_r related to each other?\n- I don’t understand line 194. So phi is the feature map of kernel (x+x’)^-r, but why does the feature inner product add squares into the kernel? This feels surprising. Why do we even care about the non-square kernel, line 190 doesn’t have that either.\n- It’s not clear what u or i mean conceptually. What is their role?\n- The story seems to restart at 4.2. Suddenly we are talking about data matrix, it’s linear kernel and its rank. I can’t follow what’s happening. Is a paragraph missing between 4.2. and end of 4.1.? Can you explain? It’s also strange that we redefine data as A, but later keep using datapoints as x, and also X appears (whatever it is, maybe data again?).\n- What does min_i |x_i|^2 mean? I can’t follow the notation. So it seems that we take the norm of the data vector, which gives a scalar, and then its minimum is just the scalar itself. Or is this the smallest element in x_i? Is x_i a vector or an element of a vector? It would be helpful to use boldfacing on vectors.\n- What does zeta mean? In usual fingerprints we have always some zero values, which would indicate that zeta is zero. But then the eq 204 seems to get screwed up. Does this mean that the theory of sec 4 requires the fingerprints to be strictly non-zero or strictly positive? That would be a significant limitation.\n- What is _op?  What is poly(n)? What is tildeO? What is tilde\\Omega?\n- It seems that the paper hides the actual method description to the appendix, or borrows the TreeSketch/SRHT/OSNAP methods to do the heavy lifting. This makes it a bit difficult to follow what’s going on. Could you summarise the main ideas of these, and how does your method extend them?\n- In 5.2. I wonder what is the computational cost between T_DP exact and approx? The exact is just few matrix products, while the approx is something more involved. Why is there such a big running time difference?\n- In table 1 are the methods approx or exact?", "limitations": "No issues", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "The idea of featurising Tanimoto is clever and novel. The two feature maps are seemingly well derived, and theoretical analysis is provided. The paper shows how these ideas open the door for feature-based learning in GP setting, which is potentially significant.  \n\nThe results demonstrate the efficiency of the feature maps, and show good results when applied to BO/GP problems.", "weaknesses": "The paper clarity is frustratingly low, and it lacks a lot of rigor and its claims or contributions are not very clear.\n\nThe paper is vague about the cost or running times of the new feature maps, and how they relate to exact kernel computation. \n\nComputation of Tanimoto is just few matrix products, and should be fast. I don’t think the paper has demonstrated sufficiently what was the bottleneck earlier in using Tanimoto, or why do we need these new feature maps in the first place. Furthermore, if we have relatively small fingerprints and linear kernels, why do we even need GPs at all: wouldn’t a Bayesian linear regression do the same thing and avoid any need to construct feature maps? Why don’t we just use the datapoints as feature vectors directly? It seems that the paper is motivated by plugging Tanimoto’s into existing (and complex) GP frameworks, without first dismissing simpler direct approaches.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your thorough review and thoughtful comments. We appreciated that you found our work to be “clever and novel” and recognized that it “open[s] the door for feature-based learning in [the] GP setting, which is potentially significant.” Due to the 6000 character limit for the rebuttal we will try to answer all the _general_ questions/concerns below, then post specific answers to all questions as official comments.\n\n**Concern 1:** clarity of notation\n\nWe agree that we could make some improvements here. Based on the reviewer’s questions, we propose to make the following notation changes:\n\n1. We will use superscripts to denote vectors in a set and subscripts to denote indices of a vector. Previously there was some ambiguity because they were both denoted with subscripts.\n2. We will make the definition of random variables and expectations more explicit throughout.\n3. Keep computational complexity terms like poly(n), but define them in the main text rather than referring the reader to the appendix.\n\n**Concern 2:** necessity of random features\n\nDue to the page limit we gave only a brief introduction to random features in sections 1-2. We realize this may not be enough to help readers unfamiliar with random features to understand their advantages and disadvantages. The key points here are:\n\n1. It is an approximate method which sacrifices accuracy to improve the time / memory complexity of model training/inference for large datasets.\n2. Given a dataset of size N, it approximates the kernel matrix $K\\approx\\Phi \\Phi^T$, where $K$ has shape $N\\times N$ and $\\Phi$ has shape $N\\times M$. $M$ controls the quality of the approximation (larger is better, but more expensive). If $M < N$ then it will be faster than an exact kernel method.\n3. Computing $\\Phi$ depends only linearly on $N$ (but could have any dependence on M and the input dimension D).\n4. This generally reduces the complexity of operations to be linear in N instead of quadratic or cubic. For example, GP inference is reduced from O(N^3) to O(NM^2) (assuming $M < N$).\n\nTherefore, even though the computation of Tanimoto is “just a few matrix products” as the reviewer points out, the bottleneck is high time/memory complexity when N is large.\n\n**Concern 3:** why not use linear models\n\nAlthough using a linear kernel and doing Bayesian linear regression is another way to avoid high computational complexity, this model is less expressive than Tanimoto. For example, with D-dimensional inputs a GP with a linear kernel can only interpolate through D data points (regardless of their spacing) while a GP with a Tanimoto kernel could interpolate through an arbitrary number of data points provided they are spaced sufficiently far apart.\n\n**Concern 4:** running times of the new features maps\n\nWe did not state this because it depends on specific choices of hyperparameters. For $T_{MM}$ the cost is $O(NMD)$ if the hash of Ioffe (2010) is used. For $T_{DP}$, the cost is at most $O(NM^2 D)$. However, please note that:\n- This cost is _additive_, not multiplicative with the inference cost above (so the overall cost is still linear in N)\n- Unlike kernel inference, computing $\\Phi$ can be done in parallel for different data points, so a high computational cost does not necessarily imply a long runtime.\nWe will add this information to the main text.\n\n**Concern 5:** novelty / improvement of random features over previous work\n\nTo clarify, we are not aware of any prior works proposing random features for the Tanimoto kernel. We view them as the “first of their kind” rather than an improvement of prior work. Of course, to develop these features we leveraged innovations and developments of prior papers, notably random hashes and TreeSketch.\n\n**Concern 6:** general spectral random features for non-stationary kernels\n\nWe were not aware of the works by Kom-Samo et al and Shen et al: thank you for pointing them out. These works are based on a generalization of Bochner’s theorem for harmonizable, non-stationary kernels. Although the approach does seem broadly applicable, we noticed the following limitations:\n\n1. Although the theorems in these papers  guarantee the existence of a spectral representation  for harmonizable kernels, this may not be easy to derive or indeed take a simple form for a given kernel.\n2. The representation does not lead in general to a random features approximation (See the discussion in Section 2.2.1 of Shen’s dissertation: Spectral Kernels for Gaussian Processes).\n3. One could attempt to approximate the Tanimoto kernel with the sparse mixture kernel by Kom-Samo and Roberts or the harmonizable mixture kernel of Shen et al. These approximations sometimes lead to random features. However, even though these mixture kernels are dense in the space of non-stationary kernels, there are no results quantifying how many mixture components would be necessary to get a reasonable approximation.\n\nDespite these limitations, we think this is a very relevant suggestion which is worth investigating in the future. We will add citations to this literature in the final version of the paper.\n\n**Summary:** We hope that this response and the associated comments answering your questions have adequately addressed all of the reviewer’s concerns. Please let us know if you have any other questions or concerns."}, {"review_id": "o5LmwrS5iN", "reviewer": "Reviewer_b6mP", "summary": "This paper presents two methods for approximating a Tanimoto kernel using random features, shows that the approximations are good, and validates the utility of the proposed approach by applying it to three different settings using molecular fingerprint data.  The writing is clear, and the contributions are clearly stated.  Scaling the Tanimoto kernel to larger datasets is of practical importance.\n", "questions": "I was confused by the status of Theorem 4.7.  Is this a novel contribution, or just a restatement of something previously shown by (Cohen et al. 2015)?\n\nWhat is the intuition for why T_{DP} achieves better  accuracy than T_{MM} in Section 5.3?\n\n", "limitations": "N/A", "rating": 8, "confidence": 2, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "Clarity of exposition.\n\nA nice combination of theoretical and empirical results.\n\nThe experiments are extensive and show convincingly that the method works well on a variety of real datasets.\n\nI particularly liked the background in Sections 2 and 3, which set up the contribution of this work very nicely.\n\nThe construction of the data-oblivious sketch (Section 4.1) is elegant.\n", "weaknesses": "One of the main theoretical result is really just to show that Equation (6) is positive definite not just on binary inputs (which was shown in 2005) but for real vectors.  This seems like a relatively straightforward result.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your kind review and very positive rating. We will respond to the negative points and questions from your review below.\n\n> One of the main theoretical result is really just to show that Equation (6) is positive definite not just on binary inputs (which was shown in 2005) but for real vectors. This seems like a relatively straightforward result.\n\nThe proof that $T_{DP}$ is positive-definite is indeed straightforward, although we would highlight that our proof technique is novel, and in fact much more concise than previous proofs of positive definiteness for binary inputs. We do however agree that it is not the most technically involved result of our paper (that would be theorem 4.7).\n\n> I was confused by the status of Theorem 4.7. Is this a novel contribution, or just a restatement of something previously shown by (Cohen et al. 2015)?\n\nThis was a novel theorem. We mention Cohen et al (2015) because that paper pioneered the study of approximate matrix multiplication in terms of stable rank, and because it gives a general bound for the error of approximate matrix multiplication, which sets the standard of what kind of error bound is considered “good” for this type of problem. However, the specific problem studied in Cohen et al is matrix multiplication, whereas our paper focuses on _creating_ random feature matrices for the Tanimoto kernel.\n\n> What is the intuition for why $T_{DP}$ achieves better accuracy than $T_{MM}$ in Section 5.3?\n\nWe don’t really know. It is hard to say why one kernel performs better than another kernel on any specific dataset. It could also be related to the hyperparameter tuning for the inducing point optimization.\n"}, {"review_id": "cQtIDp8BOc", "reviewer": "Reviewer_A9mq", "summary": "This work develops random feature maps to approximate the Tanimoto coefficient among molecular datapoints. Approximating the Tanimoto coefficients allows this measure to be used as a kernel in classical machine learning approaches on larger datasets such as clustering and Bayesian optimization.", "questions": "Please explain L156 more: $x,x' \\in \\mathbb{R}^d = \\vec{0} \\Rightarrow T_{DP}(\\vec{0}, \\vec{0}) = 1$? This does not appear to be a continuous function and therefore not differentiable at 0 (which is at odds with the prior paragraph L149 motivating the need for differentiable random feature maps for Bayesian optimization). Perhaps reformulate (6) as a sequence of functions (i.e. \\begin{cases}..)  to include 0 and then argue for continuity at 0 ($T_{DP} \\rightarrow 1$ as $x,x' \\rightarrow 0$).", "limitations": "The authors have adequately described the limitations of their approach.", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "The presentation and technical precision of this work is excellent. The authors provide sufficient background to their approach and logic such that even a reader with a rudimentary understanding of analysis and algebra can comprehend their work as well as understand the chain of logical events constituting the work's contribution to research. While much of my research focuses on deep learning, I am deeply appreciative of the authors' work.\n\nIt is my opinion that the primary contribution of this work is to develop approximations to the Tanimoto kernel. This enables its usage in classical non-parameteric machine learning algorithms such as bayesian optimization and clustering. The authors present two different approaches---one based on hashing and another on power-series expansion---that each convey disparate benefits and drawbacks. While outside my primary research focus, the solution is quite technically interesting, and I believe it would be relevant for non-deep learning application of molecular science. \n\nThe empirical evaluation on binding-affinity prediction is also well-received as it illusatrates how approximating the Tanimoto kernel can be translated into a real-world application. \n\nA final point of appreciation relates to the clarity and honesty of the writing style and evaluation --- the authors do not seek to claim sota or surpass the performance of the most recent deep learning methods, but rather make a fundamental contribution that would be relevant to scaling non-parametric algorithms operating on molecular data to larger datasets. It is also refreshing to read a paper that is very precise with mathematical jargon, using accurate terminology to leverage existing structure in analysis and not proving theorems that are only tangentially related to the thrust of the work.", "weaknesses": "With the continued march of progress in deep learning algorithms---and their application to domains of molecular science---it is likely this approach will not see significant stand-alone use in applications. The counterpoint to this estimation is that there's exciting work integrating deep learning with classical non-parametric algorithms, and such approaches may leverage non-parametric approaches as inductive biases to assist representational learning. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your kind review and very positive rating. We will try to respond to your questions / concerns below.\n\n> With the continued march of progress in deep learning algorithms---and their application to domains of molecular science---it is likely this approach will not see significant stand-alone use in applications. The counterpoint to this estimation is that there's exciting work integrating deep learning with classical non-parametric algorithms, and such approaches may leverage non-parametric approaches as inductive biases to assist representational learning.\n\nWhile deep learning has many strengths, particularly on large datasets, we think the approaches in this paper do have their niche:\n\n1. Kernel methods are generally more useful on small datasets (where neural networks typically overfit). Small datasets of size <1000 are very common in chemistry. Although small datasets generally don’t require approximate kernel methods like the one we present in our paper, one exception to this is approximate Thompson sampling from Gaussian processes which scales cubically in the number of _test_ points (this is desirable for Bayesian optimization). Even with a small training dataset, making predictions on a large test set is expensive. In section 5.2 we show that the random features in our paper are useful for this.\n2. In instances where reliability or interpretability are important, kernel methods have a clear advantage over neural networks, which often make spurious predictions (e.g. adversarial example)\n\n> Please explain L156 more: ? This does not appear to be a continuous function and therefore not differentiable at 0\n\n$T_{DP}$ is defined in equation 6, but this equation is ill-defined at (0,0). Line 156 gives  the definition of $T_{DP}$ at this point. The choice $T_{DP}(0,0)=1$ ensures that the kernel is positive definite over all of R^n.The kernel is continuous and differentiable everywhere except at (0,0).\n"}, {"review_id": "mcxOErXamn", "reviewer": "Reviewer_jyA1", "summary": "The paper investigates two low-rank approximations for the tanimoto similarity in the context of molecules. The two approximations are both random features newly proposed in the paper using random hashes and sketches. As theory, the paper contributes a proof of optimality for a certain way of constructing one of the approximations and an error bound for the operator norm of the other. They perform four experiments answering the immediate questions about the random features and investigating the claimed advantages of the methods.\n\nThank you for any time you take to answer my questions!", "questions": "1. Has the trick with random hashes for constructing random features based on an alphabet of random vectors and indexing that with the hashes really never been done before?\n2. (Only if you can find the time to answer) Regarding data oblivious sketches: What does it mean for a sketch to be data dependent? Just that it for instance makes assumptions on the sparsity of the kernel matrix? What would be an example of a data-dependent sketch?\n3. What is Attentive FP? A half-sentence explanation in the paper might be useful as well.\n4. Corollary 4.2 and Theorem 4.1 imply that TDP is an extension of set-valued Tanimoto similarity to real vectors. Is there any inherent value to that or why is the Tanimoto similarity put on such a high pedestal?", "limitations": "The paper seems very transparent with weaknesses and does not shy away from e.g. showing results of better-performing methods on certain tasks. They point out limitations in a dedicated section. These are not just conjured up little limitations or remarks about possible small extensions as we see them all too often in ML papers. Instead, they are significant such as computation difficulties that result in speed impairments. The authors could have easily hidden this to make their method seem better but they do not. Thank you!\n\nMaybe the paper should discuss the limitations compared to deep learned fingerprints more and e.g. the weakness compared to attentive FP instead of just showing the results.", "rating": 7, "confidence": 2, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "1. They propose a scheme to use random hashes to index an alphabet of random vectors to construct random features as a low rank approximation of the tanimoto similarity. It has several desirable properties in terms of their variance and they prove what the optimal choice is for the alphabet of random vectors.\n2. They prove that a previously known kernel is for binary inputs is also a kernel for real vectors. For real vectors their first method of constructing random features does not work and they propose another novel method for constructing random features thus arriving at random features with desirable properties like differentiability which allows for using them in e.g., bayesian optimization. For this they hold significant promise since they are fast to compute as well.\n3. The paper provides an error bound on the operator norm of the constructed T_DP random features which apparently is important for sketches to work well according to a citation that I did not check.\n4. All claimed advantages are investigated in 3 experiments and the authors do not shy away from showing better performing methods for some tasks and the weaknesses of the random features/kernel methods in general. The experiments are insightful and well chosen.", "weaknesses": "1. Main concern: How relevant are hand-crafted kernels, and what are the advantages compared to deep learning embeddings? It seems that also in the DOCKSTRING experiments, deep-learned features fare better. Are there any efficiency and scalability advantages? If the answer is that there are no real advantages, I would see it as a large drawback even if we can say that, in practice hand crafted kernels are still often commonly used. If the answer is otherwise, I think it would be important to discuss this in detail.\n2. A more extensive investigation of deep learned fingerprints would be valuable (maybe I am wrong, and you can explain why not): There already is an investigation of deep learned embeddings in terms of Attentive FP in Table 1, but how do deep learning FPs compare in the other experiments, could the f(lambda) computation be extended to deep learned fingerprints? \n3. A more extensive investigation of deep learned fingerprints would be valuable (maybe I am wrong, and you can explain why not): If other embeddings from SSL methods, such as those from “Self-supervised Graph-level Representation Learning with Local and Global Structure” would be tried as well that would be interesting or the pretrained features from “Uni-Mol: A Universal 3D Molecular Representation Learning Framework”.\n4. Are the baselines in the PMO benchmark actually strong, or are they just easily implemented baselines, and that is why they already were in the benchmark, and one should instead, e.g., compare against other kernels?\n5. How do we know that the DOCKSTRING scores reflect/have the same ranking as actual binding affinities, and is the DockString task actually relevant as a real-world benchmark? I think you are making a mistake in saying they are binding affinity tasks - am I wrong in that?\n6. Why can we assume that T_MM without random features is the best and should be a baseline to compare to in, e.g., Figure 1? Should a comparison with other kernels not also be considered?\n7. Is it fair to say that we now have justification for using T_DP for nonbinary vectors through your Theorem 4.1, but no new capability is unlocked in the sense of the kernel being known and used previously?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your kind review and very positive rating. We will try to respond to your questions / concerns below.\n\n> How relevant are hand-crafted kernels?\n\nAlthough hand-crafted kernels are less relevant than a decade ago, they are still very relevant in applications with small datasets (where neural networks overfit) or where interpretability/robustness are important (because their behavior can be theoretically characterized). An example where both of these considerations apply is Bayesian optimization of molecules, wherein predictions and uncertainty estimates of a model are used to choose promising molecules for experiments. In contrast, the conditions of the DOCKSTRING experiment in section 5.3 are very favorable to neural networks (large dataset, only assessing mean prediction accuracy). It is not surprising that the deep learning approach did better here.\n\nWe think that “hand-crafted kernels” vs “deep learning” is perhaps the wrong dichotomy, because most kernels (including the ones described in our paper) could be applied to learned fingerprints from neural networks instead of hand-engineered features. We did not explore this angle in our paper since our focus was scalable estimation of these kernels rather than choosing the best features for these kernels. We used fingerprints because this is the most common molecular feature used with Tanimoto similarity. \n\n> A more extensive investigation of deep learned fingerprints would be valuable (maybe I am wrong, and you can explain why not)\n\nWe aimed to compare with deep learning methods in situations where it was appropriate. In particular:\n- The experiment in section 5.1 aims to reconstruct matrices of Tanimoto coefficients using random features. The key claim of our work is that the random features we propose can do this, so the goal in this section was to demonstrate this experimentally. Deep learning fingerprints are not trained to do this, and there is no reason to expect that the inner product of deep learning fingerprints would approximate a Tanimoto coefficient (in fact the inner product of deep learning features is usually unbounded). Therefore we did not think this was a reasonable comparison. We think adding it would distract from the main point of the section, which is experimentally verifying the claims of sections 3-4.\n- Deep learning methods were used by Gao et al (2022) as part of the PMO benchmark, and therefore _we did (implicitly) compare against them_ (and Tanimoto GP performed better). We can change the text to emphasize that this comparison was made. The Thompson sampling experiment is a unique capability of methods which make correlated predictions (like GPs), which deep learning methods generally do not. Therefore we could not perform a comparison here.\n- In section 5.3 we compared a scalable kernel method operating on hand-engineered features with an end-to-end deep learning system. We thought this was a sensible comparison because both of these approaches might be used by practitioners. We chose not to include more deep learning methods because it was already clear that the kernel method performs worse (which was the expected outcome).\n\nTherefore we think that it would be difficult to include a thorough investigation of deep learning fingerprints without diluting the overall narrative or significantly expanding the scope of the work. We are however excited to explore applications of scalable Tanimoto kernels in future work, which will involve comparing to deep learning methods on other tasks.\n\n> Are the baselines in the PMO benchmark actually strong?\n\nWe think that they are reasonably strong: they included a wide variety of state-of-the-art methods and performed hyperparameter tuning.\n\n> How do we know that the DOCKSTRING scores reflect actual binding affinities, and is the DockString task actually relevant as a real-world benchmark?\n\nWe don’t know exactly, but DOCKSTRING targets were chosen to have a reasonably strong correlation between docking score and empirically measured binding affinity. However, we see this package as simply a tool to benchmark our methods against moderately complex molecular properties, and it is one of the few large and publicly-available labeled datasets of small molecules.\n\n> Should a comparison with other kernels not also be considered [in Figure 1]?\n\nThe aim of Figure 1 is to evaluate the reconstruction error of various random features approximations of the Tanimoto kernel matrix. This figure does not evaluate the statistical performance of kernel methods for a given molecular prediction task. We would not argue that T_MM should be the definitive baseline in such a task, and we include a limited experimental evaluation of this sort in Sections 5.3.\n\n> Is it fair to say that we now have justification for using T_DP for nonbinary vectors through your Theorem 4.1, but no new capability is unlocked in the sense of the kernel being known and used previously?\n\nI think a new capability is unlocked because even though it was used, it doesn’t look like it was used on non-binary vectors.\n\n> Has the trick with random hashes for constructing random features based on an alphabet of random vectors and indexing that with the hashes really never been done before?\n\nAs far as we are aware, no. This is indeed a very natural idea, but we were unable to find a reference to it in the literature.\n\n> What does it mean for a sketch to be data dependent?\n\nA data dependent sketch is one which applies to a specific dataset only. Two examples are Nystrom approximations or inducing point methods.\n\n> What is Attentive FP?\n\nIt is a type of graph neural network.\n\n> [...] why is the Tanimoto similarity put on such a high pedestal?\n\nEven though the Tanimoto kernel is a predominant metric in cheminformatics, we see no reason why it should be put on a pedestal. It is simply the kernel which our work focuses on approximating.\n\nPlease let us know if this answers your questions and addresses your concerns."}, {"review_id": "z0rGm9PDO4", "reviewer": "Reviewer_pMsa", "summary": "Random feature map can overcome the limitation of kernel methods, and reduces its computational cost to linear in the number of training data points. Recent advances in random feature map broaden its application to not only to stationary kernels, but also to non-stationary kernels such as product kernels. The authors propose two kinds of random features for Tanimoto coefficient such that kernel methods based on Tanimoto coefficient can be applied to large datasets.", "questions": "“2. Uncertainty on the correctness of the error bound”\nThis reviewer finds this a serious point. What is the reason behind this ? Bound not tight enough, or bug in experiments ?\nPlease describe what the authors can do to resolve it.", "limitations": "The authors clearly mentioned the limitations of their work.\n", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. Proposal of two random features\n2. Theoretical analysis on the variance and approximation error for the proposed random features.\n3. Demonstrated efficiency through computational experiments.\n4. Proposal of positive semidefinite Tanimoto product kernel for real-valued vectors.", "weaknesses": "As mentioned by the authors in Discussion,\n1. Unexpectedly slow computation of T_MM \n2. Uncertainty on the correctness of the error bound", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for reviewing our paper. We will try to respond to your questions / concerns below.\n\n> Uncertainty on the correctness of the error bound” This reviewer finds this a serious point. What is the reason behind this ? Bound not tight enough, or bug in experiments ? Please describe what the authors can do to resolve it.\n\nThe reviewer may have misinterpreted our comment: we are not uncertain about the correctness of the bound. What we intended to say is:\n\n1. In practice, the error could be much lower than the bound suggests (the bound might not be tight).\n2. The bound used a specific set of hyperparameters. However, it is possible that tuning the hyperparameters differently could result in a lower error in practice.\n\nTo be clear: we are not uncertain about the correctness of the error bound and we do not believe there is a bug in any of our experiments.\n\n> Unexpectedly slow computation of T_MM \n\nThis is not fundamental, it was just because we could not find an easy way to vectorize its computation in python. A good C++ implementation could be fast; we just didn’t do this. The time complexity still scales linearly with N and D.\n\nPlease let us know whether this addresses your concerns and whether you have any outstanding questions or concerns."}]

Author rebuttal:
We would like to thank all reviewers for their constructive comments. We have left a response to each reviewer answering their questions, but would like to provide a summary here.

First, we appreciate that many reviewers complimented various aspects of our work including:

- Theoretical analysis and technical contribution (pMsa, jyA1, A9mq, b6mP)
- Clarity (b6mP, A9mq, jyA1)
- Honesty about limitations (A9mq, jyA1)

The main concerns raised by the reviewers were:
- Notation (stq5): there was some ambiguity in subscripts and instances where we defined probability distributions implicitly. _We will change this in the next version of the manuscript._
- Status of theorem 4.7 (b6mP, pMsa): _this is a novel theorem, and we are certain about its correctness_. Cohen et al (2015) pioneered the analysis technique and set a reference for what a “good” error bound is in these types of problems, and our comment in the discussion was about the bound not being tight or optimal, _not_ a statement that we are unsure about its correctness.
- Relationship to deep learning methods (A9mq, jyA1): we believe and clearly acknowledge in the paper that there are circumstances where deep learning methods perform much better than kernel methods: for example on problems with large datasets like in section 5.3. However, kernel methods have distinct advantages on small datasets or in instances where robustness / interpretability are important. The Thompson sampling experiment in 5.2 is a typical setting where we envision our method will be used: even with a small number of training data points, sampling scales poorly with the number of _evaluation_ points but can be greatly accelerated using our random features (at the cost of some approximation error).

  Furthermore, because the kernels studied in this work are differentiable and valid for continuous inputs, they could be combined with learned features in the future, for example via “Deep Kernel Learning” (Wilson et al 2016). We leave this to future work.

We hope we have adequately addressed the concerns of all reviewers and appreciate the time that you have put into the reviewing process. We would be happy to answer any further questions.


Meta-review: The paper focuses on computing efficient approximations to the Tanimoto coefficient, an important similarity measure in many applications, especially in applications on molecules. The paper puts forward two different random feature based approximations of the Tanimoto kernel, and analyze their theoretical properties. The methods are tested empirically in extensive experiments that showed the practical relevance of the method. The reviewers found the paper technically solid. The minor concerns of the reviewers were adequately addressed in the rebuttals and the following discussion.
Decision: Accept (poster)