マテリアルズインフォマティクス（MI）入門⑦【多目的最適化とパレートフロント】

この記事から学べること

• マルチターゲット回帰モデルの構築: 複数の物性値を同時に予測するモデルを、matminer の steel_strength データセットを使って構築する方法を学びます。
• トレードオフの可視化: 予測結果を2軸のグラフにプロットし、2つの強さ指標の間にどのような関係性があるのかを視覚的に理解します。
• パレートフロントの探求: 複数の性能を天秤にかけた際の「最も効率的な解の集合」であるパレートフロントとは何かを理論的に学び、それを既存の予測データから特定し、グラフ上に描画するハンズオンを実践します。

関連理論の解説

1. トレードオフ：材料科学における現実的な課題

二つの異なる特性を同時に最大化しようとすると、多くの場合、一方を向上させることがもう一方の向上を制約するトレードオフの関係に直面します。

• 降伏強度 (Yield Strength): 材料が塑性変形（元の形に戻らない変形）を起こし始める応力。材料の弾性的な限界を示します。
• 引張強さ (Tensile Strength): 材料が耐えることができる最大の応力。これを超えると、材料は「くびれ」始め、やがて破断に至ります。

降伏強度と引張強さは強い正の相関がありますが、その比率（降伏比）は材料によって異なり、両者を独立して向上させることはできません。MIにおける多目的最適化は、このような複雑な関係性をデータから学習し、与えられた組成の範囲内で達成可能な特性の限界（パレートフロント）を明らかにすることを目的とします。

2. マルチターゲット回帰モデル

これまでは y （ターゲット）として単一の物性値を予測してきましたが、今回は y が「降伏強度」と「引張強さ」の2つの値を持つことになります。このような複数のターゲットを同時に予測するタスクを マルチターゲット回帰 （Multi-target Regression）と呼びます。

今回は、高性能な勾配ブースティングモデルである CatBoost を使用します。CatBoostは、損失関数に MultiRMSE などを指定することで、複数のターゲット変数を自然に扱うことができるため、複雑なラッパーを使わずにマルチターゲットモデルをシンプルに実装できます。

3. 多目的最適化とパレートフロント

単一の目的（例：強度を最大化）であれば、「最も高い値」が唯一の解となります。しかし、目的が複数（例：降伏強度を最大化 かつ 引張強さを最大化）になると、最適な解は一つに決まりません。そこで登場するのがパレート最適という考え方です。

• パレート支配: ある材料Aが、材料Bに対して全ての特性で同等以上であり、かつ、少なくとも一つの特性で明確に優れている場合、「AはBをパレート支配する」と言います。
• パレートフロント（Pareto Front）: どの材料にも支配されない、優秀な材料たちの集まりのことです。強度-強度グラフで言えば、この線上にある点は、それよりも右（高UTS）かつ上（高降伏強度）に他の点が存在しない、究極のトレードオフ曲線を描きます。

材料設計者は、このパレートフロントという「最高の選択肢をまとめたメニュー」の中から、製品の要求仕様（「降伏強度を重視」「バランスを重視」など）に応じて、最も適した材料候補を選ぶという、合理的な意思決定が可能になるのです。

実装方法

ワークフロー

matminer の steel_strength データセットを使い、2つの強度物性を同時に予測し、パレートフロントを可視化します。

# ===================================================================
# 0. 環境構築：必要なライブラリのインストール
# 1. 必要なライブラリのインポート
# 2. matminerから実データセットの読み込みと確認
# 3. 特徴量と複数ターゲットの定義
# 4. データを訓練用とテスト用に分割
# 5. マルチターゲット回帰モデルの訓練
# 6. 予測の実行と性能評価
# 7. パレートフロントの可視化
# 8. パレートフロント上の有望な材料組成の特定
# ===================================================================

実行手順

• 以下のコードブロック全体をコピーします。
• Google Colaboratoryの新しいセルに貼り付けます。
• セルが選択されていることを確認し、Shift キーと Enter キーを同時に押してコードを実行します。

