from sklearn import datasets, model_selection
import numpy as np
    
PIX_SIZE = 28  # 画像1辺のピクセルサイズ
NUM_FEATURES = PIX_SIZE * PIX_SIZE  # 特徴量の総数 (=784)
    
    
mnist_data, mnist_label = datasets.fetch_openml(
    "mnist_784", return_X_y=True, parser="auto"
)
mnist_data /= 255  # 輝度値を [0, 1] にスケーリング
train_data, test_data, train_label, test_label = model_selection.train_test_split(
    mnist_data, mnist_label, train_size=1000, test_size=1000, random_state=0
)

import pandas as pd
    
    
def apply_mask(data: pd.DataFrame, mask: np.ndarray) -> pd.DataFrame:
    """データに対してマスクを適用 (マスクの True である要素だけを返却)"""
    return data.loc[:, mask]

from sklearn import svm, metrics
    
    
def train_and_evaluate(mask: np.ndarray) -> float:
    """マスク適用後のデータで SVM モデルを学習・評価し、誤答率を返す"""
    masked_train_data = apply_mask(train_data, mask)
    masked_test_data = apply_mask(test_data, mask)
    
    model = svm.SVC()
    model.fit(masked_train_data, train_label)
    prediction = model.predict(masked_test_data)
    return 1.0 - float(metrics.accuracy_score(test_label, prediction))

mask_all_feat = np.array([True] * NUM_FEATURES)
error_all_feat = train_and_evaluate(mask_all_feat)
print(f"{error_all_feat=:.3f}")

from amplify_bbopt import blackbox, BinaryVariable
    
# バイナリ決定変数配列を作成 (要素数はマスクの長さ = 784 に一致させる)
variables = [BinaryVariable() for _ in range(NUM_FEATURES)]
    
    
@blackbox
def my_blackbox_func(x=variables):
    # x は要素数 784 の一次元配列であることを確認
    assert len(x) == NUM_FEATURES
    # 特徴量割合の計算 (マスク中の True の割合)
    feature_ratio = sum(x) / NUM_FEATURES
    # マスク x による学習・評価を実行し、誤答率を計算
    error = train_and_evaluate(np.array(x, dtype=bool))
    # 多目的最適化 (MultiOptimizer) を実行するので list で返却 (いずれも最小化)
    return [error, feature_ratio]

def random_mask(feature_ratio: float, rng: np.random.Generator) -> np.ndarray:
    """feature_ratio の割合だけ特徴量を残すマスクをランダムに生成する"""
    num_kept = int(NUM_FEATURES * feature_ratio)
    # num_kept 個の特徴量を残すマスクを生成
    mask = np.array([True] * num_kept + [False] * (NUM_FEATURES - num_kept))
    rng.shuffle(mask)
    return mask

from amplify_bbopt import Dataset
    
# 初期学習データのサイズ
num_init_data = 10
    
rng = np.random.default_rng(0)
    
x = []  # 学習データのマスク配列
y = []  # マスクの誤答率と特徴量割合のタプルのリスト
for i in range(num_init_data):
    feature_ratio = rng.uniform(1 / NUM_FEATURES, 1.0)
    mask = random_mask(feature_ratio, rng)
    x.append(mask)
    y.append(my_blackbox_func(mask))
    
# 初期学習データの準備
init_data = Dataset(np.array(x), np.array(y))

from amplify import AmplifyAEClient
    
# ソルバークライアントを Amplify AE に設定
client = AmplifyAEClient()
client.parameters.time_limit_ms = 1000  # 1 sec
# client.token = "API トークンを入力してください"

from amplify_bbopt import MultiOptimizer, KMTrainer, ExpScaler
    
# MultiOptimizer のインスタンス化
optimizer = MultiOptimizer(
    blackbox=my_blackbox_func,
    trainer=[KMTrainer(), KMTrainer()],
    client=client,
    training_data=init_data,
    surrogate_data_transformer=[ExpScaler(), ExpScaler()],
)

import matplotlib.pyplot as plt
    
num_iterations = 5  # 最適化サイクル数
    
# 動的重みスケジュールを生成
weight_schedule = np.geomspace(3e-1, 3e1, num_iterations)
    
_, ax = plt.subplots(figsize=(10, 4))
ax.plot(range(1, num_iterations + 1), weight_schedule)
ax.set_xlabel("Number of iterations", fontsize=14)
ax.set_ylabel("Weight of feature", fontsize=14)
ax.set_yscale("log")  # y軸を対数スケールに設定
ax.grid(True)

import logging
    
from amplify_bbopt import AMPLIFY_BBOPT_LOGGER_NAME
    
logger = logging.getLogger(AMPLIFY_BBOPT_LOGGER_NAME)
    
# 最適化サイクルで得られる解を記録するリスト
objectives_opt = []
    
