Ercsey-Ravász & Toroczkai, Nature Physics 2011 — reproduced

SAT in Phase Space

Replace a formula's truth values with real numbers, and satisfiability becomes a deterministic flow. It finds every solution — but the path there is a chaotic transient, and the basins of the solutions have fractal boundaries.

n = 16, α = 4.6, 74 clauses, 6 solutions · 1.8M integrated trajectories

View code & data — the chaotic flow and the Kantz/Grassberger test

chaos_sat.py — SAT as a continuous dynamical system

"""
SAT als kontinuierliches dynamisches System (Ercsey-Ravasz/Toroczkai 2011).

DAS SYSTEM. Variablen werden entspannt: s_i in [-1,1] statt {-1,+1}.
Die Verletzung einer Klausel m ist

    K_m(s) = 2^-k * prod_{i in m} (1 - c_mi s_i)      in [0,1]

K_m = 0 genau dann, wenn ein Literal voll erfuellt ist. Jede Klausel bekommt
ein Gewicht a_m, das exponentiell waechst, solange sie verletzt ist:

    ds_i/dt = 2 * sum_m a_m K_m K_mi c_mi        (K_mi = K_m ohne den Faktor i)
    da_m/dt = a_m K_m

Der s-Teil ist Gradientenabstieg auf V = sum_m a_m K_m^2; der a-Teil hebt
genau die Klauseln an, die haengen bleiben, und macht damit jedes lokale
Minimum instabil. Folge: **es gibt keine Attraktoren ausser den Loesungen.**
Der Wuerfel [-1,1]^n ist invariant.

Der Preis: der Weg dorthin ist ein CHAOTISCHER TRANSIENT. Die Zeit bis zur
Loesung ist unbeschraenkt, die Einzugsgebiete der Loesungen haben FRAKTALE
Raender, und die Rate, mit der Trajektorien das chaotische Gebiet verlassen
(die Fluchtrate kappa), geht an der SAT/UNSAT-Schwelle gegen null. Haerte
wird damit zu einer dynamischen Groesse: 1/kappa.

Alles hier ist ueber das Gitter vektorisiert -- ein RK4-Schritt bewegt
hunderttausend Trajektorien gleichzeitig.
"""
import struct
import zlib

import numpy as np


# ------------------------------------------------------------------ Aufbau

def bau(klauseln, n):
    """lit[m,k] Variablenindex, sgn[m,k] Vorzeichen, streu[3m,n] Streumatrix."""
    k = len(klauseln[0])
    m = len(klauseln)
    lit = np.zeros((m, k), np.int32)
    sgn = np.zeros((m, k), np.float32)
    for j, c in enumerate(klauseln):
        for l, x in enumerate(c):
            lit[j, l] = abs(x) - 1
            sgn[j, l] = 1.0 if x > 0 else -1.0
    streu = np.zeros((m * k, n), np.float32)
    streu[np.arange(m * k), lit.ravel()] = 1.0
    return lit, sgn, streu


def ableitung(s, a, lit, sgn, streu):
    """ds/dt und dK. s: (P,n), a: (P,M). Gibt (ds, K) zurueck."""
    P = s.shape[0]
    M, k = lit.shape
    f = 1.0 - sgn[None, :, :] * s[:, lit]          # (P,M,k), jeder Faktor in [0,2]
    K = f.prod(axis=2) * (0.5 ** k)                # (P,M)
    # K_mi = K ohne den Faktor i  ->  Produkt der uebrigen
    Kmi = np.empty_like(f)
    for l in range(k):
        andere = [x for x in range(k) if x != l]
        p = f[:, :, andere[0]]
        for x in andere[1:]:
            p = p * f[:, :, x]
        Kmi[:, :, l] = p * (0.5 ** (k - 1))
    beitrag = (2.0 * a * K)[:, :, None] * Kmi * sgn[None, :, :]
    ds = beitrag.reshape(P, M * k) @ streu
    return ds, K


def erfuellt(s, lit, sgn):
    """Ist sign(s) eine Loesung? (P,) bool"""
    z = np.sign(s)
    z[z == 0] = 1.0
    ok = (sgn[None, :, :] * z[:, lit]) > 0        # (P,M,k)
    return ok.any(axis=2).all(axis=1)