# ===================================================================
# 0. 環境構築：必要なライブラリのインストール
# ===================================================================
!pip install matminer==0.9.3 scikit-learn==1.6.1
!pip install matplotlib==3.10.0 catboost==1.2.8
# ===================================================================
# 1. 必要なライブラリのインポート
# ===================================================================
import matplotlib.pyplot as plt
import numpy as np
from catboost import CatBoostRegressor
from matminer.datasets import load_dataset
from sklearn.metrics import r2_score
from sklearn.model_selection import train_test_split
# ===================================================================
# 2. matminerから実データセットの読み込みと確認
# ===================================================================
print("ステップ2: matminerから実データセットを読み込み、前処理します...")
df = load_dataset("steel_strength")
# 欠損値の確認と処理
print(f"元のデータ数: {len(df)}")
df_clean = df.dropna()
print(f"欠損値除去後のデータ数: {len(df_clean)}")
print("完了しました。")
# ===================================================================
# 3. 特徴量と複数ターゲットの定義
# ===================================================================
# 'composition', 'formula' とターゲット以外の列を特徴量とする
features = [ col for col in df_clean.columns if col not in [ "composition", "formula", "yield strength", "tensile strength", "elongation", ]
]
targets = ["yield strength", "tensile strength"]
X = df_clean[features]
y = df_clean[targets]
print("\n特徴量:", features)
print("ターゲット:", targets)
# ===================================================================
# 4. データを訓練用とテスト用に分割
# ===================================================================
print("\nステップ4: データを訓練データとテストデータに分割します...")
X_train, X_test, y_train, y_test = train_test_split( X, y, test_size=0.2, random_state=42
)
print("データの分割が完了しました。")
# ===================================================================
# 5. マルチターゲット回帰モデルの訓練
# ===================================================================
print("\nステップ5: マルチターゲット回帰モデルの学習を開始します...")
# CatBoostは 'MultiRMSE' を指定することでマルチターゲットをネイティブにサポート
model = CatBoostRegressor(random_seed=42, verbose=0, loss_function="MultiRMSE")
model.fit(X_train, y_train)
print("モデルの学習が完了しました。")
# ===================================================================
# 6. 予測の実行と性能評価
# ===================================================================
print("\nステップ6: 予測を実行し、モデル性能を評価します...")
y_test_pred = model.predict(X_test)
r2_yield = r2_score(y_test["yield strength"], y_test_pred[:, 0])
r2_uts = r2_score(y_test["tensile strength"], y_test_pred[:, 1])
print("\n" + "=" * 50)
print("【マルチターゲットモデル性能評価】")
print(f" 降伏強度のR2スコア (Test) : {r2_yield:.3f}")
print(f" 引張強さのR2スコア (Test) : {r2_uts:.3f}")
print("=" * 50 + "\n")
# ===================================================================
# 7. パレートフロントの可視化
# ===================================================================
print("ステップ7: 予測結果からパレートフロントを可視化します...")
# パレートフロントを特定する関数
def find_pareto_front_indices(points): pareto_indices = [] for i, p1 in enumerate(points): is_dominated = False for j, p2 in enumerate(points): if i == j: continue # p2がp1を支配しているかチェック (p2が両方の値で優れているか) if (p2[0] >= p1[0] and p2[1] > p1[1]) or (p2[0] > p1[0] and p2[1] >= p1[1]): is_dominated = True break if not is_dominated: pareto_indices.append(i) return pareto_indices
# テストデータでの予測結果からパレートフロントを計算
pareto_indices = find_pareto_front_indices(y_test_pred)
pareto_front_points = y_test_pred[pareto_indices]
# ソートして線で結べるようにする
pareto_front_points = pareto_front_points[np.argsort(pareto_front_points[:, 0])]
plt.figure(figsize=(10, 8))
plt.style.use("seaborn-v0_8-whitegrid")
# 全予測結果をプロット
plt.scatter( y_test_pred[:, 1], y_test_pred[:, 0], alpha=0.5, c="lightblue", label="Predicted Candidates",
)
# 実際のテストデータをプロット
plt.scatter( y_test["tensile strength"], y_test["yield strength"], alpha=0.5, c="gray", label="Actual Test Data",
)
# パレートフロントをプロット
plt.plot( pareto_front_points[:, 1], pareto_front_points[:, 0], "r-o", markersize=10, lw=3, label="Pareto Front",
)
plt.xlabel("Predicted Tensile Strength (MPa)", fontsize=14)
plt.ylabel("Predicted Yield Strength (MPa)", fontsize=14)
plt.title("Pareto Front for Yield Strength vs. Tensile Strength", fontsize=16)
plt.legend(fontsize=12)
plt.grid(True)
plt.show()
print("可視化処理が完了しました。")
# ===================================================================
# 8. パレートフロント上の有望な材料組成の特定
# ===================================================================
print("\nステップ8: パレートフロント上の有望な材料候補を特定します...")
# パレートフロントに対応する組成情報を表示
pareto_compositions_df = X_test.iloc[pareto_indices].copy()
pareto_compositions_df["predicted_yield_strength"] = y_test_pred[pareto_indices, 0]
pareto_compositions_df["predicted_tensile_strength"] = y_test_pred[pareto_indices, 1]
# 表示用にソートして抜粋
display_df = pareto_compositions_df.sort_values( "predicted_tensile_strength", ascending=False
)
print("\n【パレートフロント上の有望な材料候補】")
print( display_df[ [ "predicted_yield_strength", "predicted_tensile_strength", "c", "mn", "si", "cr", "ni", "mo", ] ].round(3)
)
print("\n処理はすべて完了しました。")

