加工実績データを速く、正しく、再現可能に意思決定へつなぐ

CuPy・JAX・Polars・高速化・Notebook実践：100本ノック No.091〜No.100

製造現場のデータ活用では、計算を速くするだけでは不十分です。設備・日・シフトが増えても処理を継続でき、同じ入力から同じ結果を再現でき、数値誤差と処理時間を説明できて初めて、日々の点検判断に組み込めます。本記事では、架空の加工ラインを題材に、計算基盤の選択から高速化、可視化、検証、点検優先順位の作成までを一つの流れで扱います。

[!NOTE] 本資料は、数理工房 (もしくは代表である和山個人) が過去に企業研修において使用した notebook を企業様の許可を得て再構成・編集のうえ公開しています。
掲載データはすべて架空のものであり、実在する企業・工場・数値とは一切関係ありません。

はじめに：この記事で扱う製造業の実務課題

架空の工場では、加工サイクルごとに温度、振動、消費電力、加工時間、品質スコアを記録しています。朝会までに前日のデータを集計して点検候補を出したい一方、設備増設により処理時間が延び、担当者ごとにNotebookの結果が異なり、計算方式を変えた後の妥当性確認も不足しています。

そこで、処理量と計算資源を見積もり、表処理と行列計算を適材適所で高速化し、再現可能な分析記録と検証指標を残します。最終目的はライブラリの採用ではなく、限られた点検工数を、停止・不良リスクの高い設備へ配分することです。

現場でよくある状況

• 試作時は数秒だった集計が、全工場展開で朝会に間に合わない
• GPUや並列化を導入したが、転送・起動コストを含めると速くならない
• Notebookを上から実行せず、古い変数が残ったまま報告値を作る
• 高速化前後で丸め誤差や欠損処理が変わっても検知できない
• 平均処理時間だけを測り、日々のばらつきや最悪時間を見落とす
• グラフが多い一方、点検対象と判断基準が明示されていない

なぜこの問題は判断が難しいのか

処理時間はデータ件数だけでなく、行列形状、データ型、メモリアクセス、転送、コンパイル、並列粒度に左右されます。理論上の演算性能が高い基盤でも、小さな処理では準備コストが支配します。また高速化により演算順序が変わると、浮動小数点の結果は完全には一致しません。

したがって、速度を単一の値で評価せず、総所要時間

$$T_{\mathrm{total}}=T_{\mathrm{read}}+T_{\mathrm{prepare}}+T_{\mathrm{transfer}}+T_{\mathrm{compute}}+T_{\mathrm{report}}$$

と、許容誤差、再現性、保守性を同時に評価する必要があります。

今回扱うノックの全体像

No.	テーマ	製造業での判断
091	CuPy	GPU転送を含めて採用効果を見積もる
092	JAX	純粋関数・一括計算・微分をモデル検証へ使う
093	Polars	大量の加工実績を遅延・列指向で集計する
094	行列計算高速化	ループをベクトル化し、不要な中間配列を減らす
095	並列化	分割可能な設備別処理の粒度を決める
096	Notebook活用	パラメータ、検証、実行証跡を残す
097	可視化	異常度を点検優先順位へ翻訳する
098	数値誤差	桁落ち・悪条件を検知し誤判断を防ぐ
099	ベンチマーク	候補実装を公平に反復測定する
100	行列計算の世界	技術選択を統合して運用案を作る

Python 環境の準備

NumPyで行列計算、pandasとPolarsで表処理、Matplotlibで可視化を行います。CuPyとJAXはGPU・アクセラレータ環境への依存が大きいため必須にせず、利用可能なら検出する構成です。乱数生成器は np.random.default_rng(91) で固定します。

%matplotlib inline
%config InlineBackend.figure_format = 'svg'

import hashlib
import importlib.util
import platform
import sys
import time
from concurrent.futures import ThreadPoolExecutor

import japanize_matplotlib
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import polars as pl
from IPython.display import display

rng = np.random.default_rng(91)
plt.rcParams["figure.figsize"] = (8, 4.5)

print(f"Python     : {sys.version.split()[0]}")
print(f"NumPy      : {np.__version__}")
print(f"pandas     : {pd.__version__}")
print(f"Polars     : {pl.__version__}")
print(f"Matplotlib : {matplotlib.__version__}")
print(f"OS          : {platform.system()} {platform.machine()}")
print("CuPy       :", "利用可能" if importlib.util.find_spec("cupy") else "未導入（CPUで再現）")
print("JAX        :", "利用可能" if importlib.util.find_spec("jax") else "未導入（NumPyで原理を再現）")