def kode(s):
    """Vorzeichenvektor als Ganzzahl -- Kennung der erreichten Loesung."""
    b = (s > 0).astype(np.int64)
    g = np.zeros(s.shape[0], np.int64)
    for j in range(s.shape[1]):
        g |= b[:, j] << j
    return g


# ------------------------------------------------------------------ Integration

def laufe(s0, klauseln, n, tmax=200.0, eta=0.05, dtmax=0.5, pruef=1,
          maxschritt=3000, adeckel=1e8):
    """Integriert alle Startpunkte gleichzeitig (RK4, punktweise Schrittweite).

    ZWEI Abbruchbedingungen, und beide werden gebraucht: die analoge Zeit t
    UND die Schrittzahl. Denn die Schrittweite ist ~ 1/|ds| ~ 1/a, und a
    waechst exponentiell -- im chaotischen Transienten friert t praktisch ein,
    waehrend die Trajektorie in s weiterwandert. Ohne die Schrittgrenze laeuft
    so ein Punkt endlos. Die Schrittzahl misst die BOGENLAENGE der Bahn und ist
    damit das ehrlichere Mass fuer den Rechenaufwand.

    Gibt (zeit, loesung, fertig) zurueck.
    """
    lit, sgn, streu = bau(klauseln, n)
    M = len(klauseln)
    P = s0.shape[0]
    s = s0.astype(np.float32).copy()
    a = np.ones((P, M), np.float32)
    t = np.zeros(P, np.float32)
    bogen = np.zeros(P, np.int32)
    zeit = np.full(P, np.nan, np.float32)
    schritte = np.full(P, -1, np.int32)
    loes = np.full(P, -1, np.int64)
    offen = np.ones(P, bool)
    runde = 0
    while offen.any() and runde < maxschritt:
        runde += 1
        idx = np.flatnonzero(offen)
        ss, aa = s[idx], a[idx]
        d1, K1 = ableitung(ss, aa, lit, sgn, streu)
        dt = np.clip(eta / (1e-6 + np.abs(d1).max(axis=1)), 1e-5, dtmax)[:, None]
        d2, K2 = ableitung(ss + 0.5 * dt * d1, aa * np.exp(0.5 * dt * K1), lit, sgn, streu)
        d3, K3 = ableitung(ss + 0.5 * dt * d2, aa * np.exp(0.5 * dt * K2), lit, sgn, streu)
        d4, K4 = ableitung(ss + dt * d3, aa * np.exp(dt * K3), lit, sgn, streu)
        ss = np.clip(ss + (dt / 6.0) * (d1 + 2 * d2 + 2 * d3 + d4), -1.0, 1.0)
        aa = np.minimum(aa * np.exp((dt / 6.0) * (K1 + 2 * K2 + 2 * K3 + K4)), adeckel)
        s[idx], a[idx] = ss, aa
        t[idx] += dt[:, 0]
        bogen[idx] += 1
        if runde % pruef == 0:
            fertig = erfuellt(ss, lit, sgn)
            if fertig.any():
                g = idx[fertig]
                zeit[g] = t[g]
                schritte[g] = bogen[g]
                loes[g] = kode(s[g])
                offen[g] = False
        raus = idx[t[idx] > tmax]
        offen[raus] = False
    return zeit, loes, np.isfinite(zeit), schritte


# ------------------------------------------------------------------ PNG

def png(pfad, bild):
    H, W, _ = bild.shape
    roh = b"".join(b"\x00" + bild[y].tobytes() for y in range(H))

    def brocken(typ, daten):
        return (struct.pack(">I", len(daten)) + typ + daten
                + struct.pack(">I", zlib.crc32(typ + daten) & 0xFFFFFFFF))

    d = b"\x89PNG\r\n\x1a\n"
    d += brocken(b"IHDR", struct.pack(">IIBBBBB", W, H, 8, 2, 0, 0, 0))
    d += brocken(b"IDAT", zlib.compress(roh, 6))
    d += brocken(b"IEND", b"")
    with open(pfad, "wb") as f:
        f.write(d)
    return pfad