実行結果と考察

マルチターゲットモデルの性能

==================================================
【マルチターゲットモデル性能評価】 降伏強度のR2スコア (Test) : 0.839 引張強さのR2スコア (Test) : 0.862
==================================================

matminer の steel_strength データセットを用いて学習したCatBoostモデルは、未知のテストデータに対し、降伏強度（R²=0.839）と引張強さ（R²=0.862）の両方を高い精度で予測できています。決定係数（R²スコア）が1に近いほど良いモデルとされ、0.8を超える値は非常に良好な予測性能を示していると言えます。

パレートフロントの可視化と解釈

コードを実行して得られる上のグラフは、横軸に「予測された引張強さ」、縦軸に「予測された降伏強度」を取ります。

• 青い点群 (Predicted Candidates): モデルがテストデータの組成から予測した物性値の分布です。
• 灰色の点群 (Actual Test Data): 実際のテストデータの物性値の分布です。青い点群がこの灰色の点群の傾向をうまく捉えられていることが、モデルの精度の高さを視覚的に裏付けています。
• 赤い点（Pareto Front）: 予測結果の中で、他のどの点にも「支配」されていない、最も優れた性能を持つ点を示します。

なぜパレートフロントは「線」ではなく「1点」なのか？

ここで重要な疑問が生まれます。「なぜパレートフロントは線ではなく、たった1つの点なのか？」

これはエラーではなく、データセットの特性を反映した、あり得る正しい結果です。パレートフロントの定義は「他のどの点にも支配されていない、優れた点の集まり」です。今回の結果では、テストデータの中に、他のすべての点を圧倒的に「支配」してしまう、たった一つの非常に優秀な材料候補が存在していました。

グラフ上の赤い点は、他のどの青い点よりも右（高引張強さ）かつ上（高降伏強度）に位置しています。そのため、他の全ての点はこの1点に「支配」されてしまい、結果として「支配されていない点」はこの1点のみとなったのです。

もし、テストデータの中に、

• A材料：引張強さは一番高いが、降伏強度は2番手

• B材料：降伏強度は一番高いが、引張強さは2番手

  というような、異なる強みを持つ優秀な候補が複数存在していれば、それらが互いに「支配」しあわないため、複数の点がパレートフロントを形成し、「線」として描画されます。今回の結果は、このデータセットのテストサンプルの中に、たまたま一つだけ突出して優秀な「チャンピオンデータ」が含まれていたことを示しており、MIがその最強の1点を正しく見つけ出した、と解釈できます。

パレートフロント上の有望な材料組成

