100本ノック / 線形代数 / 線形代数100本ノック

製造業の固有値分析入門|工程変動の増幅・収束をPythonで見抜く

工程変動の「増幅・収束」を見抜く固有値分析 — 製造業の線形代数100本ノック No.041〜No.050

製造ラインでは、前工程の寸法偏差、設備振動、温度ずれなどが互いに影響し、翌日の品質状態へ伝播します。本稿では、3工程からなる架空ラインを題材に、固有値・固有ベクトル・対角化・重複度・多項式・ジョルダン標準形・スペクトル分解・正定値行列を、設備調整と品質判断につながる形で確認します。

到達点は、数式を計算することではなく、(1) どの変動モードが残るか、(2) 影響が長期的に収束するか、(3) 評価指標が安全に使えるか、を説明できることです。掲載データは Python で生成した架空データで、外部データには依存しません。

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

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

ある精密部品ラインで、加工・熱処理・検査の標準化偏差を日次で記録しているとします。品質偏差ベクトルを xt\mathbf{x}_t、工程間の伝播を表す行列を AA とすると、単純化した状態更新は

xt+1=Axt+εt\mathbf{x}_{t+1}=A\mathbf{x}_t+\boldsymbol{\varepsilon}_t

です。εt\boldsymbol{\varepsilon}_t はその日に加わる小さな外乱です。知りたいのは個別係数だけではありません。偏差が何倍で残るのか、どの工程の組合せとして現れるのか、調整後に収束するのかを判断する必要があります。

現場でよくある状況

  • 各工程の管理図は正常でも、工程をまたぐ緩やかな変動が残る
  • 調整対象を一工程に絞った結果、別工程の偏差が増える
  • シミュレーションはできるが、長期安定性を説明できない
  • 品質スコアの重み行列が妥当か確認せず、最適化を実施している

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

行列 AA の各要素は局所的な影響です。一方、現場が必要とする「何日後にどうなるか」は AkA^k で決まり、要素を眺めるだけでは読み取れません。固有値は変動モードの倍率、固有ベクトルはその工程構成を与えます。ただし重複固有値や対角化不能な行列では、単純な解釈が破綻するため、重複度や最小多項式まで確認する必要があります。

今回扱うノックの全体像

No.テーマ実務上の問い
041固有値変動は増幅するか、減衰するか
042固有ベクトルどの工程の組合せが支配的か
043対角化多日先の伝播を簡潔に計算できるか
044幾何学的重複度独立な変動方向はいくつあるか
045代数的重複度同じ倍率が式の上で何回現れるか
046特性多項式固有値を定義式から確認できるか
047最小多項式行列を消す最小の関係式は何か
048ジョルダン標準形対角化不能時の一時的増幅をどう捉えるか
049スペクトル分解変動を直交モード別に分解できるか
050正定値行列品質スコアが常に非負で一意な評価になるか

Python 環境の準備

NumPy で行列計算、pandas で表、Matplotlib で可視化します。再現性確保のため乱数生成器の seed を固定します。表示上の丸めと計算精度は分け、判定は丸め前の値で行います。

import sys
import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
from IPython.display import display

np.set_printoptions(precision=4, suppress=True)
rng = np.random.default_rng(20260712)
print(f"Python {sys.version.split()[0]} / NumPy {np.__version__} / pandas {pd.__version__} / Matplotlib {matplotlib.__version__}")
Python 3.13.1 / NumPy 2.5.1 / pandas 3.0.3 / Matplotlib 3.11.0

架空データの作成

3工程の標準化偏差を120日分生成します。伝播行列 AA は対称行列とし、実務上ありがちな「隣接工程ほど影響が強い」構造を設定します。最大固有値が1未満なので、外乱がなければ偏差は長期的に減衰します。ここでは、この結論を先取りせず各ノックで検証します。

processes = ["Machining", "HeatTreat", "Inspection"]
A = np.array([[0.72, 0.12, 0.04],
              [0.12, 0.58, 0.10],
              [0.04, 0.10, 0.46]])
