数理工房

100本ノック / 数値計算 / 数値計算100本ノック

製造業の固有値問題入門｜設備振動・PCA・工程ネットワークをPythonで分析

設備振動と工程ネットワークを「固有値」で読み解く：製造業の意思決定に効く数値計算10本ノック（No.041〜No.050）

本記事では、固有値問題、反復法、行列分解、グラフ解析を、架空の精密部品工場における設備保全・センサー監視・工程改善へ結び付けます。数式を計算すること自体ではなく、「どの変動を優先して調べるか」「大規模データで何を近似するか」「停止リスクが集中する工程はどこか」を判断できることが目的です。

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

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

対象工場では、複数設備の振動センサーが相関して変動し、工程間では仕掛品、再加工、検査情報が循環しています。個別の平均や上限値だけでは、設備群に共通する異常モードや、ネットワーク全体で重要な工程を見落とします。本記事では行列の「主要な方向」を抽出し、有限な保全工数をどこへ振り向けるべきかを考えます。

現場でよくある状況

  • センサー点数が増え、アラーム件数だけが増えている
  • 相関した振動を設備ごとの単独閾値で監視している
  • 全固有値を厳密計算しているが、意思決定に使うのは上位数個だけである
  • 工程表はあるものの、再加工ループや情報依存を含む全体影響を評価できていない

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

行列の各要素は局所的な関係ですが、固有値・固有ベクトルは関係が繰り返し伝播した結果を表します。一方、行列が大きいほど、精度・計算時間・メモリの間にトレードオフが生じます。さらに、数学的に大きな成分が、そのまま因果関係や投資効果を意味するわけではありません。現場知識、データ品質、停止損失を合わせて解釈する必要があります。

今回扱うノックの全体像

No.テーマ製造業での判断
041固有値問題とは支配的な設備変動を特定
042べき乗法最大変動モードを高速推定
043逆反復法注目周波数付近のモードを抽出
044QR法小中規模行列の固有値を一括計算
045ランチョス法大規模対称問題を低次元化
046アーノルディ法非対称な工程伝播を近似
047特異値分解センサー情報を圧縮・再構成
048PCAとの関係監視指標と寄与率を設計
049PageRankとの関係工程影響度を再帰的に評価
050グラフ解析への応用分断・ボトルネック候補を発見

Python 環境の準備

NumPy で密行列、SciPy で疎行列と反復法、pandas で表、matplotlib でグラフを扱います。乱数シードを固定して再現性を確保します。

%matplotlib inline
import platform
import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
from scipy import linalg, sparse
from scipy.sparse import linalg as spla

SEED = 42
rng = np.random.default_rng(SEED)
pd.set_option("display.precision", 4)
print("Python     :", platform.python_version())
print("NumPy      :", np.__version__)
print("pandas     :", pd.__version__)
print("SciPy      :", __import__("scipy").__version__)
print("matplotlib :", matplotlib.__version__)
Python     : 3.13.1
NumPy      : 2.5.1
pandas     : 3.0.3
SciPy      : 1.18.0
matplotlib : 3.11.0

架空データの作成

120時間分、6点の振動センサーを作ります。共通の「回転アンバランス」と「支持部の緩み」に相当する潜在変動を混ぜ、後半20時間には駆動側の振幅増加を加えます。また、8工程の有向ネットワークを、通常搬送と再加工の比率から作ります。値は説明用の架空データです。

n_time = 120
sensor_names = ["Drive-X", "Drive-Y", "Bearing-A", "Bearing-B", "Frame", "Outlet"]
t = np.arange(n_time)
mode_1 = np.sin(2 * np.pi * t / 18) + 0.18 * rng.normal(size=n_time)
mode_2 = 0.65 * np.sin(2 * np.pi * t / 7 + 0.8) + 0.18 * rng.normal(size=n_time)
loadings = np.array([
    [1.00, 0.15], [0.88, 0.10], [0.58, 0.72],
    [0.52, 0.78], [0.30, 0.45], [0.20, 0.28]
])
X = np.column_stack([mode_1, mode_2]) @ loadings.T + 0.16 * rng.normal(size=(n_time, 6))
X[-20:, :2] *= 1.35
sensor_df = pd.DataFrame(X, columns=sensor_names)