最終ステップでは、この「チャンピオンデータ」の組成を確認します。この材料は「降伏強度 約2426 MPa」「引張強さ 約2488 MPa」という、他を圧倒する強度特性を持つと予測されています。その組成（C: 0.0%, Ni: 18.1%, Mo: 3.6%など）は、マルエージング鋼のような高強度鋼の特徴を示唆しており、MIがデータの中から物理的に意味のある有望な候補を自動的に発見した好例と言えます。

コードの詳細解説

ステップ3: 特徴量と複数ターゲットの定義

# ===================================================================
# 3. 特徴量と複数ターゲットの定義
# ===================================================================
# 'composition', 'formula' とターゲット以外の列を特徴量とする
features = [ col for col in df_clean.columns if col not in [ "composition", "formula", "yield strength", "tensile strength", "elongation", ]
]
targets = ["yield strength", "tensile strength"]
X = df_clean[features]
y = df_clean[targets]
print("\n特徴量:", features)
print("ターゲット:", targets)

• featuresの定義: ここでは、データフレームの全列から、ターゲット（予測したい物性値）と予測に不要な情報（composition, formula, elongation）を除いたものを、特徴量（予測の根拠となる変数）としています。この書き方は、列が増減してもコードを修正する必要がない、柔軟性の高い方法です。
• targetsの定義: 今回の予測対象である「降伏強度」と「引張強さ」をリストとして定義します。これにより、**y**には2つの列を持つデータフレームが格納され、マルチターゲット学習の準備が整います。

ステップ5: マルチターゲット回帰モデルの訓練

# ===================================================================
# 5. マルチターゲット回帰モデルの訓練
# ===================================================================
print("\nステップ5: マルチターゲット回帰モデルの学習を開始します...")
# CatBoostは 'MultiRMSE' を指定することでマルチターゲットをネイティブにサポート
model = CatBoostRegressor(random_seed=42, verbose=0, loss_function="MultiRMSE")
model.fit(X_train, y_train)
print("モデルの学習が完了しました。")

• loss_function='MultiRMSE': ここがマルチターゲット学習の鍵です。CatBoostモデルを定義する際に、損失関数として MultiRMSE を指定しています。これは、複数のターゲット（この場合は2つ）のRMSE（二乗平均平方根誤差）を総合的に最小化するように学習を進めるための設定です。これにより、CatBoostは追加のラッパーなどを使うことなく、ネイティブに複数の物性を同時に予測するモデルを構築できます。
• model.fit(X_train, y_train): 学習コマンド自体はこれまでと同じですが、y_train が複数の列を持つデータフレームである点をCatBoostが認識し、マルチターゲット回帰を実行します。

ステップ7: パレートフロントの可視化

# ===================================================================
# 7. パレートフロントの可視化
# ===================================================================
print("ステップ7: 予測結果からパレートフロントを可視化します...")
# パレートフロントを特定する関数
def find_pareto_front_indices(points): pareto_indices = [] for i, p1 in enumerate(points): is_dominated = False for j, p2 in enumerate(points): if i == j: continue # p2がp1を支配しているかチェック (p2が両方の値で優れているか) if (p2[0] >= p1[0] and p2[1] > p1[1]) or (p2[0] > p1[0] and p2[1] >= p1[1]): is_dominated = True break if not is_dominated: pareto_indices.append(i) return pareto_indices
# テストデータでの予測結果からパレートフロントを計算
pareto_indices = find_pareto_front_indices(y_test_pred)
pareto_front_points = y_test_pred[pareto_indices]
# ソートして線で結べるようにする
pareto_front_points = pareto_front_points[np.argsort(pareto_front_points[:, 0])]
plt.figure(figsize=(10, 8))
plt.style.use("seaborn-v0_8-whitegrid")
# 全予測結果をプロット
plt.scatter( y_test_pred[:, 1], y_test_pred[:, 0], alpha=0.5, c="lightblue", label="Predicted Candidates",
)
# 実際のテストデータをプロット
plt.scatter( y_test["tensile strength"], y_test["yield strength"], alpha=0.5, c="gray", label="Actual Test Data",
)
# パレートフロントをプロット
plt.plot( pareto_front_points[:, 1], pareto_front_points[:, 0], "r-o", markersize=10, lw=3, label="Pareto Front",
)
plt.xlabel("Predicted Tensile Strength (MPa)", fontsize=14)
plt.ylabel("Predicted Yield Strength (MPa)", fontsize=14)
plt.title("Pareto Front for Yield Strength vs. Tensile Strength", fontsize=16)
plt.legend(fontsize=12)
plt.grid(True)
plt.show()
print("可視化処理が完了しました。")

