Himmelblau's System
SPICE指向型解析法とデータサイエンスの融合による全解探索の軌跡
Himmelblau's System of Nonlinear Equations
Roots: 4 Points in XY Plane
1. 序論：ヒンメルブラウの方程式と「勇者」
今回は、非線形最適化のベンチマークとして世界中で使われている「ヒンメルブラウの連立方程式（Himmelblau's system）」に挑戦します。XY平面上に4つの解（谷底）を持つ、非常に美しくも厄介な方程式です。
余談ですが、私は『葬送のフリーレン』というアニメに出てくる「勇者ヒンメル」というキャラクターがすごく好きです。普段から自分の娘に「結婚するならヒンメルのような男性にしなさい」と言い聞かせているくらいなので、この方程式の名前を見たとき、個人的にとても素晴らしい名前だなと親近感が湧きました。（ちなみに『鬼滅の刃』の炭治郎も大好きなので、彼もオッケーと言っています。笑）
2. ホモトピー法からの脱却：全く新しいアプローチ
通常、SPICEを用いて未知の非線形方程式の全解を探索する場合、パラメータを連続的に変化させて解の軌跡を追う「ホモトピー法」などを検討します。 しかし、ニュートン法（SPICEのOP解析）単体では、初期値に依存して1つの解にしか辿り着けないという古典的な弱点があります。
そこで今回は、SPICE指向型解析法の新たな可能性を探るため、ホモトピー法には頼らず、「乱数探索（モンテカルロ法）」と「機械学習的クラスタリング」を融合させた別のアプローチで、4つの解をすべて一網打尽に求めてみることにしました。
3. モンテカルロ法による初期値候補の抽出（Filter & Dump）
まずはSPICEの制御構文（Nutmeg）を使って、広大な空間に1万回の乱数（ダーツ）を投げるモンテカルロ探索を行いました。 ただ闇雲に探すのではなく、ファームウェアのデバッグログのように「誤差が一定以下の有望な座標に落ちた時だけ、ターミナルにログを出力して捨てる」というFilter & Dumpの手法を実装しました。
以下が、SPICEが吐き出したログデータの一部です。
= = = Starting Monte Carlo search (Filter & Dump) = = =
Hit! X = 3.695 , Y = -2.767 , Err = 0.136394
Hit! X = 2.804 , Y = 2.839 , Err = 0.10767
Hit! X = -2.648 , Y = 3.681 , Err = 0.103959
Hit! X = -3.896 , Y = -3.879 , Err = 0.112582
...4. Pythonによるクラスタリング（解の絞り込み）
SPICEが吐き出した数十行のログには、真の解の周辺に集まった「データの塊（クラスター）」が含まれています。これがいくつあるのか、そしてそれぞれの塊の「一番誤差が小さい中心点」はどこかを自動で絞り出します。
今回は外部の重いライブラリ（scikit-learn等）を一切使わず、純粋な距離計算（ピタゴラスの定理）だけでクラスタリングを行う「ベアメタル」なPythonコードを自作しました。
import re
import math
# SPICEのログデータ
log_data = """
Hit! X = 3.695 , Y = -2.767 , Err = 0.136394
Hit! X = 2.804 , Y = 2.839 , Err = 0.10767
Hit! X = -2.648 , Y = 3.681 , Err = 0.103959
Hit! X = -3.896 , Y = -3.879 , Err = 0.112582
...
"""
# 1. 正規表現でパース
points = []
for line in log_data.strip().split('\n'):
m = re.search(r'X =\s*([-\d.]+)\s*,\s*Y =\s*([-\d.]+)\s*,\s*Err =\s*([-\d.]+)', line)
if m:
points.append({'x': float(m.group(1)), 'y': float(m.group(2)), 'err': float(m.group(3))})
# 2. 独自クラスタリング実装（距離 eps=1.0 でグループ化）
eps = 1.0
clusters = []
for p in points:
found_cluster = False
for cluster in clusters:
for cp in cluster:
dist = math.hypot(p['x'] - cp['x'], p['y'] - cp['y'])
if dist < eps:
cluster.append(p)
found_cluster = True
break
if found_cluster:
break
if not found_cluster:
clusters.append([p])
# 3. 各クラスターの中で最もエラー(Err)が小さい代表点を抽出
print(f"自動検出された解の数: {len(clusters)} 個\n")
for i, cluster in enumerate(clusters):
best_p = min(cluster, key=lambda item: item['err'])
print(f"* Root {i+1} (Err: {best_p['err']:.5f})")
print(f".NODESET V(x{i+1})={best_p['x']:.5f} V(y{i+1})={best_p['y']:.5f}\n")
5. 並列ニュートン法による全解の一括確定
ここからが最大のポイントです。抽出した4つの初期値を使ってニュートン法を起動するわけですが、SPICEのOP解析（DC解析）は「実行するたびに以前の計算状態を忘れて白紙に戻る」という厄介な仕様があります。そのため、ソフトウェア的なループ制御で順番に解こうとすると、すべて同じ解に吸い込まれてしまいます。
そこで、ループ制御に頼るのをやめ、「見つかった初期値の数（今回は4つ）だけ、方程式を解く回路ブロック自体をネットリスト上にコピペして並列に配置する」という空間的な力技を採用しました。
各回路ブロック（Core 1 〜 Core 4）に .NODESET で個別の初期値を与え、たった1回のOP解析を走らせることで、4つの独立した微分エンジンが同時に作動し、すべての解を一気に確定させます。
6. 解析結果：一撃での全解確定
上記の並列化アプローチを実行したMacSpiceの出力ログです。
= = = Starting 4-Core Parallel Newton Solver = = =
* Root 1 Exact Solution:
v(x1) = 3.58443e+00
v(y1) = -1.84813e+00
* Root 2 Exact Solution:
v(x2) = 3.00000e+00
v(y2) = 2.00001e+00
* Root 3 Exact Solution:
v(x3) = -2.80512e+00
v(y3) = 3.13131e+00
* Root 4 Exact Solution:
v(x4) = -3.77931e+00
v(y4) = -3.28319e+00
= = = All Roots Found = = =第1象限の美しい整数解（$x=3, y=2$）を含め、4つの真の解が完璧な精度で同時に求められました。
ソフトウェアの制御文で詰まったら、空間（回路構成）をパラレルにして解決する。複雑な判定はシミュレータ内で無理にやらず、ログだけ吐かせて外のPythonに任せる。 シミュレータの制約を、泥臭いエンジニアリングと「適材適所のツールチェーン構築」で突破し、未知の方程式を完全に制圧することができました。
7. 解析環境
- 解析エンジン： Mac SPICE3 (SPICE-oriented approach)
- PC： MacBook
- OS： macOS Monterey 12.7.6
- CPU： 1.2GHz デュアルコア Intel Core m5
- メモリ： 8GB