n_days = 120
states = np.zeros((n_days, 3))
states[0] = [1.8, -0.8, 0.5]
for t in range(n_days - 1):
    states[t + 1] = A @ states[t] + rng.normal(0, 0.12, 3)
quality_df = pd.DataFrame(states, columns=processes).assign(Day=np.arange(1, n_days + 1))
display(pd.DataFrame(A, index=processes, columns=processes).round(2))
display(quality_df.head().round(3))

fig, ax = plt.subplots(figsize=(9, 4))
for col in processes:
    ax.plot(quality_df["Day"], quality_df[col], label=col, linewidth=1.4)
ax.set_title("Simulated standardized quality deviations")
ax.set_xlabel("Day"); ax.set_ylabel("Standardized deviation")
ax.grid(True, alpha=0.3); ax.legend(); fig.tight_layout(); plt.show()
Machining HeatTreat Inspection
Machining 0.72 0.12 0.04
HeatTreat 0.12 0.58 0.10
Inspection 0.04 0.10 0.46
Machining HeatTreat Inspection Day
0 1.800 -0.800 0.500 1
1 1.309 -0.116 0.272 2
2 0.988 0.167 0.221 3
3 0.713 0.299 0.347 4
4 0.657 0.304 0.234 5

png

No.041:固有値とは何か

実務での意味

固有値は、工程間の影響を繰り返したときに、特定の変動モードが1期間で何倍になるかを表します。絶対値が1を超える固有値があれば、そのモードは外乱がなくても増幅するため、調整ルールの見直しが必要です。

分析・モデル化の考え方

ゼロでないベクトル v\mathbf{v} に対し、行列を掛けても向きが変わらず倍率だけが λ\lambda になるとき、λ\lambda を固有値と呼びます。安定性の一次判定にはスペクトル半径 ρ(A)=maxiλi\rho(A)=\max_i|\lambda_i| を使います。

Av=λv,ρ(A)<1Akx00A\mathbf{v}=\lambda\mathbf{v},\qquad \rho(A)<1\Rightarrow A^k\mathbf{x}_0\to\mathbf{0}

Pythonで確認する

eigvals = np.linalg.eigvalsh(A)
eigenvalue_table = pd.DataFrame({"Eigenvalue": eigvals[::-1], "Magnitude": np.abs(eigvals[::-1])})
display(eigenvalue_table.round(4))
rho = np.max(np.abs(eigvals))
print(f"Spectral radius = {rho:.4f}; asymptotically stable = {rho < 1}")

fig, ax = plt.subplots(figsize=(7, 3.5))
ax.bar(range(1, 4), eigvals[::-1], color="#2878B5")
ax.axhline(1, color="crimson", linestyle="--", label="stability boundary")
ax.set_title("Eigenvalues of the process propagation matrix")
ax.set_xlabel("Mode rank"); ax.set_ylabel("Eigenvalue")
ax.set_xticks([1, 2, 3]); ax.grid(True, axis="y", alpha=0.3); ax.legend(); fig.tight_layout(); plt.show()
Eigenvalue Magnitude
0 0.8102 0.8102
1 0.5487 0.5487
2 0.4010 0.4010
Spectral radius = 0.8102; asymptotically stable = True


png

結果の読み取り

最大固有値(スペクトル半径)は約0.82で1未満です。したがって設定した伝播だけなら全モードが減衰します。ただし0.82のモードは最も減衰が遅く、日々の外乱が重なると残存しやすいため、最優先の監視対象です。

No.042:固有ベクトルとは何か

実務での意味

固有ベクトルは、同じ倍率で一体として動く工程の組合せです。最大固有値に対応するベクトルを見れば、長く残る品質偏差がどの工程を同時に含むかを特定できます。

分析・モデル化の考え方

固有ベクトルの符号全体は反転しても同じ方向を表します。そのため正負そのものより、成分の絶対値と工程間の相対的な符号を読みます。対称行列では固有ベクトルを互いに直交する長さ1のベクトルとして選べます。

AV=VΛ,VTV=IA V=V\Lambda,\qquad V^\mathsf{T}V=I