lyapunov.py — escape rate, stretching rate, the Kantz/Grassberger test

"""
Woher der Fraktal kommt: die Streckrate lambda, und die Probe kappa/lambda.

DIE KETTE. Der Rand zwischen zwei Einzugsgebieten ist unter dem Fluss
INVARIANT -- wer genau auf ihm startet, konvergiert nie, denn er muesste den
Rand verlassen, und Raender gehen unter einem Fluss wieder in Raender ueber.
Auf dem Rand lebt also eine invariante Menge ohne Konvergenz: der chaotische
SATTEL. Der Rand selbst ist dessen STABILE MANNIGFALTIGKEIT.

Ein Sattel hat streckende und stauchende Richtungen. Weil der Wuerfel
[-1,1]^n beschraenkt und invariant ist, kann die Streckung nicht davonlaufen
-- die Bahn muss zurueckgefaltet werden. Strecken und Falten, wiederholt, ist
das Hufeisen: quer zur stabilen Richtung entsteht eine Cantormenge. Genau die
sieht man im 500-fachen Zoom als Lamination.

DIE PROBE. Fuer einen chaotischen Sattel verknuepft die Beziehung von
Kantz und Grassberger (1985) die drei Groessen:

    alpha  =  kappa / lambda          und   D_Rand = D_Raum - alpha

  kappa   Fluchtrate -- wie schnell Bahnen den Sattel verlassen
  lambda  Lyapunov-Exponent -- wie schnell Nachbarn auseinanderlaufen

Anschaulich: die Streckung erzeugt in jeder Zeiteinheit neue Unsicherheit,
die Flucht raeumt sie ab. Das Verhaeltnis ist die Fraktalitaet. Ist die
Streckung viel schneller als die Flucht, ist alpha klein und der Rand fast
flaechenfuellend.

Gemessen wird lambda nach Benettin: Paare mit Anfangsabstand delta0, gleicher
Schrittweite, regelmaessig auf delta0 zurueckskaliert; lambda ist das Mittel
von ln(Zuwachs) je Zeiteinheit -- gezaehlt nur, solange beide Bahnen noch im
Transienten sind.

VORHERSAGE, VOR DEM LAUF NOTIERT. Fuer diese Instanz wurden alpha = 0,273 und
(fuer alpha_Dichte = 4,6) kappa ungefaehr 0,67 gemessen. Dann muss lambda bei
etwa 0,67/0,273 = 2,5 liegen. Trifft das, ist die Herkunft des Fraktals nicht
nur erzaehlt, sondern geschlossen.
"""
import argparse

import numpy as np

import chaos_sat as C
import gruppen as G


def paarlauf(s0, klauseln, n, delta0=1e-6, eta=0.05, schritte=4000,
             normiere=20, rng=None):
    """Benettin: P Paare, gemeinsame Schrittweite je Paar."""
    lit, sgn, streu = C.bau(klauseln, n)
    M = len(klauseln)
    P = s0.shape[0]
    rng = rng or np.random.default_rng(0)

    d0 = rng.normal(size=(P, n)).astype(np.float32)
    d0 /= np.linalg.norm(d0, axis=1, keepdims=True)
    s = np.repeat(s0.astype(np.float32), 2, axis=0)
    s[1::2] = np.clip(s[1::2] + delta0 * d0, -1, 1)
    a = np.ones((2 * P, M), np.float32)

    summe = np.zeros(P)          # aufsummiertes ln(Zuwachs)
    zeit = np.zeros(P)           # zugehoerige analoge Zeit
    lebt = np.ones(P, bool)

    for k in range(schritte):
        if not lebt.any():
            break
        d1, K1 = C.ableitung(s, a, lit, sgn, streu)
        # Schrittweite aus der Referenzbahn, fuer beide Partner dieselbe
        dtr = np.clip(eta / (1e-6 + np.abs(d1[0::2]).max(axis=1)), 1e-5, 0.5)
        dt = np.repeat(dtr, 2)[:, None]
        d2, K2 = C.ableitung(s + 0.5 * dt * d1, a * np.exp(0.5 * dt * K1), lit, sgn, streu)
        d3, K3 = C.ableitung(s + 0.5 * dt * d2, a * np.exp(0.5 * dt * K2), lit, sgn, streu)
        d4, K4 = C.ableitung(s + dt * d3, a * np.exp(dt * K3), lit, sgn, streu)
        s = np.clip(s + (dt / 6.0) * (d1 + 2 * d2 + 2 * d3 + d4), -1.0, 1.0)
        a = np.minimum(a * np.exp((dt / 6.0) * (K1 + 2 * K2 + 2 * K3 + K4)), 1e8)

        fertig = C.erfuellt(s[0::2], lit, sgn) | C.erfuellt(s[1::2], lit, sgn)
        lebt &= ~fertig
        zeit += np.where(lebt, dtr, 0.0)

        if (k + 1) % normiere == 0:
            diff = s[1::2] - s[0::2]
            d = np.linalg.norm(diff, axis=1)
            gut = lebt & (d > 0)
            summe[gut] += np.log(d[gut] / delta0)
            skal = np.where(d > 0, delta0 / np.maximum(d, 1e-30), 1.0)[:, None]
            s[1::2] = np.clip(s[0::2] + diff * skal, -1, 1)

    brauchbar = zeit > 1.0
    lam = summe[brauchbar] / zeit[brauchbar]
    return lam, zeit[brauchbar]


