Proceedings of the Twenty-Ninth AAAI Conference on Artificial Intelligence

An SVD and Derivative Kernel Approachto Learning from Geometric Data
Eric Wong

J. Zico Kolter

School of Computer Science
Carnegie Mellon University
ericwong0@cmu.edu

School of Computer Science
Carnegie Mellon University
zkolter@cs.cmu.edu

Abstract

the energy of new configurations, or within a Bayesian optimization procedure for finding the minimum energy configuration for a molecule.
A main contribution of this paper is the development of a
(improper) kernel approach to dealing with this data, based
upon the singular value decomposition. Geometric learning
problems of the type we consider typically exhibit certain
invariances, such as translational and rotational invariance(i.e., rotating or translating all the atoms on a molecule do
not actually change any qualitative information about themolecule, since it remains the same molecule), and we derive a kernel formulation that respects these invariances. Furthermore, because many physical simulations that serve as
the input to such geometric methods (notably those based
upon adjoint methods) can also provide derivatives of theoutput quantity with respect to all point locations, we derive an approach to including derivative information into
our kernel learning. Finally, we highlight the method on examples from the molecular energy prediction task, showing
good performance in modeling previously unseen configurations. Incorporating the method into a Bayesian optimization
framework, we show significant improvement over existingstate of the art in molecular energy optimization.

Motivated by problems such as molecular energy prediction, we derive an (improper) kernel between geometric inputs, that is able to capture the relevant rotational and translation invariances in geometric data.
Since many physical simulations based upon geometric data produce derivatives of the output quantity withrespect to the input positions, we derive an approachthat incorporates derivative information into our kernel learning. We further show how to exploit the low
rank structure of the resulting kernel matrices to speed
up learning. Finally, we evaluated the method in thecontext of molecular energy prediction, showing good
performance for modeling previously unseen molecularconfigurations. Integrating the approach into a Bayesian
optimization, we show substantial improvement overthe state of the art in molecular energy optimization.

Introduction
This paper focuses on learning from geometric data, where
each input to the learning problem consists of a set ofpoints, typically in two or three dimensional space. While
many physical problems take this form, we focus here on
the task of molecular energy prediction from input data
describing the 3D configuration of atoms. This problem
has been studied before in the machine learning literature,
though in a slightly different context (Rupp et al. 2012;
Montavon et al. 2012), and has numerous applications ranging from energy materials to drug design. The basic task hereis, given a set of atoms and their respective 3D positions,to predict the energy of the resulting configuration; sincemolecules in nature will typically settle in minimum energy
configurations, in finding minima of these functions we canassist in the design of new materials. Physically-based simulation to determine these configurations, notably methods
based upon density functional theory (DFT) approaches, area very mature technology, but are a very computationally intensive process. The canonical problem we consider in thispaper, therefore, is the task of learning a faster surrogate
function for this expensive process, using machine learningtechniques. This can be used simply as a means of predicting

Related work
This work builds upon a number of different approaches
in kernel methods, computer vision, Bayesian optimization,
and applications to material science. Kernel methods havelong been applied to spatial data in general (see, e.g. (Pozdnoukhov 2009)), which shares a high-level similarity to thework we present here in that both methods use 2D or 3D coordinates as input to a learning system. However, the work
presented here differs substantially from most work in spatial learning, as the goal there is typically to predict a quantity at a single spatial point (and thus the notions of translational and rotational invariance are specifically inapplicable). In contrast, the work here looks at developing a kernelthat measures similarity between configurations of multiple
points, while respecting certain invariances in the data.
The basic singular value decomposition (SVD) approachthat we use, for finding optimal translation and rotations between molecules, has been previously studied for some time
(Horn, Hilden, and Negahdaripour 1988), and has been usedin the vision community as a subroutine for fast geometric

Copyright c 2015, Association for the Advancement of Artificial
Intelligence (www.aaai.org). All rights reserved.

1889