Pythonで確認する

eigvals, eigvecs = np.linalg.eigh(A)
order = np.argsort(eigvals)[::-1]
eigvals, eigvecs = eigvals[order], eigvecs[:, order]
# 読みやすさのため、各ベクトルで最大絶対値の成分を正にそろえる
eigvecs *= np.sign(eigvecs[np.argmax(np.abs(eigvecs), axis=0), range(3)])
loading_df = pd.DataFrame(eigvecs, index=processes, columns=[f"Mode {i+1}" for i in range(3)])
display(loading_df.round(3))

fig, ax = plt.subplots(figsize=(7, 3.5))
ax.bar(processes, np.abs(eigvecs[:, 0]), color="#F28E2B")
ax.set_title("Absolute loadings of the slowest-decaying mode")
ax.set_xlabel("Process"); ax.set_ylabel("Absolute eigenvector component")
ax.grid(True, axis="y", alpha=0.3); fig.tight_layout(); plt.show()
Mode 1 Mode 2 Mode 3
Machining 0.812 -0.575 0.095
HeatTreat 0.529 0.660 -0.534
Inspection 0.244 0.484 0.840

png

結果の読み取り

最も減衰の遅い Mode 1 は加工と熱処理の寄与が大きく、検査も同方向に含みます。単一工程の異常として扱うより、加工条件と熱履歴をセットで確認する横断的な原因調査が適切です。固有ベクトルは因果関係そのものではなく、調査の優先順位を示す仮説です。

No.043:対角化

実務での意味

対角化により、多日先の影響 AkA^k を、独立なモードごとの倍率 λik\lambda_i^k として計算できます。長期シミュレーションの説明性と計算効率が上がります。

分析・モデル化の考え方

独立な固有ベクトルを列に持つ行列 VV が正則なら A=VΛV1A=V\Lambda V^{-1} と分解できます。対称な AA では V1=VTV^{-1}=V^\mathsf{T} です。初期偏差をモード座標へ移し、倍率を掛け、工程座標へ戻します。

Ak=VΛkV1,xk=VΛkV1x0A^k=V\Lambda^kV^{-1},\qquad \mathbf{x}_k=V\Lambda^kV^{-1}\mathbf{x}_0

Pythonで確認する

k = 12
Lambda = np.diag(eigvals)
A_k_direct = np.linalg.matrix_power(A, k)
A_k_diag = eigvecs @ np.linalg.matrix_power(Lambda, k) @ eigvecs.T
x0 = states[0]
comparison = pd.DataFrame({
    "Direct": A_k_direct @ x0,
    "Diagonalized": A_k_diag @ x0,
    "Absolute error": np.abs((A_k_direct - A_k_diag) @ x0)
}, index=processes)
display(comparison.round(10))
print(f"matrix reconstruction error = {np.linalg.norm(A - eigvecs @ Lambda @ eigvecs.T):.2e}")
Direct Diagonalized Absolute error
Machining 0.076042 0.076042 0.0
HeatTreat 0.048527 0.048527 0.0
Inspection 0.022203 0.022203 0.0
matrix reconstruction error = 4.75e-16

結果の読み取り

12日後の偏差は直接計算と対角化で一致し、誤差は浮動小数点の丸め程度です。最大固有値の12乗は小さくなるため、外乱なしの偏差は大幅に減衰します。実務では AA が期間中一定という仮定が妥当か、設備改造や品種切替ごとに確認します。

No.044:幾何学的重複度

実務での意味

同じ固有値に対して独立な固有ベクトルが何本あるかは、同じ減衰率を持つ独立した変動パターンの数を意味します。不足すると対角化できず、予測挙動に多項式的な項が混ざります。

分析・モデル化の考え方

固有値 λ\lambda の幾何学的重複度は、AλIA-\lambda I の零空間の次元です。数値計算では特異値に許容誤差を設け、ランクを判定します。ここでは説明用に重複固有値2を持つ対角行列を使います。