Python     : 3.13.1
NumPy      : 2.5.1
pandas     : 3.0.3
Polars     : 1.42.1
Matplotlib : 3.11.0
OS          : Darwin arm64
CuPy       : 未導入（CPUで再現）
JAX        : 未導入（NumPyで原理を再現）

架空データの作成

4ライン・12設備について、60日、各設備20サイクルの架空データを生成します。経年傾向、夜勤負荷、設備固有差を温度・振動・消費電力・サイクル時間へ反映し、そこから不良フラグを作ります。以降の全ノックで同じデータを使い、技術比較の前提を揃えます。

実務では、設備ID・時刻・単位・校正履歴・欠損理由をデータ契約として管理し、「未測定」と正常値の0を区別します。

lines = [f"ライン{i}" for i in range(1, 5)]
equipment = [f"設備{i:02d}" for i in range(1, 13)]
n_days, cycles_per_equipment = 60, 20
n = n_days * len(equipment) * cycles_per_equipment

day = np.repeat(np.arange(n_days), len(equipment) * cycles_per_equipment)
equipment_id = np.tile(np.repeat(np.arange(len(equipment)), cycles_per_equipment), n_days)
shift = np.tile(np.arange(cycles_per_equipment) % 2, n_days * len(equipment))
line_id = equipment_id // 3
age = np.array([2, 5, 8, 3, 7, 11, 4, 9, 6, 12, 3, 10])[equipment_id]
load = np.clip(rng.normal(0.72 + 0.05 * shift, 0.09, n), 0.35, 1.0)
temperature = 43 + 13 * load + 0.18 * age + 0.025 * day + rng.normal(0, 1.7, n)
vibration = 1.3 + 1.4 * load + 0.09 * age + 0.012 * day + rng.normal(0, 0.28, n)
power = 14 + 22 * load + 0.35 * line_id + rng.normal(0, 1.4, n)
cycle_sec = 48 + 13 * load + 0.25 * age + rng.normal(0, 2.2, n)
quality_score = 100 - 0.30 * (temperature - 50) - 2.4 * (vibration - 2.5) - rng.normal(0, 1.0, n)
logit = -5.2 + 0.075 * (temperature - 50) + 0.85 * (vibration - 2.5) + 0.55 * shift
defect_prob = 1 / (1 + np.exp(-logit))
defect = rng.binomial(1, defect_prob)

production_df = pd.DataFrame({
    "日": day + 1, "ライン": np.array(lines)[line_id], "設備": np.array(equipment)[equipment_id],
    "シフト": np.where(shift == 0, "昼", "夜"), "負荷率": load, "温度_C": temperature,
    "振動_mm_s": vibration, "消費電力_kW": power, "サイクル秒": cycle_sec,
    "品質スコア": quality_score, "不良": defect,
})
features = ["負荷率", "温度_C", "振動_mm_s", "消費電力_kW", "サイクル秒"]
X = production_df[features].to_numpy(dtype=np.float64)
Xz = (X - X.mean(axis=0)) / X.std(axis=0)

display(production_df.head().style.format({c: "{:.2f}" for c in features + ["品質スコア"]}))
print(f"レコード数: {len(production_df):,} / 行列形状: {Xz.shape} / 不良率: {defect.mean():.2%}")

	日	ライン	設備	シフト	負荷率	温度_C	振動_mm_s	消費電力_kW	サイクル秒	品質スコア	不良
0	1	ライン1	設備01	昼	0.73	54.11	2.50	31.38	63.65	99.56	0
1	1	ライン1	設備01	夜	0.84	54.58	2.99	32.46	58.95	96.96	0
2	1	ライン1	設備01	昼	0.60	51.38	2.02	27.15	58.40	100.27	0
3	1	ライン1	設備01	夜	0.72	52.40	1.92	31.06	54.44	99.91	0
4	1	ライン1	設備01	昼	0.77	52.37	2.72	30.02	59.08	97.88	0

レコード数: 14,400 / 行列形状: (14400, 5) / 不良率: 2.28%

---

No.091：CuPy

実務での意味

CuPyはNumPyに近い記法でNVIDIA GPU上の配列計算を行います。大量の波形特徴量、画像、シミュレーションを同じ行列演算で繰り返す場合に候補になります。ただし、GPUへ送ってすぐ戻す処理では転送時間が便益を打ち消します。