for n_iter in range(num_iterations):
    current_weight = weight_schedule[n_iter]
    
    logger.info(f"=== Iteration: {n_iter + 1}/{num_iterations} ===")
    logger.info(f"current weight: {current_weight:.2f}")
    
    # Step 1: サロゲートモデル関数の構築
    optimizer.train_surrogate()
    
    # Step 2: サロゲートモデルの重みを設定 (2つ目のサロゲートモデル (特徴量の最小化) を徐々に重視する)
    optimizer.surrogate_model[1].weight = current_weight
    
    # Step 3: イジングマシンによるサロゲートモデル関数の最適化 (最小化)
    solutions, _ = optimizer.minimize_surrogate()
    if len(solutions) == 0:
        raise RuntimeError(f"No feasible solution was found in iteration {n_iter}")
    
    # Step 4: 得られた最良解をデータセットに追加
    new_solution = solutions[0]
    new_objective = optimizer.evaluate_objective(new_solution)
    optimizer.add_solution(new_solution, new_objective)
    
    # 得られた解の記録 (可視化用)
    objectives_opt.append(new_objective)
    
    # 最適化履歴の出力
    optimizer._log_objective()

# 初期学習データ
objectives_init = np.array(optimizer.training_data.y[:num_init_data])
    
# 最適化サイクルで得られた解の list を numpy 配列に変換
objectives_opt = np.array(objectives_opt)
    
# 特徴量割合を特徴量数に変換
objectives_init[:, 1] *= NUM_FEATURES
objectives_opt[:, 1] *= NUM_FEATURES

# 各目的関数値のプロット (idx_sm: 何番目の目的関数か (sm: surrogate model))
def plot_objectives(obj_init, obj_opt, idx_sm, title=None):
    obj_init = obj_init[:, idx_sm]
    obj_opt = obj_opt[:, idx_sm]
    num_init_data = len(obj_init)
    
    # 累積最小値の計算
    best_opt = [
        min((obj_init.tolist() + obj_opt.tolist())[: i + 1])
        for i in range(0, num_init_data + len(obj_opt))
    ]
    
    _, ax = plt.subplots(figsize=(10, 6))
    
    # 初期学習データのプロット (最適化サイクルの前なので、x軸は負の値を使用)
    ax.plot(range(-num_init_data + 1, 1), obj_init, "b", marker="o", linestyle="-")
    # 最適化サイクルで得られた解のプロット
    ax.plot(range(1, len(obj_opt) + 1), obj_opt, "r", marker="o", linestyle="-")
    # 累積最小値のプロット
    ax.plot(range(-num_init_data + 1, len(obj_opt) + 1), best_opt, "k")
    
    ax.set_xlabel("Number of iterations", fontsize=18)
    if title is not None:
        ax.set_ylabel(title, fontsize=18)
    ax.tick_params(labelsize=18)
    
    ax.grid(True)
    
    
plot_objectives(objectives_init, objectives_opt, 0, title="Error")
plot_objectives(objectives_init, objectives_opt, 1, title="Number of features")

# パレート最適集合の計算
def get_pareto_front(objs):
    cand = np.ones(objs.shape[0], dtype=bool)
    for i, c in enumerate(objs):
        if cand[i]:
            cand[cand] = np.any(objs[cand] < c, axis=1)
            cand[i] = True
    return cand
    
    
# 全ての解 (初期データ + 探索データ) からパレート最適集合を計算
objectives_combined = np.concatenate([objectives_init, objectives_opt])
pareto_mask = get_pareto_front(objectives_combined)
pareto_front = objectives_combined[pareto_mask]
pareto_front = pareto_front[pareto_front[:, 0].argsort()]

from matplotlib.colors import LinearSegmentedColormap
    
_, ax = plt.subplots(figsize=(10, 6))
# 初期学習データのプロット
ax.scatter(
    objectives_init[:, 0], objectives_init[:, 1], color="blue", label="Initial data"
)
# 最適化サイクルで得られた解のプロット (探索順にグラデーション)
cm = LinearSegmentedColormap.from_list(
    "custom_cm", [(0, "orangered"), (1, "limegreen")]
)
cols = np.linspace(0, 1, len(objectives_opt))
ax.scatter(
    objectives_opt[:, 0],
    objectives_opt[:, 1],
    c=cols,
    cmap=cm,
    edgecolors="k",
    label="Optimized solutions",
)
# パレートフロントのプロット
ax.plot(pareto_front[:, 0], pareto_front[:, 1], "--k", label="Pareto front")
# 全特徴量を使った場合の誤答率を基準線として描画
ax.axvline(
    error_all_feat,
    color="gray",
    linestyle=":",
    label=f"All-features error ({error_all_feat:.3f})",
)
# 軸ラベルの設定
ax.set_xlabel("Error", fontsize=18)
ax.set_ylabel("Number of features", fontsize=18)
    
ax.set_title("Pareto front", fontsize=18)
ax.legend(fontsize=12)
ax.grid(True)

ブラックボックス最適化による機械学習のための特徴量選択¶

はじめに¶

目次¶

1. MNIST と特徴量選択¶

1.1. MNIST とは¶

1.2. 特徴量選択のフロー¶

1.3. 2 つの目的とパレートフロント¶

1.4. データの読み込みと前処理¶

1.5. マスク適用・学習・評価のための関数¶

1.6. 全特徴量を考慮した場合のモデル精度¶

2. MultiOptimizer による多目的ブラックボックス最適化¶

2.1. ブラックボックス関数¶

2.2. 初期学習データの生成¶

2.3. ソルバーの設定¶

2.4. MultiOptimizer のインスタンス化¶

2.5. 重みスケジュールの設計¶

2.6. 最適化サイクルの実行¶

3. 最適化履歴の可視化¶

3.1. 探索された解に対する目的関数値の推移¶

3.2. パレートフロント¶

3.3. まとめ¶