gλ=dimker(AλI)=nrank(AλI)g_\lambda=\dim\ker(A-\lambda I)=n-\operatorname{rank}(A-\lambda I)

Pythonで確認する

B = np.diag([2.0, 2.0, 0.5])
lam = 2.0
rank = np.linalg.matrix_rank(B - lam * np.eye(3))
geom_mult = B.shape[0] - rank
print("Example matrix B:\n", B)
print(f"rank(B - {lam}I) = {rank}")
print(f"Geometric multiplicity of lambda={lam} = {geom_mult}")
print("Independent directions: e1 and e2")
Example matrix B:
 [[2.  0.  0. ]
 [0.  2.  0. ]
 [0.  0.  0.5]]
rank(B - 2.0I) = 1
Geometric multiplicity of lambda=2.0 = 2
Independent directions: e1 and e2

結果の読み取り

固有値2には独立方向が2本あり、同じ増幅率を持つ二つの独立現象を区別できます。たとえば二台の独立設備が同じ減衰率を持つ状況に対応します。数値データでは「ほぼ同じ固有値」を重複と断定せず、測定誤差と許容差を明示すべきです。

No.045:代数的重複度

実務での意味

代数的重複度は、特性多項式の根として同じ固有値が何回現れるかです。幾何学的重複度と比較することで、行列が対角化可能かを判断できます。

分析・モデル化の考え方

常に 1gλaλ1\le g_\lambda\le a_\lambda が成立します。すべての固有値で幾何学的重複度と代数的重複度が一致し、その合計が次元に達すれば対角化可能です。ここでは、重複度2なのに独立固有ベクトルが1本しかない例を確認します。

aλ: det(μIA) における根 λ の重複回数a_\lambda:\ \det(\mu I-A)\text{ における根 }\lambda\text{ の重複回数}

Pythonで確認する

C = np.array([[0.8, 1.0, 0.0],
              [0.0, 0.8, 0.0],
              [0.0, 0.0, 0.3]])
roots = np.roots(np.poly(C))
alg_mult = int(np.sum(np.isclose(roots, 0.8)))
geom_mult = C.shape[0] - np.linalg.matrix_rank(C - 0.8 * np.eye(3))
display(pd.DataFrame({"Measure": ["Algebraic multiplicity", "Geometric multiplicity"],
                      "Value": [alg_mult, geom_mult]}))
print("Eigenvalues:", roots)
print("Diagonalizable:", alg_mult == geom_mult)
Measure Value
0 Algebraic multiplicity 2
1 Geometric multiplicity 1
Eigenvalues: [0.8 0.8 0.3]
Diagonalizable: False

結果の読み取り

固有値0.8の代数的重複度は2、幾何学的重複度は1です。独立固有ベクトルが不足するため、この行列は対角化できません。固有値がすべて1未満でも、非対角成分の影響で一時的な増幅が起こり得る点が設備制御上の注意点です。

No.046:特性多項式

実務での意味

特性多項式は固有値を決める方程式で、工程モデルの安定性を理論的に確認する入口です。行列のトレースや行列式とも結びつき、実装結果の整合性チェックにも使えます。

分析・モデル化の考え方

NumPy の polydet(λIA)\det(\lambda I-A) の係数を返します。根が固有値であり、根の和はトレース、積は行列式に一致します。3次なら p(λ)=λ3+c2λ2+c1λ+c0p(\lambda)=\lambda^3+c_2\lambda^2+c_1\lambda+c_0 です。

pA(λ)=det(λIA)=i=1n(λλi)p_A(\lambda)=\det(\lambda I-A)=\prod_{i=1}^n(\lambda-\lambda_i)

Pythonで確認する

coeff = np.poly(A)
roots = np.roots(coeff)
poly_df = pd.DataFrame({"Power": [3, 2, 1, 0], "Coefficient": coeff})
display(poly_df.round(6))
print("Roots:", np.sort(roots)[::-1])
print(f"sum(roots)={roots.sum():.6f}, trace(A)={np.trace(A):.6f}")
print(f"product(roots)={roots.prod():.6f}, det(A)={np.linalg.det(A):.6f}")
Power Coefficient
0 3 1.000000
1 2 -1.760000
2 1 0.989600
3 0 -0.178304
Roots: [0.8102 0.5487 0.401 ]
sum(roots)=1.760000, trace(A)=1.760000
product(roots)=0.178304, det(A)=0.178304