分析・モデル化の考え方

CPUとGPUの境界をまたぐデータ量を $B$、転送帯域を $W$ とすると、往復転送時間の下限は概ね $2B/W$ です。ここでは xp という配列APIを切り替える形で同じリスクスコアを計算し、GPU実機がない環境ではNumPyへ安全にフォールバックします。採用時はウォームアップ後のエンドツーエンド時間を実測します。

Pythonで確認する

try:
    import cupy as cp
    xp, backend = cp, "CuPy (GPU)"
except ImportError:
    xp, backend = np, "NumPy (CPU fallback)"

weights = xp.asarray([0.25, 0.20, 0.30, 0.10, 0.15])
X_device = xp.asarray(Xz)
risk_device = X_device @ weights
risk_091 = cp.asnumpy(risk_device) if backend.startswith("CuPy") else np.asarray(risk_device)

sizes = np.array([10_000, 100_000, 1_000_000, 10_000_000])
bytes_roundtrip = sizes * len(features) * 8 * 2
transfer_plan = pd.DataFrame({
    "行数": sizes, "往復データ量_MB": bytes_roundtrip / 1e6,
    "転送下限_ms_帯域12GB毎秒": bytes_roundtrip / 12e9 * 1e3,
})
print("実行バックエンド:", backend)
display(transfer_plan.style.format({"往復データ量_MB": "{:.1f}", "転送下限_ms_帯域12GB毎秒": "{:.2f}"}))

plt.plot(sizes, transfer_plan["転送下限_ms_帯域12GB毎秒"], marker="o")
plt.xscale("log")
plt.title("GPU往復転送時間の下限試算（計算時間を含まない）")
plt.xlabel("レコード数")
plt.ylabel("転送時間の下限 [ms]")
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

実行バックエンド: NumPy (CPU fallback)

	行数	往復データ量_MB	転送下限_ms_帯域12GB毎秒
0	10000	0.8	0.07
1	100000	8.0	0.67
2	1000000	80.0	6.67
3	10000000	800.0	66.67

結果の読み取り

今回の実行バックエンドが明示され、GPUがなくても同じ計算結果を再現できます。試算はGPU優位を証明する値ではなく、転送だけで必要な最低時間です。実務では配列をGPU上に保持して複数処理をまとめ、CPU前処理と帳票作成まで含む総時間、GPUメモリ、運用費用を比較します。

---

No.092：JAX

実務での意味

JAXは、NumPyに近い関数をJITコンパイルし、自動微分や一括ベクトル化へ変換できます。品質予測のパラメータ調整や多数条件のシミュレーションで、同じ関数を繰り返し評価する場合に有効です。

分析・モデル化の考え方

変換しやすい計算は、外部状態を書き換えない純粋関数として記述します。ここではロジスティック損失 $L(w)$ の解析勾配 $\nabla L=X^\mathsf{T}(\sigma(Xw)-y)/n$ を有限差分と照合します。JAX導入時はこの関数を jit・grad・vmap の対象にし、初回コンパイルと2回目以降を分けて測定します。

Pythonで確認する

X_jax = np.column_stack([np.ones(len(Xz)), Xz])
y_jax = defect.astype(float)

def logistic_loss(w, X_input=X_jax, y_input=y_jax):
    z = X_input @ w
    return np.mean(np.logaddexp(0, z) - y_input * z)

def analytic_gradient(w):
    z = X_jax @ w
    p = 1 / (1 + np.exp(-z))
    return X_jax.T @ (p - y_jax) / len(y_jax)

w0 = np.zeros(X_jax.shape[1])
eps = 1e-6
grad_fd = np.array([(logistic_loss(w0 + eps * np.eye(len(w0))[j]) -
                     logistic_loss(w0 - eps * np.eye(len(w0))[j])) / (2 * eps)
                    for j in range(len(w0))])
grad_an = analytic_gradient(w0)
gradient_check = pd.DataFrame({"係数": ["切片"] + features, "解析勾配": grad_an, "有限差分": grad_fd,
                               "絶対差": np.abs(grad_an - grad_fd)})
display(gradient_check.style.format({"解析勾配": "{:.6f}", "有限差分": "{:.6f}", "絶対差": "{:.2e}"}))
print("勾配チェック最大誤差:", f"{np.max(np.abs(grad_an-grad_fd)):.2e}")

	係数	解析勾配	有限差分	絶対差
