100本ノック / 行列 / 行列100本ノック
製造業の予知保全をPythonで学ぶ｜Transformer・GNN・行列分解の実務活用
保全履歴と設備ネットワークから次の点検を設計する
Transformer・埋め込み・GNN・大規模行列計算：100本ノック No.061〜No.070
製造現場では、アラームの発生順序、作業記録、設備と部品の接続関係など、表形式だけでは捉えにくいデータが増えています。本記事では、架空の加工ラインを題材に、履歴の文脈を捉え、似た設備・部品を見つけ、ネットワークからリスクを伝播し、計算規模を見積もるまでを一つの意思決定ストーリーとして扱います。
[!NOTE] 本資料は、数理工房 (もしくは代表である和山個人) が過去に企業研修において使用した notebook を企業様の許可を得て再構成・編集のうえ公開しています。
掲載データはすべて架空のものであり、実在する企業・工場・数値とは一切関係ありません。
はじめに：この記事で扱う製造業の実務課題
架空の工場では、設備ごとにアラーム、点検、部品交換の履歴を記録しています。しかし、同じ「振動」アラームでも、その前に過負荷があったのか、給油後なのかで意味が変わります。また、新設設備は故障実績が少なく、単体データだけでは点検優先順位を決めにくい状況です。
そこで、時系列の文脈、設備・部品の類似性、設備ネットワークの接続を行列で表現し、次に確認すべき設備を絞ります。目的はAI用語を並べることではなく、停止損失と点検工数を踏まえた保全判断へ変換することです。
現場でよくある状況
- アラーム単発の閾値監視では、発生順序や同時発生の意味を扱えない
- 作業記録の表記揺れや設備ごとのデータ量の差が大きい
- 類似設備の知見を新設設備へ横展開できていない
- 部品の共用・工程の接続による影響波及を見落とす
- PoCでは動いても、全工場へ広げるとメモリや計算時間が急増する
なぜこの問題は判断が難しいのか
履歴は「行の集合」ではなく順序を持ち、設備ネットワークは互いに独立ではありません。さらに、埋め込みや深層学習のスコアは便利でも、学習データの偏り、未知設備、計算資源、説明責任を無視すると運用できません。
本記事では小さな行列で原理を可視化します。TransformerやGNNの本番実装を再現するのではなく、各手法の中心となる行列演算と、保全担当者が確認すべき出力・限界を切り分けます。
今回扱うノックの全体像
| No. | テーマ | 製造業での判断 |
|---|---|---|
| 061 | Transformer | 履歴のどの事象を参照して現在リスクを判断するか |
| 062 | 埋め込み | 設備・部品・アラームの類似性を数値化する |
| 063 | 行列分解 | 履歴行列を少数の潜在要因へ要約する |
| 064 | Matrix Factorization | 未観測の設備―部品リスクを補完する |
| 065 | LightGCN | 接続関係から類似設備へ知見を伝播する |
| 066 | Graph Neural Network | 接続と設備特徴を同時に集約する |
| 067 | 深層学習 | 非線形な故障リスクを表現する |
| 068 | 自動微分 | 損失から各パラメータの更新方向を求める |
| 069 | GPU計算 | 行列計算を高速化すべき規模か見積もる |
| 070 | 大規模行列計算 | 疎行列で全工場規模へ拡張する |
Python 環境の準備
NumPyで密行列計算、pandasで表、SciPyで疎行列、Matplotlibで可視化を行います。外部データやGPU専用ライブラリには依存しません。乱数生成器は np.random.default_rng(61) で固定します。
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
import platform
import sys
import time
import japanize_matplotlib
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scipy
from IPython.display import display
from scipy import sparse
rng = np.random.default_rng(61)
pd.set_option("display.precision", 4)
print("Python:", sys.version.split()[0])
print("OS:", platform.platform())
print("NumPy:", np.__version__)
print("pandas:", pd.__version__)
print("SciPy:", scipy.__version__)
print("Matplotlib:", matplotlib.__version__)Python: 3.13.1
OS: macOS-26.3-arm64-arm-64bit-Mach-O
NumPy: 2.5.1
pandas: 3.0.3
SciPy: 1.18.0
Matplotlib: 3.11.0架空データの作成
6台の設備、8種類の保全イベント、7種類の部品について架空データを作ります。設備には「回転系」「油圧系」「搬送系」という潜在的な特性があり、それに応じてイベントと部品の利用傾向が変わる設定です。
実務ではイベント辞書、設備階層、時刻粒度、未記録と「発生なし」の区別が分析品質を左右します。ここでは原理を見通せる規模に限定します。
equipment = ["旋盤A", "旋盤B", "研削盤A", "研削盤B", "搬送機A", "搬送機B"]
events = ["正常", "過負荷", "温度上昇", "振動", "圧力低下", "給油", "点検", "部品交換"]
parts = ["主軸軸受", "砥石", "油圧ポンプ", "シール", "搬送ベルト", "モーター", "潤滑油"]
# 行：設備、列：潜在特性（回転・油圧・搬送）
equipment_type = np.array([
[1.0, 0.4, 0.1], [0.9, 0.5, 0.1], [1.0, 0.2, 0.1],
[0.9, 0.3, 0.1], [0.2, 0.1, 1.0], [0.2, 0.2, 0.9],
])
event_affinity = np.array([
[0.2, 0.2, 0.1], [0.8, 0.3, 0.4], [0.7, 0.4, 0.2], [1.0, 0.1, 0.2],
[0.1, 1.0, 0.1], [0.4, 0.8, 0.3], [0.4, 0.4, 0.4], [0.8, 0.5, 0.7],
])
counts = rng.poisson(2 + 7 * (equipment_type @ event_affinity.T))
history_df = pd.DataFrame(counts, index=equipment, columns=events)
sequences = [
["正常", "過負荷", "温度上昇", "振動", "点検"],
["正常", "給油", "正常", "過負荷", "温度上昇"],
["正常", "振動", "振動", "部品交換", "正常"],
["正常", "温度上昇", "給油", "正常", "点検"],
["正常", "圧力低下", "過負荷", "点検", "正常"],
["正常", "過負荷", "圧力低下", "温度上昇", "点検"],
]
display(history_df)
display(pd.DataFrame(sequences, index=equipment, columns=[f"時点{i}" for i in range(1, 6)]))
fig, ax = plt.subplots(figsize=(9, 4.2))
im = ax.imshow(history_df, cmap="YlOrRd", aspect="auto")
ax.set_title("設備別の保全イベント回数（架空データ）")
ax.set_xlabel("イベント")
ax.set_ylabel("設備")
ax.set_xticks(range(len(events)), events, rotation=30, ha="right")
ax.set_yticks(range(len(equipment)), equipment)
ax.grid(True, color="white", linewidth=0.4, alpha=0.6)
fig.colorbar(im, ax=ax, label="記録回数")
plt.tight_layout()
plt.show()| 正常 | 過負荷 | 温度上昇 | 振動 | 圧力低下 | 給油 | 点検 | 部品交換 | |
|---|---|---|---|---|---|---|---|---|
| 旋盤A | 8 | 9 | 9 | 8 | 11 | 7 | 10 | 6 |
| 旋盤B | 4 | 7 | 4 | 9 | 11 | 6 | 6 | 9 |
| 研削盤A | 6 | 3 | 5 | 10 | 2 | 4 | 2 | 7 |
| 研削盤B | 3 | 11 | 6 | 10 | 5 | 4 | 4 | 8 |
| 搬送機A | 5 | 7 | 6 | 9 | 2 | 8 | 9 | 10 |
| 搬送機B | 3 | 4 | 2 | 7 | 7 | 4 | 8 | 7 |
| 時点1 | 時点2 | 時点3 | 時点4 | 時点5 | |
|---|---|---|---|---|---|
| 旋盤A | 正常 | 過負荷 | 温度上昇 | 振動 | 点検 |
| 旋盤B | 正常 | 給油 | 正常 | 過負荷 | 温度上昇 |
| 研削盤A | 正常 | 振動 | 振動 | 部品交換 | 正常 |
| 研削盤B | 正常 | 温度上昇 | 給油 | 正常 | 点検 |
| 搬送機A | 正常 | 圧力低下 | 過負荷 | 点検 | 正常 |
| 搬送機B | 正常 | 過負荷 | 圧力低下 | 温度上昇 | 点検 |
No.061：Transformer
実務での意味
Transformerは、現在の状態を判断するときに過去のどの記録をどれだけ参照するかをAttentionで表します。アラームの順序や遠く離れた保全記録を扱えるため、故障予兆、作業記録分類、保全記録の要約に使われます。
分析・モデル化の考え方
各時点の表現からQuery、Key、Valueを作り、Scaled Dot-Product Attentionを
で計算します。未来を見ない予測ではマスク を使います。ここでは直近時点が過去のどこを参照するかを可視化し、学習済みモデルの性能を主張するのではなく仕組みを確認します。
Pythonで確認する
event_to_id = {name: i for i, name in enumerate(events)}
d_model = 6
token_embedding = rng.normal(0, 0.7, size=(len(events), d_model))
position = np.arange(5)[:, None]
position_encoding = np.column_stack([
np.sin(position[:, 0] / (10000 ** (2 * k / d_model))) if k % 2 == 0
else np.cos(position[:, 0] / (10000 ** (2 * (k - 1) / d_model)))
for k in range(d_model)
])
seq = sequences[0]
X_seq = token_embedding[[event_to_id[x] for x in seq]] + position_encoding
Q = K = V = X_seq
scores = Q @ K.T / np.sqrt(d_model)
causal_mask = np.triu(np.full_like(scores, -np.inf), k=1)
masked_scores = scores + causal_mask
weights = np.exp(masked_scores - np.max(masked_scores, axis=1, keepdims=True))
weights = weights / weights.sum(axis=1, keepdims=True)
context = weights @ V
display(pd.DataFrame(weights, index=[f"判断:{x}" for x in seq], columns=[f"参照:{x}" for x in seq]).style.format("{:.2f}"))
fig, ax = plt.subplots(figsize=(7.4, 4.4))
im = ax.imshow(weights, cmap="Blues", vmin=0, vmax=1)
ax.set_title("因果マスク付きAttentionの参照重み")
ax.set_xlabel("参照する過去イベント")
ax.set_ylabel("判断対象イベント")
ax.set_xticks(range(5), seq, rotation=25, ha="right")
ax.set_yticks(range(5), seq)
ax.grid(True, color="white", linewidth=0.4)
fig.colorbar(im, ax=ax, label="Attention重み")
plt.tight_layout()
plt.show()| 参照:正常 | 参照:過負荷 | 参照:温度上昇 | 参照:振動 | 参照:点検 | |
|---|---|---|---|---|---|
| 判断:正常 | 1.00 | 0.00 | 0.00 | 0.00 | 0.00 |
| 判断:過負荷 | 0.33 | 0.67 | 0.00 | 0.00 | 0.00 |
| 判断:温度上昇 | 0.21 | 0.39 | 0.40 | 0.00 | 0.00 |
| 判断:振動 | 0.14 | 0.26 | 0.28 | 0.31 | 0.00 |
| 判断:点検 | 0.09 | 0.07 | 0.06 | 0.10 | 0.67 |
結果の読み取り
上三角が0なので、各時点は未来の記録を参照していません。直近の「点検」の文脈表現には過去の過負荷・温度上昇・振動が重み付きで統合されます。実務では重みを因果根拠と断定せず、故障前リードタイム、見逃し率、設備別性能、データ漏洩を検証します。
No.062：埋め込み
実務での意味
埋め込みは、設備やイベントを少数次元のベクトルに変換し、「似た履歴を持つ対象」を距離や内積で探せるようにします。表記が異なる記録の統合、新設設備の類似設備検索、保全部品の候補提示に利用できます。
分析・モデル化の考え方
設備―イベント行列を標準化し、SVDから設備埋め込み を作ります。コサイン類似度
で方向の近さを測ります。類似は因果や代替可能性を意味しないため、設備仕様との照合が必要です。
Pythonで確認する
H = np.log1p(history_df.to_numpy())
H_centered = H - H.mean(axis=0, keepdims=True)
U_h, s_h, Vt_h = np.linalg.svd(H_centered, full_matrices=False)
equipment_embedding = U_h[:, :2] * s_h[:2]
norms = np.linalg.norm(equipment_embedding, axis=1, keepdims=True)
cosine = (equipment_embedding @ equipment_embedding.T) / np.maximum(norms @ norms.T, 1e-12)
embedding_df = pd.DataFrame(equipment_embedding, index=equipment, columns=["埋め込み1", "埋め込み2"])
display(embedding_df.style.format("{:+.3f}"))
display(pd.DataFrame(cosine, index=equipment, columns=equipment).style.format("{:.2f}"))
fig, ax = plt.subplots(figsize=(7.2, 4.8))
ax.scatter(equipment_embedding[:, 0], equipment_embedding[:, 1], s=80, color="teal")
for i, name in enumerate(equipment):
ax.annotate(name, equipment_embedding[i], xytext=(5, 5), textcoords="offset points")
ax.set_title("イベント履歴から得た設備埋め込み")
ax.set_xlabel("埋め込み次元1")
ax.set_ylabel("埋め込み次元2")
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()| 埋め込み1 | 埋め込み2 | |
|---|---|---|
| 旋盤A | +0.965 | +0.479 |
| 旋盤B | +0.512 | -0.349 |
| 研削盤A | -1.232 | +0.039 |
| 研削盤B | -0.065 | +0.078 |
| 搬送機A | -0.257 | +0.635 |
| 搬送機B | +0.076 | -0.882 |
| 旋盤A | 旋盤B | 研削盤A | 研削盤B | 搬送機A | 搬送機B | |
|---|---|---|---|---|---|---|
| 旋盤A | 1.00 | 0.49 | -0.88 | -0.23 | 0.08 | -0.37 |
| 旋盤B | 0.49 | 1.00 | -0.84 | -0.96 | -0.83 | 0.63 |
| 研削盤A | -0.88 | -0.84 | 1.00 | 0.66 | 0.40 | -0.12 |
| 研削盤B | -0.23 | -0.96 | 0.66 | 1.00 | 0.95 | -0.82 |
| 搬送機A | 0.08 | -0.83 | 0.40 | 0.95 | 1.00 | -0.96 |
| 搬送機B | -0.37 | 0.63 | -0.12 | -0.82 | -0.96 | 1.00 |
結果の読み取り
同系統の設備が近い方向に配置され、履歴パターンの類似候補を作れます。距離が近い設備の保全手順をそのまま転用せず、型式、負荷、稼働時間、故障モードを確認します。埋め込みは候補探索の入口であり、設備台帳と現場知識による絞り込みが必要です。
No.063：行列分解
実務での意味
行列分解は、多数の設備×イベントの関係を「回転系」「油圧系」など少数の潜在要因へ要約します。監視KPIの整理、データ圧縮、ノイズ除去、設備群の特徴把握に有効です。
分析・モデル化の考え方
SVDにより と分解し、上位 成分だけの
で近似します。Frobeniusノルムによる相対再構成誤差と累積寄与率を見て、簡素化と情報損失を比較します。
Pythonで確認する
rank_rows = []
for rank in range(1, min(H.shape) + 1):
H_rank = (U_h[:, :rank] * s_h[:rank]) @ Vt_h[:rank]
rank_rows.append({
"ランク": rank,
"累積寄与率": np.sum(s_h[:rank] ** 2) / np.sum(s_h ** 2),
"相対再構成誤差": np.linalg.norm(H_centered - H_rank) / np.linalg.norm(H_centered),
})
rank_df = pd.DataFrame(rank_rows)
display(rank_df.style.format({"累積寄与率": "{:.1%}", "相対再構成誤差": "{:.1%}"}))
fig, ax = plt.subplots(figsize=(7.2, 4.0))
ax.plot(rank_df["ランク"], rank_df["累積寄与率"] * 100, marker="o", label="累積寄与率")
ax.plot(rank_df["ランク"], rank_df["相対再構成誤差"] * 100, marker="s", label="相対再構成誤差")
ax.set_title("採用ランクと情報保持・再構成誤差")
ax.set_xlabel("採用ランク")
ax.set_ylabel("割合（%）")
ax.set_ylim(0, 105)
ax.grid(True, alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()| ランク | 累積寄与率 | 相対再構成誤差 | |
|---|---|---|---|
| 0 | 1 | 46.5% | 73.1% |
| 1 | 2 | 72.2% | 52.7% |
| 2 | 3 | 86.7% | 36.5% |
| 3 | 4 | 97.7% | 15.1% |
| 4 | 5 | 100.0% | 0.0% |
| 5 | 6 | 100.0% | 0.0% |
結果の読み取り
ランクを増やすほど情報保持率は上がり、再構成誤差は下がります。少数成分で概要を捉えられても、希少故障は小さい成分に現れる可能性があります。圧縮率だけでランクを決めず、故障モード別の再現率と残差監視を併用します。
No.064：Matrix Factorization
実務での意味
設備と部品の組合せには、交換実績がないだけでリスクが低いとは限らない未観測箇所があります。Matrix Factorizationは既知の関係から潜在因子を学び、未観測の組合せを優先調査候補として補完します。
分析・モデル化の考え方
観測集合 について設備因子 と部品因子 を求め、
を確率的勾配降下法で最小化します。未観測を0と扱わず、評価からも分離することが重要です。
Pythonで確認する
part_factor = np.array([
[1.0, 0.1, 0.1], [0.9, 0.1, 0.1], [0.2, 1.0, 0.1], [0.1, 0.9, 0.2],
[0.1, 0.1, 1.0], [0.7, 0.2, 0.7], [0.4, 0.8, 0.3],
])
true_risk = 1 + 4 * (equipment_type @ part_factor.T) / 1.5
true_risk = np.clip(true_risk + rng.normal(0, 0.15, true_risk.shape), 1, 5)
observed = rng.random(true_risk.shape) < 0.68
R = np.where(observed, true_risk, np.nan)
k = 3
P = rng.normal(0, 0.3, size=(len(equipment), k))
Q = rng.normal(0, 0.3, size=(len(parts), k))
lr, reg = 0.025, 0.02
losses = []
obs_pairs = np.argwhere(observed)
for epoch in range(600):
for i, j in obs_pairs[rng.permutation(len(obs_pairs))]:
err = R[i, j] - P[i] @ Q[j]
p_old = P[i].copy()
P[i] += lr * (err * Q[j] - reg * P[i])
Q[j] += lr * (err * p_old - reg * Q[j])
pred = P @ Q.T
losses.append(np.mean((R[observed] - pred[observed]) ** 2))
prediction = np.clip(P @ Q.T, 1, 5)
candidates = pd.DataFrame([
{"設備": equipment[i], "部品": parts[j], "予測リスク": prediction[i, j]}
for i, j in np.argwhere(~observed)
]).sort_values("予測リスク", ascending=False)
display(candidates.head(8).style.format({"予測リスク": "{:.2f}"}))
fig, ax = plt.subplots(figsize=(7.2, 4.0))
ax.plot(losses, color="tab:purple")
ax.set_title("Matrix Factorizationの学習損失")
ax.set_xlabel("エポック")
ax.set_ylabel("観測箇所の平均二乗誤差")
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()| 設備 | 部品 | 予測リスク | |
|---|---|---|---|
| 0 | 旋盤A | 主軸軸受 | 3.66 |
| 2 | 旋盤B | 砥石 | 3.04 |
| 1 | 旋盤A | 潤滑油 | 2.89 |
| 3 | 旋盤B | 油圧ポンプ | 2.43 |
| 4 | 旋盤B | シール | 2.31 |
| 6 | 搬送機A | 潤滑油 | 2.11 |
| 7 | 搬送機B | 主軸軸受 | 1.71 |
| 5 | 研削盤B | 搬送ベルト | 1.33 |
結果の読み取り
未観測の設備―部品組合せに予測値が入り、追加点検や台帳確認の候補を順位付けできます。ただし、観測済み箇所だけの学習誤差は未知箇所の精度を保証しません。時系列で分けた検証、交換されやすさの偏り、必須部品・非搭載部品のマスクを設計します。
No.065：LightGCN
実務での意味
LightGCNは、設備―部品の二部グラフ上で近隣の埋め込みを伝播します。同じ部品を使う設備や、似た設備が使う部品から情報を借り、履歴の少ない対象にも候補を提示できます。
分析・モデル化の考え方
自己ループを加えず、正規化隣接行列
で埋め込みを伝播し、各層の平均を最終表現とします。層を重ねすぎると表現が均される過平滑化に注意します。
Pythonで確認する
interaction = (true_risk >= 3.2).astype(float)
n_e, n_p = len(equipment), len(parts)
A = np.block([
[np.zeros((n_e, n_e)), interaction],
[interaction.T, np.zeros((n_p, n_p))],
])
degree = A.sum(axis=1)
D_inv_sqrt = np.diag(1 / np.sqrt(np.maximum(degree, 1)))
A_norm = D_inv_sqrt @ A @ D_inv_sqrt
E0 = rng.normal(0, 0.3, size=(n_e + n_p, 4))
layers = [E0]
for _ in range(3):
layers.append(A_norm @ layers[-1])
E_light = np.mean(layers, axis=0)
light_scores = E_light[:n_e] @ E_light[n_e:].T
layer_change = [np.mean(np.linalg.norm(layers[i + 1] - layers[i], axis=1)) for i in range(3)]
display(pd.DataFrame(light_scores, index=equipment, columns=parts).style.format("{:+.2f}"))
fig, ax = plt.subplots(figsize=(7.2, 4.0))
ax.bar(["0→1層", "1→2層", "2→3層"], layer_change, color="steelblue")
ax.set_title("LightGCNの伝播層ごとの埋め込み変化")
ax.set_xlabel("伝播ステップ")
ax.set_ylabel("ノード当たり平均変化量")
ax.grid(True, axis="y", alpha=0.3)
plt.tight_layout()
plt.show()| 主軸軸受 | 砥石 | 油圧ポンプ | シール | 搬送ベルト | モーター | 潤滑油 | |
|---|---|---|---|---|---|---|---|
| 旋盤A | +0.11 | +0.09 | -0.03 | -0.01 | +0.03 | +0.06 | +0.01 |
| 旋盤B | +0.05 | +0.04 | -0.02 | -0.01 | +0.00 | +0.07 | -0.01 |
| 研削盤A | +0.11 | +0.11 | -0.03 | +0.00 | +0.03 | -0.02 | +0.04 |
| 研削盤B | +0.09 | +0.09 | -0.03 | +0.00 | +0.06 | +0.01 | +0.02 |
| 搬送機A | +0.03 | +0.05 | -0.02 | -0.02 | +0.03 | +0.07 | +0.02 |
| 搬送機B | +0.01 | +0.03 | -0.01 | -0.01 | +0.07 | +0.05 | +0.00 |
結果の読み取り
部品との直接関係だけでなく、同じ部品へ接続する設備を介した情報が表現へ入ります。層が進むほど変化量が小さくなり、過度な伝播では設備ごとの差が失われます。実務では層数、負例の定義、新設ノード、グラフ更新頻度を検証します。
No.066：Graph Neural Network
実務での意味
GNNは接続関係だけでなく、稼働率、経過年数、異常回数など各設備の特徴も近隣から集約します。工程の前後関係や部品共用ネットワークを考慮したリスク評価に使えます。
分析・モデル化の考え方
自己ループを加えた正規化隣接行列 を使い、1層のGraph Convolutionを
とします。 は設備特徴、 は学習パラメータです。ここでは固定した重みで近隣集約の効果だけを比較します。
Pythonで確認する
# 工程順と共用部品を簡略化した設備グラフ
equipment_adj = np.array([
[0, 1, 1, 0, 0, 0], [1, 0, 0, 1, 0, 0], [1, 0, 0, 1, 1, 0],
[0, 1, 1, 0, 0, 1], [0, 0, 1, 0, 0, 1], [0, 0, 0, 1, 1, 0],
], dtype=float)
X_equipment = np.column_stack([
np.array([0.82, 0.77, 0.91, 0.88, 0.70, 0.74]), # 稼働率
np.array([8, 4, 11, 6, 5, 9]) / 12, # 経過年数（縮尺済み）
history_df[["振動", "温度上昇", "圧力低下"]].sum(axis=1).to_numpy() / 40,
])
A_self = equipment_adj + np.eye(n_e)
D_self = np.diag(1 / np.sqrt(A_self.sum(axis=1)))
A_hat = D_self @ A_self @ D_self
W_gnn = np.array([[0.7, -0.3], [0.4, 0.8], [1.1, 0.5]])
H_gnn = np.maximum(0, A_hat @ X_equipment @ W_gnn)
local_only = np.maximum(0, X_equipment @ W_gnn)
gnn_df = pd.DataFrame({
"局所表現1": local_only[:, 0], "近隣集約後1": H_gnn[:, 0],
"局所表現2": local_only[:, 1], "近隣集約後2": H_gnn[:, 1],
}, index=equipment)
display(gnn_df.style.format("{:.3f}"))
fig, ax = plt.subplots(figsize=(8.0, 4.2))
x = np.arange(n_e)
ax.bar(x - 0.18, local_only[:, 0], width=0.36, label="自設備のみ")
ax.bar(x + 0.18, H_gnn[:, 0], width=0.36, label="近隣集約後")
ax.set_title("GNNによる設備特徴の近隣集約")
ax.set_xlabel("設備")
ax.set_ylabel("表現値（第1次元）")
ax.set_xticks(x, equipment)
ax.grid(True, axis="y", alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()| 局所表現1 | 近隣集約後1 | 局所表現2 | 近隣集約後2 | |
|---|---|---|---|---|
| 旋盤A | 1.611 | 1.406 | 0.637 | 0.519 |
| 旋盤B | 1.332 | 1.383 | 0.336 | 0.439 |
| 研削盤A | 1.471 | 1.506 | 0.673 | 0.549 |
| 研削盤B | 1.394 | 1.464 | 0.399 | 0.532 |
| 搬送機A | 1.124 | 1.219 | 0.336 | 0.499 |
| 搬送機B | 1.258 | 1.196 | 0.578 | 0.420 |
結果の読み取り
自設備だけの表現と、接続設備を含めた表現に差が生じます。上流設備の異常が下流へ影響する現場では有用ですが、接続を無向とするか、時間遅れや流量をどう重みにするかで意味が変わります。設備エンジニアとグラフ定義を合意し、誤った影響伝播を防ぎます。
No.067：深層学習
実務での意味
深層学習は、温度と負荷の組合せなど、単純な加算では表しにくい非線形な故障リスクを学習できます。画像、波形、履歴を統合する基盤にもなります。
分析・モデル化の考え方
2層ニューラルネットワークを
とします。ここでは架空の設備状態から二値ラベルを作り、勾配降下法で学習します。訓練精度だけでなく、設備・期間を分けた外部検証が必要です。
Pythonで確認する
n_samples = 500
X_dl = rng.normal(size=(n_samples, 3))
nonlinear_risk = 1.3 * X_dl[:, 0] * X_dl[:, 1] + 0.9 * X_dl[:, 2] ** 2 - 0.5
y_dl = (nonlinear_risk + rng.normal(0, 0.35, n_samples) > 0).astype(float)[:, None]
split = 400
mean_dl, std_dl = X_dl[:split].mean(axis=0), X_dl[:split].std(axis=0)
X_dl = (X_dl - mean_dl) / std_dl
W1 = rng.normal(0, 0.3, size=(3, 8)); b1 = np.zeros((1, 8))
W2 = rng.normal(0, 0.3, size=(8, 1)); b2 = np.zeros((1, 1))
lr_dl = 0.08
dl_losses = []
for epoch in range(800):
Xb, yb = X_dl[:split], y_dl[:split]
z1 = Xb @ W1 + b1
h = np.maximum(0, z1)
logits = h @ W2 + b2
prob = 1 / (1 + np.exp(-np.clip(logits, -30, 30)))
loss = -np.mean(yb * np.log(prob + 1e-9) + (1 - yb) * np.log(1 - prob + 1e-9))
dl_losses.append(loss)
dlogits = (prob - yb) / split
dW2, db2 = h.T @ dlogits, dlogits.sum(axis=0, keepdims=True)
dh = dlogits @ W2.T
dz1 = dh * (z1 > 0)
dW1, db1 = Xb.T @ dz1, dz1.sum(axis=0, keepdims=True)
W1 -= lr_dl * dW1; b1 -= lr_dl * db1
W2 -= lr_dl * dW2; b2 -= lr_dl * db2
def mlp_predict(X):
hidden = np.maximum(0, X @ W1 + b1)
return 1 / (1 + np.exp(-np.clip(hidden @ W2 + b2, -30, 30)))
train_acc = np.mean((mlp_predict(X_dl[:split]) >= 0.5) == y_dl[:split])
test_acc = np.mean((mlp_predict(X_dl[split:]) >= 0.5) == y_dl[split:])
display(pd.DataFrame({"区分": ["学習", "評価"], "正解率": [train_acc, test_acc]}).style.format({"正解率": "{:.1%}"}))
fig, ax = plt.subplots(figsize=(7.2, 4.0))
ax.plot(dl_losses, color="tab:red")
ax.set_title("2層ニューラルネットワークの学習曲線")
ax.set_xlabel("エポック")
ax.set_ylabel("交差エントロピー損失")
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()| 区分 | 正解率 | |
|---|---|---|
| 0 | 学習 | 90.0% |
| 1 | 評価 | 89.0% |
結果の読み取り
損失が低下し、非線形な関係を持つ架空データで一定の分類性能を得ています。ただし正解率だけでは、まれな故障の見逃しや誤停止コストを評価できません。適合率、再現率、校正、リードタイムを確認し、単純モデルに対する改善幅が運用コストに見合う場合に採用します。
No.068：自動微分
実務での意味
深層学習は多数のパラメータを持つため、損失が各パラメータへどう影響するかを手計算できません。自動微分は計算グラフへ連鎖律を適用し、勾配を正確かつ効率的に求めます。
分析・モデル化の考え方
単純な予測 、損失 でも、連鎖律により
です。小さな逆向き自動微分クラスを実装し、有限差分と照合します。本番では十分に検証されたフレームワークを利用します。
Pythonで確認する
class Value:
def __init__(self, data, parents=(), grads=()):
self.data = float(data)
self.grad = 0.0
self.parents = parents
self.grads = grads
def __add__(self, other):
other = other if isinstance(other, Value) else Value(other)
return Value(self.data + other.data, (self, other), (1.0, 1.0))
__radd__ = __add__
def __mul__(self, other):
other = other if isinstance(other, Value) else Value(other)
return Value(self.data * other.data, (self, other), (other.data, self.data))
__rmul__ = __mul__
def __neg__(self):
return self * -1
def __sub__(self, other):
return self + (-other)
def __pow__(self, power):
return Value(self.data ** power, (self,), (power * self.data ** (power - 1),))
def backward(self):
topo, seen = [], set()
def build(v):
if id(v) not in seen:
seen.add(id(v))
for p in v.parents: build(p)
topo.append(v)
build(self)
self.grad = 1.0
for v in reversed(topo):
for p, local_grad in zip(v.parents, v.grads):
p.grad += v.grad * local_grad
x_ad, y_ad = 1.7, 4.2
w_ad, b_ad = Value(1.1), Value(0.3)
loss_ad = (w_ad * x_ad + b_ad - y_ad) ** 2
loss_ad.backward()
def scalar_loss(w, b):
return (w * x_ad + b - y_ad) ** 2
eps = 1e-6
fd_w = (scalar_loss(1.1 + eps, 0.3) - scalar_loss(1.1 - eps, 0.3)) / (2 * eps)
fd_b = (scalar_loss(1.1, 0.3 + eps) - scalar_loss(1.1, 0.3 - eps)) / (2 * eps)
grad_df = pd.DataFrame({
"パラメータ": ["w", "b"], "自動微分": [w_ad.grad, b_ad.grad],
"有限差分": [fd_w, fd_b], "絶対差": [abs(w_ad.grad - fd_w), abs(b_ad.grad - fd_b)],
})
display(grad_df.style.format({"自動微分": "{:.6f}", "有限差分": "{:.6f}", "絶対差": "{:.2e}"}))
fig, ax = plt.subplots(figsize=(7.2, 4.0))
ax.bar(np.arange(2) - 0.17, grad_df["自動微分"], width=0.34, label="自動微分")
ax.bar(np.arange(2) + 0.17, grad_df["有限差分"], width=0.34, label="有限差分")
ax.set_title("自動微分と有限差分による勾配の照合")
ax.set_xlabel("パラメータ")
ax.set_ylabel("損失の勾配")
ax.set_xticks(range(2), grad_df["パラメータ"])
ax.grid(True, axis="y", alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()| パラメータ | 自動微分 | 有限差分 | 絶対差 | |
|---|---|---|---|---|
| 0 | w | -6.902000 | -6.902000 | 4.69e-10 |
| 1 | b | -4.060000 | -4.060000 | 2.50e-10 |
結果の読み取り
自動微分の勾配は有限差分と数値誤差の範囲で一致します。実務では代表的な小規模入力で勾配チェックを行い、NaN、勾配消失・爆発、損失の実装ミスを監視します。勾配が計算できることと、目的関数が現場の損失を正しく表すことは別です。
No.069：GPU計算
実務での意味
GPUは多数の積和演算を並列実行でき、深層学習や大規模な埋め込み計算を高速化します。一方、小さな表データでは転送・起動のオーバーヘッドが支配し、CPUで十分なことがあります。
分析・モデル化の考え方
と の行列積には概ね FLOPが必要です。必要メモリは入力・出力・中間値・勾配を含めて見積もります。ここではGPUを要求せず、行列サイズごとの演算量と最低限の配列メモリを算出し、ベクトル化したCPU計算の時間も確認します。
Pythonで確認する
sizes = [128, 512, 1024, 2048, 4096]
gpu_plan = pd.DataFrame({
"正方行列サイズ": sizes,
"行列積の演算量_GFLOP": [2 * n**3 / 1e9 for n in sizes],
"3配列の最低メモリ_GB_float32": [3 * n**2 * 4 / 1e9 for n in sizes],
"3配列の最低メモリ_GB_float64": [3 * n**2 * 8 / 1e9 for n in sizes],
})
display(gpu_plan.style.format({
"行列積の演算量_GFLOP": "{:.2f}",
"3配列の最低メモリ_GB_float32": "{:.3f}",
"3配列の最低メモリ_GB_float64": "{:.3f}",
}))
bench_sizes = [128, 256, 512]
bench_times = []
for n in bench_sizes:
a = rng.normal(size=(n, n)).astype(np.float32)
b = rng.normal(size=(n, n)).astype(np.float32)
start = time.perf_counter()
_ = a @ b
bench_times.append(time.perf_counter() - start)
fig, ax = plt.subplots(figsize=(7.2, 4.0))
ax.plot(bench_sizes, np.array(bench_times) * 1000, marker="o", color="tab:orange")
ax.set_title("CPUでのベクトル化行列積の実測時間")
ax.set_xlabel("正方行列サイズ")
ax.set_ylabel("実行時間（ms、この環境での参考値）")
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()| 正方行列サイズ | 行列積の演算量_GFLOP | 3配列の最低メモリ_GB_float32 | 3配列の最低メモリ_GB_float64 | |
|---|---|---|---|---|
| 0 | 128 | 0.00 | 0.000 | 0.000 |
| 1 | 512 | 0.27 | 0.003 | 0.006 |
| 2 | 1024 | 2.15 | 0.013 | 0.025 |
| 3 | 2048 | 17.18 | 0.050 | 0.101 |
| 4 | 4096 | 137.44 | 0.201 | 0.403 |
結果の読み取り
行列サイズが2倍になると演算量は約8倍となり、大規模化でGPU候補になります。ただし表のメモリは最低限で、学習時は勾配・最適化状態・一時領域が追加されます。実機データで前処理、転送、学習、推論を含むエンドツーエンド時間を測り、納期・費用・消費電力を含めて導入を判断します。
No.070：大規模行列計算
実務での意味
全工場の設備―部品、設備―イベント行列は大部分が0です。密行列のまま保持すると、データがない箇所にもメモリを使い、計算規模を無駄に増やします。疎行列は実際に存在する関係だけを保存します。
分析・モデル化の考え方
行列の非ゼロ数を とすると、密行列の保存量は 、CSR形式は概ね です。行列ベクトル積も非ゼロ要素中心に計算できます。疎であること、アクセスパターン、更新頻度に応じて形式を選びます。
Pythonで確認する
n_rows, n_cols, density = 20_000, 8_000, 0.0005
large_sparse = sparse.random(
n_rows, n_cols, density=density, format="csr", dtype=np.float32, random_state=61,
)
vector = rng.normal(size=n_cols).astype(np.float32)
start = time.perf_counter()
result_sparse = large_sparse @ vector
sparse_time = time.perf_counter() - start
dense_bytes = n_rows * n_cols * np.dtype(np.float32).itemsize
csr_bytes = large_sparse.data.nbytes + large_sparse.indices.nbytes + large_sparse.indptr.nbytes
storage_df = pd.DataFrame({
"表現": ["密行列（推定）", "CSR疎行列（実測）"],
"メモリ_MB": [dense_bytes / 1e6, csr_bytes / 1e6],
"非ゼロ要素数": [n_rows * n_cols, large_sparse.nnz],
})
display(storage_df.style.format({"メモリ_MB": "{:,.2f}", "非ゼロ要素数": "{:,.0f}"}))
print("疎行列ベクトル積:", f"{sparse_time * 1000:.2f} ms", "/ 出力形状:", result_sparse.shape)
fig, ax = plt.subplots(figsize=(7.2, 4.0))
ax.bar(storage_df["表現"], storage_df["メモリ_MB"], color=["tab:red", "steelblue"])
ax.set_yscale("log")
ax.set_title("密行列とCSR疎行列の保存量比較")
ax.set_xlabel("行列表現")
ax.set_ylabel("メモリ（MB、対数目盛）")
ax.grid(True, axis="y", which="both", alpha=0.3)
plt.tight_layout()
plt.show()| 表現 | メモリ_MB | 非ゼロ要素数 | |
|---|---|---|---|
| 0 | 密行列（推定） | 640.00 | 160,000,000 |
| 1 | CSR疎行列（実測） | 0.72 | 80,000 |
疎行列ベクトル積: 0.33 ms / 出力形状: (20000,)
結果の読み取り
0.05%しか値を持たない行列では、CSR形式が密行列より大幅に少ないメモリで保持できます。全工場展開では、疎行列化に加えて、分割処理、ミニバッチ、近似探索、増分更新、監視対象期間の制限を組み合わせます。実行時間だけでなく、再計算時間、障害復旧、データ更新の整合性も設計対象です。
対象ノックを通して見える実務上の示唆
- 履歴は順序と文脈を含む：TransformerのAttentionにより、単発アラームでは見えない事象の組合せを扱えます。
- 埋め込みは候補探索を助ける：設備・部品の類似性を数値化できますが、仕様確認と因果検証は別途必要です。
- 未観測はゼロではない：Matrix Factorizationでは、記録なしとリスクなしを区別して評価します。
- ネットワーク定義がモデルの意味を決める：LightGCNやGNNでは、接続方向・重み・時間遅れを現場と合意します。
- 高度なモデルほど基準モデルとの比較が重要である：精度改善が説明・運用・計算コストに見合うか判断します。
- 計算資源はモデル設計の一部である：GPU、疎行列、分割方法をデータ規模と更新期限から選びます。
実務導入する場合に必要なこと
1. 判断単位と正解を定義する
「24時間以内の停止」「7日以内の部品交換」「次回定修での点検」など、予測対象と時間窓を決めます。アラーム、作業指示、交換、原因確定を区別し、ラベルが現場判断の結果に偏ることも記録します。
2. データとグラフを管理する
設備ID、部品ID、イベント辞書、時刻、設備構成の版を管理します。設備増設や工程変更をグラフへ反映し、未来情報の混入、記録漏れ、保全方針変更を監視します。
3. 段階的に検証する
単純なルール・回帰を基準に、設備・期間を分けた評価、故障モード別の再現率、誤警報、先行時間を比較します。Attentionや埋め込みの可視化は説明補助として使い、因果説明とは区別します。
4. 運用・計算・責任分界を設計する
誰が候補を確認し、何を追加測定し、停止・継続・経過観察を決めるかを手順化します。再学習周期、処理期限、GPU費用、障害時の代替ルール、モデルとデータの版管理を含めます。
まとめ
No.061〜No.070では、架空の設備履歴を題材に、TransformerのAttention、埋め込み、行列分解、Matrix Factorization、LightGCN、GNN、深層学習、自動微分、GPU計算の規模見積もり、疎行列による大規模計算を確認しました。
製造業でこれらを活かす鍵は、モデルを複雑にすること自体ではありません。履歴・設備構成・現場損失を正しく表現し、候補提示から確認・措置までの運用へ接続することです。
法人向けのご相談
数理工房では、製造業の保全履歴分析、設備・部品ネットワーク分析、異常検知、予知保全、推薦・グラフ機械学習、大規模行列計算について、データ整理からPoC、評価設計、現場運用までをご支援しています。
「アラーム履歴を保全判断へつなげたい」「類似設備の知見を横展開したい」「GNNやTransformerを使う妥当性を基準モデルから検証したい」「全工場展開に向け計算基盤を見直したい」といった課題をご相談いただけます。
📩 お問い合わせ: surikobo.co.jp/contact
まずはお気軽にご相談ください。