結果の読み取り

特性多項式の根は No.041 の固有値と一致し、根の和と積もトレース・行列式に一致します。このような複数の恒等関係を検算に使うと、係数の転記ミスや行列の向きの取り違えを早期に発見できます。高次行列では多項式の根を直接求める方法は数値的に不安定なため、実務計算は専用の固有値アルゴリズムを使います。

No.047:最小多項式

実務での意味

最小多項式は、行列が満たす最短の多項式関係です。状態更新の高いべき乗を低次数に置き換え、モデルの本質的な動的複雑さを把握できます。

分析・モデル化の考え方

異なる3固有値を持つ今回の対称行列では、最小多項式は各固有値を1回ずつ根に持ち、特性多項式と一致します。Cayley–Hamilton の定理により特性多項式を行列へ代入するとゼロ行列になります。重複固有値がある場合、最小多項式の次数はより低くなることがあります。

mA(A)=0,mA はこの性質を持つモニック多項式のうち最小次数m_A(A)=0,\qquad m_A\text{ はこの性質を持つモニック多項式のうち最小次数}

Pythonで確認する

# p(A) = A^3 + c2 A^2 + c1 A + c0 I を評価
pA = (np.linalg.matrix_power(A, 3) + coeff[1] * np.linalg.matrix_power(A, 2)
      + coeff[2] * A + coeff[3] * np.eye(3))
print("p(A) =\n", pA)
print(f"Frobenius norm of p(A) = {np.linalg.norm(pA):.2e}")

B_min = (B - 2 * np.eye(3)) @ (B - 0.5 * np.eye(3))
print(f"For B, ||(B-2I)(B-0.5I)|| = {np.linalg.norm(B_min):.2e}")
p(A) =
 [[ 0. -0. -0.]
 [-0. -0. -0.]
 [-0. -0. -0.]]
Frobenius norm of p(A) = 1.55e-16
For B, ||(B-2I)(B-0.5I)|| = 0.00e+00

結果の読み取り

今回の p(A)p(A) は数値誤差の範囲でゼロです。また、No.044 の行列 BB は3次元でも異なる固有値が2種類なので、2次の最小多項式で消えます。最小多項式は単なる次元ではなく、独立な動的モードとジョルダンブロックの大きさを反映します。

No.048:ジョルダン標準形

実務での意味

対角化できないモデルでも、ジョルダン標準形を使えば挙動を「固有値による倍率」と「連鎖による一時的増幅」に分けられます。制御モデルで固有値だけを見て安心しないための重要な視点です。

分析・モデル化の考え方

No.045 の 2×22\times2 ブロックは J=λI+NJ=\lambda I+NN2=0N^2=0 です。そのため Jk=λkI+kλk1NJ^k=\lambda^k I+k\lambda^{k-1}N となり、kk を含む項が現れます。λ<1|\lambda|<1 なら最終的には収束しますが、初期には偏差が増える場合があります。

Jk=(λI+N)k=λkI+kλk1N(N2=0)J^k=(\lambda I+N)^k=\lambda^k I+k\lambda^{k-1}N\quad(N^2=0)

Pythonで確認する

J = C[:2, :2]
x0_j = np.array([0.0, 1.0])
steps = np.arange(0, 26)
trajectory = np.array([np.linalg.matrix_power(J, int(k)) @ x0_j for k in steps])
peak_k = int(steps[np.argmax(np.abs(trajectory[:, 0]))])
print(f"First component peaks at step {peak_k}: {trajectory[peak_k, 0]:.4f}")