0	切片	0.477153	0.477153	9.97e-11
1	負荷率	-0.008225	-0.008225	4.26e-11
2	温度_C	-0.009080	-0.009080	4.42e-12
3	振動_mm_s	-0.010062	-0.010062	6.76e-11
4	消費電力_kW	-0.008843	-0.008843	2.26e-11
5	サイクル秒	-0.006579	-0.006579	1.38e-11

勾配チェック最大誤差: 9.97e-11

結果の読み取り

解析勾配と有限差分が十分小さい誤差で一致し、変換対象となる損失関数の実装を点検できました。JAXではさらに自動微分で保守性を高められますが、JITの初回コスト、配列形状変更による再コンパイル、乱数キー、64bit設定を管理し、学習結果だけでなく勾配チェックもテストに残します。

---

No.093：Polars

実務での意味

Polarsは列指向・並列実行・遅延評価を備え、大量の加工実績の抽出と集計に向きます。朝会用KPIの前処理を短くし、分析者がモデル検討へ使える時間を増やせます。

分析・モデル化の考え方

必要な列の選択、行の絞り込み、設備別集計を一つの遅延クエリとして表し、最後の collect() で実行します。実務ではParquetの列・行グループ削減が効くため、CSV読込だけの比較で結論を出しません。集計定義がpandas版と一致することを先に検証します。

Pythonで確認する

pl_df = pl.DataFrame(production_df.to_dict(orient="list"))
summary_pl = (
    pl_df.lazy()
    .filter(pl.col("日") > 30)
    .group_by(["ライン", "設備"])
    .agg([
        pl.len().alias("サイクル数"),
        pl.col("振動_mm_s").mean().alias("平均振動"),
        pl.col("不良").mean().alias("不良率"),
    ])
    .sort("不良率", descending=True)
    .collect()
)
summary_pd = (production_df.query("日 > 30").groupby(["ライン", "設備"], as_index=False)
              .agg(サイクル数=("不良", "size"), 平均振動=("振動_mm_s", "mean"), 不良率=("不良", "mean"))
              .sort_values("不良率", ascending=False))
check = np.allclose(summary_pl.sort(["ライン", "設備"])["不良率"].to_numpy(),
                    summary_pd.sort_values(["ライン", "設備"])["不良率"].to_numpy())
display(summary_pl.head(8))
print("pandas版との不良率一致:", check)

shape: (8, 5)

ライン	設備	サイクル数	平均振動	不良率
str	str	u32	f64	f64
”ライン4"	"設備10”	600	3.957084	0.046667
”ライン3"	"設備08”	600	3.701666	0.038333
”ライン2"	"設備05”	600	3.520651	0.036667
”ライン2"	"設備06”	600	3.855538	0.03
”ライン4"	"設備12”	600	3.766253	0.03
”ライン1"	"設備03”	600	3.598863	0.026667
”ライン3"	"設備07”	600	3.234954	0.023333
”ライン3"	"設備09”	600	3.436864	0.021667

pandas版との不良率一致: True

結果の読み取り

直近30日の設備別KPIが得られ、pandas版との数値一致も確認できました。速度比較の前に業務定義の同値性を確認することが重要です。本番ではスキーマを固定し、日時・欠損・カテゴリ型の解釈、遅延クエリの実行計画、入力ファイル数を監視します。

---

No.094：行列計算高速化

実務での意味

センサーごとの重み付き異常度を全サイクルへ計算する処理は、Pythonループより行列積へまとめる方が簡潔で高速です。朝会締切やオンライン監視周期に対する余裕を作れます。

分析・モデル化の考え方

各行のリスクスコアを $r_i=\sum_j X_{ij}w_j$ とします。ベクトル化により、ループ制御をPythonから最適化された行列演算へ移します。ただし巨大な中間配列を作る式や、データ型の暗黙変換はメモリ帯域を圧迫します。

Pythonで確認する

w_fast = np.array([0.25, 0.20, 0.30, 0.10, 0.15])

def score_loop(matrix, weights):
    out = np.empty(matrix.shape[0])
    for i in range(matrix.shape[0]):
        out[i] = sum(matrix[i, j] * weights[j] for j in range(matrix.shape[1]))
    return out

def score_vectorized(matrix, weights):
    return matrix @ weights

def median_time(func, *args, repeat=5):
    values = []
    for _ in range(repeat):
        start = time.perf_counter()
        func(*args)
        values.append(time.perf_counter() - start)
    return float(np.median(values))