processes = ["Material", "Machining", "HeatTreat", "Grinding", "Wash", "Inspection", "Rework", "Shipping"]
edges = [
    (0,1,1.00),(1,2,.82),(1,6,.18),(2,3,.92),(2,6,.08),(3,4,.88),
    (3,6,.12),(4,5,1.00),(5,7,.86),(5,6,.14),(6,1,.55),(6,3,.45)
]
A = np.zeros((8, 8))
for i, j, w in edges: A[i, j] = w

display(sensor_df.head().round(3))
print(f"sensor data shape: {sensor_df.shape}, process matrix shape: {A.shape}")
Drive-X Drive-Y Bearing-A Bearing-B Frame Outlet
0 -0.043 0.061 -0.046 0.014 0.484 -0.116
1 0.081 0.498 1.042 0.421 0.292 0.275
2 1.111 0.564 0.684 0.823 0.473 0.201
3 1.017 0.693 0.578 0.512 0.357 0.124
4 0.639 0.673 0.057 0.022 -0.011 -0.004 
sensor data shape: (120, 6), process matrix shape: (8, 8)

No.041：固有値問題とは

実務での意味

共分散行列の固有ベクトルは、複数センサーが同時に動く代表パターンです。固有値はそのパターンに沿った変動の大きさを表します。設備群のどの組合せが支配的かを整理でき、点検の優先順位付けに使えます。

分析・モデル化の考え方

正方行列 CC に対し、ゼロでないベクトル vv とスカラー λ\lambda が

Cv=λvCv=\lambda v

を満たすとき、λ\lambda を固有値、vv を固有ベクトルと呼びます。ここでは標準化後の相関行列を使うため、センサー単位の違いに左右されにくくなります。ただし固有ベクトルの符号は反転しても同じ解です。

Pythonで確認する

Z = (sensor_df - sensor_df.mean()) / sensor_df.std(ddof=1)
C = Z.corr().to_numpy()
eigvals, eigvecs = np.linalg.eigh(C)
order = np.argsort(eigvals)[::-1]
eigvals, eigvecs = eigvals[order], eigvecs[:, order]
eig_table = pd.DataFrame({"eigenvalue": eigvals, "variance_ratio": eigvals / eigvals.sum()})
display(eig_table.round(4))

plt.figure(figsize=(7, 3.5))
plt.bar(np.arange(1, 7), eigvals, color="#3568a8")
plt.title("Eigenvalues of sensor correlation matrix")
plt.xlabel("Component")
plt.ylabel("Eigenvalue")
plt.grid(True, axis="y", alpha=0.3)
plt.tight_layout()
plt.show()
eigenvalue variance_ratio
0 4.6865 0.7811
1 0.6686 0.1114
2 0.3543 0.0590
3 0.1664 0.0277
4 0.0638 0.0106
5 0.0604 0.0101

png

結果の読み取り

第1・第2固有値が残りより大きければ、6本のセンサー信号を少数の共通変動で説明できる可能性があります。個別アラームを増やす前に、上位固有ベクトルで負荷の大きいセンサー群を同時点検する判断が有効です。固有値だけで故障原因は確定できないため、回転数、負荷、整備履歴との照合が必要です。

No.042：べき乗法

実務での意味

最大固有値と対応ベクトルだけが欲しい場合、全固有値計算は過剰です。べき乗法は支配的な振動モードを、行列とベクトルの積の反復だけで推定します。

分析・モデル化の考え方

xk+1=Cxk/∥Cxk∥2x_{k+1}=Cx_k/\|Cx_k\|_2 を反復すると、初期値が最大固有ベクトルと直交せず、最大固有値の絶対値が一意なら、その方向へ収束します。固有値はレイリー商 λk=xkTCxk/(xkTxk)\lambda_k=x_k^TCx_k/(x_k^Tx_k) で見積もります。収束速度は ∣λ2/λ1∣|\lambda_2/\lambda_1| に左右されます。

Pythonで確認する

x = np.ones(C.shape[0]) / np.sqrt(C.shape[0])
history = []
for k in range(30):
    x_new = C @ x
    x_new /= np.linalg.norm(x_new)
    lam = float(x_new @ C @ x_new)
    history.append(lam)
    if np.linalg.norm(x_new - x) < 1e-10 or np.linalg.norm(x_new + x) < 1e-10:
        x = x_new
        break
    x = x_new