we define the distance between two molecules X ∈ R3×nand Z ∈ R3×n as the minimum distance over all rotations
and translations:d(X, Z) =min||X − (RZ + t1T )||2F (1)transforms (Lu, Hager, and Mjolsness 2000). However, themain contributions over this basic approach in the current
work are: 1) the introduction of a kernel based upon thisSVD-based distance metric, and 2) the derivation of functional derivatives of this kernel in a compact form, which
is of particular importance for the domain of molecular energy prediction. This last item also exploits work in kernel
methods for incorporating derivative information in Gaussian processes (Solak et al. 2003); we use exactly this approach in order to incorporate observations of the molecular forces (derivatives) into our model, though the challenging aspect here is computing the proper kernel derivatives,
which is one of the primary contributions of this paper.
Our work also uses recent approaches to Bayesian optimization (Brochu, Cora, and De Freitas 2010). From thecore algorithmic perspective, we are mainly using existing
approaches from the literature, specifically the lower confidence bound approach (Cox and John 1997) that is widelyused in the machine learning community. However, most
Bayesian optimization approaches have been applied using
relatively simple forms of kernels such as the exponential
or Matern kernels (Snoek, Larochelle, and Adams 2012). In
contrast, the use of the geometric kernel that we derive here
lets us use Bayesian optimization as a drop in replacement
for existing molecular energy optimizers, based upon traditional optimization approaches such as gradient descent or
LBFGS. In this context, the work here can be viewed as integrating Bayesian optimization in this approach to numericalscientific computing.
Finally, in the domain of material science, our approachbuilds upon a number of past works. Both neural networks
and Gaussian processes have been used to predict the potential energy of a single molecule (Behler and Parrinello 2007;
Bartók et al. 2010), but these used either black box or noninvariant kernel approaches, and do not integrate derivative
observations. More recent approaches learn across the entire chemical compound space by using data from multiple
molecules. These approaches represent molecules by theirso-called Coulomb matrix and their corresponding vectorsof eigenvalues (Rupp et al. 2012; Montavon et al. 2013).
The major difference presented here is the use of the singularvalue decomposition to define the distance between two inputs and the incorporation of the available molecular energy
derivative information.

R∈R3×3 ,t∈R3 ,RT R=I3

Although this is a non-convex problem, it can be solved efficiently using the singular value decomposition. This is awell-known property (see, e.g. (Horn, Hilden, and Negahdaripour 1988)), but we include a brief derivation for completeness.
Theorem 1. The solution to the minimization in (1) is givenbyt? = (X − R? Z)1/n, R? = V U T(2)where U SV T = XBB T Z T is a singular value decomposition and B is the centering matrix defined as B =In −11T /n. Furthermore, the distance takes the closed formd(X, Z) = ||XB||2F + ||Y B||2F − 2 Tr S
(3)where S is the diagonal matrix of singular values as above.
Proof. Taking the gradient of (1) with respect to t
∇t kX − RZ − t1T k2F = (X − RZ − t1T )1(4)
gives the the first part of the solutiont? = (X − RZ)1/n.
(5)
Substituting this back into the objective transforms (1) into
the equivalent optimization problemd(X, Z) = min ||XB − RZB||2F(6)
RT R=I3

and note that||XB − RZB||2F= kXBk2F + kRZBk2F − 2 Tr B T Z T RT XB(7)= kXBk2F + kZBk2F − 2 Tr RT XBB T Z T .
We next take the singular value decomposition U SV T =XBB T Z T
Tr RT U SV T = Tr V T RT U S = Tr R̃S(8)T
where R̃ R̃ = I. Since S is diagonal, this problem is optimized with R̃? = I, so that R? = U V T . The final form of
the distance follows from the fact that
Tr V U T U SV T = Tr U T U SV T V = Tr S.
(9)
We use this distance in a corresponding exponential kernel given by


1
2
2
K(X, Z) = exp − γ(||XB||F + ||ZB||F − 2 Tr S) .
2
(10)Note that this is not a proper (semidefinite) kernel, due to thefact that this distance function is not a true metric (it does notobey the triangle inequality); thus the kernel matrix above
may have negative eigenvalues. In this case, we can simply
treat the kernel as a set of features describing the input (i.e,using the matrix K T K to compute predictions); we can also
exploit the fact that, for large enough γ, the kernel will be
guaranteed to be positive definite, as all non-diagonal entries will fall off exponentially. In practice, we observe that
this is rarely a problem, and all values of γ chosen by crossvalidation, for instance, results in positive definite kernels.