fig, ax = plt.subplots(figsize=(8, 3.8))
ax.plot(steps, trajectory[:, 0], marker="o", markersize=3, label="component 1")
ax.plot(steps, trajectory[:, 1], marker="s", markersize=3, label="component 2")
ax.set_title("Transient amplification in a Jordan block")
ax.set_xlabel("Step k"); ax.set_ylabel("State value")
ax.grid(True, alpha=0.3); ax.legend(); fig.tight_layout(); plt.show()
First component peaks at step 5: 2.0480


png

結果の読み取り

固有値は0.8で安定範囲内ですが、第1成分はいったん増えてから減衰します。これは第2成分から第1成分への連鎖が kλk1k\lambda^{k-1} として効くためです。設備立上げ時の許容上限は最終安定性だけでなく、この過渡ピークを含めて設計する必要があります。

No.049:スペクトル分解

実務での意味

対称行列のスペクトル分解は、工程間影響を直交したモードの足し合わせとして表します。各モードの寄与を個別に可視化でき、改善施策がどの変動構造を抑えるか説明しやすくなります。

分析・モデル化の考え方

対称行列では正規直交固有ベクトル vi\mathbf{v}_i を用い、A=iλiviviTA=\sum_i\lambda_i\mathbf{v}_i\mathbf{v}_i^\mathsf{T} と分解できます。viviT\mathbf{v}_i\mathbf{v}_i^\mathsf{T} はそのモードへの射影です。上位モードだけを残せば低ランク近似になります。

A=VΛVT=i=1nλiviviTA=V\Lambda V^\mathsf{T}=\sum_{i=1}^n\lambda_i\mathbf{v}_i\mathbf{v}_i^\mathsf{T}

Pythonで確認する

components = [eigvals[i] * np.outer(eigvecs[:, i], eigvecs[:, i]) for i in range(3)]
A_rank1 = components[0]
relative_error = np.linalg.norm(A - A_rank1, "fro") / np.linalg.norm(A, "fro")
display(pd.DataFrame(A_rank1, index=processes, columns=processes).round(3))
print(f"Rank-1 relative Frobenius error = {relative_error:.3f}")

fig, axes = plt.subplots(1, 3, figsize=(10, 3))
vmax = max(np.max(np.abs(c)) for c in components)
for i, (ax, comp) in enumerate(zip(axes, components), 1):
    im = ax.imshow(comp, cmap="coolwarm", vmin=-vmax, vmax=vmax)
    ax.set_title(f"Mode {i}"); ax.set_xlabel("Source process"); ax.set_ylabel("Target process")
    ax.set_xticks(range(3), ["M", "H", "I"]); ax.set_yticks(range(3), ["M", "H", "I"])
    ax.grid(False)
fig.colorbar(im, ax=axes, shrink=0.75, label="Contribution")
fig.suptitle("Spectral components of the propagation matrix")
fig.tight_layout(); plt.show()
Machining HeatTreat Inspection
Machining 0.535 0.349 0.161
HeatTreat 0.349 0.227 0.105
Inspection 0.161 0.105 0.048
Rank-1 relative Frobenius error = 0.643


/var/folders/3y/fmw40k0x78xblvb3gkcyvy1h0000gn/T/ipykernel_35217/1914797165.py:16: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
  fig.tight_layout(); plt.show()


png

結果の読み取り

Mode 1 は全工程に同方向の寄与を持ち、共通的な変動として残ります。ただしランク1近似の相対誤差も無視できないため、実運用の予測に上位1モードだけを使うのは粗すぎます。可視化・説明には上位モード、制御計算には必要精度を満たすモード数、という使い分けが有効です。

No.050:正定値行列

実務での意味

正定値行列は、複数の品質偏差を一つの損失スコアへまとめるとき、ゼロ以外の偏差に必ず正の罰則を与えます。また、二次最適化の解を一意にし、Cholesky 分解など安定した計算を可能にします。

分析・モデル化の考え方

対称行列 QQ がすべての非ゼロ x\mathbf{x} に対して xTQx>0\mathbf{x}^\mathsf{T}Q\mathbf{x}>0 を満たすとき正定値です。対称行列では、全固有値が正であることと同値です。ここでは品質偏差の相関を考慮した損失行列を構成します。

