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

製造業の行列分解入門｜SVD・PCA・NMF・行列補完をPythonで実践

2020-07-10

欠損のある多品種品質データを意思決定へつなぐ行列分解 — 製造業の線形代数100本ノック No.061〜No.070

製造現場では、設備条件、検査特性、品種、時間帯が組み合わさり、品質データはすぐに大きな表になります。本稿では、架空の精密部品工場を題材に、連立方程式を安定して解く分解、品質変動を少数の要因へ圧縮する分解、欠測した検査値を補う分解を、現場判断につながる形で確認します。

到達点は手法名を覚えることではありません。(1) 何を速く・安定に計算したいのか、(2) どこまで情報を圧縮できるのか、(3) 欠損値をどの条件で補ってよいのかを説明し、目的に合う分解を選べることです。データは Python で生成した架空データで、外部データには依存しません。

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

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

12品種を8台の設備で加工し、寸法、真円度、粗さ、硬度、振動、温度の6特性を検査する工場を考えます。全組合せを毎回測定できるとは限らず、測定行列には欠損が生じます。また、特性同士は加工負荷や熱影響という共通要因を通じて相関します。

この状況で必要なのは、単に平均値を出すことではなく、校正計算を安定に解き、変動の主要因を抽出し、圧縮誤差を管理しながら未測定値を推定することです。行列分解は、大きな計算を意味のある小さな計算へ分ける共通言語になります。

現場でよくある状況

同じ係数行列で、日ごとに異なる測定値を何度も解いている
検査特性が増え、どの指標を会議で優先すべきか分からない
多品種データを可視化したいが、6次元のままでは比較できない
検査省略や通信断で欠損が生じ、単純な平均補完では品種差が消える
圧縮モデルの精度は高いが、現場が潜在因子を説明できない

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

行列分解には複数の選択肢があり、前提と目的が異なります。Cholesky 分解は対称正定値行列に強く、QR 分解は最小二乗に安定、SVD はランクと情報量の診断に強い一方で計算量が大きくなります。PCA は中心化した変動、NMF は非負の加法的パターンを説明します。欠損補完では、低ランク性が存在しても欠損が特定品種へ偏れば推定は危険です。

今回扱うノックの全体像

No.	テーマ	製造業での判断
061	LU分解	同じ設備モデルで複数の校正計算を効率化できるか
062	Cholesky分解	共分散・正規方程式を小さな計算量で解けるか
063	QR分解	多重共線性を含む校正回帰を安定に解けるか
064	SVD	品質行列の独立な変動方向と有効ランクは何か
065	低ランク近似	少数因子でどこまで情報を保存できるか
066	PCA	品種を主要な品質変動軸で比較できるか
067	NMF	非負の不良寄与を加法的な原因パターンへ分けられるか
068	Tensor分解	品種×設備×特性の三方向構造を保って圧縮できるか
069	CUR分解	実在する代表品種・代表特性で近似を説明できるか
070	行列補完	未測定値を低ランク構造からどこまで補えるか

Python 環境の準備

NumPy で行列計算、pandas で表、Matplotlib で可視化します。SciPy や機械学習ライブラリに依存せず、分解の中身が追える短い実装にします。乱数生成器の 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

架空データの作成

12品種、8設備、6検査特性のテンソルを作ります。背後には「加工負荷」「熱影響」「表面状態」という3つの潜在因子があり、観測値には小さな測定ノイズが加わる設定です。分析用の品質行列 $X$ は設備方向を平均した品種×特性の表とし、各列を標準化します。

products = [f"P{i:02d}" for i in range(1, 13)]
machines = [f"M{i:02d}" for i in range(1, 9)]
features = ["Dimension", "Roundness", "Roughness", "Hardness", "Vibration", "Temperature"]

product_factor = rng.normal(size=(12, 3))
machine_factor = rng.normal(scale=0.45, size=(8, 3))
feature_loading = np.array([[0.90, 0.15, 0.10], [0.75, 0.10, 0.30],
                            [0.15, 0.05, 0.95], [0.05, 0.90, 0.10],
                            [0.55, 0.25, 0.45], [0.10, 0.95, 0.05]])