display(pd.DataFrame({"sensor": sensor_names, "power_loading": x}).round(4))
print(f"iterations={len(history)}, estimate={lam:.6f}, exact={eigvals[0]:.6f}")
plt.figure(figsize=(7, 3.5))
plt.plot(range(1, len(history)+1), history, marker="o", ms=3)
plt.axhline(eigvals[0], color="crimson", linestyle="--", label="exact")
plt.title("Convergence of power iteration")
plt.xlabel("Iteration")
plt.ylabel("Estimated dominant eigenvalue")
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()
sensor power_loading
0 Drive-X 0.4026
1 Drive-Y 0.3938
2 Bearing-A 0.4377
3 Bearing-B 0.4253
4 Frame 0.4090
5 Outlet 0.3782 
iterations=11, estimate=4.686548, exact=4.686548

png

結果の読み取り

推定値が厳密計算値へ速やかに近づけば、日次監視で最大モードだけを更新する用途に適します。上位2固有値が近い設備では収束が遅くなるため、反復上限と残差 ∥Cx−λx∥\|Cx-\lambda x\| を監視し、必要ならランチョス法へ切り替えます。

No.043：逆反復法

実務での意味

共振試験では最大モードではなく、設計上注意する周波数に近いモードが重要です。逆反復法は指定値（シフト）付近の固有値を選択的に抽出します。

分析・モデル化の考え方

シフト μ\mu に対して (C−μI)yk=xk(C-\mu I)y_k=x_k を解き、xk+1=yk/∥yk∥x_{k+1}=y_k/\|y_k\| とします。(C−μI)−1(C-\mu I)^{-1} では μ\mu に最も近い固有値の成分が支配的になります。シフトが固有値そのものだと行列が特異になるため、数値的な条件を確認します。

Pythonで確認する

shift = 0.50
x_inv = np.ones(C.shape[0]) / np.sqrt(C.shape[0])
inv_history = []
for k in range(20):
    y = np.linalg.solve(C - shift * np.eye(C.shape[0]), x_inv)
    x_inv = y / np.linalg.norm(y)
    lam_inv = float(x_inv @ C @ x_inv)
    residual = np.linalg.norm(C @ x_inv - lam_inv * x_inv)
    inv_history.append((k + 1, lam_inv, residual))
    if residual < 1e-10: break
display(pd.DataFrame(inv_history, columns=["iteration", "eigenvalue", "residual"]).tail().round(8))
print("nearest exact eigenvalue:", eigvals[np.argmin(np.abs(eigvals - shift))].round(6))
iteration eigenvalue residual
15 16 0.3547 0.0110
16 17 0.3546 0.0095
17 18 0.3545 0.0082
18 19 0.3545 0.0071
19 20 0.3544 0.0062 
nearest exact eigenvalue: 0.354291

結果の読み取り

シフトに最も近い固有値へ収束するため、対象帯域を決めたモード探索に向きます。実務では「狙った帯域に固有値がある」だけでなく、対応ベクトルのセンサー位置と、回転数由来の加振周波数が重なるかを確認します。

No.044：QR法

実務での意味

小中規模のモデルで全固有値を把握したいとき、QR法は基本となる手法です。設備モデルの全モードを並べ、危険帯域との接近を確認できます。

分析・モデル化の考え方

Ak=QkRkA_k=Q_kR_k とQR分解し、Ak+1=RkQk=QkTAkQkA_{k+1}=R_kQ_k=Q_k^TA_kQ_k と更新します。相似変換なので固有値は保存され、対称行列では反復により対角成分が固有値へ収束します。実用ライブラリはシフトやHessenberg化で高速化しています。

Pythonで確認する

A_qr = C.copy()
offdiag_history = []
for k in range(80):
    Q, R = np.linalg.qr(A_qr)
    A_qr = R @ Q
    offdiag_history.append(np.linalg.norm(A_qr - np.diag(np.diag(A_qr))))