L(x)=xTQx,Q0λmin(Q)>0L(\mathbf{x})=\mathbf{x}^\mathsf{T}Q\mathbf{x},\qquad Q\succ0\Longleftrightarrow\lambda_{\min}(Q)>0

Pythonで確認する

Q = np.array([[2.0, 0.35, 0.10],
              [0.35, 1.5, 0.25],
              [0.10, 0.25, 1.0]])
q_eigs = np.linalg.eigvalsh(Q)
loss = np.einsum("ij,jk,ik->i", states, Q, states)
print("Eigenvalues of Q:", q_eigs)
print("Positive definite:", np.all(q_eigs > 0))
print(f"Minimum observed loss = {loss.min():.6f}")
display(pd.DataFrame({"Day": quality_df["Day"].head(8), "Quality loss": loss[:8]}).round(4))

fig, ax = plt.subplots(figsize=(8, 3.6))
ax.plot(quality_df["Day"], loss, color="#59A14F")
ax.set_title("Quadratic quality loss over time")
ax.set_xlabel("Day"); ax.set_ylabel("Quality loss x^T Q x")
ax.grid(True, alpha=0.3); fig.tight_layout(); plt.show()
Eigenvalues of Q: [0.8946 1.3903 2.2151]
Positive definite: True
Minimum observed loss = 0.002225
Day Quality loss
0 1 6.6620
1 2 3.4683
2 3 2.2214
3 4 1.5221
4 5 1.2650
5 6 0.6545
6 7 0.4485
7 8 0.5211

png

結果の読み取り

QQ の固有値はすべて正で、観測された損失も非負です。したがって偏差ゼロが唯一の最小点となり、改善方向を曖昧にしません。ただし QQ の重みは数学だけでは決められません。不良コスト、安全影響、顧客仕様を反映し、部門間で合意したうえで版管理する必要があります。

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

今回の行列では最大固有値が1未満で、外乱がなければ品質偏差は収束します。しかし、支配モードは加工・熱処理を中心に複数工程へ広がるため、個別工程だけの対策では再発しやすいと読めます。また、ジョルダンブロックの例が示すように、固有値が安定範囲でも一時的な増幅はあり得ます。

実務上は次の順序が有効です。

  1. スペクトル半径で長期的な収束性を確認する
  2. 支配固有ベクトルで横断的な調査対象を絞る
  3. 重複度と対角化可能性を確認し、過渡応答を評価する
  4. 対称行列ならスペクトル分解でモード別の寄与を説明する
  5. 品質損失の重み行列が正定値か確認し、最適化の一意性を担保する

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

  • データ定義:測定単位、サンプリング時刻、欠測・外れ値処理、品種切替を統一する
  • モデル同定AA を固定値として置くのではなく、履歴データから推定し、ホールドアウト期間で検証する
  • 不確実性評価:固有値・固有ベクトルの点推定だけでなく、ブートストラップ等で変動幅を確認する
  • 非定常性への対応:設備保全、金型交換、季節、材料ロットで行列が変わらないか監視する
  • 因果検証:固有モードは相関構造の要約であり、対策前には現場知識、実験計画、変更履歴と照合する
  • 運用設計:再推定頻度、警報閾値、承認者、ロールバック手順、モデル版を定める

まとめ

No.041〜No.050 では、行列の反復作用を「倍率」と「方向」に分ける固有値分析から、対角化不能な場合のジョルダン標準形、二次評価を支える正定値行列までを一続きで確認しました。固有値は安定性、固有ベクトルは変動の工程構成、重複度と最小多項式は動的構造、スペクトル分解は説明可能性、正定値性は評価・最適化の妥当性に対応します。

実データへ適用するときは、計算結果を「原因」と即断せず、工程知識と検証実験を組み合わせることが重要です。

法人向けのご相談

数理工房では、製造データを用いた品質変動の構造分析、設備状態モデル、異常検知、数理最適化、現場向け研修の設計をご支援しています。「データはあるが工程横断で説明できない」「モデルを現場判断へつなげたい」といった段階からご相談いただけます。

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