quality_tensor = np.einsum("pr,fr->pf", product_factor, feature_loading)[:, None, :] + np.einsum("mr,fr->mf", machine_factor, feature_loading)[None, :, :]
quality_tensor += rng.normal(scale=0.10, size=quality_tensor.shape)
raw_matrix = quality_tensor.mean(axis=1)
X = (raw_matrix - raw_matrix.mean(axis=0)) / raw_matrix.std(axis=0, ddof=1)
quality_df = pd.DataFrame(X, index=products, columns=features)
display(quality_df.round(2))

fig, ax = plt.subplots(figsize=(9, 4))
im = ax.imshow(X, cmap="coolwarm", aspect="auto", vmin=-2.5, vmax=2.5)
ax.set_title("Standardized product-by-quality matrix")
ax.set_xlabel("Quality feature"); ax.set_ylabel("Product")
ax.set_xticks(range(len(features)), features, rotation=30, ha="right")
ax.set_yticks(range(len(products)), products); ax.grid(False)
fig.colorbar(im, ax=ax, label="Standardized value"); fig.tight_layout(); plt.show()

	Dimension	Roundness	Roughness	Hardness	Vibration	Temperature
P01	0.64	0.70	0.64	0.78	1.07	0.72
P02	0.17	0.29	0.57	0.42	0.47	0.39
P03	-0.54	0.09	2.07	0.67	0.87	0.53
P04	0.44	0.54	0.31	0.19	0.49	0.19
P05	-1.49	-1.41	0.02	0.11	-1.05	-0.03
P06	-0.03	-0.56	-2.27	-0.81	-1.34	-0.73
P07	1.08	1.21	0.13	-1.40	0.43	-1.33
P08	-0.18	-0.51	-0.72	1.52	-0.05	1.61
P09	1.70	1.62	-0.02	-0.55	0.93	-0.50
P10	-1.79	-1.75	-0.24	-1.45	-1.94	-1.53
P11	-0.53	-0.58	-0.20	-0.81	-0.63	-0.74
P12	0.54	0.36	-0.30	1.32	0.75	1.41

png

No.061：LU分解

実務での意味

校正係数や工程バランスを $mathbf{A}mathbf{x}=mathbf{b}$ として解くとき、設備構造を表す $A$ が同じで、日々の観測 $b$ だけが変わることがあります。LU 分解を一度行えば、以後は前進代入と後退代入で複数ケースを効率よく解けます。

分析・モデル化の考え方

行列を下三角行列 $L$ と上三角行列 $U$ に $A=LU$ と分け、 $Lmathbf{y}=mathbf{b}$ 、 $Umathbf{x}=mathbf{y}$ の順に解きます。実務の汎用実装では、ゼロまたは小さすぎるピボットによる不安定化を避けるため行交換 $PA=LU$ が必要です。ここでは原理を見やすくするため、行交換不要の優対角な例を使います。

Pythonで確認する

def lu_decompose(A):
    A = A.astype(float); n = len(A); L = np.eye(n); U = np.zeros_like(A)
    for i in range(n):
        U[i, i:] = A[i, i:] - L[i, :i] @ U[:i, i:]
        if abs(U[i, i]) < 1e-12: raise ValueError("Pivoting is required")
        L[i+1:, i] = (A[i+1:, i] - L[i+1:, :i] @ U[:i, i]) / U[i, i]
    return L, U

A_cal = np.array([[4., 1., 0.5], [1., 3., 0.4], [0.5, 0.4, 2.5]])
B_cal = np.array([[12., 10.], [8., 9.], [6., 7.]])  # two operating days
L, U = lu_decompose(A_cal)
Y = np.linalg.solve(L, B_cal); coef_lu = np.linalg.solve(U, Y)
display(pd.DataFrame(L).round(3)); display(pd.DataFrame(U).round(3))
display(pd.DataFrame(coef_lu, index=["Feed", "Speed", "Cooling"], columns=["Day A", "Day B"]).round(3))
print(f"reconstruction error = {np.linalg.norm(A_cal - L @ U):.2e}")
print(f"maximum residual = {np.max(np.abs(A_cal @ coef_lu - B_cal)):.2e}")

	0	1	2
