Himmelblau's System

This time, we challenge "Himmelblau's system," a multi-modal mathematical function widely used around the world as a benchmark for non-linear optimization. It is a beautiful yet troublesome system that possesses four distinct roots (valleys) on the XY plane.

On a personal note, I am a huge fan of the character "Hero Himmel" from the anime Frieren: Beyond Journey's End. I always tell my daughter, "If you ever get married, find a man like Himmel." So, when I saw the name of this equation, I felt a wonderful sense of affinity. (By the way, I also love Tanjiro from Demon Slayer, so he has my blessing as well! Haha.)

Typically, when using SPICE to explore all roots of an unknown non-linear equation, we consider methods like the Homotopy method, tracking the trajectory of roots by continuously varying parameters. However, Newton's method (SPICE's OP analysis) alone has a classic weakness: it only converges to a single root depending entirely on the initial value provided.

To explore the new possibilities of the SPICE-oriented approach, we decided not to rely on the Homotopy method. Instead, we adopted an entirely different strategy fusing Random Search (Monte Carlo method) with Machine-Learning Clustering to capture all four roots in one fell swoop.

First, we used SPICE's control structures (Nutmeg) to throw 10,000 random "darts" across a vast search space. Instead of blindly logging everything, we implemented a Filter & Dump technique—much like a firmware debug log—where the terminal only outputs coordinates when the error drops below a certain threshold.

Below is a snippet of the log data spat out by 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
...

The dozens of log lines generated by SPICE contain "data clusters" gathered around the true roots. How many clusters are there, and where is the center point (the one with the smallest error) for each?

We wrote a "bare-metal" Python script that performs clustering purely through distance calculations (the Pythagorean theorem) without relying on heavy external libraries like scikit-learn.

import re
import math

# SPICE log data
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. Parse using Regex
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. Custom Clustering (Grouping by distance 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. Extract the representative point with the smallest error from each cluster
print(f"Number of automatically detected roots: {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")

Here lies the greatest breakthrough. We now trigger Newton's method using the four extracted initial values. However, SPICE's OP analysis (DC analysis) has a frustrating trait: "it suffers from amnesia, wiping the slate clean and forgetting previous calculation states every time it executes." Therefore, sequentially solving them via a software loop will just suck you into the same single root over and over.

To bypass this software loop limitation, we utilized brute-force spatial engineering: "Copy and paste the equation-solving circuit block in the netlist as many times as there are discovered initial values (in this case, four), arranging them in parallel."

By assigning individual initial values to each circuit block (Core 1 to Core 4) using the .NODESET command, a single OP analysis simultaneously fires up four independent derivative engines, locking down all roots instantly.

Here is the MacSpice output log executing the parallelized approach.

= = = 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 = = =

• Analysis Engine: Mac SPICE3 (SPICE-oriented approach)
• PC: MacBook
• OS: macOS Monterey 12.7.6
• CPU: 1.2GHz Dual-Core Intel Core m5
• Memory: 8GB