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

1import numpy as np 

2import numpy.linalg as la 

3 

4 

5class Kernel(): 

6 def __init__(self): 

7 pass 

8 

9 def set_params(self, params): 

10 pass 

11 

12 def kernel(self, x1, x2): 

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

14 

15 def K(self, X1, X2): 

16 """Compute the kernel matrix """ 

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

18 

19 

20class SE_kernel(Kernel): 

21 """Squared exponential kernel without derivatives""" 

22 

23 def __init__(self): 

24 Kernel.__init__(self) 

25 

26 def set_params(self, params): 

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

28 

29 Parameters: 

30 

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

32 weight: prefactor of the exponential 

33 l : scale of the kernel 

34 """ 

35 self.weight = params[0] 

36 self.l = params[1] 

37 

38 def squared_distance(self, x1, x2): 

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

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

41 

42 def kernel(self, x1, x2): 

43 """ This is the squared exponential function""" 

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

45 

46 def dK_dweight(self, x1, x2): 

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

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

49 

50 def dK_dl(self, x1, x2): 

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

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

53 

54 

55class SquaredExponential(SE_kernel): 

56 """Squared exponential kernel with derivatives. 

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

58 Jonsson. 

59 Nudged elastic band calculations accelerated with Gaussian process 

60 regression. Section 3. 

61 

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

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

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

65 the second is the scale (l). 

66 

67 Parameters: 

68 

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

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

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

72 

73 Attributes: 

74 ---------------- 

75 D: int. Dimensionality of the problem to optimize 

76 weight: float. Multiplicative constant to the exponenetial kernel 

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

78 

79 Relevant Methods: 

80 ---------------- 

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

82 attributes 

83 kernel_function: Squared exponential covariance function 

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

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

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

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

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

89 

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

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

92 process. 

93 """ 

94 

95 def __init__(self, dimensionality=None): 

96 self.D = dimensionality 

97 SE_kernel.__init__(self) 

98 

99 def kernel_function(self, x1, x2): 

100 """ This is the squared exponential function""" 

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

102 

103 def kernel_function_gradient(self, x1, x2): 

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

105 x1: first data point 

106 x2: second data point 

107 """ 

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

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

110 return prefactor 

111 

112 def kernel_function_hessian(self, x1, x2): 

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

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

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

116 return prefactor 

117 

118 def kernel(self, x1, x2): 

119 """Squared exponential kernel including derivatives. 

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

121 the manifold. 

122 """ 

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

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

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

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

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

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

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

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

131 

132 def kernel_matrix(self, X): 

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

134 is then symmetric. 

135 """ 

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

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

138 self.D = D 

139 D1 = D + 1 

140 

141 # fill upper triangular: 

142 for i in range(n): 

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

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

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

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

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

148 X[i], X[i]) 

149 return K 

150 

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

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

153 

154 # ---------Derivatives-------- 

155 def dK_dweight(self, X): 

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

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

158 

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

160 def dK_dl_k(self, x1, x2): 

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

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

163 

164 def dK_dl_j(self, x1, x2): 

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

166 respect to l 

167 """ 

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

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

170 

171 def dK_dl_h(self, x1, x2): 

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

173 to l 

174 """ 

175 I = np.identity(self.D) 

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

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

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

179 

180 def dK_dl_matrix(self, x1, x2): 

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

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

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

184 h = self.dK_dl_h(x1, x2) 

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

186 

187 def dK_dl(self, X): 

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

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

190 

191 def gradient(self, X): 

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

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

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

195 """ 

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