0	1.000	0.0	0.0
1	0.250	1.0	0.0
2	0.125	0.1	1.0

	0	1	2
0	4.0	1.00	0.500
1	0.0	2.75	0.275
2	0.0	0.00	2.410

	Day A	Day B
Feed	2.379	1.697
Speed	1.652	2.152
Cooling	1.660	2.116

reconstruction error = 0.00e+00
maximum residual = 1.78e-15

結果の読み取り

再構成誤差と残差は丸め誤差の水準で、2日分の右辺を同じ分解から同時に解けました。係数行列を再利用する日次校正に向きます。ただし実運用では自作関数ではなく、ピボット選択を備えた検証済みライブラリを使い、条件数も監視します。

No.062：Cholesky分解

実務での意味

品質特性の共分散行列や、正則化した最小二乗の係数行列は対称正定値になりやすく、Cholesky 分解を使うと一般の分解より少ない計算量とメモリで解けます。

分析・モデル化の考え方

対称正定値行列 $A$ は $A=LL^{\mathsf T}$ と一意に分解できます。正定値とは任意のゼロでない $z$ に対して $z^{\mathsf T}Az>0$ となる性質です。対称性だけでは不十分で、最小固有値が正かを確認します。

Pythonで確認する

cov = np.cov(X, rowvar=False) + 0.05 * np.eye(X.shape[1])
chol = np.linalg.cholesky(cov)
target = np.ones(len(features))
weights = np.linalg.solve(chol.T, np.linalg.solve(chol, target))
weights /= target @ weights
chol_table = pd.DataFrame({"Feature": features, "Portfolio weight": weights})
display(chol_table.round(3))
print(f"minimum eigenvalue = {np.linalg.eigvalsh(cov).min():.4f}")
print(f"factorization error = {np.linalg.norm(cov - chol @ chol.T):.2e}")

	Feature	Portfolio weight
0	Dimension	0.791
1	Roundness	0.171
2	Roughness	0.653
3	Hardness	0.244
4	Vibration	-1.139
5	Temperature	0.280

minimum eigenvalue = 0.0507
factorization error = 3.85e-16

結果の読み取り

最小固有値が正で、分解誤差も十分小さいため、正定値という前提を満たします。ここでの重みは相関を考慮した合成品質指標の例ですが、負の重みが業務上許されるかは別の制約です。数値的に解けることと KPI として妥当であることを分けて判断します。

No.063：QR分解

実務での意味

温度、速度、負荷のように相関の強い説明変数から寸法偏差を推定すると、正規方程式は誤差を増幅しやすくなります。QR 分解は設計行列を直接分解し、校正回帰を比較的安定に解きます。

分析・モデル化の考え方

$D=QR$ で $Q^{\mathsf T}Q=I$ 、 $R$ は上三角です。最小二乗問題 $\min_\beta\|D\beta-y\|_2$ は $R\beta=Q^{\mathsf T}y$ に置き換わります。正規方程式 $D^{\mathsf T}D$ は条件数をほぼ二乗するため、QR の方が丸め誤差に強い選択です。

Pythonで確認する

n = 80
temperature = rng.normal(180, 4, n)
load = 0.85 * temperature + rng.normal(0, 2, n)
speed = rng.normal(1200, 60, n)
D = np.column_stack([np.ones(n), temperature - 180, load - load.mean(), (speed - 1200) / 100])
y = 0.20 + 0.035 * D[:, 1] - 0.020 * D[:, 2] + 0.08 * D[:, 3] + rng.normal(0, 0.08, n)
Q, R = np.linalg.qr(D, mode="reduced")
beta_qr = np.linalg.solve(R, Q.T @ y)
beta_lstsq = np.linalg.lstsq(D, y, rcond=None)[0]
display(pd.DataFrame({"Term": ["Intercept", "Temperature", "Load", "Speed/100"], "QR": beta_qr, "lstsq": beta_lstsq}).round(4))
print(f"orthogonality error = {np.linalg.norm(Q.T @ Q - np.eye(Q.shape[1])):.2e}")
print(f"condition(D)={np.linalg.cond(D):.1f}, condition(D.T@D)={np.linalg.cond(D.T @ D):.1f}")

	Term	QR	lstsq