qr_vals = np.sort(np.diag(A_qr))[::-1]
display(pd.DataFrame({"QR_iteration": qr_vals, "numpy_eigh": eigvals, "abs_error": np.abs(qr_vals-eigvals)}).round(8))
plt.figure(figsize=(7, 3.5))
plt.semilogy(range(1, 81), offdiag_history)
plt.title("QR iteration: off-diagonal norm")
plt.xlabel("Iteration")
plt.ylabel("Off-diagonal Frobenius norm")
plt.grid(True, which="both", alpha=0.3)
plt.tight_layout()
plt.show()
QR_iteration numpy_eigh abs_error
0 4.6865 4.6865 0.0000e+00
1 0.6686 0.6686 0.0000e+00
2 0.3543 0.3543 0.0000e+00
3 0.1664 0.1664 0.0000e+00
4 0.0638 0.0638 3.0000e-08
5 0.0604 0.0604 3.0000e-08

png

結果の読み取り

非対角成分のノルムが低下し、対角値が既製関数の固有値と一致することを確認できます。教材としての単純QR法は仕組みの確認用です。本番では自作反復ではなく、十分に検証された numpy.linalg.eigh などを使用します。

No.045：ランチョス法

実務での意味

有限要素モデルや多数センサーの共分散行列は大規模でも、必要なのは低次の数モードだけということがあります。ランチョス法は対称行列を小さな三重対角行列へ射影します。

分析・モデル化の考え方

Krylov部分空間 Km(A,q)=span(q,Aq,…,Am−1q)\mathcal{K}_m(A,q)=\mathrm{span}(q,Aq,\ldots,A^{m-1}q) に基底を作り、Tm=QmTAQmT_m=Q_m^TAQ_m の固有値（Ritz値）で元の固有値を近似します。有限精度では直交性が失われるため、実装では再直交化が重要です。

Pythonで確認する

n = 240
main = 2.2 + 0.15 * rng.random(n)
off = -1.0 * np.ones(n-1)
K = sparse.diags([off, main, off], [-1, 0, 1], format="csr")

def lanczos(A, m, seed=42):
    local_rng = np.random.default_rng(seed)
    q = local_rng.normal(size=A.shape[0]); q /= np.linalg.norm(q)
    Q, alphas, betas = [], [], []
    q_prev = np.zeros_like(q); beta = 0.0
    for j in range(m):
        Q.append(q.copy())
        z = A @ q - beta * q_prev
        alpha = float(q @ z); z -= alpha * q
        # 教材でも数値直交性を保つため、既存基底へ再直交化する
        for qi in Q: z -= (qi @ z) * qi
        alphas.append(alpha)
        beta_new = np.linalg.norm(z)
        if beta_new < 1e-12: break
        if j < m-1: betas.append(beta_new)
        q_prev, q, beta = q, z / beta_new, beta_new
    return np.array(alphas), np.array(betas)

rows = []
exact_large = spla.eigsh(K, k=1, which="LA", return_eigenvectors=False)[0]
for m in [5, 10, 20, 30]:
    alpha, beta = lanczos(K, m)
    T = np.diag(alpha) + np.diag(beta, 1) + np.diag(beta, -1)
    ritz = np.linalg.eigvalsh(T)[-1]
    rows.append((m, ritz, abs(ritz-exact_large)))
display(pd.DataFrame(rows, columns=["Krylov_dim", "largest_Ritz", "abs_error"]).round(8))
Krylov_dim largest_Ritz abs_error
0 5 4.1556 0.1362
1 10 4.2555 0.0363
2 20 4.2772 0.0145
3 30 4.2846 0.0071

結果の読み取り

240次元の行列でも、はるかに小さいKrylov部分空間から最大固有値を近似できます。設備モデルでは、必要モード数、残差、計算時間を受入基準にします。対称性が崩れるモデルには、そのまま適用せずアーノルディ法などを選びます。

No.046：アーノルディ法

実務での意味

工程間の流れは上流・下流があり、行列は一般に非対称です。アーノルディ法は、再加工や遅延が方向を持って伝播するモデルの主要固有値を近似できます。

分析・モデル化の考え方

アーノルディ法はKrylov部分空間の正規直交基底 QmQ_m を作り、AQm≈QmHmAQ_m\approx Q_mH_m となる上Hessenberg行列 HmH_m へ射影します。ランチョス法を非対称行列へ拡張した位置付けです。複素固有値が現れる場合は、絶対値と偏角の両方を解釈します。

