-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathexp_1_4_c.py
155 lines (123 loc) · 5.89 KB
/
exp_1_4_c.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
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
import numpy as np
import scipy.linalg as la
from sklearn.gaussian_process.kernels import RBF, Matern
from sklearn.gaussian_process import GaussianProcessRegressor
from utilities import dataset_generation_uniform_normal, aposteriori_scaling, aposteriori_scalings_generator, aposteriori_rescalings_generator, check_bounds_on_grid
from utilities import sample_rkhs_se_onb
from experiments import run_learning_instance_experiment
from joblib import Parallel, delayed
noise_level = 0.5
rkhs_norm = 2
eps = 0.1
# Implementation of the robust uncertainty bound for unstructured misspecification of the kernel
def add_scaling(eps, N, Kr, ks_norm, scaling):
return (1/scaling + la.norm(Kr, ord=2))*(ks_norm + np.sqrt(N)*eps) + la.norm(Kr, ord=2)*np.sqrt(N)*eps
def add_var(eps, N, Kr, ks, Cs):
return eps + np.sqrt(N)*eps*la.norm(Kr @ ks, ord=2, axis=0) \
+ (np.sqrt(N)*eps + la.norm(ks, ord=2, axis=0)).flatten()*Cs
def calc_robust_width(eps, N, K, ks, post_mean, post_var, ys_train, B, R, scaling, delta):
Kr = la.inv(K + scaling*np.identity(N))
alpha_bar = np.max([1, scaling]).item()
detval =np.log(la.det(alpha_bar/scaling*K + alpha_bar*np.identity(K.shape[0])))
beta = B + R/scaling*np.sqrt(detval - 2*np.log(delta))
Cs = add_scaling(eps, N, Kr, la.norm(ks, ord=2, axis=0), scaling)
return Cs*la.norm(ys_train, ord=2) + beta*np.sqrt(post_var.flatten() + add_var(eps, N, Kr, ks, Cs))
def run_robust_learning_instance_experiment(xs, ys, config):
"""Run one learning instance on a given target function
A learning instance consists of generating a dataset (details specified by
the 'dataset_generation' entry of config), training with Gaussian Process Regression
(settings specified by the 'training' entry of config) and testing the uniform
validity of various bounds (assuming tube form with widths given by scaled posterior SD
with the scalings generated by the function in the 'scalings_generator' entry in config).
Arguments
xs 1d numpy array containing the inputs for evaluation
ys 1d numpy array containing the target function evaluated on xs
config Dictionary containing the settings for this experimental instance
Returns
A 1d numpy array of length n_bounds, where n_bounds is the number of different widths
generated by the scalings_generator, with 0 if the bound holds and 1 if the
corresponding bound is violated at least once in the evaluation grid.
"""
# Generate training data
n_samples = config['dataset_generation']['n_samples']
dataset_generator = config['dataset_generation']['dataset_generator']
(xs_train, ys_train) = dataset_generator(xs, ys, n_samples)
# Learn with GPR
kernel = config['training']['kernel']
gpr = GaussianProcessRegressor(
kernel=kernel,
alpha=config['training']['noise_level_train'],
optimizer=None).fit(xs_train,ys_train)
(post_mean, post_sd) = gpr.predict(xs, return_std=True)
K = kernel(xs_train)
# Build uncertainty tube
deltas = np.array([0.1, 0.01, 0.001, 0.0001])
robust_widths = np.zeros([len(deltas), len(xs)])
ks = kernel(xs_train, xs)
for i in range(len(deltas)):
robust_widths[i,:] = calc_robust_width(eps,
n_samples,
K,
ks,
post_mean,
post_sd**2,
ys_train,
rkhs_norm,
noise_level,
1,
deltas[i])
# Check bounds
upper_bounds = post_mean.reshape([1,-1]) + robust_widths
lower_bounds = post_mean.reshape([1,-1]) - robust_widths
return (robust_widths, check_bounds_on_grid(ys, upper_bounds, lower_bounds))
def run_function_instance(config, func_id):
# Generate function
func_config = config['target_function']
xs = func_config['xs'].reshape([-1,1])
kernel = func_config['kernel']
rkhs_norm = func_config['rkhs_norm']
n_kernels = func_config['n_kernels']
ys = sample_rkhs_se_onb(xs, rkhs_norm, 0.2)
# Run learning instances
instance_config = {
'dataset_generation': config['dataset_generation'],
'training': config['training'],
}
instance_experiment = lambda: run_robust_learning_instance_experiment(xs, ys, instance_config)
results = Parallel(n_jobs=config['n_jobs'])(delayed(instance_experiment)()
for _ in range(config['n_rep_training']))
scalings_list = [result[0] for result in results]
failures_list = [result[1] for result in results]
scalings_array = np.vstack(scalings_list)
failures_array = np.vstack(failures_list)
# Save
with open('outputs/' + config['experiment_prefix'] + '_func_' + str(func_id) + '.npz', 'wb') as f:
np.savez(f, xs=xs, ys=ys, scalings=scalings_array, failures=failures_array)
# Config
kernel_gpr = RBF(length_scale=0.5)
dataset_generation_config = {
'n_samples': 50,
'dataset_generator': lambda xs, ys, n_samples: dataset_generation_uniform_normal(xs, ys, 50, noise_level)
}
training_config = {
'kernel': kernel_gpr,
'noise_level_train': noise_level
}
func_config = {
'xs': np.linspace(-1, 1, 1000),
'kernel': kernel_gpr,
'rkhs_norm': rkhs_norm,
'n_kernels': 100
}
config = {
'target_function': func_config,
'dataset_generation': dataset_generation_config,
'training': training_config,
'n_jobs': 14,
'n_rep_training': 10000,
'n_rep_funcs': 50,
'experiment_prefix': 'exp_1_4_c'
}
# Run and store
for i in range(config['n_rep_funcs']):
run_function_instance(config, i)