t_loop = median_time(score_loop, Xz, w_fast)
t_vec = median_time(score_vectorized, Xz, w_fast)
risk_fast = score_vectorized(Xz, w_fast)
speed_df = pd.DataFrame({"実装": ["Pythonループ", "行列積"], "中央値_ms": [t_loop*1e3, t_vec*1e3]})
display(speed_df.style.format({"中央値_ms": "{:.3f}"}))
print(f"結果一致: {np.allclose(score_loop(Xz, w_fast), risk_fast)} / 高速化倍率: {t_loop/t_vec:.1f}倍")

	実装	中央値_ms
0	Pythonループ	15.081
1	行列積	0.097

結果一致: True / 高速化倍率: 155.1倍

結果の読み取り

この環境では行列積がループより高速で、結果も許容誤差内で一致しました。倍率はハードウェア、配列サイズ、BLAS、同時負荷で変わるため固定値として転用しません。まずアルゴリズムとデータ配置を改善し、それでもSLAを満たさない箇所だけをコンパイルやGPUの候補にします。

---

No.095：並列化

実務での意味

設備ごとの独立した集計やシミュレーションは並列化できます。ただし小さな仕事を細かく分けすぎると、分割・起動・結合の時間が増え、運用も複雑になります。

分析・モデル化の考え方

総時間を $T_p\approx T_s+T_{parallel}/p+T_{overhead}$ と捉えます。ここではNumPyが内部計算でGILを解放し得る設備別のSVD処理を、逐次とスレッド並列で比較します。結果順序を設備IDで固定し、値の一致も検査します。

Pythonで確認する

matrices = {name: rng.normal(size=(220, 80)) for name in equipment[:8]}

def equipment_job(item):
    name, matrix = item
    singular_values = np.linalg.svd(matrix, compute_uv=False)
    return name, singular_values[:3]

def run_serial(items):
    return [equipment_job(item) for item in items]

def run_parallel(items):
    with ThreadPoolExecutor(max_workers=4) as executor:
        return list(executor.map(equipment_job, items))

items = list(matrices.items())
t_serial = median_time(run_serial, items, repeat=3)
t_parallel = median_time(run_parallel, items, repeat=3)
serial_result, parallel_result = run_serial(items), run_parallel(items)
parallel_df = pd.DataFrame({"方式": ["逐次", "4スレッド"], "中央値_ms": [t_serial*1e3, t_parallel*1e3]})
display(parallel_df.style.format({"中央値_ms": "{:.2f}"}))
print("設備順と結果の一致:", all(a[0] == b[0] and np.allclose(a[1], b[1]) for a, b in zip(serial_result, parallel_result)))

	方式	中央値_ms
0	逐次	4.48
1	4スレッド	4.92

設備順と結果の一致: True

結果の読み取り

並列化の効果は実行環境のCPU数やBLAS内部スレッドと競合するため、必ずしも4倍にはなりません。今回の表も採否を決める実測例です。実務では設備単位など意味のある粒度で分割し、タイムアウト、部分失敗、再実行、メモリ上限、出力順序を設計します。

---

No.096：Notebook活用

実務での意味

Notebookはコード、説明、表、グラフを一体化でき、製造条件の検討記録に向きます。一方、セル順序への依存を放置すると、同じファイルでも報告値が再現しません。

分析・モデル化の考え方

入力パラメータを一か所へ集約し、データ指紋、環境、検証結果を実行証跡として残します。セルを上から全実行できること、入力行数・範囲・欠損率・KPIを assert で検査することを最低条件にします。

Pythonで確認する

PARAMS = {"analysis_day_from": 46, "risk_threshold": 0.85, "seed": 91, "model_version": "risk-v1"}
fingerprint_cols = ["日", "設備", "温度_C", "振動_mm_s", "不良"]
data_fingerprint = hashlib.sha256(
    pd.util.hash_pandas_object(production_df[fingerprint_cols], index=True).values.tobytes()
).hexdigest()[:16]

validation = {
    "行数が期待値": len(production_df) == n,
    "主要列に欠損なし": not production_df[features + ["不良"]].isna().any().any(),
    "負荷率が0〜1": production_df["負荷率"].between(0, 1).all(),
    "不良が0/1": production_df["不良"].isin([0, 1]).all(),
}
assert all(validation.values()), validation
run_log = pd.DataFrame({
    "項目": ["モデル版", "乱数seed", "分析開始日", "入力行数", "データ指紋", "検証"],
    "値": [PARAMS["model_version"], PARAMS["seed"], PARAMS["analysis_day_from"], len(production_df),
           data_fingerprint, f"{sum(validation.values())}/{len(validation)} 合格"],
})
display(run_log)

	項目	値
