Coverage for /builds/ase/ase/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 params: [weight, l] Parameters of the kernel: 

34 weight: prefactor of the exponential 

35 l : scale of the kernel 

36 """ 

37 self.weight = params[0] 

38 self.l = params[1] 

39 

40 def squared_distance(self, x1, x2): 

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

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

43 

44 def kernel(self, x1, x2): 

45 """ This is the squared exponential function""" 

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

47 

48 def dK_dweight(self, x1, x2): 

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

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

51 

52 def dK_dl(self, x1, x2): 

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

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

55 

56 

57class SquaredExponential(SE_kernel): 

58 """Squared exponential kernel with derivatives. 

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

60 Jonsson. 

61 Nudged elastic band calculations accelerated with Gaussian process 

62 regression. Section 3. 

63 

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

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

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

67 the second is the scale (l). 

68 

69 Parameters: 

70 

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

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

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

74 

75 Attributes: 

76 ---------------- 

77 D: int. Dimensionality of the problem to optimize 

78 weight: float. Multiplicative constant to the exponenetial kernel 

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

80 

81 Relevant Methods: 

82 ---------------- 

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

84 attributes 

85 kernel_function: Squared exponential covariance function 

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

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

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

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

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

91 

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

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

94 process. 

95 """ 

96 

97 def __init__(self, dimensionality=None): 

98 self.D = dimensionality 

99 SE_kernel.__init__(self) 

100 

101 def kernel_function(self, x1, x2): 

102 """ This is the squared exponential function""" 

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

104 

105 def kernel_function_gradient(self, x1, x2): 

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

107 x1: first data point 

108 x2: second data point 

109 """ 

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

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

112 return prefactor 

113 

114 def kernel_function_hessian(self, x1, x2): 

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

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

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

118 return prefactor 

119 

120 def kernel(self, x1, x2): 

121 """Squared exponential kernel including derivatives. 

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

123 the manifold. 

124 """ 

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

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

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

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

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

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

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

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

133 

134 def kernel_matrix(self, X): 

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

136 is then symmetric. 

137 """ 

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

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

140 self.D = D 

141 D1 = D + 1 

142 

143 # fill upper triangular: 

144 for i in range(n): 

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

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

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

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

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

150 X[i], X[i]) 

151 return K 

152 

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

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

155 

156 # ---------Derivatives-------- 

157 def dK_dweight(self, X): 

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

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

160 

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

162 def dK_dl_k(self, x1, x2): 

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

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

165 

166 def dK_dl_j(self, x1, x2): 

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

168 respect to l 

169 """ 

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

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

172 

173 def dK_dl_h(self, x1, x2): 

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

175 to l 

176 """ 

177 I = np.identity(self.D) 

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

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

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

181 

182 def dK_dl_matrix(self, x1, x2): 

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

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

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

186 h = self.dK_dl_h(x1, x2) 

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

188 

189 def dK_dl(self, X): 

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

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

192 

193 def gradient(self, X): 

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

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

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

197 """ 

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