Pythonで確認する

def arnoldi(A, m):
    n = A.shape[0]
    Q = np.zeros((n, m+1)); H = np.zeros((m+1, m))
    Q[:, 0] = np.ones(n) / np.sqrt(n)
    for k in range(m):
        v = A @ Q[:, k]
        for j in range(k+1):
            H[j, k] = Q[:, j] @ v
            v -= H[j, k] * Q[:, j]
        H[k+1, k] = np.linalg.norm(v)
        if H[k+1, k] < 1e-12: return H[:k+1, :k+1]
        Q[:, k+1] = v / H[k+1, k]
    return H[:m, :m]

H = arnoldi(A.T, 6)
ritz = np.linalg.eigvals(H)
exact = np.linalg.eigvals(A.T)
arnoldi_table = pd.DataFrame({
    "method": ["Arnoldi Ritz", "Exact"],
    "largest_abs_eigenvalue": [np.max(np.abs(ritz)), np.max(np.abs(exact))]
})
display(arnoldi_table.round(6))

plt.figure(figsize=(6, 4))
plt.scatter(exact.real, exact.imag, label="Exact", s=55)
plt.scatter(ritz.real, ritz.imag, marker="x", s=70, label="Arnoldi Ritz")
plt.title("Eigenvalues of directed process matrix")
plt.xlabel("Real part")
plt.ylabel("Imaginary part")
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()
method largest_abs_eigenvalue
0 Arnoldi Ritz 0.8004
1 Exact 0.7457

png

結果の読み取り

小さなHessenberg行列のRitz値が、元行列の外側の固有値を捉える様子を確認できます。工程モデルで絶対値1に近い固有値があると、影響が循環して減衰しにくい可能性があります。ただし行列の正規化規則で意味が変わるため、重量が確率か数量かを明示します。

No.047：特異値分解

実務での意味

センサーデータは長方形行列です。特異値分解（SVD）は、時間方向とセンサー方向の代表パターンへ分け、保存容量の削減、ノイズ除去、代表波形の抽出に使えます。

分析・モデル化の考え方

Z=UΣVTZ=U\Sigma V^T と分解します。特異値 σi\sigma_i は非負で、ZTZZ^TZ の固有値との間に σi2=λi\sigma_i^2=\lambda_i があります。上位 kk 成分による Zk=UkΣkVkTZ_k=U_k\Sigma_kV_k^T は、Frobeniusノルムの意味で最良の階数 kk 近似です。

Pythonで確認する

U, s, Vt = np.linalg.svd(Z.to_numpy(), full_matrices=False)
svd_rows = []
for k in range(1, 7):
    Zk = U[:, :k] @ np.diag(s[:k]) @ Vt[:k, :]
    rel_error = np.linalg.norm(Z.to_numpy()-Zk, "fro") / np.linalg.norm(Z.to_numpy(), "fro")
    svd_rows.append((k, 1-rel_error, rel_error))
display(pd.DataFrame(svd_rows, columns=["rank", "reconstruction_score", "relative_error"]).round(4))

plt.figure(figsize=(7, 3.5))
plt.plot(range(1, 7), s, marker="o")
plt.title("Singular values of standardized sensor data")
plt.xlabel("Component")
plt.ylabel("Singular value")
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
rank reconstruction_score relative_error
0 1 0.5321 0.4679
1 2 0.6722 0.3278
2 3 0.7799 0.2201
3 4 0.8562 0.1438
4 5 0.8997 0.1003
5 6 1.0000 0.0000

png

結果の読み取り

上位少数の特異値が大きければ、少数成分で主要変動を保持できます。圧縮率だけでなく、故障直前の小さな兆候まで消していないかを検証してください。再構成誤差は異常スコアにもできますが、正常期間だけで基準を作る運用が必要です。

No.048：PCAとの関係

実務での意味

主成分分析（PCA）は、相関した多数の監視値を、互いに無相関な少数のKPIへまとめます。監視画面の指標削減や、設備状態の二次元表示に利用できます。

分析・モデル化の考え方