def lauf(args):
    n, alpha, seed = args.n, args.alpha, args.seed
    kl = G.zufalls_cnf(n, int(round(alpha * n)), 3, np.random.default_rng(1000 * n + seed))
    rng = np.random.default_rng(args.rng)
    print(f"\n  Instanz der Bilder: n = {n}, alpha = {alpha}, seed = {seed}\n")

    # kappa fuer GENAU diese Instanz
    s0 = rng.uniform(-1, 1, (args.punkte, n))
    z, l, f, sch = C.laufe(s0, kl, n, tmax=1e9, eta=args.eta, pruef=1, maxschritt=6000)
    zs = np.sort(z[f])
    p = 1.0 - np.arange(1, len(zs) + 1) / len(f)
    m = (p < 0.5) & (p > 0.02) & (zs > 0)
    k, b = np.polyfit(zs[m], np.log(p[m]), 1)
    kappa = -k
    r = np.log(p[m]) - (k * zs[m] + b)
    r2 = 1 - r.var() / np.log(p[m]).var()
    print(f"  Fluchtrate      kappa  = {kappa:.4f}   (R^2 = {r2:.3f}, "
          f"{len(zs)} von {len(f)} konvergiert)")

    # lambda auf Startpunkten im Randgebiet (lange Transienten)
    schwelle = np.percentile(sch[f], args.perzentil)
    kand = s0[f & (sch >= schwelle)]
    if len(kand) > args.paare:
        kand = kand[rng.choice(len(kand), args.paare, replace=False)]
    lam, zt = paarlauf(kand.astype(np.float32), kl, n, eta=args.eta,
                       schritte=args.schritte, rng=rng)
    L = float(np.median(lam))
    print(f"  Streckrate      lambda = {L:.4f}   (Median ueber {len(lam)} Paare, "
          f"Quartile {np.percentile(lam,25):.3f} / {np.percentile(lam,75):.3f})")
    print(f"\n  Vorhersage      alpha  = kappa/lambda = {kappa/L:.4f}")
    print(f"  Gemessen (Bild) alpha  = 0,273           ->  D_Rand = {2 - kappa/L:.3f} "
          f"gegen gemessene 1,727")


if __name__ == "__main__":
    p = argparse.ArgumentParser()
    p.add_argument("--n", type=int, default=16)
    p.add_argument("--alpha", type=float, default=4.6)
    p.add_argument("--seed", type=int, default=3)
    p.add_argument("--punkte", type=int, default=20000)
    p.add_argument("--paare", type=int, default=500)
    p.add_argument("--schritte", type=int, default=3000)
    p.add_argument("--perzentil", type=float, default=90.0)
    p.add_argument("--eta", type=float, default=0.05)
    p.add_argument("--rng", type=int, default=99)
    lauf(p.parse_args())

