Surrogate model-based Derivative-Free Optimizations

1
import numpy as np
2

3
# -------------------------
4
# 목적 함수 (예시)
5
# -------------------------
6
def objective(x):
7
    """노이즈가 있는 quadratic"""
8
    return np.sum(np.array(x)**2) + np.random.randn()*0.1
9

10
# -------------------------
11
# RBF Surrogate
12
# -------------------------
13
class RBF_Surrogate:
14
    def __init__(self, phi="cubic"):
15
        if phi == "gaussian":
16
            self.phi = lambda r: np.exp(-(r**2))
17
        elif phi == "linear":
18
            self.phi = lambda r: r
19
        else:  # cubic
20
            self.phi = lambda r: r**3
21

22
    def fit(self, X, y):
23
        self.X = np.array(X)
24
        self.y = np.array(y)
25
        N, d = self.X.shape
26

27
        Phi = np.zeros((N, N))
28
        for i in range(N):
29
            for j in range(N):
30
                Phi[i, j] = self.phi(np.linalg.norm(self.X[i] - self.X[j]))
31

32
        P = np.hstack([self.X, np.ones((N, 1))])
33
        A = np.block([
34
            [Phi, P],
35
            [P.T, np.zeros((d+1, d+1))]
36
        ])
37
        b = np.concatenate([self.y, np.zeros(d+1)])
38

39
        sol = np.linalg.solve(A, b)
40
        self.lmbda = sol[:N]
41
        self.beta = sol[N:]
42

43
    def predict(self, X_new):
44
        X_new = np.atleast_2d(X_new)
45
        y_pred = np.zeros(len(X_new))
46

47
        for k, xk in enumerate(X_new):
48
            rbf_sum = 0
49
            for i in range(len(self.X)):
50
                r = np.linalg.norm(xk - self.X[i])
51
                rbf_sum += self.lmbda[i] * self.phi(r)
52
            trend = np.dot(self.beta[:-1], xk) + self.beta[-1]
53
            y_pred[k] = rbf_sum + trend
54
        return y_pred
55

56
# -------------------------
57
# Data 관리 클래스
58
# -------------------------
59
class DataStore:
60
    def __init__(self, dim, m_init=5, bounds=(-2, 2)):
61
        self.dim = dim
62
        self.bounds = bounds
63
        self.m = m_init
64

65
        self.S = np.random.uniform(bounds[0], bounds[1], size=(m_init, dim))
66
        self.Y = np.array([objective(s) for s in self.S])
67
        self.rbf = RBF_Surrogate()
68
        self.update_rbf()
69

70
    def update_rbf(self):
71
        self.rbf.fit(self.S, self.Y)
72

73
    def add_point(self, x):
74
        y = objective(x)
75
        self.S = np.vstack([self.S, x])
76
        self.Y = np.append(self.Y, y)
77
        self.m += 1
78
        self.update_rbf()
79

80
# -------------------------
81
# 후보 점수 계산 (RBF 예측 + 거리)
82
# -------------------------
83
def RBF_score(candidates, data, omega=0.5):
84
    s = data.rbf.predict(candidates)
85
    dist = np.array([np.min(np.linalg.norm(data.S - c, axis=1)) for c in candidates])
86

87
    # 0~1 scaling
88
    s_scaled = (s - s.min()) / (s.max() - s.min() + 1e-8)
89
    d_scaled = (dist.max() - dist) / (dist.max() - dist.min() + 1e-8)
90

91
    score = omega*d_scaled + (1-omega)*s_scaled
92
    return score
93

94
# -------------------------
95
# 최적화 루프
96
# -------------------------
97
def run_rbf(dim=2, n_iter=10):
98
    np.random.seed(0)
99
    data = DataStore(dim=dim)
100

101
    for it in range(n_iter):
102
        best_idx = np.argmin(data.Y)
103
        best_x = data.S[best_idx]
104

105
        # 후보 생성
106
        local_candidates = best_x + np.random.uniform(-0.5, 0.5, size=(20, dim))
107
        global_candidates = np.random.uniform(data.bounds[0], data.bounds[1], size=(20, dim))
108
        candidates = np.vstack([local_candidates, global_candidates])
109

110
        # 점수 계산
111
        score = RBF_score(candidates, data, omega=0.5)
112

113
        # 다음 점 선택
114
        next_x = candidates[np.argmin(score)]
115
        data.add_point(next_x)
116

117
        print(f"Iter {it}: best_y={data.Y.min():.4f}, next_x={next_x}")
118

119
    return data
120

121
# 실행
122
if __name__ == "__main__":
123
    result = run_rbf(dim=2, n_iter=10)