An SVD kernel for geometric data
In this section we present the main algorithmic contributionof this paper, an improper kernel for geometric input data.
In the sequel, we will use this kernel in the context of Gaussian process estimation and Bayesian optimization, but the
key contribution here relates to the derivation of the kernel
itself and its corresponding derivatives that can enable us toincorporate derivative information into the kernel.
Let X ∈ R3×n represent molecules of n atoms, where
each column of the matrix denotes the 3D coordinate of thecorresponding atom. Since translating or rotating a molecule
results in the an unchanged molecule, a measure of similarity between two molecules should be invariant to these properties. To construct a Gaussian process model for such data,

1890

Kernel derivatives

with ⊗ as the Kronecker product, B̄ = (B ⊗ I3 ), and


0 0
0
0 0 0
0
0 0
1
−1
 0 ŝ12
0
0 0
0
0 0 
s12


−1
1
 0 0
0
0
0
0
0 
s13
s13


1
 0 −1
0
0 0
0
0 0 
s12
s12


0
0 0 0
0
0 0 
A=
.
 0 0
1
−1

 0 0
0
0
0
0
0
s
s


23
23
1
−1
 0 0
0 0 0
0 0 


s13
s13
1
 0 0
0
0 0 s−1
0
0 
s23
23
0 0
0
0 0 0
0
0 0
(15)
where sij = Sii + Sjj .

A common theme in many learning methods that use geometric data, specifically those where the geometric data isitself generated through some simulation process, also produce derivatives of the output quantity with respect to all the
input coordinates. For instance, the general class of methods known as adjoint methods produce these derivatives using approximately twice as much computation as it takes tocompute the output originally. Such methods have therefore
become extremely widely used in simulation and optimization, and we want to enable our kernel approach to use suchinformation. As shown in previous work (Solak et al. 2003),
a Gaussian process can incorporate derivative observations
by defining an extended kernel matrix that also models thecovariance between the inputs and their derivatives, and between the derivatives themselves. In particular, since the covariance is a bilinear operator:

Proof. To derive these expressions, we determine the derivatives of the singular value decomposition using a methodology similar to the work of (Papadopoulo and Lourakis
2000), beginning with the following identity:
U SV T = XBBZ T(16)
We take the derivative of both sides with respect to someT
input coordinate Xi,j , multiply on the left and right
by

U



Cov y (s) , y (t) = K X (s) , X (t)!
∂y (s) (t)∂K(X (s) , X (t) )
Cov
,
y
=
(s)
(s)
∂Xij
∂Xij
!
∂y (s) ∂y (t)∂ 2 K(X (s) , X (t) )
Cov
,
=
.
(s)
(t)
(s)
(t)
∂Xij ∂Xk`∂Xij ∂Xk`

and V respectively, and take the trace. Since U T
 T∂V∂Xij V are antisymmetric, we get:


∂S
tr
= tr(J ij BZ T V U T )
∂Xij

(11)







∂K(X,Z)∂X11∂K(X,Z)∂Z11∂ 2 K(X,Z)∂X11 ∂Z11

..
.

···
..
.

∂K(X,Z)∂X3n∂ 2 K(X,Z)∂X3n ∂Z11···

K(X, Z)..
.

···∂K(X,Z)∂Z3n∂ 2 K(X,Z)∂X11 ∂Z3n..
.

∂ 2 K(X,Z)∂X3n ∂Z3n(12)

Theorem 2. For the rotationally invariant distance kernelin (10), the block of derivatives in (12) has the following
form:


1
K̄11
K̄(X, Z) = exp − γd(X, Z)
K̄21
2

K̄12
K̄22


(13)

where
K̄11 =1
K̄12 = − γ vec(Z − V U T X)T B̄ T
K̄21 = − γ B̄ vec(X − U V T Z)
K̄22 =γ B̄ (In ⊗ U V T )(17)∂d(X, Z)= 2(bj ⊗ ei )T vec(X − U V T Z)
(18)
∂Xij
By varying i and j, we can thus construct a row of all thefirst derivatives of X:∂d(X, Z)= 2(B ⊗ I3 )T vec(X − U V T Z)
(19)
∂ vec(X)We substitute into the derivative of the kernel to get:∂K(X, Z)1
= exp − γd(X, Z) K̄21
(20)
∂ vec(X)2
An almost identical process results in the respective expression for the derivative with respect to the second input:
∂K(X, Z)1
= exp − γd(X, Z) K̄12(21)
∂ vec(Z)2
For the second derivatives, we begin by directly taking thesecond derivative of the kernel:∂ 2 K(X, Z)= K(X, Z) (α − β)(22)∂Zkl ∂Xij
where