0	モデル版	risk-v1
1	乱数seed	91
2	分析開始日	46
3	入力行数	14400
4	データ指紋	7731a861f84f2ce0
5	検証	4/4 合格

結果の読み取り

パラメータとデータ指紋、4件の検証結果が一つの証跡になりました。指紋はデータ内容の同一性確認であり、品質保証そのものではありません。本番ではGitのコミット、依存関係ロック、実行日時、承認者、成果物保存先も記録し、Notebookから定期ジョブへ移す境界を決めます。

---

No.097：可視化

実務での意味

可視化の目的はきれいな図ではなく、どの設備を、なぜ、いつ点検するかを共通認識にすることです。平均だけでなく振動・温度・不良率を並べ、複数兆候の重なりを見ます。

分析・モデル化の考え方

直近15日の設備別KPIを標準化してヒートマップにします。色は設備間の相対比較であり、保全限界ではありません。絶対閾値、過去推移、サンプル数を併記して、色だけで判断しない設計にします。

Pythonで確認する

analysis_day_from = PARAMS["analysis_day_from"]
recent = production_df.query("日 >= @analysis_day_from")
kpi = recent.groupby("設備").agg(
    平均温度=("温度_C", "mean"), 最大振動=("振動_mm_s", "max"),
    平均サイクル秒=("サイクル秒", "mean"), 不良率=("不良", "mean"),
)
kpi_z = (kpi - kpi.mean()) / kpi.std(ddof=0)

fig, ax = plt.subplots(figsize=(8.5, 5.0))
im = ax.imshow(kpi_z.to_numpy(), cmap="RdYlBu_r", aspect="auto", vmin=-2, vmax=2)
ax.set_xticks(range(len(kpi_z.columns)), kpi_z.columns, rotation=25, ha="right")
ax.set_yticks(range(len(kpi_z.index)), kpi_z.index)
ax.set_title("直近15日の設備別KPI（列ごとの標準化値）")
ax.set_xlabel("KPI")
ax.set_ylabel("設備")
ax.grid(False)
fig.colorbar(im, ax=ax, label="標準化値")
plt.tight_layout()
plt.show()
display(kpi.sort_values("不良率", ascending=False).head(5).style.format("{:.3f}"))

	平均温度	最大振動	平均サイクル秒	不良率
設備
設備08	55.818	4.975	59.931	0.050
設備10	56.251	4.962	60.728	0.050
設備03	55.448	4.695	59.869	0.047
設備05	55.247	4.530	59.407	0.040
設備11	54.651	4.318	58.522	0.030

結果の読み取り

赤いセルが複数KPIで重なる設備は、点検候補を絞る入口になります。ただし標準化値は集団が変わると動きます。設備仕様に基づく警報値、前日差、長期傾向、測定品質を確認し、グラフから自動的に原因を断定しません。

---

No.098：数値誤差

実務での意味

センサー補正、微小差の累積、連立方程式では、表示上は同じ入力でも計算結果が大きく変わることがあります。誤差を故障兆候と取り違えると、不要点検や見逃しにつながります。

分析・モデル化の考え方

浮動小数点では結合法則が厳密には成立しません。また $Ax=b$ で条件数 $\kappa(A)$ が大きいと、入力の小さな誤差が解で増幅されます。逆行列を明示的に作らず solve を使い、条件数と残差 $\lVert Ax-b\rVert/\lVert b\rVert$ を併記します。

Pythonで確認する

cancel_values = np.array([1e8, 1.0, -1e8], dtype=np.float32)
sum_forward = cancel_values.sum(dtype=np.float32)
sum_reordered = cancel_values[[0, 2, 1]].sum(dtype=np.float32)

A_good = np.array([[1.0, 0.2], [0.2, 1.0]])
A_bad = np.array([[1.0, 1.0], [1.0, 1.0 + 1e-10]])
b = np.array([2.0, 2.0 + 1e-10])
rows = []
for label, A in [("安定した校正行列", A_good), ("ほぼ重複した校正行列", A_bad)]:
    solution = np.linalg.solve(A, b)
    residual = np.linalg.norm(A @ solution - b) / np.linalg.norm(b)
    rows.append({"行列": label, "条件数": np.linalg.cond(A), "相対残差": residual,
                 "解1": solution[0], "解2": solution[1]})
