Coverage for ase / optimize / gpmin / kernel.py: 71.25%

1# fmt: off 

2 

3import numpy as np 

4import numpy.linalg as la 

5 

6 

7class Kernel(): 

8 def __init__(self): 

9 pass 

10 

11 def set_params(self, params): 

12 pass 

13 

14 def kernel(self, x1, x2): 

15 """Kernel function to be fed to the Kernel matrix""" 

16 

17 def K(self, X1, X2): 

18 """Compute the kernel matrix """ 

19 return np.block([[self.kernel(x1, x2) for x2 in X2] for x1 in X1]) 

20 

21 

22class SE_kernel(Kernel): 

23 """Squared exponential kernel without derivatives""" 

24 

25 def __init__(self): 

26 Kernel.__init__(self) 

27 

28 def set_params(self, params): 

29 """Set the parameters of the squared exponential kernel. 

30 

31 Parameters 

32 ---------- 

33 

34 params: [weight, l] Parameters of the kernel: 

35 weight: prefactor of the exponential 

36 l : scale of the kernel 

37 """ 

38 self.weight = params[0] 

39 self.l = params[1] 

40 

41 def squared_distance(self, x1, x2): 

42 """Returns the norm of x1-x2 using diag(l) as metric """ 

43 return np.sum((x1 - x2) * (x1 - x2)) / self.l**2 

44 

45 def kernel(self, x1, x2): 

46 """ This is the squared exponential function""" 

47 return self.weight**2 * np.exp(-0.5 * self.squared_distance(x1, x2)) 

48 

49 def dK_dweight(self, x1, x2): 

50 """Derivative of the kernel respect to the weight """ 

51 return 2 * self.weight * np.exp(-0.5 * self.squared_distance(x1, x2)) 

52 

53 def dK_dl(self, x1, x2): 

54 """Derivative of the kernel respect to the scale""" 

55 return self.kernel * la.norm(x1 - x2)**2 / self.l**3 

56 

57 

58class SquaredExponential(SE_kernel): 

59 """Squared exponential kernel with derivatives. 

60 For the formulas see Koistinen, Dagbjartsdottir, Asgeirsson, Vehtari, 

61 Jonsson. 

62 Nudged elastic band calculations accelerated with Gaussian process 

63 regression. Section 3. 

64 

65 Before making any predictions, the parameters need to be set using the 

66 method SquaredExponential.set_params(params) where the parameters are a 

67 list whose first entry is the weight (prefactor of the exponential) and 

68 the second is the scale (l). 

69 

70 Parameters 

71 ---------- 

72 

73 dimensionality: The dimensionality of the problem to optimize, typically 

74 3*N where N is the number of atoms. If dimensionality is 

75 None, it is computed when the kernel method is called. 

76 

77 Attributes 

78 ---------------- 

79 D: int. Dimensionality of the problem to optimize 

80 weight: float. Multiplicative constant to the exponenetial kernel 

81 l : float. Length scale of the squared exponential kernel 

82 

83 Relevant Methods: 

84 ---------------- 

85 set_params: Set the parameters of the Kernel, i.e. change the 

86 attributes 

87 kernel_function: Squared exponential covariance function 

88 kernel: covariance matrix between two points in the manifold. 

89 Note that the inputs are arrays of shape (D,) 

90 kernel_matrix: Kernel matrix of a data set to itself, K(X,X) 

91 Note the input is an array of shape (nsamples, D) 

92 kernel_vector Kernel matrix of a point x to a dataset X, K(x,X). 

93 

94 gradient: Gradient of K(X,X) with respect to the parameters of 

95 the kernel i.e. the hyperparameters of the Gaussian 

96 process. 

97 """ 

98 

99 def __init__(self, dimensionality=None): 

100 self.D = dimensionality 

101 SE_kernel.__init__(self) 

102 

103 def kernel_function(self, x1, x2): 

104 """ This is the squared exponential function""" 

105 return self.weight**2 * np.exp(-0.5 * self.squared_distance(x1, x2)) 

106 

107 def kernel_function_gradient(self, x1, x2): 

108 """Gradient of kernel_function respect to the second entry. 

109 x1: first data point 

110 x2: second data point 

111 """ 

112 prefactor = (x1 - x2) / self.l**2 

113 # return prefactor * self.kernel_function(x1,x2) 