∂d(X, Z)1 2 ∂d(X, Z)
α= γ
4
∂Zkl
∂Xij
(23)
2
1 ∂ d(X, Z)
β= γ2 ∂Zkl ∂Xij








and

with (J ij )kl = δik δjl . Let bj be the jth column of B. Substituting the value in equation 17 for the derivative of the singular values, the first derivative of the distance metric withrespect to Xij can be expressed as the following:

In the remainder of this section, we will derive a closed form
expression for these derivatives for the SVD kernel. For convenience, it is useful to combine all the kernel derivatives
between two molecules X, Z ∈ R3×n in the block form:


∂U
∂Xij

(14)


−(Z T V ⊗ U )A(U T X ⊗ V T ) B̄ T
+ K̄21 K̄12

1891

By varying the indices i, j, k, l we construct a block of all
the second derivatives of the kernel, resulting in the desired
form for K̄22 This completes the expression for the second
derivative of the kernel with respect to Xij and Zkl :The first part is expressed from the first derivative results:
α = [K̄21 ](ij) [K̄12 ](kl)(24)

where the subscripts (ij), (kl) refer to the indices of thederivative with respect to Xij , Zkl . For the second part, taking the second derivative of the distance metric shows that
we need the second derivative of the singular values:∂2S∂ 2 d(X, Z)= −2 tr(25)∂Zkl ∂Xij∂Zkl ∂Xij∂ 2 K(X, Z)= K(X, Z)K̄22∂ vec(X)∂ vec(Z)T

Low rank structure
One potential issue with the methodology proposed here isthat the size of the kernel matrix, especially with derivatives, can grow very quickly: for m molecules with n atoms
each, the full K matrix will be m(3n + 1) × m(3n + 1).
Even though the goal of our approach is to use relatively
few molecules for training, this can quickly grow infeasible
for large n: even multiplying a vector by the full K matrix
will naively have cost O(m2 n2 ).
However, there is also substantial structure in the full K
matrix that lets us reduce this time substantially. Namely,
each 22 block of the collection of derivatives is the sum of arank-one matrix, and two matrices

To get this quantity, we first take the second derivative of
the identity U SV T = XBB T Z T and multiply on the left
and right by U T and V respectively. Taking the trace of theresulting equation, we get



∂2S
tr
= tr U T J ij BJ lk V
∂Zkl ∂Xij



ij
kl
kl
− tr Ωij
+
Ω
S
Ω
+
Ω
U
V
U
V
(26)
where
T
Ωij
U =U

∂V
Ωij
V =

∂U
∂Xij

T
Ωkl
U =U

T

∂Xij

V

Ωkl
V =

∂U
∂Zkl

∂V
V
∂Zkl

B̄(In ⊗ U V T )B̄ T = B ⊗ U V T
(27)

(32)

and
B̄(Z T V ⊗ U )A(U T X ⊗ V T )B̄ T= (BZ T V ⊗ U )A(BU T X ⊗ V T )ij

and (J )kl = δik δjl . For the latter term, we can solve for an
ij
explicit form of Ωij
U and ΩV by solving the system of equations given by the off-diagonal entries of the first derivativeof equation 16. Using these, we can compute

ij
ij 
ij
ij
W12
−W21
W13
−W31
0
s1 +s2s1 +s3
 W ij −W ij
ij
ij 
ij
W23
−W32

 − 12 21
Ωij
+
Ω
=
0
U
V


s1 +s2s2 +s3
ij
ij
ij
ij
W −W
W −W
0
− s131 +s3 31 − s232 +s3 32(28)
where W ij = U T J ij BZ T V . We can get an identical exklkl
= U T XBJ lk V in placepression for Ωkl
U + ΩV , using W
ij
of W . Substituting these into the trace in equation 26 andusing A, we derive the following expression:
tr

(31)

(33)