python3 lyapunov.py — live output (n=16, α=4.6, seed 3)


  Instanz der Bilder: n = 16, alpha = 4.6, seed = 3

  Fluchtrate      kappa  = 0.1308   (R^2 = 0.937, 20000 von 20000 konvergiert)
  Streckrate      lambda = 0.6126   (Median ueber 500 Paare, Quartile 0.455 / 0.822)

  Vorhersage      alpha  = kappa/lambda = 0.2135
  Gemessen (Bild) alpha  = 0,273           ->  D_Rand = 1.786 gegen gemessene 1,727

A separate sample reproduces the same order of magnitude (α = κ/λ = 0.21 here versus 0.27 in the text) and the same conclusion: κ/λ predicts D_boundary close to the measured 1.727 — the origin of the fractal is closed, not just narrated.

1.71

Boundary dimension

500×

Scale range tested

0.28

Uncertainty exponent

4,700×

Precision per decade

The System

The variables are relaxed: s_i ∈ [−1,1] instead of {−1,+1}. The violation of a clause is a product that vanishes exactly when a literal is satisfied. Every clause carries a weight a_m, which grows exponentially as long as it remains violated.

K_m(s) = 2^−k ∏_i∈m (1 − c_mi s_i) ∈ [0,1] ds_i/dt = 2 ∑_m a_m K_m K_mi c_mi da_m/dt = a_m K_m

The first part is gradient descent on V = ∑ a_m K_m². The second is the trick: it raises exactly the clauses on which the trajectory gets stuck, and thereby makes every local minimum unstable. Consequence — apart from the solutions there are no attractors. The price shows up in the images.

The Basins

A two-dimensional slice through the cube [−1,1]¹⁶: two coordinates are varied over the grid, the remaining fourteen stay fixed. Every pixel is a starting value. Hue shows which of the six solutions was reached, brightness shows how long the path was.

1×

Window: 2.0
Pixel: 2.8·10⁻³
Grid: 720²

Sechs farbige Einzugsgebiete mit einem dunklen gesprenkelten Band dazwischen — The whole slice. The cores of the basins are smooth and large — starting there, the formula is easy. Between them lies a dark, speckled band: long paths, and the color changes from pixel to pixel.

25×

Window: 0.08
Pixel: 1.3·10⁻⁴
Grid: 640²

Ineinandergreifende Zungen dreier Einzugsgebiete mit selbstähnlichen Zacken — Zoomed into the band. What looked like noise from outside is structure: interlocking tongues, each jagged again, and the remaining basins sit in the jags. **The boundary measure rises from 0.062 to 0.116** — for a smooth boundary it would have to fall by a factor of 25.

500×

Window: 0.004
Pixel: 6.3·10⁻⁶
Grid: 640²

Feine cantorartige Lamination aller sechs Einzugsgebiete — Five hundredfold. The boundary is not a line but a lamination: all six basins lie interleaved in ever finer layers, Cantor-like. This magnification reveals the same construction as two levels up.

From the Image to a Number

A point is called ε-uncertain if a perturbation of size ε can move it into a different basin. Its fraction scales as f(ε) ~ ε^α, and in the planar slice D_boundary = 2 − α holds (Grebogi, McDonald, Ott, Yorke 1983). Box counting runs independently of that.

Boundary dimension, three scales, two methods
Slice	Window	α	R²	D from α	D from box counting
Overall	2.0	0.293	0.871	1.707	1.691
Zoom 25×	0.08	0.263	0.999	1.737	1.714
Zoom 500×	0.004	0.273	0.999	1.727	1.776

Two independent methods give the same number, and it holds steady over a scale factor of 500: D ≈ 1.72 in a two-dimensional slice. The boundary is not a curve, it is nearly space-filling.

In practice this means: to reduce the uncertainty in the outcome by one order of magnitude, the starting value must be known 10^1/0.28 ≈ 4,700 times more precisely. Computational precision buys almost nothing.

The Chaos, Measured

A chaotic saddle attracts and repels at the same time: trajectories are captured, wander around, and eventually escape. Its signature is an exponential decay, p(t) ~ e^−κt. The dwell time 1/κ is the compute time of the analog system — hardness becomes a dynamical quantity.

