100本ノック / 行列 / 行列100本ノック
製造業の品質予測をPythonで学ぶ｜Ridge・Logistic回帰・Attention入門
加工条件から品質リスクを予測し、改善条件を決める
最適化・回帰・Attentionで学ぶ製造業の行列分析：100本ノック No.051〜No.060
製造現場では、加工速度、送り量、温度、振動、工具摩耗などが品質へ同時に影響します。本記事では、架空の精密部品ラインを題材に、品質誤差を予測し、不良リスクを定量化し、計算を安定させ、重点確認すべき工程を絞るまでを行列の視点で扱います。
[!NOTE] 本資料は、数理工房 (もしくは代表である和山個人) が過去に企業研修において使用した notebook を企業様の許可を得て再構成・編集のうえ公開しています。
掲載データはすべて架空のものであり、実在する企業・工場・数値とは一切関係ありません。
はじめに：この記事で扱う製造業の実務課題
架空の工場では、精密シャフトをロット単位で加工し、寸法誤差と最終検査の合否を記録しています。製造技術者は、加工条件と設備状態から寸法誤差を予測し、不良確率の高いロットを検査へ優先的に回し、条件見直しの根拠を作りたいと考えています。
重要なのは、予測値を作ることだけではありません。どの計算法が安定しているか、非線形な品質変化を捉えられるか、どの工程時点に注意を向けるべきかまで説明できて、初めて現場の意思決定に使えます。
現場でよくある状況
- 加工速度と送り量が連動し、説明変数同士の相関が強い
- 寸法誤差のような連続値と、不良・良品のような二値判断が混在する
- 表計算ソフトでは係数が得られても、反復計算の収束や数値安定性が見えにくい
- 温度と振動が同時に高いときだけ品質が悪化するなど、線形モデルでは表しにくい関係がある
- 時系列センサーが長く、担当者が全時点を同じ密度で確認できない
なぜこの問題は判断が難しいのか
品質モデルは、予測精度だけでなく、説明変数の重複、損失関数の曲率、反復法の収束、閾値による見逃しと過検出のトレードオフに左右されます。さらに、高精度な非線形モデルほど、適用範囲や判断根拠を現場へ説明する設計が必要です。
そこで本記事では、同じ架空データに複数の行列計算法を適用します。単に手法を並べるのではなく、線形予測、不良確率、安定な大規模計算、非線形補正、工程時系列の重点化という判断の流れとして比較します。
今回扱うノックの全体像
| No. | テーマ | 製造業での判断 |
|---|---|---|
| 051 | 勾配降下法 | 品質誤差を小さくする係数を反復的に学習する |
| 052 | 最小二乗法 | 観測誤差を最小にする品質予測式を求める |
| 053 | 正規方程式 | 行列式から回帰係数を直接計算し、前提を確認する |
| 054 | Ridge回帰 | 相関の強い加工条件でも係数を安定させる |
| 055 | Logistic回帰 | ロットごとの不良確率を推定する |
| 056 | Newton法 | 曲率を利用して不良モデルを効率よく更新する |
| 057 | ヘッセ行列 | 学習の不安定方向と識別しにくい条件を診断する |
| 058 | 共役勾配法 | 大規模な連立方程式を逆行列なしで解く |
| 059 | カーネル法 | 温度と振動の非線形な品質関係を捉える |
| 060 | Attention | 工程時系列から重点確認すべき時点を可視化する |
Python 環境の準備
NumPyで行列計算、pandasで表の確認、Matplotlibで可視化を行います。Logistic回帰の評価指標だけscikit-learnを利用します。外部データには依存せず、乱数生成器は np.random.default_rng(42) で固定します。
import sys
import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
import japanize_matplotlib
from sklearn.metrics import confusion_matrix, roc_auc_score
from IPython.display import display
np.set_printoptions(precision=4, suppress=True)
pd.set_option("display.precision", 4)
print("Python :", sys.version.split()[0])
print("NumPy :", np.__version__)
print("pandas :", pd.__version__)
print("Matplotlib :", matplotlib.__version__)Python : 3.13.1
NumPy : 2.5.1
pandas : 3.0.3
Matplotlib : 3.11.0架空データの作成
240ロットについて、加工速度、送り量、加工温度、振動、油圧、工具摩耗度を生成します。寸法誤差、表面粗さ、不良フラグを品質結果とし、表面粗さには温度・振動の非線形効果を含めます。不良フラグは設備状態から決まる確率に従って生成します。
先頭180ロットを学習、残り60ロットを評価に使います。実務では時系列順の分割、設備別の外部検証、保全前後の分離が必要ですが、ここでは各行列手法の違いを読みやすくするため固定分割とします。
rng = np.random.default_rng(42)
n = 240
speed = rng.normal(1800, 120, n)
feed = 0.18 + 0.00028 * (speed - 1800) + rng.normal(0, 0.018, n)
temperature = 61 + 0.006 * (speed - 1800) + 16 * (feed - 0.18) + rng.normal(0, 1.4, n)
vibration = 1.5 + 0.0018 * (speed - 1800) + 2.2 * (feed - 0.18) + rng.normal(0, 0.22, n)
pressure = 5.1 - 0.9 * (feed - 0.18) + rng.normal(0, 0.13, n)
tool_wear = rng.uniform(0, 1, n)
dim_error = (
0.4 + 0.0015 * (speed - 1800) + 4.0 * (feed - 0.18)
+ 0.12 * (temperature - 61) + 0.55 * (vibration - 1.5)
+ 1.25 * tool_wear**2
+ 0.10 * np.maximum(temperature - 62, 0) * np.maximum(vibration - 1.55, 0)
+ rng.normal(0, 0.38, n)
)
z_temp = (temperature - temperature.mean()) / temperature.std()
z_vib = (vibration - vibration.mean()) / vibration.std()
logit_true = -2.5 + 0.65 * z_temp + 0.95 * z_vib + 1.55 * tool_wear + 0.45 * z_temp * z_vib
defect_prob_true = 1 / (1 + np.exp(-logit_true))
defect = rng.binomial(1, defect_prob_true)
surface_roughness = (
0.75 + 0.28 * z_temp**2 + 0.35 * z_vib**2
+ 0.25 * np.maximum(z_temp, 0) * np.maximum(z_vib, 0)
+ rng.normal(0, 0.12, n)
)
df = pd.DataFrame({
"ロット": [f"L{i:03d}" for i in range(1, n + 1)],
"加工速度_rpm": speed, "送り量_mm_rev": feed, "加工温度_C": temperature,
"振動_mm_s": vibration, "油圧_MPa": pressure, "工具摩耗度": tool_wear,
"寸法誤差_um": dim_error, "表面粗さ_Ra": surface_roughness, "不良": defect,
})
display(df.head().round(3))
print(f"全体不良率: {df['不良'].mean():.1%}")| ロット | 加工速度_rpm | 送り量_mm_rev | 加工温度_C | 振動_mm_s | 油圧_MPa | 工具摩耗度 | 寸法誤差_um | 表面粗さ_Ra | 不良 | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | L001 | 1836.566 | 0.174 | 61.443 | 1.532 | 5.315 | 0.207 | 0.390 | 0.873 | 0 |
| 1 | L002 | 1675.202 | 0.143 | 61.671 | 1.443 | 5.150 | 0.423 | -0.784 | 0.475 | 0 |
| 2 | L003 | 1890.054 | 0.174 | 61.566 | 1.146 | 4.975 | 0.176 | -0.053 | 1.287 | 0 |
| 3 | L004 | 1912.868 | 0.185 | 62.573 | 1.385 | 5.081 | 0.135 | 1.164 | 0.700 | 0 |
| 4 | L005 | 1565.876 | 0.153 | 59.080 | 0.816 | 5.120 | 0.860 | 0.006 | 2.269 | 0 |
全体不良率: 20.4%feature_cols = ["加工速度_rpm", "送り量_mm_rev", "加工温度_C", "振動_mm_s", "油圧_MPa", "工具摩耗度"]
train_idx = np.arange(180)
test_idx = np.arange(180, 240)
X_raw = df[feature_cols].to_numpy()
y = df["寸法誤差_um"].to_numpy()
y_cls = df["不良"].to_numpy()
mean_x = X_raw[train_idx].mean(axis=0)
std_x = X_raw[train_idx].std(axis=0)
X_std = (X_raw - mean_x) / std_x
X = np.column_stack([np.ones(n), X_std])
X_train, X_test = X[train_idx], X[test_idx]
y_train, y_test = y[train_idx], y[test_idx]
summary = df.loc[train_idx, feature_cols + ["寸法誤差_um", "表面粗さ_Ra", "不良"]].describe().T[["mean", "std", "min", "max"]]
display(summary.round(3))| mean | std | min | max | |
|---|---|---|---|---|
| 加工速度_rpm | 1793.188 | 104.082 | 1544.154 | 2149.663 |
| 送り量_mm_rev | 0.179 | 0.033 | 0.080 | 0.276 |
| 加工温度_C | 60.848 | 1.788 | 56.053 | 66.787 |
| 振動_mm_s | 1.475 | 0.343 | 0.470 | 2.594 |
| 油圧_MPa | 5.109 | 0.113 | 4.805 | 5.460 |
| 工具摩耗度 | 0.504 | 0.301 | 0.005 | 0.997 |
| 寸法誤差_um | 0.738 | 0.873 | -1.438 | 3.218 |
| 表面粗さ_Ra | 1.454 | 1.002 | 0.475 | 8.842 |
| 不良 | 0.189 | 0.393 | 0.000 | 1.000 |
No.051：勾配降下法
実務での意味
勾配降下法は、予測誤差が小さくなる方向へ係数を少しずつ更新する方法です。データ量や変数数が多く、全データを一度に直接計算しにくい品質予測や機械学習の基礎になります。
分析・モデル化の考え方
寸法誤差の予測を 、平均二乗誤差の半分を
とすると、勾配は です。更新式 を繰り返します。学習率 が大きすぎると発散し、小さすぎると収束が遅くなります。変数の標準化は収束速度を揃えるためにも重要です。
Pythonで確認する
beta_gd = np.zeros(X_train.shape[1])
learning_rate = 0.08
loss_history = []
for _ in range(600):
residual = X_train @ beta_gd - y_train
loss_history.append(np.mean(residual**2) / 2)
gradient = X_train.T @ residual / len(y_train)
beta_gd -= learning_rate * gradient
pred_gd = X_test @ beta_gd
rmse_gd = np.sqrt(np.mean((y_test - pred_gd) ** 2))
print(f"最終損失: {loss_history[-1]:.5f}")
print(f"評価データRMSE: {rmse_gd:.3f} μm")
plt.figure(figsize=(7, 4))
plt.plot(loss_history)
plt.title("勾配降下法の収束")
plt.xlabel("反復回数")
plt.ylabel("損失（MSE / 2）")
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()最終損失: 0.08004
評価データRMSE: 0.414 μm
結果の読み取り
損失は反復とともに低下し、一定値へ収束します。現場導入時は最終精度だけでなく、損失曲線が振動・発散していないか、反復上限に達する前に改善が止まったかをログへ残します。ライン追加後に収束挙動が変わった場合は、単位、外れ値、データ分布の変化を疑う材料になります。
No.052：最小二乗法
実務での意味
最小二乗法は、実測値と予測値のずれの二乗和を最小にする係数を求めます。加工条件から寸法、サイクルタイム、電力使用量などの連続KPIを説明・予測する基本手法です。
分析・モデル化の考え方
目的は
です。np.linalg.lstsq はSVDを用いて解くため、明示的に逆行列を作るより数値的に扱いやすい方法です。RMSEは寸法誤差と同じ単位で読めますが、大きな誤差を強く罰する指標である点も理解して使います。
Pythonで確認する
beta_ls, residuals, rank, singular_values = np.linalg.lstsq(X_train, y_train, rcond=None)
pred_ls = X_test @ beta_ls
rmse_ls = np.sqrt(np.mean((y_test - pred_ls) ** 2))
coef_table = pd.DataFrame({
"項目": ["切片"] + feature_cols,
"標準化係数": beta_ls,
})
display(coef_table.round(4))
print(f"行列ランク: {rank} / {X_train.shape[1]}")
print(f"評価データRMSE: {rmse_ls:.3f} μm")| 項目 | 標準化係数 | |
|---|---|---|
| 0 | 切片 | 0.7377 |
| 1 | 加工速度_rpm | 0.1792 |
| 2 | 送り量_mm_rev | 0.1658 |
| 3 | 加工温度_C | 0.2595 |
| 4 | 振動_mm_s | 0.1257 |
| 5 | 油圧_MPa | -0.0331 |
| 6 | 工具摩耗度 | 0.3758 |
行列ランク: 7 / 7
評価データRMSE: 0.414 μm結果の読み取り
係数の絶対値は、標準化後の各条件が寸法誤差とどの程度連動するかを比較する目安になります。ただし、観測データ上の関連であり、加工条件を変更したときの因果効果を保証しません。評価RMSEは測定器の分解能や公差幅と比較し、「精度が高い」ではなく、検査省略や条件調整に足りる誤差かを判断します。
No.053：正規方程式
実務での意味
正規方程式は、最小二乗法が行列の連立方程式として解けることを示します。係数の算出根拠を理解し、変数の重複によって解が不安定になる理由を把握するために重要です。
分析・モデル化の考え方
損失の勾配を0と置くと
を得ます。 が正則なら と書けますが、実装では逆行列を明示せず solve を使います。正規方程式は条件数を悪化させるため、大規模・悪条件の問題ではSVDやQR分解を優先します。
Pythonで確認する
gram = X_train.T @ X_train
rhs = X_train.T @ y_train
beta_normal = np.linalg.solve(gram, rhs)
comparison = pd.DataFrame({
"項目": ["切片"] + feature_cols,
"lstsq": beta_ls,
"正規方程式": beta_normal,
"差": beta_normal - beta_ls,
})
display(comparison.round(8))
print(f"X の条件数: {np.linalg.cond(X_train):.1f}")
print(f"X^T X の条件数: {np.linalg.cond(gram):.1f}")| 項目 | lstsq | 正規方程式 | 差 | |
|---|---|---|---|---|
| 0 | 切片 | 0.7377 | 0.7377 | 0.0 |
| 1 | 加工速度_rpm | 0.1792 | 0.1792 | 0.0 |
| 2 | 送り量_mm_rev | 0.1658 | 0.1658 | 0.0 |
| 3 | 加工温度_C | 0.2595 | 0.2595 | -0.0 |
| 4 | 振動_mm_s | 0.1257 | 0.1257 | -0.0 |
| 5 | 油圧_MPa | -0.0331 | -0.0331 | 0.0 |
| 6 | 工具摩耗度 | 0.3758 | 0.3758 | 0.0 |
X の条件数: 4.4
X^T X の条件数: 19.1結果の読み取り
このデータでは両手法の係数はほぼ一致します。一方、 の条件数は より大きくなります。実務では計算式が理論上正しいだけで採用せず、変数数、相関、精度要件を確認します。説明変数追加後に条件数が急増した場合は、そのKPIが新しい情報を持つかを見直します。
No.054：Ridge回帰
実務での意味
設備設定値と実測値、元センサーと派生KPIなど、似た変数を同時に入れると回帰係数が大きく揺れます。Ridge回帰は係数を縮小し、予測と運用を安定させます。
分析・モデル化の考え方
Ridge回帰は
を解き、 となります。ここでは送り量にほぼ同じ派生指標を追加し、相関を意図的に強くします。切片は正則化しません。 は評価データや交差検証で選びます。
Pythonで確認する
feed_proxy = X_std[:, 1] + np.random.default_rng(0).normal(0, 0.001, n)
X_col = np.column_stack([X, feed_proxy])
Xc_train, Xc_test = X_col[train_idx], X_col[test_idx]
penalty = np.eye(Xc_train.shape[1])
penalty[0, 0] = 0
ridge_rows = []
ridge_models = {}
for lam in [0, 0.01, 0.1, 1, 10, 100]:
beta = np.linalg.solve(Xc_train.T @ Xc_train + lam * penalty, Xc_train.T @ y_train)
ridge_models[lam] = beta
ridge_rows.append({
"lambda": lam,
"評価RMSE": np.sqrt(np.mean((y_test - Xc_test @ beta) ** 2)),
"係数ノルム": np.linalg.norm(beta[1:]),
"送り量係数": beta[2],
"派生送り係数": beta[-1],
})
ridge_result = pd.DataFrame(ridge_rows)
display(ridge_result.round(4))
plt.figure(figsize=(7, 4))
plt.semilogx(ridge_result.loc[1:, "lambda"], ridge_result.loc[1:, "評価RMSE"], marker="o")
plt.title("Ridge正則化強度と評価誤差")
plt.xlabel("lambda")
plt.ylabel("評価RMSE（μm）")
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()| lambda | 評価RMSE | 係数ノルム | 送り量係数 | 派生送り係数 | |
|---|---|---|---|---|---|
| 0 | 0.00 | 0.4139 | 10.4952 | -7.3286 | 7.4956 |
| 1 | 0.01 | 0.4140 | 0.5279 | 0.0223 | 0.1436 |
| 2 | 0.10 | 0.4140 | 0.5207 | 0.0769 | 0.0891 |
| 3 | 1.00 | 0.4144 | 0.5186 | 0.0833 | 0.0845 |
| 4 | 10.00 | 0.4182 | 0.5000 | 0.0906 | 0.0907 |
| 5 | 100.00 | 0.4624 | 0.3881 | 0.1064 | 0.1064 |
結果の読み取り
正則化なしでは、ほぼ同じ2つの送り指標へ係数が不安定に分配されます。適度なRidgeは係数ノルムを抑えつつ評価誤差を維持します。大きすぎる は必要な信号まで縮小します。係数を原因寄与として説明する場合は、Ridgeを使っても重複変数の意味を整理し、選んだ と評価方法を記録します。
No.055：Logistic回帰
実務での意味
不良・良品の二値を直接予測するだけでなく、不良確率として出力すると、全数検査、追加検査、通常出荷など複数の運用閾値を設計できます。
分析・モデル化の考え方
Logistic回帰は
で不良確率を表します。交差エントロピー損失を最小化し、ここでは勾配降下法で学習します。確率閾値は0.5固定ではなく、見逃しコスト、追加検査能力、不良率によって決めます。
Pythonで確認する
def sigmoid(z):
z = np.clip(z, -30, 30)
return 1 / (1 + np.exp(-z))
beta_logit = np.zeros(X_train.shape[1])
logloss_history = []
for _ in range(1500):
p = sigmoid(X_train @ beta_logit)
logloss = -np.mean(y_cls[train_idx] * np.log(p + 1e-12) + (1 - y_cls[train_idx]) * np.log(1 - p + 1e-12))
logloss_history.append(logloss)
beta_logit -= 0.08 * (X_train.T @ (p - y_cls[train_idx]) / len(train_idx))
p_test = sigmoid(X_test @ beta_logit)
threshold = 0.35
pred_cls = (p_test >= threshold).astype(int)
cm = confusion_matrix(y_cls[test_idx], pred_cls)
tn, fp, fn, tp = cm.ravel()
metrics = pd.Series({
"AUC": roc_auc_score(y_cls[test_idx], p_test),
"再現率（不良捕捉率）": tp / (tp + fn),
"適合率": tp / (tp + fp),
"追加検査率": pred_cls.mean(),
})
display(metrics.to_frame("値").round(3))
display(pd.DataFrame(cm, index=["実際_良品", "実際_不良"], columns=["予測_良品", "予測_不良"]))| 値 | |
|---|---|
| AUC | 0.776 |
| 再現率（不良捕捉率） | 0.533 |
| 適合率 | 0.615 |
| 追加検査率 | 0.217 |
| 予測_良品 | 予測_不良 | |
|---|---|---|
| 実際_良品 | 40 | 5 |
| 実際_不良 | 7 | 8 |
結果の読み取り
AUCはロットのリスク順位付け能力、再現率は実不良をどれだけ拾えたか、適合率は追加検査の効率を表します。閾値0.35では検査対象を広めに取り、見逃しを抑える設計です。実務では不良件数が少ない期間ほど指標が揺れるため、信頼区間、品種別評価、確率校正、閾値変更後の検査負荷を合わせて確認します。
No.056：Newton法
実務での意味
Newton法は、勾配だけでなく損失関数の曲率を使って更新幅を調整します。反復回数を減らせる一方、各更新で連立方程式を解くため、変数数が多い場合は計算負荷と安定性の設計が必要です。
分析・モデル化の考え方
勾配 とヘッセ行列 を用い、
と更新します。Logistic回帰では 、 です。実装では逆行列を作らず solve を使い、数値安定化のため小さな対角項を加えます。
Pythonで確認する
beta_newton = np.zeros(X_train.shape[1])
newton_history = []
for iteration in range(10):
p = sigmoid(X_train @ beta_newton)
loss = -np.mean(y_cls[train_idx] * np.log(p + 1e-12) + (1 - y_cls[train_idx]) * np.log(1 - p + 1e-12))
newton_history.append(loss)
gradient = X_train.T @ (p - y_cls[train_idx])
weights = p * (1 - p)
hessian = X_train.T @ (weights[:, None] * X_train) + 1e-6 * np.eye(X_train.shape[1])
beta_newton -= np.linalg.solve(hessian, gradient)
display(pd.DataFrame({
"反復": np.arange(1, len(newton_history) + 1),
"Newton法の損失": newton_history,
}).round(6))
print(f"勾配降下法の最終損失: {logloss_history[-1]:.6f}")
print(f"Newton法の最終損失: {newton_history[-1]:.6f}")| 反復 | Newton法の損失 | |
|---|---|---|
| 0 | 1 | 0.6931 |
| 1 | 2 | 0.3831 |
| 2 | 3 | 0.3437 |
| 3 | 4 | 0.3378 |
| 4 | 5 | 0.3376 |
| 5 | 6 | 0.3376 |
| 6 | 7 | 0.3376 |
| 7 | 8 | 0.3376 |
| 8 | 9 | 0.3376 |
| 9 | 10 | 0.3376 |
勾配降下法の最終損失: 0.337600
Newton法の最終損失: 0.337592結果の読み取り
Newton法は少ない反復で損失を下げます。ただし「反復回数が少ない」ことと「計算時間が短い」ことは同じではありません。変数が多いモデルではヘッセ行列の作成・求解が支配的になります。データが完全分離に近い、または説明変数が重複する場合は不安定になるため、正則化やステップ幅制御を組み合わせます。
No.057：ヘッセ行列
実務での意味
ヘッセ行列は、損失が各パラメータ方向にどれだけ曲がっているかを表します。小さな固有値方向は、データから係数を識別しにくく、わずかなデータ変更で推定が揺れやすい方向です。
分析・モデル化の考え方
ヘッセ行列は2階偏微分を並べた行列
です。固有値がすべて正なら局所的に凸です。最大固有値と最小固有値の比が大きいと、方向ごとの曲率差が大きく、最適化が進みにくくなります。Ridgeの対角項は小さい固有値を押し上げます。
Pythonで確認する
p_col = sigmoid(Xc_train @ np.pad(beta_newton, (0, 1)))
W_col = p_col * (1 - p_col)
H_col = Xc_train.T @ (W_col[:, None] * Xc_train)
H_ridge = H_col + 1.0 * penalty
eig_plain = np.linalg.eigvalsh(H_col)
eig_ridge = np.linalg.eigvalsh(H_ridge)
hessian_table = pd.DataFrame({
"指標": ["最小固有値", "最大固有値", "条件数"],
"正則化なし": [eig_plain.min(), eig_plain.max(), np.linalg.cond(H_col)],
"Ridgeあり": [eig_ridge.min(), eig_ridge.max(), np.linalg.cond(H_ridge)],
})
display(hessian_table.round(4))
plt.figure(figsize=(7, 4))
plt.semilogy(np.arange(1, len(eig_plain) + 1), eig_plain, marker="o", label="正則化なし")
plt.semilogy(np.arange(1, len(eig_ridge) + 1), eig_ridge, marker="s", label="Ridgeあり")
plt.title("ヘッセ行列の固有値")
plt.xlabel("固有値の順位（昇順）")
plt.ylabel("固有値（対数軸）")
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()| 指標 | 正則化なし | Ridgeあり | |
|---|---|---|---|
| 0 | 最小固有値 | 0.0000e+00 | 1.0000 |
| 1 | 最大固有値 | 7.6713e+01 | 77.5857 |
| 2 | 条件数 | 8.1267e+06 | 77.5850 |
結果の読み取り
重複に近い送り指標を含むと最小固有値が小さくなり、損失面に平らな方向ができます。Ridgeはその方向の曲率を加えて条件数を改善します。現場では、条件数の悪化を単なる計算問題とせず、「ほぼ同じ意味のKPIを重複管理していないか」「独立した試験条件が不足していないか」というデータ設計の問題として扱います。
No.058：共役勾配法
実務での意味
設備数・品種数・特徴量が増えると、連立方程式の行列を直接分解するコストが大きくなります。共役勾配法は、対称正定値行列に対し、逆行列を作らず行列ベクトル積の反復で解を求めます。
分析・モデル化の考え方
Ridge回帰の正規方程式
は正則化により対称正定値になります。共役勾配法は互いに 共役な探索方向を使い、二次形式を効率よく最小化します。収束判定には残差 を使います。
Pythonで確認する
def conjugate_gradient(A, b, tol=1e-10, max_iter=100):
x = np.zeros_like(b)
r = b - A @ x
p = r.copy()
residual_norms = [np.linalg.norm(r)]
for _ in range(max_iter):
Ap = A @ p
alpha = (r @ r) / (p @ Ap)
x = x + alpha * p
r_new = r - alpha * Ap
residual_norms.append(np.linalg.norm(r_new))
if residual_norms[-1] < tol:
break
beta = (r_new @ r_new) / (r @ r)
p = r_new + beta * p
r = r_new
return x, residual_norms
lam = 1.0
A_cg = Xc_train.T @ Xc_train + lam * penalty
b_cg = Xc_train.T @ y_train
beta_cg, cg_residuals = conjugate_gradient(A_cg, b_cg)
beta_direct = np.linalg.solve(A_cg, b_cg)
print(f"反復回数: {len(cg_residuals) - 1}")
print(f"直接法との差: {np.linalg.norm(beta_cg - beta_direct):.3e}")
plt.figure(figsize=(7, 4))
plt.semilogy(cg_residuals, marker="o")
plt.title("共役勾配法の残差収束")
plt.xlabel("反復回数")
plt.ylabel("残差ノルム")
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()反復回数: 9
直接法との差: 1.491e-14
結果の読み取り
残差は少ない反復で低下し、直接法に近い解が得られます。今回の小さな行列では直接法で十分ですが、共役勾配法は行列を明示的に保持しにくい大規模問題で有効です。実務では正定値性、許容誤差、最大反復、前処理の有無を管理し、未収束時に古い係数を配信しない運用を設計します。
No.059：カーネル法
実務での意味
温度や振動が基準から離れた領域、あるいは両方が高い領域で表面粗さが急に悪化することがあります。カーネル法は、元の変数から多数の非線形特徴を明示的に作らず、似た運転条件同士の関係から非線形な予測を行います。
分析・モデル化の考え方
RBFカーネル
で学習点同士の類似度行列 を作ります。Kernel Ridgeでは
を解き、新しい点とのカーネル値から予測します。 は影響範囲、 は滑らかさを制御し、学習範囲外への外挿は慎重に扱います。
Pythonで確認する
kernel_cols = [2, 3] # 標準化後の加工温度、振動
Z_train = X_std[train_idx][:, kernel_cols]
Z_test = X_std[test_idx][:, kernel_cols]
def rbf_kernel(A, B, gamma=0.7):
sq_dist = ((A[:, None, :] - B[None, :, :]) ** 2).sum(axis=2)
return np.exp(-gamma * sq_dist)
rough_train = df.loc[train_idx, "表面粗さ_Ra"].to_numpy()
rough_test = df.loc[test_idx, "表面粗さ_Ra"].to_numpy()
K_train = rbf_kernel(Z_train, Z_train, gamma=0.1)
alpha = np.linalg.solve(K_train + 0.01 * np.eye(len(train_idx)), rough_train)
pred_kernel = rbf_kernel(Z_test, Z_train, gamma=0.1) @ alpha
Z_linear = np.column_stack([np.ones(len(train_idx)), Z_train])
coef_linear = np.linalg.lstsq(Z_linear, rough_train, rcond=None)[0]
pred_linear_2d = np.column_stack([np.ones(len(test_idx)), Z_test]) @ coef_linear
kernel_compare = pd.DataFrame({
"モデル": ["温度・振動の線形回帰", "RBF Kernel Ridge"],
"評価RMSE": [
np.sqrt(np.mean((rough_test - pred_linear_2d) ** 2)),
np.sqrt(np.mean((rough_test - pred_kernel) ** 2)),
],
})
display(kernel_compare.round(4))
plt.figure(figsize=(6, 5))
plt.scatter(rough_test, pred_kernel, c=df.loc[test_idx, "工具摩耗度"], cmap="viridis", alpha=0.8)
limits = [min(rough_test.min(), pred_kernel.min()), max(rough_test.max(), pred_kernel.max())]
plt.plot(limits, limits, "--", color="gray", label="実測=予測")
plt.title("カーネル法による表面粗さ予測")
plt.xlabel("実測表面粗さ（Ra）")
plt.ylabel("予測表面粗さ（Ra）")
plt.grid(True, alpha=0.3)
plt.colorbar(label="工具摩耗度")
plt.legend()
plt.tight_layout()
plt.show()| モデル | 評価RMSE | |
|---|---|---|
| 0 | 温度・振動の線形回帰 | 0.8336 |
| 1 | RBF Kernel Ridge | 0.1498 |
結果の読み取り
カーネル法は温度と振動の曲線的な関係を捉え、同じ2変数の線形回帰より表面粗さの評価誤差を小さくします。ただし、この2変数だけでは工具摩耗などの影響が残ります。色で示した摩耗度に偏りが見えるなら追加候補です。精度向上だけで変数を増やさず、学習データから離れた条件を検知し、適用範囲外では予測を保留する仕組みが必要です。
No.060：Attention
実務での意味
1ロットの加工中に多数の時系列データがあると、すべての時点を同じ重みで平均するだけでは、異常兆候が強い短時間を薄めてしまいます。Attentionは、現在の判断目的に関連する時点へ大きな重みを割り当てます。
分析・モデル化の考え方
Scaled Dot-Product Attentionは
です。 は探したい状態、 は各時点の特徴、 は集約したい情報を表します。ここでは学習済みTransformerではなく、温度・振動・負荷が高い状態をクエリとする簡略例で、重みの意味を確認します。重みは因果的な故障原因を証明するものではありません。
Pythonで確認する
stages = np.arange(1, 13)
stage_temp = np.array([59.8, 60.2, 60.7, 61.0, 61.4, 62.0, 62.6, 63.7, 64.2, 63.4, 62.8, 62.1])
stage_vib = np.array([1.20, 1.24, 1.28, 1.31, 1.35, 1.42, 1.51, 1.82, 2.05, 1.72, 1.55, 1.44])
stage_load = np.array([0.40, 0.43, 0.48, 0.52, 0.56, 0.61, 0.67, 0.91, 1.00, 0.78, 0.65, 0.55])
stage_matrix = np.column_stack([stage_temp, stage_vib, stage_load])
keys = (stage_matrix - stage_matrix.mean(axis=0)) / stage_matrix.std(axis=0)
query = np.array([0.7, 1.0, 0.6])
scores = keys @ query / np.sqrt(keys.shape[1])
weights = np.exp(scores - scores.max())
weights /= weights.sum()
values = 0.4 * keys[:, 0] + 0.8 * keys[:, 1] + 0.5 * keys[:, 2]
attention_summary = weights @ values
attention_table = pd.DataFrame({
"工程時点": stages, "温度_C": stage_temp, "振動_mm_s": stage_vib,
"負荷指数": stage_load, "Attention重み": weights,
}).sort_values("Attention重み", ascending=False)
display(attention_table.head(5).round(4))
print(f"Attentionで集約したリスク表現: {attention_summary:.3f}")
plt.figure(figsize=(8, 4))
plt.bar(stages, weights, color=np.where(weights >= np.quantile(weights, 0.75), "tomato", "steelblue"))
plt.title("工程時点ごとのAttention重み")
plt.xlabel("工程時点")
plt.ylabel("Attention重み")
plt.xticks(stages)
plt.grid(True, axis="y", alpha=0.3)
plt.tight_layout()
plt.show()| 工程時点 | 温度_C | 振動_mm_s | 負荷指数 | Attention重み | |
|---|---|---|---|---|---|
| 8 | 9 | 64.2 | 2.05 | 1.00 | 0.4829 |
| 7 | 8 | 63.7 | 1.82 | 0.91 | 0.2036 |
| 9 | 10 | 63.4 | 1.72 | 0.78 | 0.1143 |
| 10 | 11 | 62.8 | 1.55 | 0.65 | 0.0498 |
| 6 | 7 | 62.6 | 1.51 | 0.67 | 0.0444 |
Attentionで集約したリスク表現: 2.309
結果の読み取り
温度・振動・負荷が同時に高まる工程時点へ大きな重みが付き、担当者が波形を重点確認する候補を絞れます。ただし、この重みは設定したクエリに依存します。実務のAttentionモデルでは、学習データの偏り、マスキング、センサー欠損、品種ごとの工程長、重みと予測寄与の違いを検証し、元波形と併記してレビューします。
対象ノックを通して見える実務上の示唆
- 目的変数で手法を分ける：寸法誤差には最小二乗、不良確率にはLogistic回帰というように、現場の意思決定形式から損失関数を選びます。
- 収束と安定性を品質管理する：勾配、損失、残差、条件数、ヘッセ行列の固有値を記録し、値が出ただけで正常終了としません。
- 相関の強いKPIは整理する：Ridgeは安定化に有効ですが、重複した指標の業務定義を放置する理由にはなりません。
- 非線形化には適用範囲管理が必要である：カーネル法は複雑な品質関係を捉えますが、学習範囲外への外挿と説明可能性に注意します。
- Attentionは確認箇所を絞る補助線である：重みを原因断定に使わず、元のセンサー、保全履歴、現物確認へ接続します。
実務導入する場合に必要なこと
1. 判断目的と損失を定義する
寸法予測、追加検査、条件推奨、設備停止では誤差のコストが異なります。見逃し、過検出、検査工数、停止損失を整理し、RMSEやAUCだけでなく業務KPIで合否を決めます。
2. データの系譜と分割を管理する
設備、品種、工具、材料ロット、作業者、保全履歴、測定器校正を結び付けます。将来情報の混入を避け、時系列・設備・品種をまたぐ検証を行います。
3. 数値計算を監視する
標準化基準、欠損処理、正則化係数、収束条件、ライブラリ版を固定します。条件数、反復回数、残差、予測範囲外フラグをモデルの監視項目に含めます。
4. 現場フローへ接続する
確率閾値ごとの検査対象数、担当部署、確認期限、解除条件を決めます。Attention重みや係数は元データとともに表示し、担当者が根拠を追跡できるようにします。
5. PoC後の責任分界を決める
データ品質、モデル更新、閾値承認、設備停止判断、監査ログの責任者を明確にします。モデル性能が低下したときの縮退運転と旧ルールへの切り戻しも準備します。
まとめ
No.051〜No.060では、勾配降下法と最小二乗法から始め、正規方程式、Ridge回帰、Logistic回帰、Newton法、ヘッセ行列、共役勾配法、カーネル法、Attentionまでを、製造品質の一つの意思決定ストーリーとして確認しました。
行列計算は、予測値を作るためだけの技術ではありません。収束、曲率、相関、正則化、類似度、時点の重みを可視化することで、どの結果をどの条件で信用し、次に何を確認するかを組織で共有するための道具になります。
法人向けのご相談
数理工房では、製造業における品質予測、不良要因分析、加工条件最適化、設備データ活用、モデルの数値安定性評価、PoCから現場運用への移行をご支援しています。データの所在や目的がまだ整理されていない段階でも、業務課題の棚卸しからご相談いただけます。
📩 お問い合わせ: surikobo.co.jp/contact
まずはお気軽にご相談ください。