0	Intercept	0.1934	0.1934
1	Temperature	0.0334	0.0334
2	Load	-0.0203	-0.0203
3	Speed/100	0.0829	0.0829

orthogonality error = 6.77e-16
condition(D)=8.9, condition(D.T@D)=78.4

結果の読み取り

QR と NumPy の最小二乗解は一致します。正規方程式の条件数が設計行列より大きいことから、相関の強い条件で正規方程式を明示的に作らない理由が分かります。ただし係数の因果解釈には実験計画や交絡管理が別途必要です。

No.064：SVD（特異値分解）

実務での意味

SVD は品質行列を「品種側のパターン」「重要度」「特性側のパターン」に分けます。特異値の減衰を見れば、6特性が実質何個の独立要因で動いているかを診断できます。

分析・モデル化の考え方

任意の $m\times n$ 行列は $X=U\Sigma V^{\mathsf T}$ と分解できます。特異値 $\sigma_i$ は非負で降順に並び、その二乗は各成分が持つエネルギーです。小さい特異値はノイズ方向や冗長性を示しますが、工程上重要な小変動を自動的に不要とみなしてはいけません。

Pythonで確認する

U, s, Vt = np.linalg.svd(X, full_matrices=False)
sv_df = pd.DataFrame({"Component": np.arange(1, len(s)+1), "Singular value": s, "Energy ratio": s**2 / np.sum(s**2)})
display(sv_df.round(4))
fig, ax = plt.subplots(figsize=(7, 3.5))
ax.plot(sv_df["Component"], sv_df["Singular value"], marker="o")
ax.set_title("Singular values of the quality matrix")
ax.set_xlabel("Component"); ax.set_ylabel("Singular value")
ax.grid(True, alpha=0.3); fig.tight_layout(); plt.show()
print(f"reconstruction error = {np.linalg.norm(X - U @ np.diag(s) @ Vt):.2e}")

	Component	Singular value	Energy ratio
0	1	6.0129	0.5478
1	2	4.3933	0.2924
2	3	3.2389	0.1589
3	4	0.1811	0.0005
4	5	0.1097	0.0002
5	6	0.0889	0.0001

png

reconstruction error = 1.07e-14

結果の読み取り

特異値は上位から急速に小さくなり、生成時に設定した少数の潜在因子が反映されています。SVD はランク診断の基準になります。保持成分数は折れ曲がりだけでなく、再構成誤差、異常検知感度、説明可能性を合わせて決めます。

No.065：低ランク近似

実務での意味

品質ダッシュボードやエッジ端末では、全データをそのまま保持せず、主要因だけで近似したいことがあります。低ランク近似は情報量と保存・計算コストのトレードオフを定量化します。

分析・モデル化の考え方

SVD の上位 $k$ 成分だけを残した $X_k=U_k\Sigma_kV_k^{\mathsf T}$ は、ランク $k$ の行列の中で Frobenius ノルム誤差を最小にします。相対誤差と累積寄与率を併記し、業務上許容できる $k$ を選びます。

Pythonで確認する

records = []
for k in range(1, len(s)+1):
    Xk = U[:, :k] @ np.diag(s[:k]) @ Vt[:k]
    records.append([k, np.linalg.norm(X-Xk)/np.linalg.norm(X), np.sum(s[:k]**2)/np.sum(s**2)])