1
import numpy as np
2
import array
3
from deap import base, creator, tools, algorithms
4
from sklearn.gaussian_process import GaussianProcessRegressor
5
from sklearn.gaussian_process.kernels import RBF
6
from scipy.stats import norm
7

8
# -------------------------
9
# 목적 함수 (black-box 예시)
10
# -------------------------
11
def objective(x):
12
    """노이즈가 있는 2차 함수"""
13
    return np.sum(np.array(x) ** 2) + np.random.randn() * 0.1
14

15
# -------------------------
16
# 데이터 저장용 클래스 흉내
17
# -------------------------
18
class DataStore:
19
    def __init__(self, dim, m_init=5, bounds=(-2, 2)):
20
        self.dim = dim
21
        self.bounds = bounds
22
        self.m = m_init  # 초기 샘플 수
23

24
        self.S = np.random.uniform(bounds[0], bounds[1], size=(m_init, dim))
25
        self.Y = np.array([objective(s) for s in self.S])
26

27
        # GP surrogate
28
        self.gpr = GaussianProcessRegressor(kernel=RBF(), random_state=0)
29
        self.update_gp()
30

31
    def update_gp(self):
32
        self.gpr.fit(self.S, self.Y.reshape(-1, 1))
33

34
    def add_point(self, x):
35
        y = objective(x)
36
        self.S = np.vstack([self.S, x])
37
        self.Y = np.append(self.Y, y)
38
        self.m += 1
39
        self.update_gp()
40

41
# -------------------------
42
# EI 함수 (스크린샷 기반)
43
# -------------------------
44
def Expected_improvement(x, data):
45
    x_to_predict = np.array(x).reshape(1, -1)
46
    mu, sigma = data.gpr.predict(x_to_predict, return_std=True)
47

48
    greater_is_better = False
49
    if greater_is_better:
50
        loss_optimum = np.max(data.Y[:data.m])
51
    else:
52
        loss_optimum = np.min(data.Y[:data.m])
53

54
    scaling_factor = (-1) ** (not greater_is_better)
55

56
    with np.errstate(divide='ignore'):
57
        Z = scaling_factor * (mu - loss_optimum) / sigma
58
        expected_improvement = scaling_factor * (mu - loss_optimum) * norm.cdf(Z) + sigma * norm.pdf(Z)
59
        expected_improvement[sigma == 0.0] = 0.0
60
    return -expected_improvement[0]  # deap은 최소화 문제
61

62
# -------------------------
63
# DEAP 설정
64
# -------------------------
65
def run_gp_ei(dim=2, n_iter=10):
66
    data = DataStore(dim=dim)
67

68
    # Fitness 정의
69
    creator.create("FitnessMin", base.Fitness, weights=(-1.0,))
70
    creator.create("Individual", array.array, typecode='d', fitness=creator.FitnessMin)
71

72
    toolbox = base.Toolbox()
73

74
    # 각 차원 범위
75
    for i in range(dim):
76
        INT_MIN, INT_MAX = data.bounds
77
        toolbox.register(f"attr_float_{i}", np.random.uniform, INT_MIN, INT_MAX)
78

79
    # Individual 생성
80
    toolbox.register("individual", tools.initCycle, creator.Individual,
81
                     (toolbox.__getattribute__(f"attr_float_{i}") for i in range(dim)), n=1)
82
    toolbox.register("population", tools.initRepeat, list, toolbox.individual)
83

84
    # Evaluate = EI
85
    toolbox.register("evaluate", Expected_improvement, data=data)
86
    toolbox.register("mate", tools.cxTwoPoint)
87
    toolbox.register("mutate", tools.mutUniformInt,
88
                     low=[data.bounds[0]]*dim, up=[data.bounds[1]]*dim, indpb=0.2)
89
    toolbox.register("select", tools.selTournament, tournsize=3)
90

91
    # -------------------------
92
    # 최적화 루프
93
    # -------------------------
94
    for it in range(n_iter):
95
        pop = toolbox.population(n=20)
96
        hof = tools.HallOfFame(1)
97

98
        algorithms.eaSimple(pop, toolbox, cxpb=0.5, mutpb=0.3,
99
                            ngen=15, halloffame=hof, verbose=False)
100

101
        next_x = np.array(hof[0])
102
        data.add_point(next_x)
103

104
        print(f"Iter {it}: best_y={data.Y.min():.4f}, next_x={next_x}")
105

106
    return data
107

108
# 실행
109
if __name__ == "__main__":
110
    result = run_gp_ei(dim=2, n_iter=10)