-
Notifications
You must be signed in to change notification settings - Fork 0
/
KTCRegularization.py
55 lines (46 loc) · 2.08 KB
/
KTCRegularization.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
import numpy as np
class SMPrior:
def __init__(self, ginv, corrlength, var, mean, covariancetype=None):
self.corrlength = corrlength
self.mean = mean
self.c = 1e-9 # default value
if covariancetype is not None:
self.covariancetype = covariancetype
else:
self.covariancetype = 'Squared Distance' # default
self.compute_L(ginv, corrlength, var)
def compute_L(self, g, corrlength, var):
ng = g.shape[0]
a = var - self.c
b = np.sqrt(-corrlength**2 / (2 * np.log(0.01)))
Gamma_pr = np.zeros((ng, ng))
for ii in range(ng):
for jj in range(ii, ng):
dist_ij = np.linalg.norm(g[ii, :] - g[jj, :])
if self.covariancetype == 'Squared Distance':
gamma_ij = a * np.exp(-dist_ij**2 / (2 * b**2))
elif self.covariancetype == 'Ornstein-Uhlenbeck':
gamma_ij = a * np.exp(-dist_ij / corrlength)
else:
raise ValueError('Unrecognized prior covariance type')
if ii == jj:
gamma_ij = gamma_ij + self.c
Gamma_pr[ii, jj] = gamma_ij
Gamma_pr[jj, ii] = gamma_ij
self.L = np.linalg.cholesky(np.linalg.inv(Gamma_pr)).T
def draw_samples(self, nsamples):
samples = self.mean + np.linalg.solve(self.L, np.random.randn(self.L.shape[0], nsamples))
return samples
def eval_fun(self, args):
sigma = args[0]
res = 0.5 * np.linalg.norm(self.L @ (sigma - self.mean))**2
return res
def compute_hess_and_grad(self, args, nparam):
sigma = args[0]
Hess = self.L.T @ self.L
grad = Hess @ (sigma - self.mean)
if nparam > len(sigma):
Hess = np.block([[Hess, np.zeros((len(sigma), nparam - len(sigma)))],
[np.zeros((nparam - len(sigma), len(sigma))), np.zeros((nparam - len(sigma), nparam - len(sigma)))]])
grad = np.concatenate([grad, np.zeros(nparam - len(sigma))])
return Hess, grad