From 44b8f551d2bddc6dee6e102e250d1cf9a32f6b1b Mon Sep 17 00:00:00 2001
From: Tizian Wenzel <tizian.wenzel@ians.uni-stuttgart.de>
Date: Thu, 9 Feb 2023 18:40:36 +0100
Subject: [PATCH] Some improvements for the case when not a list of kernels is
provided.
---
example_files/01_testfile.py | 2 +-
vkoga_2L.py | 88 ++++++++++++++++--------------------
2 files changed, 40 insertions(+), 50 deletions(-)
diff --git a/example_files/01_testfile.py b/example_files/01_testfile.py
index 50bfe80..64b944f 100644
--- a/example_files/01_testfile.py
+++ b/example_files/01_testfile.py
@@ -3,7 +3,7 @@
import numpy as np
from matplotlib import pyplot as plt
import kernels, tkernels
-from vkoga import VKOGA_2L
+from vkoga_2L import VKOGA_2L
np.random.seed(1)
diff --git a/vkoga_2L.py b/vkoga_2L.py
index 6625bb2..fbd6f3f 100644
--- a/vkoga_2L.py
+++ b/vkoga_2L.py
@@ -7,31 +7,37 @@ from sklearn.utils.validation import check_X_y, check_array, check_is_fitted
from scipy.spatial import distance_matrix
import torch
from utilities import OptimizedKernel
-
+
+
# VKOGA implementation
class VKOGA_2L(BaseEstimator):
-
+
def __init__(self, kernel=Matern(k=2), flag_2L_optimization=False,
verbose=True, n_report=10,
greedy_type='f_greedy', reg_par=0, restr_par=0,
tol_f=1e-10, tol_p=1e-10,
reg_para_optim=1e-3, learning_rate=5e-3, n_epochs_optim=10, batch_size=32):
-
+
# Set the verbosity on/off
self.verbose = verbose
-
+
# Set the frequency of report
self.n_report = n_report
-
- # Set the params defining the method
+
+ # Set the params defining the method
self.flag_2L_optimization = flag_2L_optimization
- if self.flag_2L_optimization:
- assert type(kernel) is type([]), 'If flag_2L_optimization, then two kernels need to be provided!'
+ if self.flag_2L_optimization and type(kernel) is type([]):
self.kernel_t = kernel[1]
self.kernel = kernel[0]
else:
self.kernel_t = None
- self.kernel = kernel
+ if type(kernel) is type([]):
+ print('If kernel optimization is desired, set flag_2L_optimization to True!')
+ self.kernel = kernel[0]
+ else:
+ print('If kernel optimization is desired, provide numpy and torch kernel!')
+ self.flag_2L_optimization = False
+ self.kernel = kernel
self.greedy_type = greedy_type
self.reg_par = reg_par
@@ -53,16 +59,15 @@ class VKOGA_2L(BaseEstimator):
self.Cut_ = None
self.c = None
-
# Initialize the convergence history
self.train_hist = {}
self.train_hist['n'] = []
self.train_hist['f'] = []
self.train_hist['p'] = []
- self.train_hist['p selected'] = [] # list of selected power vals
+ self.train_hist['p selected'] = [] # list of selected power vals
self.train_hist['f val'] = []
self.train_hist['p val'] = []
-
+
def selection_rule(self, f, p):
if self.restr_par > 0:
p_ = np.max(p)
@@ -85,7 +90,7 @@ class VKOGA_2L(BaseEstimator):
f_max = np.max(f)
idx = np.argmax(p)
p_max = p[idx]
- elif self.greedy_type == 'rand': # pick some random point - David Holzmüller asked me about its performance
+ elif self.greedy_type == 'rand': # pick some random point - David Holzmüller asked me about its performance
f_max = np.max(f)
p_max = np.max(p)
idx = np.random.randint(len(p))
@@ -110,14 +115,12 @@ class VKOGA_2L(BaseEstimator):
if len(y_val.shape) == 1:
y_val = y_val[:, None]
-
# Check whether already fitted - restart in case we used kernel optimization
if self.ctrs_ is None or self.flag_2L_optimization:
self.ctrs_ = np.zeros((0, X.shape[1]))
self.Cut_ = np.zeros((0, 0))
self.c = np.zeros((0, y.shape[1]))
-
# Check whether "new X" and previously chosen centers overlap
list_truly_new_X = []
if self.ctrs_.shape[0] > 0:
@@ -131,7 +134,6 @@ class VKOGA_2L(BaseEstimator):
X = X[list_truly_new_X, :]
y = y[list_truly_new_X, :]
-
# Optimize the kernel
if self.flag_2L_optimization:
model_OptimKernel = OptimizedKernel(kernel=self.kernel_t, dim=X.shape[1],
@@ -144,12 +146,11 @@ class VKOGA_2L(BaseEstimator):
self.A = model_OptimKernel.A.detach().numpy()
- X = X @ self.A # we can incorporate the kernel modification directly into the data
+ X = X @ self.A # we can incorporate the kernel modification directly into the data
else:
self.A = np.eye(X.shape[1])
-
# Initialize the residual and update the given y values by substracting the current model
y = np.array(y)
if len(y.shape) == 1:
@@ -158,45 +159,38 @@ class VKOGA_2L(BaseEstimator):
if self.flag_val:
y_val = y_val - self.predict(X_val)
-
# Get the data dimension
N, q = y.shape
if self.flag_val:
N_val = y_val.shape[0]
-
# Set maxIter_continue
if maxIter is None or maxIter > N:
self.maxIter = 100
else:
self.maxIter = maxIter
-
# Check compatibility of restriction
if self.greedy_type == 'p_greedy':
self.restr_par = 0
if not self.reg_par == 0:
self.restr_par = 0
-
# Initialize list for the chosen and non-chosen indices
indI_ = []
notIndI = list(range(N))
c = np.zeros((self.maxIter, q))
-
# Compute the Newton basis values (related to the old centers) on the new X
Vx_new_X_old_ctrs = self.kernel.eval(X, self.ctrs_) @ self.Cut_.transpose()
if self.flag_val:
Vx_val_new_X_old_ctrs = self.kernel.eval(X_val, self.ctrs_) @ self.Cut_.transpose()
-
# Initialize arrays for the Newton basis values (related to the new centers) on the new X
Vx = np.zeros((N, self.maxIter))
if self.flag_val:
Vx_val = np.zeros((N_val, self.maxIter))
-
# Compute the powervals on X and X_val
p = self.kernel.diagonal(X) + self.reg_par
p = p - np.sum(Vx_new_X_old_ctrs ** 2, axis=1)
@@ -204,13 +198,11 @@ class VKOGA_2L(BaseEstimator):
p_val = self.kernel.diagonal(X_val) + self.reg_par
p_val = p_val - np.sum(Vx_val_new_X_old_ctrs ** 2, axis=1)
-
# Extend Cut_ matrix, i.e. continue to build on old self.Cut_ matrix
N_ctrs_so_far = self.Cut_.shape[0]
Cut_ = np.zeros((N_ctrs_so_far + self.maxIter, N_ctrs_so_far + self.maxIter))
Cut_[:N_ctrs_so_far, :N_ctrs_so_far] = self.Cut_
-
# Iterative selection of new points
self.print_message('begin')
for n in range(self.maxIter):
@@ -244,19 +236,18 @@ class VKOGA_2L(BaseEstimator):
break
# compute the nth basis (including normalization)# ToDo: Also old Vx need to be substracted here!
- Vx[notIndI, n] = self.kernel.eval(X[notIndI, :], X[indI_[n], :])[:, 0]\
- - Vx_new_X_old_ctrs[notIndI, :] @ Vx_new_X_old_ctrs[indI_[n], :].transpose()\
- - Vx[notIndI, :n+1] @ Vx[indI_[n], 0:n+1].transpose()
+ Vx[notIndI, n] = self.kernel.eval(X[notIndI, :], X[indI_[n], :])[:, 0] \
+ - Vx_new_X_old_ctrs[notIndI, :] @ Vx_new_X_old_ctrs[indI_[n], :].transpose() \
+ - Vx[notIndI, :n + 1] @ Vx[indI_[n], 0:n + 1].transpose()
Vx[indI_[n], n] += self.reg_par
Vx[notIndI, n] = Vx[notIndI, n] / np.sqrt(p[indI_[n]])
if self.flag_val:
- Vx_val[:, n] = self.kernel.eval(X_val, X[indI_[n], :])[:, 0]\
- - Vx_val_new_X_old_ctrs[:, :] @ Vx_new_X_old_ctrs[indI_[n], :].transpose()\
- - Vx_val[:, :n+1] @ Vx[indI_[n], 0:n+1].transpose()
+ Vx_val[:, n] = self.kernel.eval(X_val, X[indI_[n], :])[:, 0] \
+ - Vx_val_new_X_old_ctrs[:, :] @ Vx_new_X_old_ctrs[indI_[n], :].transpose() \
+ - Vx_val[:, :n + 1] @ Vx[indI_[n], 0:n + 1].transpose()
Vx_val[:, n] = Vx_val[:, n] / np.sqrt(p[indI_[n]])
-
# update the change of basis
Cut_new_row = np.ones(N_ctrs_so_far + n + 1)
Cut_new_row[:N_ctrs_so_far + n] = \
@@ -288,16 +279,14 @@ class VKOGA_2L(BaseEstimator):
self.print_message('end')
# Define coefficients and centers
- self.c = np.concatenate((self.c, c[:n + 1]))
+ self.c = np.concatenate((self.c, c[:n + 1]))
self.Cut_ = Cut_[:N_ctrs_so_far + n + 1, :N_ctrs_so_far + n + 1]
- self.indI_ = indI_[:n + 1] # Mind: These are only the indices of the latest points
+ self.indI_ = indI_[:n + 1] # Mind: These are only the indices of the latest points
self.coef_ = self.Cut_.transpose() @ self.c
self.ctrs_ = np.concatenate((self.ctrs_, X[self.indI_, :]), axis=0)
-
return self
-
def predict(self, X):
# Check is fit has been called
# check_is_fitted(self, 'coef_') # ToDo: Remove this one!
@@ -352,10 +341,11 @@ class VKOGA_2L(BaseEstimator):
print(' |_ regularization par. : %2.2e' % self.reg_par)
print(' |_ restriction par. : %2.2e' % self.restr_par)
print('')
-
+
if self.verbose and when == 'end':
print('Training completed with')
- print(' |_ selected points : %8d / %8d' % (self.train_hist['n'][-1], self.ctrs_.shape[0] + self.maxIter))
+ print(' |_ selected points : %8d / %8d' % (
+ self.train_hist['n'][-1], self.ctrs_.shape[0] + self.maxIter))
if self.flag_val:
print(' |_ train, val residual : %2.2e / %2.2e, %2.2e' %
(self.train_hist['f'][-1], self.tol_f, self.train_hist['f val'][-1]))
@@ -364,10 +354,11 @@ class VKOGA_2L(BaseEstimator):
else:
print(' |_ train residual : %2.2e / %2.2e' % (self.train_hist['f'][-1], self.tol_f))
print(' |_ train power fun : %2.2e / %2.2e' % (self.train_hist['p'][-1], self.tol_p))
-
+
if self.verbose and when == 'track':
print('Training ongoing with')
- print(' |_ selected points : %8d / %8d' % (self.train_hist['n'][-1], self.ctrs_.shape[0] + self.maxIter))
+ print(' |_ selected points : %8d / %8d' % (
+ self.train_hist['n'][-1], self.ctrs_.shape[0] + self.maxIter))
if self.flag_val:
print(' |_ train, val residual : %2.2e / %2.2e, %2.2e' %
(self.train_hist['f'][-1], self.tol_f, self.train_hist['f val'][-1]))
@@ -378,17 +369,16 @@ class VKOGA_2L(BaseEstimator):
print(' |_ train power fun : %2.2e / %2.2e' % (self.train_hist['p'][-1], self.tol_p))
-
-#%% Utilities to
+# %% Utilities to
import pickle
+
+
def save_object(obj, filename):
with open(filename, 'wb') as output:
pickle.dump(obj, output, pickle.HIGHEST_PROTOCOL)
-def load_object(filename):
+
+def load_object(filename):
with open(filename, 'rb') as input:
- obj = pickle.load(input)
+ obj = pickle.load(input)
return obj
-
-
-
--
GitLab