rank_df = pd.DataFrame(records, columns=["Rank", "Relative error", "Cumulative energy"])
display(rank_df.round(4))
fig, ax = plt.subplots(figsize=(7, 3.5))
ax.plot(rank_df["Rank"], rank_df["Relative error"], marker="o", label="relative error")
ax.plot(rank_df["Rank"], 1-rank_df["Cumulative energy"], marker="s", label="unexplained energy")
ax.set_title("Accuracy versus retained rank")
ax.set_xlabel("Retained rank"); ax.set_ylabel("Ratio")
ax.grid(True, alpha=0.3); ax.legend(); fig.tight_layout(); plt.show()

	Rank	Relative error	Cumulative energy
0	1	0.6725	0.5478
1	2	0.3997	0.8403
2	3	0.0283	0.9992
3	4	0.0174	0.9997
4	5	0.0109	0.9999
5	6	0.0000	1.0000

png

結果の読み取り

ランクを上げるほど誤差は単調に減りますが、改善幅は次第に小さくなります。少数成分で大半の変動を保持できても、規格逸脱のような局所現象は残差側に現れる可能性があります。圧縮率だけでなく、残差管理図を併用する設計が必要です。

No.066：PCA（主成分分析）

実務での意味

PCA は相関する品質特性を少数の合成軸へまとめ、品種の似た変動や外れた品種を二次元で比較します。品質会議で「総合的にどの品種が異なるか」を可視化する入口になります。

分析・モデル化の考え方

中心化済み行列 $X$ の SVD から、スコア $Z=U\Sigma=XV$ 、負荷量 $V$ を得ます。標準化 PCA は単位の違いを除きますが、工程能力や許容差の重要度を同一にするとは限りません。符号は反転可能なので、相対的な向きで解釈します。

Pythonで確認する

scores = U[:, :2] * s[:2]
loadings = Vt[:2].T
display(pd.DataFrame(loadings, index=features, columns=["PC1", "PC2"]).round(3))
fig, ax = plt.subplots(figsize=(7, 5))
ax.scatter(scores[:, 0], scores[:, 1], color="#2878B5")
for i, name in enumerate(products): ax.annotate(name, scores[i], xytext=(4, 3), textcoords="offset points")
ax.axhline(0, color="gray", linewidth=0.8); ax.axvline(0, color="gray", linewidth=0.8)
ax.set_title("Product map on the first two principal components")
ax.set_xlabel("PC1 score"); ax.set_ylabel("PC2 score")
ax.grid(True, alpha=0.3); fig.tight_layout(); plt.show()

	PC1	PC2
Dimension	-0.413	-0.440
Roundness	-0.452	-0.430
Roughness	-0.270	0.105
Hardness	-0.355	0.561
Vibration	-0.543	-0.083
Temperature	-0.362	0.537

png

結果の読み取り

近い品種は6特性の変動パターンが似ています。負荷量から、各軸が寸法・真円度寄りか、硬度・温度寄りかを確認できます。離れた品種は調査候補ですが、不良と断定せず、品種仕様の違いを層別してから原因を確認します。

No.067：NMF（非負値行列因子分解）

実務での意味

不良件数や寄与度のような非負データでは、正負が相殺する PCA より、非負の原因パターンを足し合わせる NMF が説明しやすい場合があります。「各品種に、各不良モードがどれだけ含まれるか」と読めます。

分析・モデル化の考え方

$D\approx WH$ とし、 $D,W,H\ge0$ を課します。ここでは Frobenius 誤差を減らす乗法更新を使います。解は初期値やスケールに依存し、一意ではありません。複数 seed の安定性と現場ラベルとの整合を確認します。

Pythonで確認する

D_defect = np.maximum(raw_matrix - raw_matrix.min(axis=0) + 0.15, 0)
k_nmf = 3
W = rng.random((D_defect.shape[0], k_nmf)) + 0.1
H = rng.random((k_nmf, D_defect.shape[1])) + 0.1
eps = 1e-10
for _ in range(800):
    H *= (W.T @ D_defect) / (W.T @ W @ H + eps)
    W *= (D_defect @ H.T) / (W @ H @ H.T + eps)