error_df = pd.DataFrame(rows)
print(f"float32の加算順序: {sum_forward:.1f} と {sum_reordered:.1f}")
display(error_df.style.format({"条件数": "{:.2e}", "相対残差": "{:.2e}", "解1": "{:.4f}", "解2": "{:.4f}"}))

float32の加算順序: 0.0 と 1.0

	行列	条件数	相対残差	解1	解2
0	安定した校正行列	1.50e+00	1.57e-16	1.6667	1.6667
1	ほぼ重複した校正行列	4.00e+10	0.00e+00	1.0000	1.0000

結果の読み取り

float32では加算順序だけで結果が変わり、ほぼ重複した校正行列は非常に大きな条件数を持ちます。残差が小さくても解が信頼できるとは限りません。単位と尺度を揃え、float64、安定な分解、再校正、正則化を検討し、業務上許容できる差をテストにします。

---

No.099：ベンチマーク

実務での意味

ベンチマークは、GPUや新ライブラリを採用する根拠を作ります。一度だけの最速値ではなく、実データに近い形状で反復し、中央値とばらつき、結果一致を比較します。

分析・モデル化の考え方

初回準備を分け、同じ入力・同じ出力条件で候補を交互に反復します。ここではリスク計算のPythonループと行列積を複数回測り、中央値、四分位範囲、95パーセンタイルを報告します。小さすぎる処理では測定器の分解能にも注意します。

Pythonで確認する

def benchmark(func, *args, repeat=12):
    func(*args)  # warm-up
    samples = []
    for _ in range(repeat):
        start = time.perf_counter_ns()
        func(*args)
        samples.append((time.perf_counter_ns() - start) / 1e6)
    return np.array(samples)

bench = {"Pythonループ": benchmark(score_loop, Xz, w_fast),
         "行列積": benchmark(score_vectorized, Xz, w_fast)}
bench_df = pd.DataFrame([
    {"実装": name, "中央値_ms": np.median(v), "IQR_ms": np.percentile(v, 75)-np.percentile(v, 25),
     "p95_ms": np.percentile(v, 95)} for name, v in bench.items()
])
display(bench_df.style.format({"中央値_ms": "{:.3f}", "IQR_ms": "{:.3f}", "p95_ms": "{:.3f}"}))

plt.boxplot([bench[k] for k in bench], tick_labels=list(bench))
plt.title("リスク計算の反復ベンチマーク")
plt.xlabel("実装")
plt.ylabel("処理時間 [ms]")
plt.grid(True, axis="y", alpha=0.3)
plt.tight_layout()
plt.show()

	実装	中央値_ms	IQR_ms	p95_ms
0	Pythonループ	14.968	0.267	15.325
1	行列積	0.005	0.001	0.009

結果の読み取り

中央値に加えてIQRとp95から安定性を確認できます。この結果は現在の端末・データ形状に限定されます。本番候補は読込、前処理、計算、出力を含むSLAで評価し、CIでは性能劣化の傾向を監視します。数値一致、メモリ、費用、保守性も採用基準に含めます。

---

No.100：行列計算の世界

実務での意味

行列計算の実務価値は、手法を増やすことではなく、データを共通表現へ整理し、計算結果を具体的な行動へ接続することにあります。最後に、直近データから設備リスクを統合し、翌日の点検候補を作ります。

分析・モデル化の考え方

設備 $e$ の優先度を、標準化した平均リスク、最大振動、不良率から $P_e=0.45R_e+0.35V_e+0.20D_e$ と定義します。重みは例示であり、実務では停止損失、安全、品質影響、点検時間を含む目的関数として関係者と合意します。上位候補には寄与KPIも添えます。

Pythonで確認する

risk_series = pd.Series(risk_fast, index=production_df.index, name="行列リスク")
decision = recent.assign(行列リスク=risk_series.loc[recent.index]).groupby("設備").agg(
    平均リスク=("行列リスク", "mean"), 最大振動=("振動_mm_s", "max"), 不良率=("不良", "mean"),
    対象サイクル=("不良", "size"),
)
decision_z = (decision[["平均リスク", "最大振動", "不良率"]] -
              decision[["平均リスク", "最大振動", "不良率"]].mean()) /              decision[["平均リスク", "最大振動", "不良率"]].std(ddof=0)
