Text · original UTF-8 line locators · 15.5 KiB
"""Disordered XXZ chain: site <Z_i>(t) by exact Krylov evolution in the
conserved total-Z sector, with an MPS/TEBD fallback for chains whose
sector is too large to store.
Conventions: bit 0 = spin up = Z eigenvalue +1; site i <-> bit i of the
basis-state integer; Pauli operators; hbar = 1.
No training, no preprocessing, no use of reference traces: every number
comes from the Hamiltonian terms and initial state given in the instance.
"""
import itertools
import math
import time
import numpy as np
PRIORS = {
"family": "xxz-chain-disordered",
"parameters": {
"exact_method": "Lanczos/Krylov propagation of the state restricted to "
"its fixed total-Z sector (exact up to krylov_tol)",
"krylov_dim": 40,
"krylov_tol": 1e-12,
"max_sector_dim": 2_000_000,
"fallback_method": "second-order TEBD (Trotter-Suzuki) on an MPS, used "
"only when the sector exceeds max_sector_dim "
"(e.g. the 32-site runtime/sanity instance)",
"tebd_dt": 0.1,
"tebd_chi_max": 128,
"tebd_svd_cutoff": 1e-10,
"tebd_wall_seconds": 120.0,
},
"notes": (
"Non-adaptive in physics, hybrid in size: (1) Hamiltonian is built "
"from the instance's Pauli terms; because XX+YY, ZZ and Z terms "
"conserve total Z and the initial state is a Z product state, the "
"state is evolved exactly inside one magnetisation sector using a "
"short-iteration Lanczos propagator with a-posteriori error control. "
"(2) For sectors larger than max_sector_dim (the 32-site sanity "
"instance) a truncated-MPS TEBD is used instead; it is approximate "
"(bond-dimension and Trotter error grow with entanglement) and has a "
"wall-clock guard that holds the last computed profile if time runs "
"out. Non-sector-conserving terms fall back to full-Hilbert-space "
"Krylov for small systems. No offline training, fitting, or lookup "
"of reference traces; numpy only (scipy.sparse used if importable)."
),
}
_PAULI = {
"I": np.eye(2, dtype=complex),
"X": np.array([[0, 1], [1, 0]], dtype=complex),
"Y": np.array([[0, -1j], [1j, 0]], dtype=complex),
"Z": np.array([[1, 0], [0, -1]], dtype=complex),
}
def _parse_terms(instance):
n = int(instance["n_sites"])
terms = []
for coef, ops in instance["terms"]:
ops = {int(k): str(v).upper() for k, v in ops.items() if str(v).upper() != "I"}
terms.append((float(coef), ops))
return n, terms
def _initial_bits(instance):
st = instance["initial_state"]
bits = [int(b) for b in st["bits"]]
return bits
def _classify(n, terms):
"""Split into Z fields, ZZ couplings, and XX/YY couplings per pair.
Returns None if any term is outside this set or XX != YY on a pair."""
hz = np.zeros(n)
zz = {}
xx = {}
yy = {}
const = 0.0
for c, ops in terms:
if not ops:
const += c
continue
items = sorted(ops.items())
paulis = "".join(p for _, p in items)
if len(items) == 1 and paulis == "Z":
hz[items[0][0]] += c
elif len(items) == 2 and paulis == "ZZ":
key = (items[0][0], items[1][0])
zz[key] = zz.get(key, 0.0) + c
elif len(items) == 2 and paulis == "XX":
key = (items[0][0], items[1][0])
xx[key] = xx.get(key, 0.0) + c
elif len(items) == 2 and paulis == "YY":
key = (items[0][0], items[1][0])
yy[key] = yy.get(key, 0.0) + c
else:
return None
for key in set(xx) | set(yy):
if abs(xx.get(key, 0.0) - yy.get(key, 0.0)) > 1e-14:
return None
return hz, zz, xx, const
# ---------------------------------------------------------------------------
# sparse operator helpers
# ---------------------------------------------------------------------------
class _Op:
def __init__(self, diag, rows, cols, vals):
self.dim = diag.shape[0]
self.diag = diag
self._sp = None
try:
import scipy.sparse as sp
if rows.size:
self._sp = sp.csr_matrix((vals, (rows, cols)), shape=(self.dim, self.dim))
else:
self._sp = sp.csr_matrix((self.dim, self.dim), dtype=vals.dtype)
except Exception:
self.rows, self.cols, self.vals = rows, cols, vals
def __call__(self, v):
out = self.diag * v
if self._sp is not None:
out = out + self._sp @ v
elif self.rows.size:
w = self.vals * v[self.cols]
out = out + (np.bincount(self.rows, w.real, self.dim)
+ 1j * np.bincount(self.rows, w.imag, self.dim))
return out
def _sector_op(n, cls, bits):
hz, zz, xx, const = cls
k = sum(bits) # number of "1" bits (down spins)
combos = np.array(list(itertools.combinations(range(n), k)), dtype=np.int64)
if k == 0:
states = np.zeros(1, dtype=np.int64)
else:
states = np.sum(np.left_shift(np.int64(1), combos), axis=1)
states = np.sort(states)
dim = states.size
zs = 1 - 2 * ((states[:, None] >> np.arange(n)[None, :]) & 1) # dim x n
diag = np.full(dim, const) + zs @ hz
for (i, j), c in zz.items():
diag = diag + c * zs[:, i] * zs[:, j]
rows, cols, vals = [], [], []
for (i, j), c in xx.items():
differ = np.nonzero(zs[:, i] != zs[:, j])[0]
if differ.size == 0 or c == 0.0:
continue
flipped = states[differ] ^ ((1 << i) | (1 << j))
dst = np.searchsorted(states, flipped)
rows.append(dst)
cols.append(differ)
vals.append(np.full(differ.size, 2.0 * c)) # (XX+YY)|ud> = 2|du>
if rows:
rows = np.concatenate(rows); cols = np.concatenate(cols); vals = np.concatenate(vals)
else:
rows = cols = np.zeros(0, dtype=np.int64); vals = np.zeros(0)
op = _Op(diag.astype(float), rows, cols, vals.astype(complex))
init = np.searchsorted(states, int(sum(b << i for i, b in enumerate(bits))))
return op, zs.astype(float), init
def _full_op(n, terms, bits):
dim = 1 << n
states = np.arange(dim, dtype=np.int64)
zs = 1 - 2 * ((states[:, None] >> np.arange(n)[None, :]) & 1)
diag = np.zeros(dim)
acc = {}
for c, ops in terms:
xmask = 0
phase = np.ones(dim, dtype=complex)
for s, p in ops.items():
if p in ("X", "Y"):
xmask |= 1 << s
if p == "Z":
phase = phase * zs[:, s]
elif p == "Y":
# Y|0> = i|1>, Y|1> = -i|0> (acting on source state)
phase = phase * (1j * zs[:, s])
if xmask == 0:
diag = diag + c * phase.real
else:
acc.setdefault(xmask, np.zeros(dim, dtype=complex))
acc[xmask] = acc[xmask] + c * phase
rows, cols, vals = [], [], []
for xmask, amp in acc.items():
nz = np.nonzero(np.abs(amp) > 0)[0]
rows.append(states[nz] ^ xmask); cols.append(nz); vals.append(amp[nz])
if rows:
rows = np.concatenate(rows); cols = np.concatenate(cols); vals = np.concatenate(vals)
else:
rows = cols = np.zeros(0, dtype=np.int64); vals = np.zeros(0, dtype=complex)
op = _Op(diag.astype(complex), rows, cols, vals)
init = int(sum(b << i for i, b in enumerate(bits)))
return op, zs.astype(float), init