H_scaled = H / H.sum(axis=1, keepdims=True)
display(pd.DataFrame(H_scaled, index=["Pattern 1", "Pattern 2", "Pattern 3"], columns=features).round(3))
print(f"relative NMF error = {np.linalg.norm(D_defect-W@H)/np.linalg.norm(D_defect):.3f}")
fig, ax = plt.subplots(figsize=(8, 3.5))
ax.imshow(H_scaled, cmap="YlOrRd", aspect="auto")
ax.set_title("Quality-feature composition of NMF patterns")
ax.set_xlabel("Quality feature"); ax.set_ylabel("Latent pattern")
ax.set_xticks(range(len(features)), features, rotation=30, ha="right")
ax.set_yticks(range(3), ["Pattern 1", "Pattern 2", "Pattern 3"]); ax.grid(False)
fig.tight_layout(); plt.show()

	Dimension	Roundness	Roughness	Hardness	Vibration	Temperature
Pattern 1	0.426	0.302	0.002	0.025	0.182	0.063
Pattern 2	0.028	0.003	0.063	0.404	0.081	0.421
Pattern 3	0.009	0.108	0.684	0.033	0.165	0.001

relative NMF error = 0.036

png

結果の読み取り

各パターンは非負の特性構成として読めるため、原因仮説を現場用語へ置き換えやすくなります。ただし潜在パターンは統計的な共起であり、物理原因の証明ではありません。初期値を変えて再現するか、保全履歴や実験で裏づけます。

No.068：Tensor分解

実務での意味

品種×設備×特性を行列へ潰すと、どの設備で現れたパターンかが失われます。Tensor 分解は三方向の構造を保ち、品種群、設備群、特性群を同時に圧縮します。

分析・モデル化の考え方

ここでは各モードへ展開した行列を SVD し、上位基底へ射影する HOSVD を使います。コアテンソル $G$ と因子行列 $U_p,U_m,U_f$ により、 $T\approx G\times_1U_p\times_2U_m\times_3U_f$ と表します。ランク $(r_p,r_m,r_f)$ は軸ごとに選びます。

Pythonで確認する

T = quality_tensor
Up = np.linalg.svd(T.reshape(T.shape[0], -1), full_matrices=False)[0][:, :3]
Um = np.linalg.svd(T.transpose(1,0,2).reshape(T.shape[1], -1), full_matrices=False)[0][:, :2]
Uf = np.linalg.svd(T.transpose(2,0,1).reshape(T.shape[2], -1), full_matrices=False)[0][:, :3]
core = np.einsum("pa,mb,fc,pmf->abc", Up, Um, Uf, T)
T_hat = np.einsum("pa,mb,fc,abc->pmf", Up, Um, Uf, core)
machine_loading = pd.DataFrame(Um, index=machines, columns=["Machine mode 1", "Machine mode 2"])
display(machine_loading.round(3))
print(f"original shape={T.shape}, core shape={core.shape}")
print(f"relative HOSVD error = {np.linalg.norm(T-T_hat)/np.linalg.norm(T):.3f}")
fig, ax = plt.subplots(figsize=(7, 3.5))
ax.bar(machines, np.abs(Um[:, 0]), color="#59A14F")
ax.set_title("Absolute loadings of the leading machine mode")
ax.set_xlabel("Machine"); ax.set_ylabel("Absolute loading")
ax.grid(True, axis="y", alpha=0.3); fig.tight_layout(); plt.show()

	Machine mode 1	Machine mode 2
M01	-0.381	0.148
M02	-0.309	-0.463
M03	-0.352	-0.144
M04	-0.405	0.709
M05	-0.349	-0.034
M06	-0.318	-0.425
M07	-0.368	0.158
M08	-0.337	-0.182

original shape=(12, 8, 6), core shape=(3, 2, 3)
relative HOSVD error = 0.308

png

結果の読み取り

元の12×8×6データを3×2×3のコアと三つの因子行列へ圧縮できました。設備モードの絶対値が大きい設備は共通変動への寄与候補です。ただし符号と回転の任意性があり、「設備の悪さ」そのものではないため、残差を設備別に確認します。

