100本ノック / 線形代数 / 線形代数100本ノック
製造業の設備推薦と大規模行列計算｜Pythonで学ぶ行列分解10本ノック
設備×品種の適合度補完と大規模工程計算 — 製造業の線形代数100本ノック No.071〜No.080
多品種生産では、すべての設備と品種の組合せを試作することはできません。一方、日々の生産計画や設備条件の設計では、未評価の組合せを含む大きな行列を扱い、限られた時間で答えを得る必要があります。本稿では、架空の設備×品種データを用い、行列分解による推薦、潜在因子、Truncated / Randomized SVD、Krylov 部分空間法、Arnoldi 法、共役勾配法、前処理、疎行列ソルバー、大規模行列計算を、製造現場の意思決定に結びつけます。
[!NOTE] 本資料は、数理工房 (もしくは代表である和山個人) が過去に企業研修において使用した notebook を企業様の許可を得て再構成・編集のうえ公開しています。
掲載データはすべて架空のものであり、実在する企業・工場・数値とは一切関係ありません。
はじめに：この記事で扱う製造業の実務課題
12台の設備で18品種を加工するラインを考えます。試作済みの組合せでは、良品率・段取り時間・エネルギー原単位を統合した「適合度」が分かりますが、多くは未評価です。本稿の前半では、観測済み評価から未評価組合せの優先候補を推定します。後半では、同じ工場で発生する大規模な温度場の線形方程式を、疎行列と反復法で効率よく解きます。
現場でよくある状況
- 実績のない設備×品種を、担当者の経験だけで割り当てている
- 全組合せ試験は時間、材料、停止損失の面で実施できない
- シミュレーションのメッシュを細かくすると、計算時間とメモリが急増する
- 解が得られても、誤差・収束条件・候補選定ルールが管理されていない
なぜこの問題は判断が難しいのか
未観測は「不適合」ではなく「まだ試していない」です。また、低ランクモデルのスコアは因果効果や安全保証ではありません。計算面でも、巨大な密行列を明示的に作るとメモリが先に尽きます。精度、計算費用、説明可能性、安全制約を同時に扱う必要があります。
今回扱うノックの全体像
No.071〜074で設備候補を推薦する低ランクモデルを作り、No.075〜080で大規模線形計算を「必要な部分だけ計算する」方法へ進みます。前半は探索する組合せの優先順位、後半は制約時間内に信頼できる数値解を返す方法が意思決定の中心です。
Python 環境の準備
NumPy と pandas で計算・表作成、Matplotlib で可視化、SciPy で疎行列と反復法、scikit-learn で Randomized SVD を扱います。乱数 seed は固定し、再実行可能にします。
import sys, time
import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
from IPython.display import display
from scipy import sparse
from scipy.sparse.linalg import LinearOperator, cg, eigsh, spsolve
from sklearn.utils.extmath import randomized_svd
rng = np.random.default_rng(71080)
np.set_printoptions(precision=4, suppress=True)
pd.set_option("display.max_columns", 12)
print("Python:", sys.version.split()[0])
print("NumPy:", np.__version__, "| pandas:", pd.__version__)
print("Matplotlib:", matplotlib.__version__)Python: 3.13.1
NumPy: 2.5.1 | pandas: 3.0.3
Matplotlib: 3.11.0架空データの作成
設備能力と品種要求には3つの潜在軸（精密加工、熱安定性、量産適性）があると仮定します。真の適合度に小さな測定ノイズを加え、約55%だけを試作済みとします。欠測マスクは評価値とは独立に生成し、各設備・各品種に最低1件の観測があることを確認します。
n_machines, n_products, rank_true = 12, 18, 3
machines = [f"M{i:02d}" for i in range(1, n_machines + 1)]
products = [f"P{i:02d}" for i in range(1, n_products + 1)]
U_true = rng.normal(size=(n_machines, rank_true))
V_true = rng.normal(size=(n_products, rank_true))
raw = U_true @ V_true.T
true_score = 50 + 12 * (raw - raw.mean()) / raw.std()
true_score = np.clip(true_score, 5, 95)
observed_score = np.clip(true_score + rng.normal(0, 2.0, true_score.shape), 0, 100)
mask = rng.random(true_score.shape) < 0.55
for i in range(n_machines): mask[i, rng.integers(n_products)] = True
for j in range(n_products): mask[rng.integers(n_machines), j] = True
R = np.where(mask, observed_score, np.nan)
print(f"Matrix shape: {R.shape}, observed: {mask.sum()}/{R.size} ({mask.mean():.1%})")
display(pd.DataFrame(R, index=machines, columns=products).iloc[:6, :9].round(1))Matrix shape: (12, 18), observed: 137/216 (63.4%)| P01 | P02 | P03 | P04 | P05 | P06 | P07 | P08 | P09 | |
|---|---|---|---|---|---|---|---|---|---|
| M01 | 54.1 | NaN | 41.2 | 59.5 | 49.1 | 59.3 | 56.3 | NaN | NaN |
| M02 | 57.8 | 59.8 | NaN | 62.0 | 59.0 | NaN | 61.3 | 58.2 | 62.4 |
| M03 | 50.2 | 46.6 | 45.9 | NaN | 51.7 | NaN | NaN | 51.2 | 40.1 |
| M04 | NaN | 61.9 | 38.8 | 62.5 | 72.4 | NaN | NaN | 60.9 | NaN |
| M05 | NaN | NaN | NaN | 58.5 | 40.4 | 49.1 | 47.7 | 53.6 | 12.7 |
| M06 | 50.7 | 48.1 | NaN | 46.4 | 38.0 | 51.4 | 52.3 | 45.8 | 67.4 |
No.071：推薦システムの行列分解
実務での意味
推薦システムはEC商品だけでなく、未評価の設備×品種から「次に試作する価値が高い組合せ」を順位付けする仕組みとして使えます。安全認定を自動化するものではなく、試験計画の絞り込みに使います。
分析・モデル化の考え方
観測集合を とし、低ランク行列 が観測値を再現するよう、 を交互最小二乗（ALS）で解きます。未観測セルは損失に入れません。
Pythonで確認する
def als(R, mask, k=3, reg=3.0, epochs=30, seed=1):
g = np.random.default_rng(seed); m, n = R.shape
U = g.normal(0, .2, (m, k)); V = g.normal(0, .2, (n, k)); I = np.eye(k)
for _ in range(epochs):
for i in range(m):
js = np.where(mask[i])[0]; U[i] = np.linalg.solve(V[js].T @ V[js] + reg*I, V[js].T @ R[i, js])
for j in range(n):
ii = np.where(mask[:, j])[0]; V[j] = np.linalg.solve(U[ii].T @ U[ii] + reg*I, U[ii].T @ R[ii, j])
return U, V
U, V = als(np.nan_to_num(R), mask)
pred = U @ V.T
candidates = [(pred[i,j], machines[i], products[j]) for i,j in zip(*np.where(~mask))]
top5 = pd.DataFrame(sorted(candidates, reverse=True)[:5], columns=["Predicted score","Machine","Product"])
display(top5.round(1))| Predicted score | Machine | Product | |
|---|---|---|---|
| 0 | 90.6 | M12 | P10 |
| 1 | 84.2 | M12 | P09 |
| 2 | 72.1 | M06 | P03 |
| 3 | 71.5 | M12 | P03 |
| 4 | 69.8 | M09 | P02 |
結果の読み取り
上位5件は「適合確定」ではなく試作候補です。候補を、設備仕様、治具可否、法規・安全条件で除外した後に、少量試作へ回します。スコアの高い順だけでなく、予測の不確実性や試験費用も次段階で加味します。
No.072：Latent Factor Model
実務での意味
潜在因子は、表に直接記録されていない設備能力・品種要求の共通軸を圧縮して表します。設備や品種を似た傾向で配置でき、候補理由を技術者と議論する入口になります。
分析・モデル化の考え方
予測値 は、設備側と品種側の因子の内積です。ただし回転しても予測値は変わらないため、各軸を自動的に「剛性」などと断定せず、既知の仕様との相関で解釈します。
Pythonで確認する
factor_df = pd.DataFrame(U, index=machines, columns=["Factor 1","Factor 2","Factor 3"])
display(factor_df.round(2))
fig, ax = plt.subplots(figsize=(7, 5))
ax.scatter(U[:,0], U[:,1], s=70)
for i, name in enumerate(machines): ax.annotate(name, (U[i,0], U[i,1]), xytext=(4,4), textcoords="offset points")
ax.set_title("Machine map in latent-factor space")
ax.set_xlabel("Factor 1"); ax.set_ylabel("Factor 2"); ax.grid(True, alpha=.3); plt.tight_layout(); plt.show()| Factor 1 | Factor 2 | Factor 3 | |
|---|---|---|---|
| M01 | -2.13 | -6.01 | -7.03 |
| M02 | -2.79 | -3.09 | -9.36 |
| M03 | -2.23 | -4.10 | -7.06 |
| M04 | -6.09 | -5.42 | -6.85 |
| M05 | 0.29 | -6.76 | -5.46 |
| M06 | 0.11 | -0.83 | -10.22 |
| M07 | -2.50 | -7.18 | -5.65 |
| M08 | 3.10 | -4.49 | -6.91 |
| M09 | -6.86 | -0.43 | -9.29 |
| M10 | -3.02 | -5.29 | -6.61 |
| M11 | 0.40 | -2.80 | -8.37 |
| M12 | -3.00 | 2.20 | -11.07 |
結果の読み取り
近くに配置された設備は、観測された適合度パターンが似ています。保全方式や主軸仕様などの台帳と照合し、因子と既知属性の対応が一貫するかを確認します。因子マップは因果説明ではなく、設備群の比較と追加調査のための地図です。
No.073：Truncated SVD
実務での意味
Truncated SVD は、大きな評価行列を上位 成分だけで近似し、主要パターンを残しながらデータ量を減らします。全成分を保存・計算しないため、予測や可視化を軽量化できます。
分析・モデル化の考え方
とし、特異値の二乗比で保持情報量を評価します。SVD は本来完全行列向けなので、ここでは欠測を品種平均で仮補完したベースラインとして扱い、ALSとの違いを明記します。
Pythonで確認する
filled = R.copy(); col_means = np.nanmean(filled, axis=0)
filled[np.where(np.isnan(filled))] = np.take(col_means, np.where(np.isnan(filled))[1])
Uc, s, Vt = np.linalg.svd(filled, full_matrices=False)
errors=[]
for k in range(1, 9):
approx=(Uc[:,:k]*s[:k])@Vt[:k]
errors.append(np.linalg.norm(filled-approx,"fro")/np.linalg.norm(filled,"fro"))
display(pd.DataFrame({"rank":range(1,9),"relative_error":errors}).round(4))
fig, ax=plt.subplots(figsize=(7,4)); ax.plot(range(1,9),errors,marker="o")
ax.set_title("Truncated SVD: rank and reconstruction error"); ax.set_xlabel("Retained rank"); ax.set_ylabel("Relative Frobenius error")
ax.grid(True,alpha=.3); plt.tight_layout(); plt.show()| rank | relative_error | |
|---|---|---|
| 0 | 1 | 0.1669 |
| 1 | 2 | 0.1191 |
| 2 | 3 | 0.0944 |
| 3 | 4 | 0.0765 |
| 4 | 5 | 0.0574 |
| 5 | 6 | 0.0472 |
| 6 | 7 | 0.0378 |
| 7 | 8 | 0.0262 |
結果の読み取り
ランクを増やすほど再構成誤差は下がりますが、複雑性も増えます。折れ曲がりだけで決めず、観測データを訓練・検証に分けた予測誤差と、運用可能な計算時間でランクを選びます。平均補完は欠測機構に偏りがある場合、系統的な誤差を生みます。
No.074：Randomized SVD
実務での意味
Randomized SVD は、必要な上位特異成分を確率的な射影で高速に近似します。設備・品種が増え、毎日の再学習時間が制約になる場合に有効です。
分析・モデル化の考え方
ランダム行列で主要な列空間を抽出し、小さな行列に対してSVDを行います。seed を固定して再現性を持たせ、厳密SVDとの差を再構成誤差で測ります。
Pythonで確認する
k=3
Ur,sr,Vtr=randomized_svd(filled,n_components=k,n_iter=5,random_state=71080)
rand_approx=(Ur*sr)@Vtr; exact_approx=(Uc[:,:k]*s[:k])@Vt[:k]
comparison=pd.DataFrame({"method":["Exact truncated SVD","Randomized SVD"],
"relative_error":[np.linalg.norm(filled-exact_approx,"fro")/np.linalg.norm(filled,"fro"),np.linalg.norm(filled-rand_approx,"fro")/np.linalg.norm(filled,"fro")],
"top_singular_value":[s[0],sr[0]]})
display(comparison.round(6))| method | relative_error | top_singular_value | |
|---|---|---|---|
| 0 | Exact truncated SVD | 0.094438 | 741.215763 |
| 1 | Randomized SVD | 0.094438 | 741.215763 |
結果の読み取り
小規模例でも Randomized SVD の誤差は厳密なランク3近似に近づきます。速度差は大規模データで現れるため、本番相当サイズで時間とメモリを測ります。乱数によるばらつき、反復回数、許容誤差を運用設定として残します。
No.075：Krylov部分空間法
実務での意味
温度場や変形解析では、行列全体の分解より、行列とベクトルの積を繰り返して必要な解・固有モードだけ得る方が現実的です。Krylov 法はその基盤です。
分析・モデル化の考え方
を作ります。ここでは1次元熱伝導の離散行列について基底を直交化し、部分空間が拡大する様子を確認します。
Pythonで確認する
n=120
A_heat=sparse.diags([-np.ones(n-1),2.2*np.ones(n),-np.ones(n-1)],[-1,0,1],format="csr")
b=np.zeros(n); b[n//3:2*n//3]=1.0
Q=[]; v=b/np.linalg.norm(b)
for _ in range(8):
for q in Q: v-=q*(q@v)
v/=np.linalg.norm(v); Q.append(v.copy()); v=A_heat@v
Q=np.column_stack(Q)
orth_error=np.linalg.norm(Q.T@Q-np.eye(Q.shape[1]))
print("Krylov basis shape:",Q.shape,"| orthogonality error:",f"{orth_error:.2e}")
fig,ax=plt.subplots(figsize=(8,4));
for j in [0,1,3,7]: ax.plot(Q[:,j],label=f"q{j+1}")
ax.set_title("Selected Krylov basis vectors"); ax.set_xlabel("Grid point"); ax.set_ylabel("Basis value")
ax.grid(True,alpha=.3); ax.legend(); plt.tight_layout(); plt.show()Krylov basis shape:
(120, 8) | orthogonality error: 2.89e-14
結果の読み取り
8本の直交基底だけで、熱源から周辺へ伝わる代表的な方向を表現しています。行列を逆行列へ変換せず、行列ベクトル積を中心に計算できる点が大規模問題で重要です。基底数は固定せず、残差が要求精度を満たすまで増やします。
No.076：Arnoldi法
実務での意味
Arnoldi 法は非対称な工程伝播行列の支配的固有値を近似し、変動が増幅するモードを少ない次元で監視する方法です。ここでは対称行列にも適用できる eigsh を使い、部分空間反復の考え方を確認します。
分析・モデル化の考え方
Arnoldi 法は を満たす直交基底と小さな Hessenberg 行列を作り、その固有値（Ritz値）で元の固有値を近似します。対称問題では Lanczos 法へ簡約されます。
Pythonで確認する
ritz=eigsh(A_heat,k=3,which="LM",return_eigenvectors=False)
exact=np.linalg.eigvalsh(A_heat.toarray())[-3:]
display(pd.DataFrame({"exact":exact[::-1],"Ritz approximation":np.sort(ritz)[::-1],
"absolute_error":np.abs(exact[::-1]-np.sort(ritz)[::-1])}).round(8))| exact | Ritz approximation | absolute_error | |
|---|---|---|---|
| 0 | 4.199326 | 4.199326 | 0.0 |
| 1 | 4.197304 | 4.197304 | 0.0 |
| 2 | 4.193936 | 4.193936 | 0.0 |
結果の読み取り
最大側の3固有値が高精度に得られ、全固有値・全固有ベクトルを計算する必要がありません。実務では対象が非対称か対称か、求めたい固有値が最大・最小・特定区間のどれかでアルゴリズム設定を変えます。
No.077：共役勾配法
実務での意味
共役勾配法（CG）は、対称正定値な大規模連立方程式 を、逆行列を作らずに解きます。熱伝導・構造解析・二次最適化で頻出します。
分析・モデル化の考え方
CG は二次関数 を互いに 共役な方向で最小化します。停止判定は反復回数ではなく相対残差 で管理します。
Pythonで確認する
residuals=[]
def cb(xk): residuals.append(np.linalg.norm(b-A_heat@xk)/np.linalg.norm(b))
x_cg,info=cg(A_heat,b,rtol=1e-10,callback=cb)
print("info:",info,"| iterations:",len(residuals),"| final relative residual:",f"{residuals[-1]:.2e}")
fig,ax=plt.subplots(figsize=(7,4)); ax.semilogy(range(1,len(residuals)+1),residuals)
ax.set_title("Convergence of conjugate gradient"); ax.set_xlabel("Iteration"); ax.set_ylabel("Relative residual")
ax.grid(True,which="both",alpha=.3); plt.tight_layout(); plt.show()info: 0 | iterations: 51 | final relative residual: 7.13e-11
結果の読み取り
残差が単調に十分小さくなり、info=0 は収束を示します。ただし小さい残差がそのまま物理モデルの妥当性を保証するわけではありません。境界条件、材料定数、メッシュ依存性は別途検証します。
No.078：前処理法
実務での意味
係数スケールが大きく異なると、同じCGでも反復回数が増えます。前処理は解を変えずに「解きやすい形」へ変換し、日次バッチや制御周期に間に合わせます。
分析・モデル化の考え方
ここでは対角成分を使う Jacobi 前処理 を用います。条件数を改善するほど一般に収束が速くなりますが、前処理自体の作成・適用費用との合計で判断します。
Pythonで確認する
scale=np.geomspace(1,1e3,n); D=sparse.diags(scale); A_scaled=D@A_heat@D; b_scaled=D@b
counts={}
for label,M in [("None",None),("Jacobi",LinearOperator((n,n),matvec=lambda x:x/A_scaled.diagonal()))]:
hist=[]
x,info=cg(A_scaled,b_scaled,rtol=1e-8,M=M,callback=lambda xk,h=hist:h.append(np.linalg.norm(b_scaled-A_scaled@xk)/np.linalg.norm(b_scaled)),maxiter=2000)
counts[label]=(len(hist),hist[-1],info)
display(pd.DataFrame(counts,index=["iterations","final_residual","info"]).T)| iterations | final_residual | info | |
|---|---|---|---|
| None | 1288.0 | 7.501723e-09 | 0.0 |
| Jacobi | 45.0 | 7.642589e-09 | 0.0 |
結果の読み取り
対角スケーリングで悪化させた問題に対し、Jacobi 前処理は反復回数を大きく削減します。本番では不完全Choleskyなども候補ですが、前処理の構築時間、追加メモリ、複数右辺での再利用回数を含めて比較します。
No.079：疎行列ソルバー
実務での意味
局所的な相互作用から生じる行列の大半はゼロです。疎形式を使えば、ゼロを保存・演算せず、より細かいメッシュを同じ計算資源で扱えます。
分析・モデル化の考え方
の密行列は概ね byte を必要とします。CSR形式は非ゼロ値、列番号、行ポインタを保持し、今回の三重対角行列では非ゼロ数が約 に留まります。直接法 spsolve とCGの解も比較します。
Pythonで確認する
n_big=5000
A_big=sparse.diags([-np.ones(n_big-1),2.2*np.ones(n_big),-np.ones(n_big-1)],[-1,0,1],format="csr")
b_big=np.ones(n_big)
dense_mb=n_big*n_big*8/1024**2
sparse_mb=(A_big.data.nbytes+A_big.indices.nbytes+A_big.indptr.nbytes)/1024**2
t0=time.perf_counter(); x_direct=spsolve(A_big,b_big); elapsed=time.perf_counter()-t0
display(pd.DataFrame({"representation":["Dense (estimated)","CSR (actual)"],"memory_MB":[dense_mb,sparse_mb]}).round(3))
print("nonzeros:",A_big.nnz,"| spsolve time:",f"{elapsed:.4f}s","| relative residual:",f"{np.linalg.norm(b_big-A_big@x_direct)/np.linalg.norm(b_big):.2e}")| representation | memory_MB | |
|---|---|---|
| 0 | Dense (estimated) | 190.735 |
| 1 | CSR (actual) | 0.191 |
nonzeros: 14998 | spsolve time: 0.0018s | relative residual: 1.27e-16結果の読み取り
CSRの実メモリは密形式の推定値より桁違いに小さく、直接法でも十分小さい残差が得られます。ただし疎行列でも分解時の fill-in によりメモリが増えるため、行列サイズ・構造・右辺数に応じて直接法と反復法を選びます。
No.080：大規模行列計算
実務での意味
大規模計算では「行列を作れるか」より、必要な演算だけを定義し、精度・時間・メモリの予算内で結果を返せるかが重要です。Digital Twin や生産計画の更新頻度を左右します。
分析・モデル化の考え方
LinearOperator は行列を明示せず だけを定義します。1,000,000次元の三重対角作用を、密行列なら約7.3 TiB必要なところ、数本のベクトルで実行します。性能評価はウォームアップや反復測定を含めるのが本来ですが、ここでは構造上の差を確認します。
Pythonで確認する
n_matrix_free=1_000_000
def heat_matvec(x):
y=2.2*x.copy(); y[1:]-=x[:-1]; y[:-1]-=x[1:]; return y
op=LinearOperator((n_matrix_free,n_matrix_free),matvec=heat_matvec,dtype=float)
x=np.ones(n_matrix_free); t0=time.perf_counter(); y=op@x; elapsed=time.perf_counter()-t0
dense_tib=n_matrix_free**2*8/1024**4
print(f"dimension: {n_matrix_free:,}")
print(f"dense matrix estimate: {dense_tib:,.1f} TiB")
print(f"matrix-free matvec: {elapsed:.4f}s | vector result memory: {y.nbytes/1024**2:.1f} MiB")
print("result sample:",y[:3],y[-3:])dimension: 1,000,000
dense matrix estimate: 7.3 TiB
matrix-free matvec: 0.0026s | vector result memory: 7.6 MiB
result sample: [1.2 0.2 0.2] [0.2 0.2 1.2]結果の読み取り
100万次元でも、演算規則だけなら通常のメモリで行列ベクトル積を実行できます。これにKrylov反復法を組み合わせると、巨大な行列の明示生成を避けられます。ただし実運用では、計算ノード間通信、前処理、停止条件、障害時の再開設計まで含めた性能設計が必要です。
対象ノックを通して見える実務上の示唆
設備×品種の未評価セルは、低ランク構造を使って試作優先順位へ変換できます。一方、予測スコアは安全性・加工可否を保証しないため、ハード制約による候補除外と実機試験が不可欠です。計算面では、上位成分、部分空間、非ゼロ要素、行列ベクトル積だけを扱うことで、問題規模を大きくできます。
意思決定では、モデル精度だけでなく次の三つを同時に管理します。
- 業務価値：試作回数、立上げ期間、計算待ち時間をどれだけ減らすか
- 数値信頼性：検証誤差、残差、収束情報、再現性を記録しているか
- 安全な運用：推薦を自動承認にせず、設備制約と技術者レビューを通すか
実務導入する場合に必要なこと
- 適合度の定義を、良品率・CT・エネルギー・安全制約から部門横断で合意する
- 欠測がランダムか、成功例だけが記録される選択バイアスがないかを監査する
- 時系列分割や設備・品種を丸ごと外す検証で、未知条件への性能を測る
- 推薦候補、採否、試作結果を記録し、モデル更新へ戻す
- 疎形式、許容誤差、前処理、時間・メモリ上限を本番相当規模で負荷試験する
- モデル、データ、seed、ライブラリ、判断者を版管理し、異常時の停止条件を決める
まとめ
No.071〜080では、低ランク行列分解による設備候補の推薦から、Krylov部分空間、固有値近似、CG、前処理、疎行列、matrix-free計算までを一つの製造課題として確認しました。共通する原則は、巨大な情報をすべて計算するのではなく、意思決定に必要な構造を選んで計算することです。数学的な近似誤差と現場の判断リスクを分けて評価することで、PoCから安定運用へ進めます。
法人向けのご相談
数理工房では、設備×品種の試験計画、推薦・最適化モデル、CAEや大規模疎行列計算の高速化、現場データを用いた企業研修まで、課題整理から運用設計をご支援します。
📩 お問い合わせ: surikobo.co.jp/contact
まずはお気軽にご相談ください。