標準化データ ZZ の相関行列 C=ZTZ/(n−1)C=Z^TZ/(n-1) を固有分解する方法と、ZZ をSVDする方法は同じ主成分方向を与えます。寄与率は λj/∑iλi\lambda_j/\sum_i\lambda_i です。PCAは分散を説明する手法であり、品質不良を最もよく予測する方向とは限りません。

Pythonで確認する

scores = Z.to_numpy() @ eigvecs[:, :2]
loadings_df = pd.DataFrame(eigvecs[:, :2], index=sensor_names, columns=["PC1", "PC2"])
display(loadings_df.round(3))
print("Cumulative variance ratio (PC1+PC2):", round((eigvals[:2].sum()/eigvals.sum()), 4))

plt.figure(figsize=(7, 4))
plt.scatter(scores[:-20, 0], scores[:-20, 1], alpha=0.65, label="Earlier")
plt.scatter(scores[-20:, 0], scores[-20:, 1], color="crimson", label="Latest 20h")
plt.title("PCA score map for equipment monitoring")
plt.xlabel("PC1 score")
plt.ylabel("PC2 score")
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()
PC1 PC2
Drive-X -0.403 0.558
Drive-Y -0.394 0.600
Bearing-A -0.438 -0.185
Bearing-B -0.425 -0.302
Frame -0.409 -0.336
Outlet -0.378 -0.301 
Cumulative variance ratio (PC1+PC2): 0.8925

png

結果の読み取り

負荷量表は各主成分をどのセンサーが構成するかを示します。最新20時間の点群が過去からずれる場合、単一センサー値ではなく設備状態全体が変化した可能性があります。管理限界は同じデータに後付けせず、正常期間と検証期間を分けて設定します。

No.049：PageRankとの関係

実務での意味

工程の重要度を単純な接続数だけで測ると、重要工程から影響を受ける工程を過小評価します。PageRankは「重要な工程から参照・流入される工程も重要」という再帰的評価を行います。

分析・モデル化の考え方

行確率行列 PP を用い、定常ベクトル rr を

r=αPTr+(1−α)vr=\alpha P^Tr+(1-\alpha)v

で定義します。α\alpha はネットワーク伝播を続ける確率、vv は一様な再出発分布です。これはGoogle行列の固有値1に対応する固有ベクトル問題です。ここでの順位は工程価値ではなく、定義したリンク上の構造的重要度です。

Pythonで確認する

P = A.copy()
row_sum = P.sum(axis=1)
for i in range(len(P)):
    P[i] = P[i] / row_sum[i] if row_sum[i] > 0 else np.ones(len(P)) / len(P)
alpha = 0.85
G = alpha * P + (1-alpha) * np.ones_like(P) / len(P)
r = np.ones(len(P)) / len(P)
pr_history = []
for k in range(100):
    r_new = G.T @ r
    pr_history.append(np.linalg.norm(r_new-r, 1))
    if pr_history[-1] < 1e-12: break
    r = r_new
pagerank_df = pd.DataFrame({"process": processes, "PageRank": r_new}).sort_values("PageRank", ascending=False)
display(pagerank_df.round(4))

plt.figure(figsize=(8, 3.8))
plt.bar(pagerank_df["process"], pagerank_df["PageRank"], color="#4f8f5b")
plt.title("Process influence based on PageRank")
plt.xlabel("Process")
plt.ylabel("PageRank score")
plt.xticks(rotation=35, ha="right")
plt.grid(True, axis="y", alpha=0.3)
plt.tight_layout()
plt.show()
process PageRank
5 Inspection 0.1682
3 Grinding 0.1609
7 Shipping 0.1585
4 Wash 0.1560
2 HeatTreat 0.1130
1 Machining 0.1111
6 Rework 0.0967
0 Material 0.0356

png

結果の読み取り

再加工ループへ接続する工程が上位なら、停止時の直接損失だけでなく、循環を介した波及を重点評価する候補です。順位はリンク重量と減衰率に依存します。改善投資ではPageRankに、処理量、停止時間、代替可否、品質コストを掛け合わせます。

No.050：グラフ解析への応用

実務での意味

設備・工程をノード、関係をエッジとすることで、工場全体をグラフとして扱えます。スペクトル解析は、密に結び付く工程群、分断しやすい境界、監視単位の候補を見つけます。