114 return prefactor 

115 

116 def kernel_function_hessian(self, x1, x2): 

117 """Second derivatives matrix of the kernel function""" 

118 P = np.outer(x1 - x2, x1 - x2) / self.l**2 

119 prefactor = (np.identity(self.D) - P) / self.l**2 

120 return prefactor 

121 

122 def kernel(self, x1, x2): 

123 """Squared exponential kernel including derivatives. 

124 This function returns a D+1 x D+1 matrix, where D is the dimension of 

125 the manifold. 

126 """ 

127 K = np.identity(self.D + 1) 

128 K[0, 1:] = self.kernel_function_gradient(x1, x2) 

129 K[1:, 0] = -K[0, 1:] 

130 # K[1:,1:] = self.kernel_function_hessian(x1, x2) 

131 P = np.outer(x1 - x2, x1 - x2) / self.l**2 

132 K[1:, 1:] = (K[1:, 1:] - P) / self.l**2 

133 # return np.block([[k,j2],[j1,h]])*self.kernel_function(x1, x2) 

134 return K * self.kernel_function(x1, x2) 

135 

136 def kernel_matrix(self, X): 

137 """This is the same method than self.K for X1=X2, but using the matrix 

138 is then symmetric. 

139 """ 

140 n, D = np.atleast_2d(X).shape 

141 K = np.identity(n * (D + 1)) 

142 self.D = D 

143 D1 = D + 1 

144 

145 # fill upper triangular: 

146 for i in range(n): 

147 for j in range(i + 1, n): 

148 k = self.kernel(X[i], X[j]) 

149 K[i * D1:(i + 1) * D1, j * D1:(j + 1) * D1] = k 

150 K[j * D1:(j + 1) * D1, i * D1:(i + 1) * D1] = k.T 

151 K[i * D1:(i + 1) * D1, i * D1:(i + 1) * D1] = self.kernel( 

152 X[i], X[i]) 

153 return K 

154 

155 def kernel_vector(self, x, X, nsample): 

156 return np.hstack([self.kernel(x, x2) for x2 in X]) 

157 

158 # ---------Derivatives-------- 

159 def dK_dweight(self, X): 

160 """Return the derivative of K(X,X) respect to the weight """ 

161 return self.K(X, X) * 2 / self.weight 

162 

163 # ----Derivatives of the kernel function respect to the scale --- 

164 def dK_dl_k(self, x1, x2): 

165 """Returns the derivative of the kernel function respect to l""" 

166 return self.squared_distance(x1, x2) / self.l 

167 

168 def dK_dl_j(self, x1, x2): 

169 """Returns the derivative of the gradient of the kernel function 

170 respect to l 

171 """ 

172 prefactor = -2 * (1 - 0.5 * self.squared_distance(x1, x2)) / self.l 

173 return self.kernel_function_gradient(x1, x2) * prefactor 

174 

175 def dK_dl_h(self, x1, x2): 

176 """Returns the derivative of the hessian of the kernel function respect 

177 to l 

178 """ 

179 I = np.identity(self.D) 

180 P = np.outer(x1 - x2, x1 - x2) / self.l**2 

181 prefactor = 1 - 0.5 * self.squared_distance(x1, x2) 

182 return -2 * (prefactor * (I - P) - P) / self.l**3 

183 

184 def dK_dl_matrix(self, x1, x2): 

185 k = np.asarray(self.dK_dl_k(x1, x2)).reshape((1, 1)) 

186 j2 = self.dK_dl_j(x1, x2).reshape(1, -1) 

187 j1 = self.dK_dl_j(x2, x1).reshape(-1, 1) 

188 h = self.dK_dl_h(x1, x2) 

189 return np.block([[k, j2], [j1, h]]) * self.kernel_function(x1, x2) 

190 

191 def dK_dl(self, X): 

192 """Return the derivative of K(X,X) respect of l""" 

193 return np.block([[self.dK_dl_matrix(x1, x2) for x2 in X] for x1 in X]) 

194 

195 def gradient(self, X): 

196 """Computes the gradient of matrix K given the data respect to the 

197 hyperparameters. Note matrix K here is self.K(X,X). 

198 Returns a 2-entry list of n(D+1) x n(D+1) matrices 

199 """ 

200 return [self.dK_dweight(X), self.dK_dl(X)]