Multiplying a vector by this product just involves centeringand multiplication by 3x3 matrices, which together are O(n)operators. Thus, the cost of a matrix-vector product can easily be brought down to O(m2 n). Combined with conjugate
gradient approaches, these methods can be used to bring thecost of the SVD kernel with derivatives to approximately thesame level as the standard SVD kernel.

Application to Molecular Energy Prediction
We now come to the main applied focus of this paper: using the above kernel approach to model the free energyof molecules. As described above, the basic task is, given
a configuration of atoms in 3D space (the 3D coordinateof each atom along with it’s type), we want to predict themolecular energy of the atom. This is a task with numerous
applications in energy, physics, and biology, and such methods are widely studied in chemical engineering (see, e.g.
(Sholl and Steckel 2011)). Existing methods, such as those
based upon density functional theory (DFT), are capable ofcomputing these quantities to reasonably high accuracy, but
the methods are computationally intensive (they amount tosolving a large partial differential equation to compute theenergy). The goal of our overall approach, similar to thework of (Rupp et al. 2012), is to use a machine learning approach to “replicate” the results of an expensive simulation
procedure with a much faster surrogate function. Unlike thispast work, however, we explicitly incorporate energy derivatives (forces), which are automatically generated by DFT
methods when computing the energy. We also focus upon
the task of optimizing the molecular configuration to find theconfiguration with minimum energy; this is a canonical task



 

ij
kl
Ωij
S Ωkl
U + ΩVU + ΩV


1
ij
kl
ij
kl
ij
kl
ij
kl
=
−W12
W12
+ W12
W21
− W21
W21
+ W21
W12
s1 + s2


1
ij
kl
ij
kl
ij
kl
ij
kl
+
−W13
W13
+ W13
W31
− W31
W31
+ W31
W13
s1 + s3


1
ij
kl
ij
kl
ij
kl
ij
kl
−W23
W23
+ W23
W32
− W32
W32
+ W32
W23
+
s2 + s3TT= vec W ij
A vec W kl
(29)

We substitute into equation 26 to get an expression for the
trace of the second derivative of the singular values. Afterrearranging, we get the final form of β:β = − γ(bj ⊗ ei )T (I ⊗ U V T )(bl ⊗ ek )+ γ(bj ⊗ ei )T (Z T V ⊗ U )A(U T X ⊗ V T )(bl ⊗ ek )
(30)

1892

KernelMean predictor
Eigenspectrum representation
SVD kernel
SVD kernel with derivatives

RMSE
2.41540.42840.71820.6239KernelMean predictor
Eigenspectrum representation
SVD kernel
SVD kernel with derivatives

Table 1: Prediction error on water

Table 2: Prediction error on glycerol
5

in molecular modeling, and we show Bayesian optimization
methods, based upon our approach, can substantially outperform existing state-of-the-art solvers used in the chemicalengineering community, such as LBFGS.
To begin, we will evaluate the performance of our approach just on predicting the energy of previously unseen
configurations of a molecule. In this section and the next,
all computations were carried out using the GPAW numerical code (Mortensen, Hansen, and Jacobsen 2005), a gridbased implementation of the DFT calculator. We also use thePython Atomic Simulation Environment (ASE) (Bahn andJacobsen 2002) to set up the computations and later to perform the molecular optimization. As mentioned above obtaining the potential energies is extremely computationallyexpensive using these tools, so the algorithm must be able togeneralize from a comparatively small subset of the molecular space. To get accurate energy results, the time to perform
a DFT calculation can take anywhere from a few hours for a
simple molecule to days for complex compounds.
We applied a Gaussian process using the SVD kernel tomodel the energies of water and glycerol molecules. For
both molecules, we generated 100 data points by adding
noise to the original coordinates. Denote X as the set of all
our inputs (3D positions of the atoms) and y the vector ofcorresponding energies. We form a prediction on new examples X 0 by the standard Gaussian process equations

4.5

σ(X 0 ) = K(X 0 , X 0 ) −K(X , X )(K(X , X ) + λI)

mean squared error

SVD kernel

3.5
3
2.5
2
1.5
1
0.5
0
10

20

30

40
50
# samples

60

70

80