No.069：CUR分解

実務での意味

SVD の潜在ベクトルは複数特性の混合で、説明が難しいことがあります。CUR 分解は元データの実在する列 $C$ と行 $R$ を選ぶため、「代表特性」と「代表品種」を使った近似として説明できます。

分析・モデル化の考え方

$X\approx CUR$ とし、列・行の重要度には上位特異ベクトルの leverage score を使います。ここでは再現性を優先して上位スコアを決定的に選び、中央行列を $U=C^+XR^+$ とします。代表性は得られますが、同質な列へ偏らないか確認が必要です。

Pythonで確認する

k_cur, n_cols, n_rows = 2, 3, 5
col_score = np.sum(Vt[:k_cur]**2, axis=0)
row_score = np.sum(U[:, :k_cur]**2, axis=1)
col_idx = np.argsort(col_score)[-n_cols:]
row_idx = np.argsort(row_score)[-n_rows:]
C_cur, R_cur = X[:, col_idx], X[row_idx, :]
U_cur = np.linalg.pinv(C_cur) @ X @ np.linalg.pinv(R_cur)
X_cur = C_cur @ U_cur @ R_cur
display(pd.DataFrame({"Selected feature": np.array(features)[col_idx], "Leverage score": col_score[col_idx]}).round(3))
display(pd.DataFrame({"Selected product": np.array(products)[row_idx], "Leverage score": row_score[row_idx]}).round(3))
print(f"relative CUR error = {np.linalg.norm(X-X_cur)/np.linalg.norm(X):.3f}")

	Selected feature	Leverage score
0	Roundness	0.39
1	Temperature	0.42
2	Hardness	0.44

	Selected product	Leverage score
0	P05	0.189
1	P08	0.206
2	P09	0.295
3	P07	0.330
4	P10	0.382

relative CUR error = 0.241

結果の読み取り

選ばれた列と行は、潜在軸をよく表す実在の特性・品種です。現場説明や代表サンプル設計に使いやすい一方、選択数が少ないと誤差は SVD の最適低ランク近似より大きくなります。代表性と精度の交換条件を明示します。

No.070：行列補完

実務での意味

検査省略や一時的なセンサー停止で空いた値を補えれば、監視を継続できます。ただし補完値は観測値ではなくモデル推定であり、出荷判定を無条件に代替するものではありません。

分析・モデル化の考え方

観測集合 $\Omega$ 上の誤差を小さくしながら低ランク行列を求めます。ここでは観測値を固定し、欠損値をランク3の SVD 近似で反復更新する簡易 Soft-Impute 型手順を使います。評価は学習に使わない人工的な検証欠損で RMSE を測ります。

Pythonで確認する

mask_observed = rng.random(X.shape) > 0.22
# 各行・各列に観測値が残るよう保証
mask_observed[:, 0] = True; mask_observed[0, :] = True
X_missing = X.copy(); X_missing[~mask_observed] = np.nan
col_mean = np.nanmean(X_missing, axis=0)
X_fill = np.where(mask_observed, X_missing, col_mean)
for _ in range(150):
    u_i, s_i, vt_i = np.linalg.svd(X_fill, full_matrices=False)
    low_rank = u_i[:, :3] @ np.diag(s_i[:3]) @ vt_i[:3]
    updated = np.where(mask_observed, X, low_rank)
    if np.linalg.norm(updated-X_fill) < 1e-7: break
    X_fill = updated
missing_rmse = np.sqrt(np.mean((X_fill[~mask_observed]-X[~mask_observed])**2))
mean_fill = np.where(mask_observed, X, col_mean)
baseline_rmse = np.sqrt(np.mean((mean_fill[~mask_observed]-X[~mask_observed])**2))
print(f"missing cells = {(~mask_observed).sum()} / {X.size}")
print(f"low-rank completion RMSE = {missing_rmse:.3f}")
print(f"column-mean baseline RMSE = {baseline_rmse:.3f}")
fig, ax = plt.subplots(figsize=(6, 4))
ax.scatter(X[~mask_observed], X_fill[~mask_observed], color="#E15759")
lims = [X[~mask_observed].min()-.2, X[~mask_observed].max()+.2]
ax.plot(lims, lims, "--", color="gray", label="ideal")
ax.set_title("Held-out values versus low-rank estimates")
ax.set_xlabel("True standardized value"); ax.set_ylabel("Estimated standardized value")
ax.grid(True, alpha=0.3); ax.legend(); fig.tight_layout(); plt.show()