Escape rate over clause density · n = 16, six satisfiable instances per row, 4000 starting values per instance
α	κ	1/κ	R²	Median steps
2.00	0.843	1.19	0.995	19.5
2.50	0.701	1.43	0.991	23.8
3.00	0.795	1.26	0.991	25.5
3.40	0.745	1.34	0.973	32.3
3.80	0.553	1.81	0.959	33.7
4.10	0.213	4.69	0.900	81.2
4.30	0.304	3.28	0.898	47.6
4.60	0.672	1.49	0.938	59.2
5.00	1.076	0.93	0.943	55.2

The decay is cleanly exponential across all densities (R² = 0.90 to 0.995) — that is the proof that a genuine chaotic saddle is present and not just a wide spread. And the curve has a maximum, not a plateau: the dwell time rises to four times its value up to α ≈ 4.1 and falls again afterward. The transition from easy to hard is here not a matter of combinatorial counting, but the growth of an invariant set in phase space.

How the Fractal Arises

The boundary is invariant. Whoever starts exactly on it stays on it: the flow is a homeomorphism, so it maps interior to interior. If the trajectory were ever inside a basin, a whole neighborhood would converge there — traced back, that would mean the starting point was never on the boundary.

An invariant set that never converges therefore lives on the boundary: it attracts along the boundary and repels across it. That is a chaotic saddle. And hence the precise statement: the basin boundary is the stable manifold of this saddle. The images do not show the saddle itself — it has measure zero — but everything that runs toward it before sliding off sideways.

Where the stretching comes from

The weight dynamics da_m/dt = a_mK_m is a positive feedback with no saturation. In a frustrated neighborhood, which clause sets the tone therefore keeps switching — and every switch is a moment where two weighted forces nearly cancel. There, a difference of 10⁻⁶ decides.

Where the folding comes from

The cube is invariant: at the wall s_i = 1 exactly the terms that would push further outward vanish. Exponential divergence can therefore not run away, the trajectories get folded back. Stretching plus folding is the horseshoe — across the stable direction a Cantor set inevitably arises.

Picture a segment crossing the boundary region: it gets stretched, a sub-interval escapes early to one solution — a wide tongue. The rest gets further stretched and folded, a thinner piece escapes later to another — the next, darker generation of tongues. What remains is built like a Cantor set, and that is why every zoom reproduces the same pattern: you only see a later generation of the same construction.

And why 1.72 and not 1.05

Two rates stand against each other. The stretching λ generates uncertainty, the escape κ clears it away. Kantz and Grassberger (1985) link exactly this: α = κ/λ. For the instance behind the images, both were measured:

The balance — measured on the same instance
Quantity	Value	Source
Escape rate κ	0.131	20,000 starting values, exponential decay, R² = 0.94
Stretching rate λ	0.613	500 pairs by Benettin, quartiles 0.46 / 0.82
α = κ/λ	0.214	prediction → D = 1.79
α measured	0.273	from the image → D = 1.73

The relation holds strictly only for two-dimensional maps with a single positive exponent; this system has 16 + 74 = 90 dimensions and a whole spectrum. The spread of λ alone (quartiles 0.46 to 0.82) covers α = 0.16 to 0.29 — the measured value 0.273 falls within it. The origin of the number is thus not proven, but the balance checks out.

What This Does Not Show

n = 16 is small. At this size the satisfiability threshold lies above 4.27 due to boundary effects, and the measured maximum at 4.1 is correspondingly blurred. That κ goes to zero at the threshold does not follow from this data — only that it drops sharply there.
The analog computer is not a win. What is saved in time t is reclaimed by the amplitude: the weights a_m grow exponentially. Hardness changes currency, it does not disappear.
Not a computational artifact. Numerical noise produces uniform speckle, not ordered tongues — and it would bring neither two independent dimension measures into agreement nor yield an exponential decay with R² = 0.99.
The step count is the honest measure. The analog time practically freezes in the transient, because the step size falls as 1/a. For the images the arc length of the trajectory was therefore colored, not t.