decision["点検優先度"] = decision_z @ np.array([0.45, 0.35, 0.20])
decision["主な寄与KPI"] = decision_z.idxmax(axis=1)
priority = decision.sort_values("点検優先度", ascending=False)
display(priority.head(5).style.format({"平均リスク": "{:.3f}", "最大振動": "{:.3f}", "不良率": "{:.2%}", "点検優先度": "{:.3f}"}))

top = priority.head(6).sort_values("点検優先度")
plt.barh(top.index, top["点検優先度"], color="#D55E00")
plt.title("翌日の点検優先候補（架空データ）")
plt.xlabel("点検優先度")
plt.ylabel("設備")
plt.grid(True, axis="x", alpha=0.3)
plt.tight_layout()
plt.show()

	平均リスク	最大振動	不良率	対象サイクル	点検優先度	主な寄与KPI
設備
設備10	0.748	4.962	5.00%	300	1.582	平均リスク
設備08	0.487	4.975	5.00%	300	1.195	不良率
設備06	0.585	4.782	2.67%	300	0.777	平均リスク
設備12	0.518	4.842	2.67%	300	0.740	最大振動
設備03	0.340	4.695	4.67%	300	0.612	不良率

結果の読み取り

上位設備と主な寄与KPIが表示され、点検班が確認すべき候補を絞れます。優先度は故障確率ではなく、架空の相対スコアです。現場導入では安全上の必須点検を別ルールで優先し、点検結果を記録して重みと閾値を更新します。行列計算は意思決定を代替するのではなく、確認順序を一貫させる支援として使います。

---

対象ノックを通して見える実務上の示唆

1. 基盤は処理全体で選ぶ：CuPyやJAXの演算性能だけでなく、転送、コンパイル、前後処理を含めます。
2. 表処理と行列計算を分担する：Polarsで必要な行・列へ絞り、NumPy等で数値計算する境界を明確にします。
3. 先にアルゴリズムを改善する：ベクトル化と不要コピー削減を確認してから、GPU・並列化を検討します。
4. 再現性は成果物の一部である：入力指紋、パラメータ、環境、検証結果をNotebookに残します。
5. 速さと正しさを同時に測る：ベンチマークには結果一致、許容誤差、中央値、ばらつきを含めます。
6. 可視化を行動へ接続する：点検候補、寄与KPI、確認期限、担当を明確にします。
7. 高度な技術を使わない判断も価値がある：データ量とSLAが小さければ、保守しやすいCPU処理が合理的です。

実務導入する場合に必要なこと

1. SLAと判断単位を定義する

「毎朝7時まで」「設備停止後30秒以内」など期限を定め、対象設備、更新頻度、許容遅延、同時利用者を数値化します。高速化率ではなく、意思決定に間に合うかで評価します。

2. 正解・許容誤差・性能基準を固定する

基準実装と代表データを保存し、KPI、順位、欠損処理、数値許容差を自動検査します。中央値とp95、ピークメモリ、失敗率、再実行時間を受入基準にします。

3. 実行環境とデータを版管理する

Python・ライブラリ・ドライバ、乱数seed、入力スキーマ、設備台帳、モデル重みを記録します。GPUを使う場合はドライバ互換性とCPUフォールバック、Notebookは全セル再実行を確認します。

4. 小さく運用検証する

一つのラインで、候補提示、現場確認、措置、結果記録まで試します。見逃し・不要点検・短縮時間・停止損失を評価し、全工場展開の計算資源と支援体制を見積もります。

まとめ

No.091〜No.100では、架空の加工実績を使い、CuPyのGPU転送見積もり、JAXにつながる純粋関数と勾配検証、Polarsの遅延集計、行列計算高速化、並列化、Notebookの再現性、意思決定向け可視化、数値誤差、反復ベンチマーク、点検優先順位への統合を確認しました。

製造業で重要なのは、最速のライブラリを選ぶことではありません。期限内に、許容誤差で、再現可能な結果を出し、その根拠と限界を示して現場の行動へつなぐことです。

法人向けのご相談

数理工房では、製造業の加工・品質・設備データについて、データ基盤整理、Python処理の高速化、GPU・並列計算の採否評価、Notebookの標準化、ベンチマーク設計、点検・品質判断への可視化までをご支援しています。

「全工場展開で処理が間に合わない」「GPU導入効果を投資前に検証したい」「Notebook分析を再現可能な業務へ移したい」「速度改善後の数値差をどう検証すべきか」といった課題をご相談いただけます。

📩 お問い合わせ: surikobo.co.jp/contact
まずはお気軽にご相談ください。