Figure 1: Learning curves for the SVD kernel on glycerolwith and without derivatives. The presence of kernel derivatives results in faster convergence to a lower error.
are able to accurately model the system, and all other methods perform virtually identical to simple mean predictions.
This is also born out in the learning curves in Figure 1. The
SVD kernel with derivatives is able to achieve low error using approximately 40 simulations, whereas the SVD-only
approach (and similar approaches), never achieve good performance (the fluctuations in the learning performance for
the SVD kernel in Figure 1 are a product of random folds,
and are all higher than simply predicting the mean energy inthe data set).

(34)
−1

SVD kernel with derivatives

4µ(X 0 ) = k(X 0 , X )(K(X , X ) + λI)−1 y0RMSE0.74180.73690.76090.22760

K(X , X )Bayesian optimization

The kernel hyperparameters, namely exponential parameterγ and regularization parameter λ, were chosen by a gridsearch with an inner 4-fold cross validation and optimized
for the root mean squared error. We introduce an additional
regularization parameter λderiv to account for the differencein magnitude between the energies and derivatives.
In the case of the simple molecule water, the SVD kernel is able to capture much of the structure of the potential energy function, with the derivative information giving asmall improvement in performance. Since water only has 3
atoms, the SVD kernel is already able to model the molecule
quite well without the derivatives. Furthermore, the eigenvalue method from (Rupp et al. 2012) performs marginallybetter in this application, likely because it exploits additional
information about the energy of single atoms (such information could be included in our setting as well, though we do
not pursue this approach here).
In contrast, on the larger glycerol molecule with 14 atoms,
none of the methods except the SVD kernel with derivatives

Finally, the ultimate goal of a fast DFT approximation is
not just to predict the energy, but to find the minimum energy configuration. Since the target function has a high cost,
we take a Bayesian optimization approach to minimize thepotential energy of a molecule. We use a Gaussian processmodel with the SVD derivative kernel as an approximationof the target energy function, and we use the lower confidence bound (LCB) as our acquisition function to sample
new points:
LCB(X ) = µ(X ) − κσ(X )(35)

where µ and σ are the mean and variance from the Gaussian process. We optimized this surrogate function itself using LBFGS, though importantly, this optimization procedure
is very fast, as we can compute the GP predictions very
quickly.
For the testing environment, we optimize three molecules
whose atoms have had noise added to their initial coordinates. We treat the DFT energy calculation as the function

1893

SVD kernel with derivatives
65

0

10

0

5

10
# iterations

15

60
55
50
45
40

20

L−BFGS
65
Minimum energy found

SVD derivative kernel
L−BFGS

Minimum energy foundlog average difference to optimallog average difference to optimallog average difference to optimal

H2O

0

10
# iterations

20

60
55
50
45
40

0

10
# iterations

20

C5H12

2

10

Figure 3: A plot of all 40 individual optimizations for eachalgorithm on water. The SVD kernel converges rapidly inalmost all cases regardless of starting position. The LBFGS
method depends on the starting configuration and has greater
variance in convergence rate and final energy value.

SVD derivative kernel
L−BFGS
1

10

0

10

kernel with derivatives uses fewer iterations to get closer tothe optimal value, converging more rapidly than LBFGS in
Figure 2 across multiple molecules. Another benefit that the
SVD method has is its ability to converge quickly to thesame optimal value regardless of the starting position. For
example, after 5 iterations, nearly all optimizations of waterusing the SVD method are close to the optimal value as seenin Figure 3. In contrast, the LBFGS method depends significantly on the starting configuration. Some LBFGS optimizations converge slower than others and reach a variety of final
energy values.

−1

10

0

5

10
# iterations

15

20

C3H8O3
SVD derivative kernel
L−BFGS
0

10

0

5

10
# iterations

15

Conclusion

20

When using physical systems or geometric data, there arecertain properties of translational and rotational invarianceand an availability of derivative information that we can exploit to improve the performance of learning in this space. In
the case of molecular energies, predicting accurate energies
is especially important since calculating the exact energy ofa compound typically takes days, and scales badly as thecomplexity of the molecule increases.
In this paper, we have developed a kernel based on the singular value decomposition that maintains the translationaland rotational invariants expected in geometric data. We incorporated derivative data into the kernel and constructed asimplified form to compute the larger covariance matrices.
Our results show that the SVD kernel has good predictive
ability, and can outperform other methods. Practically, ourwork shows that the SVD kernel can optimize molecules
faster and more reliably than current methods.

