Text · original UTF-8 line locators · 14.5 KiB
"""Site-Z dynamics for the open Ising chain with transverse and longitudinal fields.
n <= EXACT_MAX_SITES: exact state-vector propagation. H is read from the instance
`terms`; diagonal (Z-string) terms become an energy vector, single-site X terms are
applied by flipping one axis of the reshaped state. exp(-i H dt) between consecutive
output times is applied with a Chebyshev expansion whose spectral window comes from
a Lanczos estimate of the extremal eigenvalues plus a safety margin.
Larger n: second-order TEBD on an MPS (right-canonical form with singular values,
Hastings' update) with a capped bond dimension. This is an approximation used only
for chains too large for the state vector.
No training, fitting, preprocessing or lookup of reference data is used.
"""
import time
import numpy as np
from scipy.special import jv
PRIORS = {
"family": "exact state-vector Chebyshev propagation (n<=20); bond-capped 2nd-order TEBD MPS fallback (n>20)",
"parameters": {
"exact_max_sites": 20,
"chebyshev_bessel_cutoff": 1e-16,
"norm_guard": "redo window with Gershgorin bounds if |norm-1| > 1e-8",
"threads": 4,
"lanczos_iterations": 30,
"chebyshev_window": 2.0,
"spectral_margin": "5% of half-width + 1.0",
"tebd_max_bond": 64,
"tebd_dt": 0.05,
"tebd_svd_cutoff": 1e-10,
"tebd_trotter_order": 2,
},
"notes": (
"Physics-only method: the Hamiltonian is built from the instance terms (Z/ZZ diagonal, X fields) "
"and the initial product state from instance bits (bit 0 -> Z=+1). For n<=20 the dynamics is exact up "
"to Chebyshev truncation (Bessel coefficients < 1e-16) and floating point; no adaptive switching other than the fixed size "
"threshold, and a wall-clock guard on the TEBD fallback that lowers the bond cap if a run is slow. "
"For n>20 (only the 32-site runtime/sanity run) TEBD with bond dimension <=64 is an uncontrolled "
"approximation at late times in this chaotic chain. No offline training, preprocessing, or use of "
"reference traces."
),
}
EXACT_MAX_SITES = 20
_NTHREADS = 4
def _parse(instance):
n = int(instance["n_sites"])
diag_terms = [] # (coeff, [sites]) for Z strings
x_terms = [] # (coeff, site)
for coeff, ops in instance["terms"]:
ops = {int(k): v for k, v in ops.items()}
paulis = set(ops.values())
if paulis <= {"Z"}:
diag_terms.append((float(coeff), sorted(ops)))
elif len(ops) == 1 and paulis == {"X"}:
x_terms.append((float(coeff), next(iter(ops))))
else:
raise ValueError("unsupported term %r" % (ops,))
init = instance["initial_state"]
if "bits" not in init:
raise ValueError("initial state needs explicit bits")
bits = [int(b) for b in init["bits"]]
return n, diag_terms, x_terms, bits
# ---------------------------------------------------------------- exact
def _zvals(n):
# site i is axis i of the state reshaped to (2,)*n; index 0 -> Z=+1
idx = np.arange(2 ** n, dtype=np.int64)
return [(1 - 2 * ((idx >> (n - 1 - i)) & 1)).astype(np.float64) for i in range(n)]
class _ExactH:
def __init__(self, n, diag_terms, x_terms):
self.n = n
zs = _zvals(n)
d = np.zeros(2 ** n)
for c, sites in diag_terms:
v = np.full(2 ** n, c)
for s in sites:
v = v * zs[s]
d += v
self.diag = d
self.zs = zs
xc = {}
for c, s in x_terms:
xc[s] = xc.get(s, 0.0) + c
self.x = [(s, c) for s, c in sorted(xc.items()) if c != 0.0]
self.xnorm = sum(abs(c) for _, c in self.x)
def apply(self, v, out):
np.multiply(self.diag, v, out=out)
n = self.n
for s, c in self.x:
shp = (2 ** s, 2, 2 ** (n - s - 1))
o = out.reshape(shp)
w = v.reshape(shp)
if c == 1.0:
o += w[:, ::-1, :]
else:
o += c * w[:, ::-1, :]
return out
def _lanczos_bounds(H, iters):
# real start vector: H is real, so the Krylov space stays real
dim = H.diag.size
rng = np.random.default_rng(12345)
q = rng.standard_normal(dim)
q /= np.linalg.norm(q)
q_prev = np.zeros_like(q)
alphas, betas = [], []
beta = 0.0
for _ in range(iters):
w = H.apply(q, np.empty_like(q))
a = float(q @ w)
w -= a * q + beta * q_prev
alphas.append(a)
beta = float(np.linalg.norm(w))
if beta < 1e-12:
break
betas.append(beta)
q_prev, q = q, w / beta
k = len(alphas)
T = np.diag(alphas) + np.diag(betas[: k - 1], 1) + np.diag(betas[: k - 1], -1)
ev = np.linalg.eigvalsh(T)
return ev[0], ev[-1]
def _gersh(H):
return H.diag.min() - H.xnorm, H.diag.max() + H.xnorm
class _Cheb:
"""Chunk-parallel Chebyshev recurrence for Hs = (H - c)/a."""
def __init__(self, H, lo, hi, pool):
n = H.n
self.n, self.dim = n, 2 ** n
self.c = 0.5 * (hi + lo)
self.a = 0.5 * (hi - lo)
self.d2 = 2.0 * (H.diag - self.c) / self.a
groups = {}
for s, coef in H.x:
groups.setdefault(2.0 * coef / self.a, []).append(s)
self.groups = list(groups.items())
self.nb = 2 if n >= 6 else 0
self.nchunk = 2 ** self.nb
self.csz = self.dim // self.nchunk
self.bufs = [np.empty(self.csz, dtype=np.complex128) for _ in range(self.nchunk)]
self.tmps = [np.empty(self.csz, dtype=np.complex128) for _ in range(self.nchunk)]
self.pool = pool
def _chunk(self, k, u1, u0, out, accs, first):
n, nb, csz = self.n, self.nb, self.csz
sl = slice(k * csz, (k + 1) * csz)
o = out[sl]
np.multiply(self.d2[sl], u1[sl], out=o)
f = self.bufs[k]
for xc, sites in self.groups:
started = False
for s in sites:
if s < nb:
src = k ^ (1 << (nb - 1 - s))
w = u1[src * csz:(src + 1) * csz]
fv = f
else:
shp = (2 ** (s - nb), 2, 2 ** (n - s - 1))
w = u1[sl].reshape(shp)[:, ::-1, :]
fv = f.reshape(shp)
if started:
np.add(fv, w, out=fv)
else:
np.copyto(fv, w)
started = True
np.multiply(f, xc, out=f)
np.add(o, f, out=o)
if first:
np.multiply(o, 0.5, out=o)
else:
np.subtract(o, u0[sl], out=o)
t = self.tmps[k]
for acc, coef in accs:
if coef == 0.0:
continue
np.multiply(o, coef, out=t)
a_ = acc[sl]
np.add(a_, t, out=a_)
def rec(self, u1, u0, out, accs, first=False):
if self.pool is None:
for k in range(self.nchunk):
self._chunk(k, u1, u0, out, accs, first)
else:
list(self.pool.map(lambda k: self._chunk(k, u1, u0, out, accs, first), range(self.nchunk)))