• find_pareto_front_indices関数: この自作関数がパレートフロント特定の核心部です。
  • for i, p1 in enumerate(points): 予測された点（候補材料）を一つずつ取り出します（これを p1 とします）。
  • for j, p2 in enumerate(points): p1を他のすべての点（p2）と比較するためのループです。
  • 支配判定: if (p2 >= p1 and p2 > p1) or ... の条件式が「パレート支配」を判定しています。これは、「もし p2 が p1 に対して、両方の物性で同等以上、かつ、少なくとも片方の物性で明確に優れているならば」という条件を意味します。
  • is_dominated = True: もし p1 を支配する p2 が一つでも見つかれば、p1 はパレートフロント上の点ではないと判断し、内側のループを抜けます。
  • if not is_dominated: 内側のループを最後まで回りきってもどの点にも支配されなかった場合、その点 p1 はパレートフロント上の点であるとみなし、そのインデックスを pareto_indices リストに追加します。
• この関数により、「どの点にも支配されない優秀な点のインデックス」のリストが得られます。

ステップ8: パレートフロント上の有望な材料組成の特定

# ===================================================================
# 8. パレートフロント上の有望な材料組成の特定
# ===================================================================
print("\nステップ8: パレートフロント上の有望な材料候補を特定します...")
# パレートフロントに対応する組成情報を表示
pareto_compositions_df = X_test.iloc[pareto_indices].copy()
pareto_compositions_df["predicted_yield_strength"] = y_test_pred[pareto_indices, 0]
pareto_compositions_df["predicted_tensile_strength"] = y_test_pred[pareto_indices, 1]
# 表示用にソートして抜粋
display_df = pareto_compositions_df.sort_values( "predicted_tensile_strength", ascending=False
)
print("\n【パレートフロント上の有望な材料候補】")
print( display_df[ [ "predicted_yield_strength", "predicted_tensile_strength", "c", "mn", "si", "cr", "ni", "mo", ] ].round(3)
)
print("\n処理はすべて完了しました。")

• .iloc[pareto_indices]: 前のステップで得られた「優秀な点のインデックス」のリストを使い、元のテストデータ X_test から該当する行を抽出します。これにより、パレートフロント上の点が、具体的にどのような化学組成（特徴量）を持っていたのかを特定することができます。この操作こそが、予測結果を具体的な材料設計指針へと結びつける重要なステップです。
• .copy() : 抽出したデータに対して後から列を追加するため、元のデータフレームに影響を与えないよう、.copy() メソッドで明示的にコピーを作成しておくのが安全な作法です。

最後に

今回は、単一の物性予測から一歩踏み出し、matminer の実際の公開データセットを用いて「降伏強度」と「引張強さ」という複数の物性を同時に最適化する、より現実的な課題に挑戦しました。マルチターゲット回帰モデルを構築し、その予測から パレートフロント を導き出すことで、無数の選択肢の中から「最良のバランスを持つ候補群」を科学的に絞り込む手法を学びました。

また、パレートフロントが必ずしも「線」ではなく「点」になり得ることも確認し、その背景にあるデータの特性を考察しました。

次回は、今回構築した予測モデルを「探索エンジン」として活用し、未知の有望材料を発見するプロセスに挑戦します。膨大な数の仮想的な材料組成をランダムに生成し、モデルで高速に評価する「ハイスループット・バーチャルスクリーニング（HTVS）」を実践します。このプロセスを通じて、今回の「最強の1点」を超える、あるいは異なる強みを持つ新たな候補を発見し、より豊かなパレートフロントを形成していく過程をぜひご覧ください。