Figure 2: The average minimum found over all 40 runs periteration for various molecules, which corresponds to thenumber of function samples. Results for the SVD kernel
converge faster and closer than LBFGS to the optimal value.

to be minimized, and provide the negative forces as the gradient function. Each algorithm begins with only the initialstarting set of coordinates, and the optimization is repeatedfor 40 randomly perturbed starting configurations.
We compare our optimization method against the LBFGS
quasi-Newton method (directly performing LBFGS using
the derivatives provided by the DFT). We use an implementation made specifically for optimizing molecular structures from the Python ASE package, which is a state-of-theart method for performing molecular structure optimization.
Because the time used by the DFT solver far exceeds any of
the computation time used either by LBFGS or our Bayesian
optimization, we report results in terms of the number of (expensive) simulation iterations required.
We see that the SVD kernel with derivatives within theBayesian optimization framework offers a significant improvement over the LBFGS method. On average, the SVD

References
Bahn, S. R., and Jacobsen, K. W. 2002. An object-oriented
scripting interface to a legacy electronic structure code.
Computing in Science & Engineering 4(3):56–66.
Bartók, A. P.; Payne, M. C.; Kondor, R.; and Csányi, G.
2010. Gaussian approximation potentials: The accuracy of1894tronic properties in chemical compound space. New Journalof Physics 15(9):095003.
Mortensen, J. J.; Hansen, L. B.; and Jacobsen, K. W. 2005.
Real-space grid implementation of the projector augmented
wave method. Physical Review B 71(3):035109.
Papadopoulo, T., and Lourakis, M. I. 2000. Estimating
the jacobian of the singular value decomposition: Theory
and applications. In Computer Vision-ECCV 2000. Springer.
554–570.
Pozdnoukhov, A. 2009. Machine learning for spatial environmental data: theory, applications, and software. EPFL
press.
Rupp, M.; Tkatchenko, A.; Müller, K.-R.; and von Lilienfeld, O. A. 2012. Fast and accurate modeling of molecularatomization energies with machine learning. Physical review letters 108(5):058301.
Sholl, D., and Steckel, J. A. 2011. Density functional theory:a practical introduction. John Wiley & Sons.
Snoek, J.; Larochelle, H.; and Adams, R. P. 2012. Practicalbayesian optimization of machine learning algorithms. In
Advances in Neural Information Processing Systems, 2951–
2959.
Solak, E.; Murray-Smith, R.; Leithead, W. E.; Leith, D. J.;
and Rasmussen, C. E. 2003. Derivative observations in gaussian process models of dynamic systems. In Neural Information Processing Systems. MIT Press.

quantum mechanics, without the electrons. Physical review
letters 104(13):136403.
Behler, J., and Parrinello, M. 2007. Generalized neuralnetwork representation of high-dimensional potentialenergy surfaces. Physical review letters 98(14):146401.
Brochu, E.; Cora, V. M.; and De Freitas, N. 2010. A tutorial on bayesian optimization of expensive cost functions,with application to active user modeling and hierarchical reinforcement learning. arXiv preprint arXiv:1012.2599.
Cox, D. D., and John, S. 1997. Sdo: A statistical method for
global optimization. Multidisciplinary design optimization:state of the art 315–329.
Horn, B. K.; Hilden, H. M.; and Negahdaripour, S. 1988.
Closed-form solution of absolute orientation using orthonormal matrices. JOSA A 5(7):1127–1135.
Lu, C.-P.; Hager, G. D.; and Mjolsness, E. 2000. Fast andglobally convergent pose estimation from video images. Pattern Analysis and Machine Intelligence, IEEE Transactions
on 22(6):610–622.
Montavon, G.; Hansen, K.; Fazli, S.; Rupp, M.; Biegler, F.;
Ziehe, A.; Tkatchenko, A.; Lilienfeld, A. V.; and Müller, K.R. 2012. Learning invariant representations of molecules
for atomization energy prediction. In Advances in Neural
Information Processing Systems, 440–448.
Montavon, G.; Rupp, M.; Gobre, V.; Vazquez-Mayagoitia,A.; Hansen, K.; Tkatchenko, A.; Müller, K.-R.; and vonLilienfeld, O. A. 2013. Machine learning of molecular elec-

1895