missing cells = 18 / 72
low-rank completion RMSE = 0.126
column-mean baseline RMSE = 0.767

png

結果の読み取り

人工的に隠したセルでは、低ランク補完と列平均の RMSE を直接比較できます。今回のように共通因子があるデータでは低ランク構造を利用する価値があります。ただし本番で「故障時だけ欠ける」なら欠損はランダムではなく、検証結果が楽観的になります。補完値にはフラグと不確実性を付け、重要検査は再測定します。

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

行列分解は一つの万能手法ではなく、目的別の道具箱です。

方程式を解く：係数行列を再利用するなら LU、対称正定値なら Cholesky、最小二乗なら QR を候補にする
構造を知る：SVD で有効ランクを診断し、低ランク近似で精度とコストを管理する
説明する：中心化した変動には PCA、非負の加法的寄与には NMF、実在する代表行列が必要なら CUR を使う
多方向を保つ：品種×設備×特性を安易に平坦化せず、Tensor 分解を検討する
欠損を扱う：行列補完の前提、検証欠損、ベースライン比較、推定値フラグを必須にする

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

目的と損失を定義する：圧縮、可視化、予測、欠損補完のどれが主目的かを明確にする
データ粒度を統一する：品種、設備、時刻、測定単位、規格変更履歴を揃える
前提を監視する：対称正定値性、条件数、ランク、非負性、欠損メカニズムを確認する
時系列で検証する：ランダム分割だけでなく、将来期間・新規品種・設備変更後で評価する
残差を管理する：低ランクモデルが説明しない局所異常を管理図やアラートへ残す
意思決定ルールを分離する：推定値だけで出荷可否を自動決定せず、再測定や承認条件を定める
再現性を確保する：コード、seed、前処理、モデル次数、評価結果、データ版を記録する

まとめ

No.061〜No.070 では、三角分解による安定な求解から、SVD 系の圧縮、説明可能な因子分解、三次元構造の保持、欠損補完までを一つの架空品質データで確認しました。重要なのは、分解誤差が小さいことだけではありません。前提が現場データに合い、残差と不確実性を含めて意思決定者へ説明できることが導入条件です。

法人向けのご相談

数理工房では、製造データの構造診断、品質 KPI の設計、次元削減・異常検知・欠損補完の PoC、既存分析の数値安定性レビュー、現場向け研修まで支援しています。課題と利用可能なデータを伺い、意思決定に必要な最小構成から整理します。

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

製造業の行列分解入門｜SVD・PCA・NMF・行列補完をPythonで実践

欠損のある多品種品質データを意思決定へつなぐ行列分解 — 製造業の線形代数100本ノック No.061〜No.070

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

現場でよくある状況

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

今回扱うノックの全体像

Python 環境の準備

架空データの作成

No.061：LU分解

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.062：Cholesky分解

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.063：QR分解

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.064：SVD（特異値分解）

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.065：低ランク近似

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.066：PCA（主成分分析）

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.067：NMF（非負値行列因子分解）

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.068：Tensor分解

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.069：CUR分解

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.070：行列補完

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

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

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

まとめ

法人向けのご相談

Related Posts

製造業で学ぶ線形代数10選｜FEM・最適化・デジタルツイン・AIをPythonで実践

製造業で学ぶ線形代数10選｜回帰・PCA・GCN・カルマンフィルタをPython実践

製造業の設備推薦と大規模行列計算｜Pythonで学ぶ行列分解10本ノック

製造工程の安定性を固有値で見抜く｜Pythonで学ぶ線形代数 No.051〜060

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