分析・モデル化の考え方

方向をいったん無視した重み行列 WW から次数行列 DD とラプラシアン L=D−WL=D-W を作ります。正規化ラプラシアン

Lsym=I−D−1/2WD−1/2L_{\mathrm{sym}}=I-D^{-1/2}WD^{-1/2}

の第2最小固有値（代数的連結度）が小さいほど分断されやすく、対応するFiedlerベクトルの符号で2群への分割候補を得られます。

Pythonで確認する

W = (A + A.T) / 2
d = W.sum(axis=1)
D_inv_sqrt = np.diag(1 / np.sqrt(d))
Lsym = np.eye(len(W)) - D_inv_sqrt @ W @ D_inv_sqrt
lap_vals, lap_vecs = np.linalg.eigh(Lsym)
fiedler = lap_vecs[:, 1]
cluster = np.where(fiedler >= 0, "Group A", "Group B")
graph_df = pd.DataFrame({"process": processes, "Fiedler_value": fiedler, "candidate_group": cluster})
display(graph_df.sort_values("Fiedler_value").round(4))
print("normalized algebraic connectivity:", round(lap_vals[1], 6))

plt.figure(figsize=(8, 3.8))
colors = np.where(fiedler >= 0, "#3568a8", "#d1792f")
plt.bar(processes, fiedler, color=colors)
plt.axhline(0, color="black", linewidth=0.8)
plt.title("Spectral partition candidate by Fiedler vector")
plt.xlabel("Process")
plt.ylabel("Fiedler vector value")
plt.xticks(rotation=35, ha="right")
plt.grid(True, axis="y", alpha=0.3)
plt.tight_layout()
plt.show()
process Fiedler_value candidate_group
5 Inspection -0.5414 Group B
7 Shipping -0.4346 Group B
4 Wash -0.3317 Group B
3 Grinding 0.0197 Group A
6 Rework 0.1728 Group A
2 HeatTreat 0.2314 Group A
0 Material 0.3464 Group A
1 Machining 0.4519 Group A 
normalized algebraic connectivity: 0.183177

png

結果の読み取り

符号で分かれる2群は、保全担当範囲、データ収集単位、工程間バッファの検討候補です。ゼロ付近の工程や群をつなぐエッジは境界候補ですが、スペクトル分割は物理レイアウトや安全規則を知りません。現場導線と照合して採否を決めます。

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

  1. 全部を計算する必要はない：最大モードだけならべき乗法、対称大規模ならランチョス法、非対称ならアーノルディ法と、意思決定に必要な固有対へ計算を絞れます。
  2. 同じ分解が複数業務につながる：固有分解とSVDは、振動モード、PCA、圧縮、異常監視の共通基盤です。
  3. ネットワークの定義が結論を決める：PageRankやラプラシアンの順位・分割は、何をエッジとし、どの重量を与えるかで変わります。
  4. 数値結果は原因ではなく候補：大きな固有値や中心性は調査の優先順位を示しますが、故障や品質不良の因果を単独で証明しません。

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

  • センサー校正、欠測、時刻同期、設備停止区間を含むデータ品質基準
  • 正常期間、異常期間、負荷条件を分けた検証データ
  • 行列の定義、標準化、シフト、収束許容値、乱数シードの記録
  • 残差、再構成誤差、計算時間、誤報率を含む受入基準
  • 工程責任者・保全・品質保証が結果を確認するレビュー手順
  • アラームから点検、原因確認、モデル更新までの業務フロー

まとめ

No.041〜No.050では、固有値問題の基本から、べき乗法・逆反復法・QR法、Krylov部分空間法、SVD・PCA、PageRank・グラフラプラシアンまでを、製造業の一貫した架空例で確認しました。重要なのは高度な手法名ではなく、必要な情報だけを安定して計算し、現場の点検・投資・監視設計へ翻訳することです。

法人向けのご相談

数理工房では、設備データ解析、異常検知、工程ネットワーク分析、数値計算の高速化、現場担当者向け研修まで、課題定義から実装・運用設計を支援します。手元データで何が判断できるか、検証設計からご一緒できます。

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

Related Posts

関連する記事