diff --git a/docs/src/assets/mpskit.bib b/docs/src/assets/mpskit.bib index 296f4194b..1c2e16e7a 100644 --- a/docs/src/assets/mpskit.bib +++ b/docs/src/assets/mpskit.bib @@ -29,6 +29,21 @@ @article{capponi2025 abstract = {We investigate the nature of the quantum phase transition in modulated Heisenberg spin chains. In the odd- case, the transition separates a trivial nondegenerate phase to a doubly degenerate gapped chiral symmetry-protected topological (SPT) phase which breaks spontaneously the inversion symmetry. The transition is not an Ising transition associated to the breaking of the inversion symmetry, but is governed by the delocalization of the edge states of the SPT phase. In this respect, a modulated Heisenberg spin chain provides a simple example in one dimension of a non-Landau phase transition which is described by the conformal field theory. We show that the chiral SPT phase exhibits fractionalized spinon excitations, which can be confined by slightly changing the model parameters.} } +@article{ceruti2022, + title = {An Unconventional Robust Integrator for Dynamical Low-Rank Approximation}, + author = {Ceruti, Gianluca and Lubich, Christian}, + year = {2022}, + month = mar, + journal = {BIT Numerical Mathematics}, + volume = {62}, + number = {1}, + pages = {23--44}, + publisher = {Springer}, + doi = {10.1007/s10543-021-00873-0}, + url = {https://doi.org/10.1007/s10543-021-00873-0}, + abstract = {We propose and analyze a numerical integrator for computing the low-rank approximation to solutions of matrix differential equations. The proposed method is based on a variant of the projector-splitting integrator, but here the sub-steps are chosen such that the numerical integrator is robust to the presence of small singular values in the solution.} +} + @article{chepiga2017, title = {Excitation Spectrum and Density Matrix Renormalization Group Iterations}, author = {Chepiga, Natalia and Mila, Fr{\'e}d{\'e}ric}, diff --git a/docs/src/changelog.md b/docs/src/changelog.md index f7cd0741e..3fe4dbfbd 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -21,6 +21,16 @@ When releasing a new version, move the "Unreleased" changes to a new version sec ### Added +- `BUG` time-evolution algorithm: a symmetric second-order Basis-Update & Galerkin integrator for + finite MPS. Unlike `TDVP` it has no backward-in-time substep (stable for imaginary-time evolution), + and passing a truncating `trscheme` enables rank-adaptivity (the bond dimension grows and shrinks + automatically to track entanglement). +- `ParallelBUG` time-evolution algorithm (experimental): the first-order *parallel* Basis-Update & + Galerkin integrator for finite MPS, in which every local problem is solved from the same frozen + snapshot of the state — there is no sweep, so the local integrations are mutually independent. + It is intrinsically rank-adaptive (every bond is augmented and then truncated back down); the + default `notrunc()` restores the pre-step virtual spaces (fixed-rank variant). + ### Changed - `environments` now follows a single positional contract for every state and operator kind: diff --git a/docs/src/man/algorithms.md b/docs/src/man/algorithms.md index 3cf71d52a..694e46fbb 100644 --- a/docs/src/man/algorithms.md +++ b/docs/src/man/algorithms.md @@ -104,6 +104,7 @@ This procedure is commonly referred to as the [`TDVP`](@ref) algorithm, which ag ```@docs; canonical=false TDVP TDVP2 +BUG ``` ### Time evolution MPO diff --git a/research/PARALLELBUG_STATUS.md b/research/PARALLELBUG_STATUS.md new file mode 100644 index 000000000..932fdc44c --- /dev/null +++ b/research/PARALLELBUG_STATUS.md @@ -0,0 +1,98 @@ +# `ParallelBUG` implementation status + +Companion to `PARALLELBUG_design.md`. Records what is implemented and what remains. +**The integrator works** (first-order gates pass); it remains marked experimental because the +local solves are still executed serially and step rejection is not implemented. + +## What is implemented (branch `ld-parallelbug`) + +* `ParallelBUG <: Algorithm` struct + kw-constructor + docstring, registered and exported + (`src/algorithms/timestep/bug.jl`, which holds both BUG integrators). +* A `timestep!` + copying `timestep` implementing Ceruti et al. 2024 (arXiv:2412.00858) Alg. 1–4 + specialized to the caterpillar tree rooted at site `L`, in two phases: + 1. **Frozen-snapshot Galerkin evolutions** (Alg. 2/3): one frozen `ψ₀ = copy(ψ)` + `envs₀`; the + interior amplitude-weighted centers `AC[i]` (the paper's `Y_τ⁰ = U_τ⁰S_τ⁰`) and the root + center `AC[L]` are each evolved forward from that snapshot. These `L` local solves are + mutually independent — the parallel-in-time structure (currently executed serially). + 2. **Leaves→root augmentation** (Alg. 4): at bond `i` the new directions `Ũᵢ` are orthonormalized + against the zero-padded old isometry `[AL⁰ᵢ; 0]` from the evolved center *stacked with the + first-order coupling block* `C̃ᵢ` on the new rows of bond `i-1`. The couplings are one-site + effective derivatives with the **mixed** left environment `⟨Ũ-chain|H|AL⁰-chain⟩` (maintained + with two `TransferMatrix` applications per site) and the frozen old right environment. The + interior site tensors of the augmented state are the isometries `[old │ Ũᵢ]`; the root tensor + is `[C̄_L(t₁); C̃_L]` — the amplitude and all first-order content enter exactly once, at the + root. A final `SvdCut` sweep truncates; `notrunc()` restores the pre-step virtual space of + every bond (per-bond `truncspace`, fixed-rank parallel BUG). +* `test/algorithms/parallelbug.jl` (40 tests, all passing): 2-site dense exactness, energy + + eigenstate-phase conservation, TDVP agreement, imaginary-time monotone lowering + norm, + bond growth under tight/loose `trscheme`, LazySum, **convergence order ≥ 1** (measured slope + ≈ 1.95–1.98 at these bond dimensions), **accuracy improves with ϑ** (vs a ϑ→0 run of the same + integrator, isolating the `c·n·ϑ` term: 5e-4 → 2.5e-6 → 1e-15), U(1)/Z2 charge + graded-structure + preservation. + +## What closed the first-order gap (was the open problem) + +Two changes relative to the first WIP driver: + +1. **The coupling blocks participate in the orthonormalization** (the `M = Û'U₀` reconciliation of + Alg. 4): the earlier driver placed the coupling in the `(new-row, old-col)` block of each + interior site tensor, which routes it through the amplitude-carrying root block and multiplies + it by the old bond matrix (σ-suppressed) — the interior first-order terms effectively vanished + (measured slope ≈ 0). In the correct assembly the interior tensors are *pure isometries* + `[old │ Ũᵢ]` and the coupling data enters `Ũᵢ`'s span via the stacked `Ĉ¹ᵢ = [C̄¹ᵢ; C̃ᵢ]`, so + deep new directions propagate to the root, where the single coupling row + `C̃_L = dt′·⟨Ũ¹_{L-1}|H|ψ₀⟩` captures every first-order component in one exact projection. + (Rank counting is fine: the needed new subspace at bond `b` is the range of the *summed* + tangent components, ≤ r_b-dimensional, not one r-dim family per site.) +2. **Amplitude-weighted kets**: the interior Galerkin solves and couplings act on `AC⁰[i]` + (`= AL⁰[i]·C⁰ᵢ`), not the bare isometry `AL⁰[i]`. With bare-`AL` kets the spans miss the needed + directions (H_eff and the bond matrix do not commute through the parent leg; measured slope + stayed ≈ 0). The amplitude these objects carry is discarded with the R-factor in the + orthonormalization, so no phase/energy overcounting occurs — this resolves the earlier + "full-environment effective Hamiltonian trilemma" (`assembly A/B/C`): the phase only ever + enters the state through the root blocks. + +## Reconciliation with the reference (2026-07-10) + +Compared against the authors' reference `github.com/JonasKu/Publication-Parallel-BUG-for-TTNs` +(MATLAB `Section5.1`, Julia `Section5.2/5.3`). The core assembly matches on every load-bearing +point (frozen-`t₀` independent solves; old-first `[U₀│Ũ]` QR augmentation; first-order coupling with +zeroed new–new corners; amplitude carried once at the root; SVD truncation). Changes landed this +session (all 41 tests green, incl. threaded 4-thread run): + +* **Threading (done)**: phase-1's `L` solves swapped from `map` onto `tmap!` gated on + `Defaults.scheduler[]` (mirrors `tdvp.jl`). Finite envs mutate on lazy recompute, so they are + warmed serially (`_pbug_warmup_envs!`, dispatching through `MultipleEnvironments`) before the + parallel region. +* **Time-dependent coupling (fixed)**: `_pbug_coupling_hamiltonian(...)` was applied with the + one-arg form — a `MethodError` for a genuine `TimedOperator` (no one-arg apply). Now applied at the + midpoint `t + dt/2` (the anchor `integrate` freezes coefficients at). New `TimedOperator` smoke test. +* **Error estimator (done)**: `η = ‖C̃‖/dt` accumulated in `_pbug_assemble` — the coupling blocks + `C̃ᵢ` ARE the frozen derivative projected onto the new directions (Ceruti et al. 2024, eq. 6; + cf. `eta_check.m`, whose `ttm(F_root, {Ũ-projectors})` is exactly this projection). Logged `@debug`. +* **Step rejection (partial, opt-in `maxiter_rejection > 0`)**: the **rank-saturation** trigger + (`rejection_check.m`) — if a truncating step keeps the full doubled space on a bond, recompute as + two half-steps (sub-stepping lowers the per-step normal component so one doubling suffices). + Bounded recursion. Fields `c` (default 10.0) and `maxiter_rejection` (default 0) added. +* **Type stability / cleanups (done)**: `Vector{Any}`→concrete `As`/`Cevo`. + +## Remaining work (why still experimental) + +1. **η-threshold rejection trigger**: `η` is computed but only the rank-saturation trigger drives + the retry. Wiring `dt·η > c·ϑ` needs a scheme-generic way to extract the tolerance `ϑ` from an + arbitrary `trscheme` (trivial for `truncerror`, unclear in general). The current retry also uses + sub-stepping rather than the reference's **basis-enrichment** recompute (re-embed `ψ₀` in the + augmented 2r basis, re-assemble → 4r); enrichment keeps the same `dt` and is the true + rank-adaptive mechanism, but needs the matrix-paper [4] retry semantics + numerical validation + (the caterpillar has no dropped corner until `L ≥ 3`, so there is no cheap 2-site oracle). +2. **Second-order variant** (`TTN_integrator_parallel_2nd_order_nonglobal.m` + `RK_2`): the reference + is **NOT** a Strang composition — it is an **augmented-Galerkin** scheme that *keeps* the new–new + coupling corners (a `3r` core: old block, new–old coupling `C̄`, and the new–new block `Cᵢ`), + solving the enlarged core ODE. For MPS this means retaining (not discarding via the R-factor) the + new–new corner amplitude at interior nodes — effectively a parallel analogue of the sequential + `BUG`'s second order. Substantial; recommend as a focused follow-up. Same subtle corner + machinery as (1)'s enrichment, so validate them together against an `L ≥ 3` reference. +3. Minor: `timestep!`'s `envs` argument is still only used for `L == 1`; the frozen-snapshot + environments are recomputed each step. Reusing the caller's `envs` is blocked by the frozen-copy + design (`ψ₀ = copy(ψ)` ⇒ identity-keyed envs cannot match); evolving from `ψ` directly (it is not + mutated before the final overwrite) would enable reuse but trades away the copy's isolation. diff --git a/research/PARALLELBUG_design.md b/research/PARALLELBUG_design.md new file mode 100644 index 000000000..014ce33ec --- /dev/null +++ b/research/PARALLELBUG_design.md @@ -0,0 +1,225 @@ +# Design: `ParallelBUG` integrator for `FiniteMPS` in MPSKit + +Stage 3 of the BUG family (see `BUG_finiteMPS_design.md` §7, which left this "not yet pinned +down"). This note **derives** the concrete MPS-level recipe from the two 2024 parallel papers and +records the design decisions for a `ParallelBUG <: Algorithm`. + +Written on branch `ld-parallelbug` (worktree), branched off `ld-bug` (which carries the sequential +fixed-rank + rank-adaptive `BUG`). + +## 1. Source algorithm (transcribed from the papers) + +References in `research/`: +* Ceruti, Kusch, Lubich 2024, *A Parallel Rank-Adaptive Integrator for DLRA*, SISC 46(3) — the **matrix** case (the block formulas below). +* Ceruti, Kusch, Lubich, Sulz 2024, *A parallel BUG Integrator for Tree Tensor Networks*, arXiv:2412.00858 — the **tree** generalization (Alg. 1–4). **Has no MPS section**; the caterpillar specialization here is our derivation. + +### 1.1 Matrix parallel BUG (the load-bearing block structure) + +Start from `Y₀ = U₀ S₀ V₀*`, rank `r`. One step `t₀ → t₁ = t₀ + h`. **Three independent forward +ODEs, all from `t₀` data, all mutually parallel:** + +* **K-step** (evolve the left factor + core, right factor frozen): + `K̇ = F(t, K V₀*) V₀`, `K(t₀) = U₀ S₀`. Then `[U₀ | K(t₁)]` QR → `Û = [U₀ | Ũ₁]` (old-first). +* **L-step** (evolve the right factor + core, left factor frozen): + `L̇ = F(t, U₀ L*)* U₀`, `L(t₀) = V₀ S₀*`. Then `[V₀ | L(t₁)]` QR → `V̂ = [V₀ | Ṽ₁]`. +* **S-step** (evolve the core, BOTH factors frozen — this is the *parallel* twist; sequential BUG + would use the augmented bases here, forcing a `2r×2r` solve after K/L): + `Ṡ̄ = U₀* F(t, U₀ S̄ V₀*) V₀`, `S̄(t₀) = S₀`. Size stays `r×r`. + +**Augmented core (eq. 3.4)** — assembled *algebraically*, no extra ODE: +``` + old V₀ new Ṽ₁ + ┌───────────┬──────────┐ +old U₀│ S̄(t₁) │ S̃ᴸ │ S̃ᴷ = Ũ₁* K(t₁) (K coupling, bottom-left) + ├───────────┼──────────┤ S̃ᴸ = L(t₁)* Ṽ₁ (L coupling, top-right) +new Ũ₁│ S̃ᴷ │ 0 │ new–new corner = 0 ⇒ FIRST ORDER + └───────────┴──────────┘ +``` +Discarding the `0` corner costs `O(h² + hε)` local error — this is exactly why parallel BUG is +**first order** (sequential BUG keeps that block and is second order). Truncate `Ŝ₁` by SVD to +tolerance `ϑ` → new `U₁, V₁, S₁`. + +**Step rejection (§3.3):** normal-component estimator `η = ‖Ũ₁* F(t₀, Y₀) Ṽ₁‖` (= the discarded +corner / `h`). Reject & recompute the step on the *augmented* bases `Û, V̂` if either +(a) `r₁ = 2r` (truncation saturated the doubling cap), or (b) `h·η > c·ϑ` (`c ≈ 10`). + +### 1.2 Tree/MPS generalization + +Per node `τ` with children `τᵢ`: Galerkin-evolve the connecting tensor with **all surrounding bases +frozen at `t₀`** (Alg. 3); build each bond's new directions `Ũ¹_{τᵢ}` old-first (Alg. 4A); assemble +the augmented connecting tensor by **stacking coupling blocks mode-by-mode** (Alg. 4B): +``` +Ĉ ← C̄¹; for each child mode i: Ĉ ← Tenᵢ( [ Matᵢ(Ĉ) ; Matᵢ(C̃ᵢ) ] ) +C̃ᵢ = h · F_τ(Y₀) ×_{j≠i} U₀*_{τⱼ} ×ᵢ Ũ¹*_{τᵢ} (all multi-new corners left at 0) +``` +Robust first-order global bound (Thm 4.5): `‖Yₖ − A(tₖ)‖ ≤ c₁h + c₂ε + c₃δ + c₄·k·ϑ`, constants +independent of the bond singular values. Truncation term `c₄·k·ϑ` ⇒ scale `ϑ ∝ h`. + +## 2. MPS specialization (caterpillar rooted at site L) + +`FiniteMPS` = linear tree, rooted at site `L` (left-canonical: `AL[1..L-1]`, center at `L`). Physical +legs are **uncompressed leaves** (`U = 𝟙`), so they are never augmented — **only the `L-1` virtual +bonds grow.** Treat each **bond `b`** (between sites `b`, `b+1`) as the matrix-DLRA object: + +| matrix object | MPS object at bond `b` | MPSKit primitive | +|---|---|---| +| `U₀` (left factor) | left block `AL[1..b]`, locally `AL[b]` | — | +| `V₀` (right factor) | right block `AR[b+1..L]`, locally `AR[b+1]` | — | +| `S₀` (core) | bond tensor `C[b]` | — | +| K-step `K(t₁)` | evolve site `b`: `AC_hamiltonian(b)` on `AC[b]`, frozen `t₀` | `integrate` | +| L-step `L(t₁)` | evolve site `b+1`: `AC_hamiltonian(b+1)` on `AC[b+1]`, frozen `t₀` | `integrate` | +| S-step `S̄(t₁)` | evolve bond `b`: `C_hamiltonian(b)` on `C[b]`, frozen `t₀` | `integrate` | +| `Ũ₁` (new left dirs) | `_bug_augment_left(AL[b], K(t₁))` (already in `bug.jl`) | `left_null`/`catdomain` | +| `Ṽ₁` (new right dirs) | `_bug_augment_right(AR[b+1], L(t₁))` (already in `bug.jl`) | `right_null`/`catcodomain` | + +This reuses **exactly** the effective operators (`AC_hamiltonian`, `C_hamiltonian`, `integrate`) and +augmentation helpers (`_bug_augment_left/right`) already present — the only genuinely new pieces are +the frozen-`t₀` scheduling, the `2×2` core assembly, and step rejection. + +**All local ODEs read from one frozen snapshot** `(ψ₀ = copy(ψ), envs₀ = environments(ψ₀, H, ψ₀))`, +so every `AC_hamiltonian(i)`/`C_hamiltonian(b)` is well-defined simultaneously (identical to how the +`InfiniteMPS` `timestep` in `tdvp.jl` evolves all `AC`/`C` from one frozen `envs`). This is the +parallelizable structure — reuse that file's `@sync`/`Threads.@spawn` + `tmap!` gated on +`Defaults.scheduler[]`. + +## 2b. CORRECTION (implementation finding, must supersede §2's per-bond picture) + +The per-bond matrix-DLRA mapping in §2 (an S-step via `C_hamiltonian(b)` on **every** bond) is +**wrong for MPS** and was falsified in implementation: + +* **Phase overcounting.** MPSKit's `C_hamiltonian(b)` is the *full-environment* effective + Hamiltonian, so each bond's S-step carries the full energy phase `exp(-iE·dt)`. With a core on + every one of the `L-1` bonds these phases **compose**, giving `exp(-i(L-1)E·dt)` — the eigenstate + phase is overcounted by exactly `(L-1)×` (measured: ratio 2.0 at `L=3`, 3.0 at `L=4`). There is no + backward substep to cancel it (unlike TDVP). +* **No bond growth.** Augmenting from the *frozen-evolved `AC[i]`* cannot grow a bond, because + `range(AC[i]) = range(AL[i])` up to the evolution, and worse, `_bug_augment_left(AL[i], AC[i])` + extracts `g = N'·AC` which is ~0. + +**The correct MPS mapping is the literal tree recursion, not L−1 independent matrix problems:** + +* The caterpillar rooted at `L` is **one** tree. Only the **root** connecting tensor (`AC[L]`) + carries amplitude/phase in the *assembled* state; every interior site tensor of the augmented + state is an **isometry**. So there is **no per-bond `C_hamiltonian` S-step at all** — the + amplitude enters the assembled state **once**, at the root. +* First order comes from the **explicit-Euler coupling blocks with zeroed multi-new corners** + (Alg. 4B), assembled by the leaves→root recursion — *not* from transporting fully-evolved tensors + (that accidentally recovers 2nd order, as observed). + +### 2c. SECOND CORRECTION (supersedes parts of §2b; this is what landed) + +§2b's further conclusion that the interior nodes must Galerkin-evolve the bare **isometry** +`AL[i]` was **also falsified**: the paper's subflows act on the amplitude-weighted subtree objects +`Y_τ⁰ = U_τ⁰S_τ⁰` (locally `AC[i] = AL[i]·C[i]`, cf. Alg. 2's leaf K-step `K(t₀) = U⁰S⁰`), and +with bare-`AL` kets the augmented spans miss the needed first-order directions (the effective +Hamiltonian does not commute with the bond matrix on the parent leg; measured slope stayed ≈ 0). +The phase/energy overcounting §2b worried about does not occur because the evolved interior +tensors and the coupling blocks only contribute their **range** — the amplitude/phase they carry +is discarded with the R-factor in the leaves→root orthonormalization (`Q̂ᵢ = orth([Ĉ⁰ᵢ | Ĉ¹ᵢ])`, +old-first), which is precisely Alg. 4's `M = Û'U₀` reconciliation. The couplings must be stacked +into `Ĉ¹ᵢ = [C̄¹ᵢ; C̃ᵢ]` *before* orthonormalizing, so that deep new directions propagate to the +root; the root tensor `[C̄_L(t₁); C̃_L]` then carries all amplitude and first-order content. +See `PARALLELBUG_STATUS.md` for the implemented recipe and the numerical gates. + +## 2d. SECOND-ORDER variant (from the reference `TTN_integrator_parallel_2nd_order_nonglobal.m`) + +The reference achieves second order **not** by Strang composition but by an **augmented-Galerkin** +scheme that *keeps the new–new coupling corner* the first-order method drops. Exact recipe (verbatim +from the MATLAB, per node with `m` children): + +1. **Pre-step**: build the augmented child bases `U0_hat{jj}` (same QR construction as first order) + and overlap matrices `Mi{jj} = Q0ᵢ' (⊗_{i≠jj} Y0ᵢ'U0_hat_i) Q0ᵢ`. +2. **Subflow Φᵢ** (per child, parallel): K-step evolves `Y0ᵢ·Mi{i}` (leaf: `RK_4`; subtree: recurse); + augment the child basis old-first `[Y0{i} Y1_i]` via QR, then **zero-pad to `3·r`** and take + `U_tilde{i} = ` columns `r+1:end` (the ≤ `2r` new directions). +3. **Subflow Ψ** (core): `C0 = ttm(Y0{end}, Mi, 1:m)`; evolve the core `C1_bar = + RK_4_tensor(C0, func_ODE, U0_hat, …)` where `func_ODE` = `F` projected onto the **old** bases `U0` + (a Galerkin K-step of the core). +4. **Augmentation** (the second-order piece): `F{end} ← (t1-t0)·F{end}`; for each mode `jj`, + `Ci = ttm(F{end}, dum, 1:m)` with `dum{i} = Mat0Mat0(U_tilde{i}, F{i})` for `i=jj` and + `Mat0Mat0(U0_hat{i}, F{i})` otherwise — i.e. `Δt·F` projected onto the **new** directions on mode + `jj` and the augmented bases elsewhere. This `Ci` is placed into rows `r+1:3r` of the core: + `tmp(rr(jj)+1:3*rr(jj), :) = tenmat(Ci,jj,vv)`, `ss(jj)=3*ss(jj)`. **The block the first-order + method leaves at 0 is now filled by `Ci`.** (`Ci` is exactly the quantity whose norm is the + estimator `η`; keeping it is what raises the order.) +5. Truncate the `3r` bonds back down (elsewhere). + +**MPS caterpillar mapping (`m = 1` per node → the loops collapse):** each interior connecting tensor +`AC[i]` gets a **`3r`** child bond — old `r`, plus the ≤`2r` new directions from `[old | K-step | Ci]` +— with the coupling `Ci = dt·⟨Ũᵢ₋₁|H|ψ₀⟩` **kept in the core** instead of only informing `Ũ`'s range. +Concretely this means the interior tensors are no longer pure isometries: the new–new block carries +amplitude, so the assembled state must be re-gauged (`FiniteMPS` from the raw tensors handles this). + +**Implementation status / risk**: this is a genuine cross-framework port — the reference is a +Tucker-tree (cell arrays, `tenmat`/`ttm`, explicit `Mi`/`Q0_i`), MPSKit is `TensorMap` MPS. Validation +has no cheap oracle (the caterpillar 2-site step is exact; dropped corners appear only at `L ≥ 3`), +so it must be validated behaviourally by the **local one-step error order**: `‖step − exact‖ ∝ dt²` +for first order, `∝ dt³` for a genuine second order. + +**2026-07-10 attempt + NEGATIVE RESULT (not landed).** Implemented the tractable route — build the +first-order `2r` augmented basis `Û = [old | Ũ]`, then a Galerkin evolution of the *root* center on +the augmented state's own environments `⟨Û|H|Û⟩` (embedding `ψ₀` via a chain-transport `T`, evolving +with `AC_hamiltonian` from `t₀`). It is **exact at 2 sites** and consistently **~2× more accurate** +than first order in a rank-limited regime — but the measured **local slope is ≈ 2, not 3** (o1: +1.98/1.99/1.97; o2: 1.97/1.96/1.88 at `dt = 0.04→0.005`). So it is only a better-constant *first*-order +scheme, **not** second order, and was reverted rather than shipped mislabeled. Diagnosis: a Galerkin +core on the *frozen* first-order `2r` basis caps the local error at `O(dt²)` — the basis misses the +`O(dt²)` corner directions. Genuine second order needs the reference's **`3r` augmentation that adds +the `Ci` (new–new corner) directions to the basis and keeps their amplitude** (step 4 above), i.e. a +second family of new directions beyond `Ũ`, re-gauged. That is the remaining work for a true +second-order `ParallelBUG`; the sequential [`BUG`] already provides a second-order option meanwhile. + +## 3. Key design decisions (recorded) + +1. **New name, not a flag on `BUG`.** A separate `ParallelBUG <: Algorithm` struct rather than + `BUG(; parallel=true)`. Rationale: different order (1st vs 2nd), different control flow (no sweep, + frozen envs, step rejection), and different field set (needs `c`, max rejections). Mirrors how + `TDVP`/`TDVP2` are distinct structs. Keeps the well-tested `BUG` byte-for-byte. +2. **Always augment-then-truncate.** Unlike sequential `BUG`, the parallel integrator is + *intrinsically* rank-adaptive: it doubles every bond then truncates. `notrunc()` therefore means + "truncate back to the pre-step per-bond dimensions" (fixed-rank parallel BUG); any + `truncerror`/`truncrank` gives genuine rank adaptivity. Reuse `SvdCut` for the truncation sweep, + exactly as `_bug_truncate!` does. +3. **Reuse `_bug_augment_left/right`.** The old-first `[U₀ | Ũ₁]` construction with the `[𝟙;0]` + overlap is already implemented, tested (`bug_augment.jl`), and sector-correct. The coupling blocks + `S̃ᴷ = Ũ₁* K(t₁)`, `S̃ᴸ = L(t₁)* Ṽ₁` are contractions against the returned `Ũ₁`/`Ṽ₁`. +4. **Frozen-`t₀` snapshot + lazy env self-heal.** One `copy(ψ)` + its `environments`. After + installing the truncated result, `FiniteEnvironments` recomputes lazily (no `recalculate!` exists + for finite envs — confirmed). No old/new env bookkeeping (that was sequential BUG's H8 subtlety); + here every read is from `t₀`, which is *simpler* than sequential BUG. +5. **Correctness gate = behavioral, mirroring `BUG` tests.** Rather than white-box block matching, the + MPS assembly is validated against the same oracles the sequential `BUG` uses: (a) overlap vs a + dense `exp(-iH·T)` reference, (b) agreement with `TDVP` to `O(dt)`, (c) **first-order** log–log + slope ≈ 1 (vs `BUG`'s 2), (d) imaginary-time monotone energy lowering + norm, (e) symmetric-tensor + sector/charge preservation. This is what "verify correctness similar to how the regular BUG is + tested" means, and it makes the assembly detail self-correcting. +6. **Step rejection: opt-in, capped.** `maxiter_rejection` (default e.g. 3) bounds recomputes; `c` + (default 10, the paper's value) sets the `h·η > c·ϑ` threshold. With `notrunc()` (fixed rank) the + rank-saturation trigger is disabled (we deliberately cut back to the old rank). Rejection is a + genuine-adaptivity feature; keep it simple and cheap (`η` from one frozen `F₀` apply). + +## 4. Build plan (small steps, each validated) + +1. **Struct + registration.** `ParallelBUG` mirroring `BUG` fields (+`c`, `maxiter_rejection`); kw + constructor; `include` + export. No logic yet. Sanity: `ParallelBUG()` constructs. +2. **Serial driver, no rejection.** `timestep!` (+ copying `timestep`): frozen snapshot → per-site + `AC_hamiltonian` + per-bond `C_hamiltonian` local solves → per-bond `2×2` augmented assembly → + `SvdCut` truncation. Serial scheduler first (correctness before parallelism). + Gate: energy conservation on a ground state; first-order slope; TDVP agreement; dense overlap. +3. **Imaginary time + norm.** Renormalize after truncation; monotone energy-lowering test. +4. **Symmetric tensors.** U(1)/Z2/SU(2) sector + charge preservation, dynamic grading (reuse + `bug.jl`'s symmetric test scaffold). +5. **Step rejection.** `η` estimator + recompute loop; a test that a rank-1 start grows past a single + doubling within one step via rejection. +6. **Parallelism.** Swap the per-site/per-bond loops onto `@sync`/`Threads.@spawn` + `tmap!` gated on + `Defaults.scheduler[]` (mirror `tdvp.jl`'s `InfiniteMPS` branch). Gate: identical results to the + serial path (scheduler must not change the answer). + +## 5. Risk register (inherits `BUG_finiteMPS_design.md` §9 H1–H10, plus) + +| # | Risk | Mitigation | +|---|---|---| +| P1 | MPS assembly of the `2×2` block across *both* bonds of each site is subtle (each bond has one new space, from its left neighbor); easy to double-count the core (cf. §3 sequential warning) | Behavioral gate (§3.5): first-order slope + dense overlap will expose a wrong assembly; build the 2-site case first where it reduces to the exact matrix formulas | +| P2 | New-new corner must be *exactly* zero (first-order); a stray coupling makes it inconsistent | Assemble with explicit `zerovector!` block; assert bond dims double before truncation | +| P3 | Thread safety of shared frozen graded snapshot (H9) | Genuine `copy`; read-only `envs₀`; follow `tdvp.jl` infinite threading structure; no shared per-sector buffers | +| P4 | `ϑ` accumulates as `c₄·k·ϑ`; users expect TDVP-like accuracy | Document `ϑ ∝ h`; note first-order in docstring like `BUG`'s note | diff --git a/research/PARALLEL_BUG_second_order_status.md b/research/PARALLEL_BUG_second_order_status.md new file mode 100644 index 000000000..ddeea1265 --- /dev/null +++ b/research/PARALLEL_BUG_second_order_status.md @@ -0,0 +1,186 @@ +# Second-order parallel BUG (`ParallelBUG2`) — status and diagnosis + +## RESOLVED (2026-07-11): genuine second order — faithful Kusch Variant 2 + +`ParallelBUG2` clears the decisive local-order gate with a **uniform** single-step slope **≈ 3** in `dt` +(vs the first-order `ParallelBUG`'s slope ≈ 2). Measured with the shared harness +(`test/algorithms/parallelbug.jl`, TFIsing, `dts = [0.04, 0.02, 0.01, 0.005]`, `trscheme = truncerror(atol=1e-12)`): + +| L | bond | slopes | note | +|---|---|---|---| +| 2 | ℙ^4 | (exact) | machine-exact, no dropped corner | +| 3 | ℙ^8 | (exact) | full-rank enriched basis completes the block | +| 4 | ℙ^8 | ≈ 3.0 | genuine second order | +| 5 | ℙ^8 | ≈ 3.0 | genuine second order | +| 6 | ℙ^4 | ≈ 3.0 | unsaturated interior bond (a real dropped corner) | +| 7 | ℙ^4 | ≈ 3.0 | stays second order as the chain grows | + +The construction implemented in `src/algorithms/timestep/bug.jl` (`_pbug2_assemble_core` and helpers): + +1. **Enrich** both the left (`Û0`) and right (`V̂0`) bond bases to rank `2r` with one frozen `H·ψ₀` + application, old-first, with the mixed-env coupling propagating deep directions to the root; fold the + enriched environments `GLhat`/`GRhat` by **explicit transfer matrices** (no `FiniteMPS` round-trip, so + the zero-weight enriched directions never collapse — the route-1 blocker below). +2. **K-step** every center on the enriched environments, **freezing the enriched right basis `V̂0`** + (freezing the *old* right basis gives slope 2 — see route 2). Threaded like the first-order phase-1. +3. **Assemble** leaves→root: interior tensors are the `4r` isometries `[Û0 | Ũ2]` (new directions from the + K-evolved centers), and the **evolved-amplitude coupling `R = Ũ2ᵀĈ`** is transported one site at a time + through the frozen right isometry `V̂0`. This is the ingredient the earlier attempts missed: the + first-order *frozen-derivative* coupling `h·F(Y0)` projected onto `Ũ2` is ≈ 0 (because `Ũ2 ⟂ Û0` and + `Û0` already contains the `H·ψ₀` direction), so it transports no second-order content and the interior + bonds stay first order (slope 2 at `L ≥ 6`). The evolved amplitude `Ũ2ᵀĈ` is the genuine `O(h²)` term, + the caterpillar analogue of the matrix recipe's `Ũ2ᵀK(t1)` block. +4. **Root** core = the exact matrix `Ŝ1` assembly: the rank-`2r` Galerkin `S̄1 = Kevo[L]` in the `Û0` rows + stacked with the transported coupling in the `Ũ2` rows, and the **new–new corner ZERO** (at the root the + right bond is trivial, so the `LᵀṼ2` block vanishes too). This gives the published local error `O(h³)`. + +This matches Kusch (2024) Variant 2 (`4r`, old-first bases, no predicted basis / no transfer matrix `Mi`), +specialized to the caterpillar rooted at the last site. The faithfulness signature is the **uniform** slope 3 +(an earlier hybrid that replaced the root `Ŝ1` with a full `4r` Galerkin over-resolved to slope ≈ 4 at +small L — more accurate, but not the published structure; the exact zeroed-corner assembly removes it). +It clears the full shared verification (2-site exactness, energy/phase, TDVP agreement, imaginary-time +lowering, Z2 charge preservation, serial=threaded). Only plain Hamiltonians are supported so far (no +`LazySum`/time-dependent operators). It is NOT the reference-MATLAB path (that is Variant 1 with predicted +bases `U0_hat` + transfer matrices `Mi` + `3r`); the simpler Variant 2 was implemented as the handoff +recommended. `ParallelBUG2` is still marked experimental. + +The diagnosis of the earlier dead ends is retained below for the record. + +--- + +## Original status (superseded): the dead ends + +One-sentence status (historical): the **first-order** `ParallelBUG` is audited and passes the decisive +local-order slope-≈2 gate at `L = 3, 4`; a genuine **second-order** `ParallelBUG2` is **not** achieved — the +tractable implementation routes measure local slope ≈ 2 (still first order) or break, and the correct +route is a substantial faithful port of the reference tree-tensor-network algorithm that remains to be +done. `ParallelBUG2` is committed as an experimental struct that is exact for `L ≤ 2` and `throw`s for +`L ≥ 3`, so it can never be silently mislabeled as second order. + +This note records what was tried, the measured slopes, why each route fails, and the precise remaining +work, so the port can be resumed without re-deriving the dead ends. + +## The decisive gate + +Single-step local error `‖step(dt) − exp(−iH·dt)ψ₀‖` vs `dt`, log–log slope, full bonds (`ℙ^8`), +transverse-field Ising, `dts = [0.04, 0.02, 0.01, 0.005]`: + +- **First order ⇒ slope ≈ 2**; **genuine second order ⇒ slope ≈ 3**. +- `L = 2` is exact (no dropped corner) and does **not** distinguish the orders — the gate needs `L ≥ 3`. +- Empirically `ParallelBUG` (first order) gives clean slope ≈ 2.0 at `L = 3, 4, 5, 6`. This is the + committed gate in `test/algorithms/parallelbug.jl` (`@testset "local order slope L=$Lc"`). + +## What was tried for second order, and the measured slope + +The correct *mechanism* (Kusch 2024): pre-augment each bond basis with one `H·ψ₀` application **before** +evolving, run the Galerkin center at the enlarged rank on that basis, assemble with the new–new corner +kept zero (now `O(h³)`). Three MPS realizations were attempted: + +1. **Reuse the first-order assembly on an `H·ψ₀`-enriched `FiniteMPS` snapshot `ψ̂0`** (the approved + plan's central idea). **Structurally impossible.** The enriched directions carry *zero weight* + (`ψ̂0 == ψ₀` as a state — verified overlap 1.0), so `FiniteMPS` canonicalization *collapses* them: + the stored `AL[2]` keeps the enriched space (`ℙ^4`) but the mixed-gauge `AC[2]` collapses to `ℙ^2`, + and the first-order `_pbug_assemble_core` then hits a `SpaceMismatch` at `L ≥ 3`. A zero-weight + enriched basis cannot survive a canonical-form round-trip — full stop. + +2. **Explicit enriched-left environment + rank-2r root Galerkin, no 4r augmentation** (materialize a + `FiniteMPS` only from the final, nonzero-weight tensors — avoids the collapse). **Runs and is + stable, but slope ≈ 2.0 at `L = 4`** (`L = 2, 3` come out machine-exact only because the enriched + basis happens to complete the small left block). Enriching the left basis + a root Galerkin *alone* + is a better-constant first-order method — the second-order correction also needs the interior K-step + directions augmented to `4r`/`3r` and the enriched *right* basis frozen in the K-step. + +3. **Predicted-basis Galerkin (`ψhat = compress(H·ψ₀)` to rank r) + transfer matrices + old-first + augmentation.** **Broken: slope ≈ 0** (O(1) error ≈ 0.15, independent of `dt`). Evolving `ψ₀`'s + center on `ψhat`'s environments and transferring back with an ad-hoc overlap matrix does *not* + reproduce the reference's prolongation/restriction: the basis reconciliation between the + predicted-basis core `C1_bar` and the old-first augmented basis is not captured, so the assembled + state is O(1) wrong. + +Route 2's slope ≈ 2 reproduces the prior (reverted, never-committed) attempt's result; route 1 is a new +structural finding; route 3 shows an ad-hoc predicted-basis handling is worse, not better. + +## Why it is hard: the mixed-basis / basis-bookkeeping problem + +The reference `TTN_integrator_parallel_2nd_order_nonglobal.m` (Variant 1, `3r`) does **not** enrich a +single shared bond space. Per node it builds, per child: + +- a rank-`r` **predicted basis** `U0_hat` = orthonormalized range of the `H·ψ₀` image *through the old + environment* `Q0_i` (line 38), **not** the old-first `[old | new]` 2r basis, and **not** a global + SVD-compression of `H·ψ₀`; +- a **transfer matrix** `Mi = ⟨old | U0_hat⟩` (line 49) that maps the old core into the predicted basis; +- the core is evolved from `C0 = ttm(core, Mi)` on the `U0_hat` bases (`func_ODE`), giving `C1_bar` in + the *predicted* basis; +- the augmented basis is `[old | K-evolved]` zero-padded to `3r`, and `C1_bar` is placed in rows + `1:r` with the coupling `Ci` in rows `r+1:3r` (lines 135–164). + +The subtlety that defeats the shortcuts: `C1_bar` lives in the `U0_hat` (predicted) basis while the +augmented columns are the `[old | K-new]` basis; the two are reconciled *implicitly* through `Mi` and +the specific `3r` block placement. Getting local slope 3 requires reproducing this reconciliation +exactly — approximating `U0_hat` (route 3) or skipping the predicted basis (route 2) loses it. + +The left/right asymmetry compounds this: the matrix Variant-2 recipe uses *separate* enriched left +(`Û0`) and right (`V̂0`) bases whose new subspaces differ, so the interior K-step center is a genuinely +*mixed* (rectangular) object — exactly what prolongation/restriction manage and what a naive shared-bond +MPS representation cannot. + +## Remaining work (to resume the port) + +1. **Faithfully port the reference recursion** `TTN_integrator_parallel_2nd_order_nonglobal.m` to the + `FiniteMPS` caterpillar (root at `L`), reusing MPSKit's effective-Hamiltonian environments as the + prolongation/restriction. Build the rank-`r` predicted basis `U0_hat` per bond as the orthonormalized + `H·ψ₀` image *through the old right environment* (not a global compression), the transfer `Mi`, evolve + the core on `U0_hat` from `C0 = ttm(core, Mi)`, and assemble the `3r` core with `C1_bar` in the + old-block rows and the coupling in rows `r+1:3r`, new–new corner zero. Re-gauge only the final + (nonzero-weight) tensors into a `FiniteMPS`. +2. **Gate at `L ≥ 3` (use `L = 4, 5`; `L = 3` full-rank can be exact and masks the corner)** on local + slope ≈ 3 before labeling it second order. Mirror the `ParallelBUG` slope harness already in + `test/algorithms/parallelbug.jl`. +3. Cross-check against a MATLAB trace of the reference on a 3-site chain to pin down the `Mi` / `3r` + basis reconciliation numerically — this is where the shortcuts went wrong and where a from-scratch + derivation is most error-prone. + +## Faithful-port attempt (Variant 2, `4r`) — how far it got and the real blocker + +A second, deeper attempt targeted **Variant 2** (`4r`, [K24] §5.3) specifically because it uses +**old-first enriched bases** `Û0 = orth([U0, F0·V0])`, `V̂0 = orth([V0, F0ᵀ·U0])` and **no predicted +basis / no transfer matrix `Mi`** — so it sidesteps the predicted↔old reconciliation that broke route 3 +above. Progress and findings: + +- **Working building blocks** (verified to run and produce sensible spaces): + - old-first **left** enrichment `Û0[i]` (2r) via the existing `_pbug_newdirs`/`_pbug_stack_child` + machinery fed the `H·ψ₀` image (this is the retained `_pbug_preaugment`); + - old-first **right** enrichment `V̂0[i]` (2r) via `_bug_augment_right` fed the `H·ψ₀` image; + - enriched left/right **environments** by explicit `TransferMatrix` folding through `Û0` / `V̂0` + (`GLhat[i]`, `GRhat[i]`) — no `FiniteMPS` round-trip, so no zero-weight collapse; + - **old→enriched center embedding** via `isometry(enriched_bond ← old_bond)` (old-first ⇒ the + canonical `[I;0]` injection) on the front (left) and tail (right) legs. + - Confirmed the left- and right-enriched bond spaces genuinely **differ** (e.g. L=4, `ℙ^8`: bond 1 + is `ℙ^2` from the left but `ℙ^4` from the right; bond 3 the reverse), matching the matrix recipe's + distinct `Û0`/`V̂0`. The final state uses the **left**-augmented bonds; `V̂0` is only the frozen + right environment for the K-step. + +- **The real blocker (precise):** the enrichment, K-step, and `4r` augmentation **cannot be + precomputed as separate passes and then assembled** — they must be **interleaved in a single + leaves→root sweep with growing bonds**. Precomputing `Û0` at rank `2r` and then augmenting to `4r` + yields **inconsistent chained bonds**: site `i`'s right bond becomes `4r` after augmentation, but + the precomputed `Û0[i+1]` was built with a `2r` left bond, so `As[i+1].left ≠ As[i].right`. The + first-order `_pbug_assemble_core` avoids this by doing the augmentation *inside* the sweep (stacking + each site's old block onto the previous site's already-grown new bond). The second-order sweep must + do the same but additionally (a) enrich each site's old block with `H·ψ₀` on the *current* (already + grown) bond, (b) run the K-step on the mixed `(Û0-left, V̂0-right)` enriched environments, (c) + augment to `4r`, and (d) propagate the off-diagonal coupling blocks (`Ũ2ᵀK`, `Lᵀ Ṽ2`) — all in one + pass. This interleaved, growing-bond, mixed-basis sweep with coupling propagation is the substantial + remaining implementation; it is a rewrite of the assembly sweep, not an add-on. + +## What is delivered now + +- `ParallelBUG` (first order): audited against the §3.1 checklist (interior tensors are + `[AL⁰ | Ũ]` isometries; solves act on the amplitude-weighted center `AC⁰`), and given a **committed + local-order slope-≈2 gate at `L = 3, 4`** — the gate that would catch any regression or a mislabeled + order. No code change was required (the audit confirmed the existing assembly). +- `ParallelBUG2`: struct + keyword constructor + `timestep!`/`timestep` wiring + export, sharing the + `AbstractParallelBUG` supertype and the first-order truncation/rejection helpers. Exact for `L ≤ 2`; + `throw`s an informative `ArgumentError` for `L ≥ 3` (pointing here) so it is never silently first + order. The `_pbug_preaugment` pass (correct `H·ψ₀` left enrichment, verified to preserve the state + to overlap 1.0) is retained as a building block for the port. diff --git a/src/MPSKit.jl b/src/MPSKit.jl index b7387f3da..547dd4f2e 100644 --- a/src/MPSKit.jl +++ b/src/MPSKit.jl @@ -35,7 +35,7 @@ export VUMPS, VOMPS, DMRG, DMRG2, IDMRG, IDMRG2, GradientGrassmann export excitations export FiniteExcited, QuasiparticleAnsatz, ChepigaAnsatz, ChepigaAnsatz2 export time_evolve, timestep, timestep!, make_time_mpo -export TDVP, TDVP2, WI, WII, TaylorCluster +export TDVP, TDVP2, BUG, ParallelBUG, ParallelBUG2, WI, WII, TaylorCluster export changebonds, changebonds! export VUMPSSvdCut, OptimalExpand, SvdCut, RandExpand, SketchedExpand export propagator @@ -161,6 +161,7 @@ include("algorithms/changebonds/randexpand.jl") include("algorithms/changebonds/sketchedexpand.jl") include("algorithms/timestep/tdvp.jl") +include("algorithms/timestep/bug.jl") include("algorithms/timestep/taylorcluster.jl") include("algorithms/timestep/wii.jl") include("algorithms/timestep/integrators.jl") diff --git a/src/algorithms/timestep/bug.jl b/src/algorithms/timestep/bug.jl new file mode 100644 index 000000000..dd1e0e5cf --- /dev/null +++ b/src/algorithms/timestep/bug.jl @@ -0,0 +1,639 @@ +""" +$(TYPEDEF) + +Single site MPS time-evolution algorithm based on the Basis-Update & Galerkin (BUG) integrator, +an unconventional robust integrator for dynamical low-rank approximation. + +Unlike [`TDVP`](@ref), BUG advances both the basis (K-step) and the core (Galerkin C-step) tensors +*forward* in time, with no backward-in-time substep. This makes it a natural choice for +imaginary-time / dissipative evolution, where the backward core step of TDVP can become unstable. + +Passing a truncating `trscheme` (anything other than `notrunc()`) switches on **rank-adaptivity**: +each half-sweep augments every bond with the new directions discovered by the evolved connecting +tensor (old basis first, `[U₀ │ K₁]`) and then truncates back down to the tolerance of `trscheme`. +The bond dimension grows and shrinks automatically to track the entanglement; `notrunc()` recovers +the fixed-rank integrator. + +!!! note + By default the state is not renormalized, so the norm keeps useful information (the + accumulated truncation error in real time, or the decaying weight in imaginary time). + Pass `normalize = true` to `timestep`/`time_evolve` to renormalize after every half-sweep + instead. This is independent of `imaginary_evolution`. + +## Fields + +$(TYPEDFIELDS) + +## References + +* [Ceruti et al. BIT Numer. Math. 62 (2022)](@cite ceruti2022) +""" +struct BUG{A, O, T, S, F} <: Algorithm + "algorithm used in the exponential solvers" + integrator::A + + "tolerance for gauging algorithm" + tolgauge::Float64 + + "maximal amount of iterations for gauging algorithm" + gaugemaxiter::Int + + "algorithm used to re-orthonormalize the basis after each local update" + alg_orth::O + + "truncation scheme used to cut the bond back down for rank-adaptive BUG" + trscheme::T + + "algorithm used for the singular value decomposition" + alg_svd::S + + "callback function applied after each iteration, of signature `finalize(iter, ψ, H, envs) -> ψ, envs`" + finalize::F +end +function BUG(; + integrator = Defaults.alg_expsolve(), tolgauge = Defaults.tolgauge, + gaugemaxiter = Defaults.maxiter, alg_orth = Defaults.alg_orth(), + trscheme = notrunc(), alg_svd = Defaults.alg_svd(), + finalize = Defaults._finalize + ) + return BUG(integrator, tolgauge, gaugemaxiter, alg_orth, trscheme, alg_svd, finalize) +end + +function timestep!( + ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::BUG, + envs::AbstractMPSEnvironments = environments(ψ, H, ψ); + imaginary_evolution::Bool = false, normalize::Bool = false + ) + # symmetric 2nd-order: a dt/2 left→right half-sweep composed with its dt/2 mirror + L = length(ψ) + h = dt / 2 + truncates = !(alg.trscheme isa MatrixAlgebraKit.NoTruncation) + svdcut = SvdCut(; trscheme = alg.trscheme, alg_svd = alg.alg_svd) + + # left→right half-sweep (root = last site) + ψ.AC[1] # gauge center to site 1 + ψ_old = copy(ψ) # frozen bases / reprojection inputs + T = isomorphism(scalartype(ψ), left_virtualspace(ψ_old, 1) ← left_virtualspace(ψ_old, 1)) + for i in 1:(L - 1) + Ĉ = _mul_front(T, ψ_old.AC[i]) # reproject old connecting tensor + AC = integrate(AC_hamiltonian(i, ψ, H, ψ, envs), Ĉ, t, h, alg.integrator; imaginary_evolution) + U₀ = _mul_front(T, ψ_old.AL[i]) # old left isometry in the new frame + if truncates + U = _bug_augment_left(U₀, AC, alg.alg_orth) # old-first augment; cut deferred to changebonds! + C = U' * AC + else + U, C = left_gauge(AC, alg.alg_orth) + end + T = U' * U₀ # transport (new ← old) + ψ.AC[i] = (U, C) + end + AC = integrate( + AC_hamiltonian(L, ψ, H, ψ, envs), _mul_front(T, ψ_old.AC[L]), + t, h, alg.integrator; imaginary_evolution + ) + normalize && normalize!(AC) + ψ.AC[L] = AC + truncates && changebonds!(ψ, svdcut; normalize) + + # right→left half-sweep (root = first site), the mirror + ψ.AC[L] # gauge center to site L + ψ_old = copy(ψ) + T = isomorphism(scalartype(ψ), right_virtualspace(ψ_old, L) ← right_virtualspace(ψ_old, L)) + for i in L:-1:2 + Ĉ = ψ_old.AC[i] * T # reproject old connecting tensor + AC = integrate(AC_hamiltonian(i, ψ, H, ψ, envs), Ĉ, t + h, h, alg.integrator; imaginary_evolution) + U₀ = ψ_old.AR[i] * T # old right isometry in the new frame + if truncates + U = _bug_augment_right(U₀, AC, alg.alg_orth) + C = _transpose_tail(AC) * _transpose_tail(U)' + else + C, U = right_gauge(AC, alg.alg_orth) + end + T = _transpose_tail(U₀) * _transpose_tail(U)' # transport (old ← new) + ψ.AC[i] = (C, U) + end + AC = integrate( + AC_hamiltonian(1, ψ, H, ψ, envs), ψ_old.AC[1] * T, + t + h, h, alg.integrator; imaginary_evolution + ) + normalize && normalize!(AC) + ψ.AC[1] = AC + truncates && changebonds!(ψ, svdcut; normalize) + + return ψ, envs +end + +# augment the RIGHT bond (left→right sweep): orthonormalize the stacked `[U₀ │ K₁]` (old isometry +# first) so `U₀` stays the leading per-sector block and only the new directions of `K₁` are appended. +function _bug_augment_left(U₀, K₁, alg_orth = Defaults.alg_orth()) + Û, _ = left_orth(catdomain(U₀, K₁); alg = alg_orth) + return Û +end + +# mirror of `_bug_augment_left` for the right→left sweep, on the `_transpose_tail` form (right-isometry +# with orthonormal rows): an LQ orthonormalizes the stacked rows `[U₀; K₁]`, keeping `U₀` leading. +function _bug_augment_right(U₀, K₁, alg_orth = Defaults.alg_orth()) + stacked = catcodomain(_transpose_tail(U₀), _transpose_tail(K₁)) + _, Û = right_orth(stacked; alg = alg_orth) + return _transpose_front(Û) +end + +# Parallel BUG +# ------------ + +# Shared supertype for the parallel-BUG integrators (`ParallelBUG` first order, `ParallelBUG2` second +# order); the assembly / truncation / rejection helpers below dispatch on it. +abstract type AbstractParallelBUG <: Algorithm end + +""" +$(TYPEDEF) + +Single site MPS time-evolution algorithm based on the *parallel* Basis-Update & Galerkin (BUG) +integrator for tree tensor networks, specialized to the linear (`FiniteMPS`) tree. + +Unlike the sequential [`BUG`](@ref), every local center `AC[i]` is evolved forward from the **same +frozen `t₀` snapshot**: there is no sweep and no sequential dependency, so the local integrations are +mutually independent and parallelizable. A cheap leaves→root pass then augments every bond with the +new directions the frozen evolutions discovered (old basis first, `[U₀ │ Ũ₁]`), and a final SVD sweep +truncates the (at most doubled) bonds back down. Like [`BUG`](@ref) it advances every tensor *forward* +in time (no backward substep), which suits imaginary-time / dissipative evolution. + +The integrator is **first-order** in time. Any truncating `trscheme` makes the bond dimension grow +and shrink to track the entanglement; the default `notrunc()` restores the pre-step virtual spaces +(fixed-rank parallel BUG). Setting `maxiter_rejection > 0` enables **step rejection**: a step that +saturates the doubling on some bond is recomputed as two half-steps (cf. Ceruti et al. 2024). + +!!! warning "Experimental" + This integrator is **work in progress**; the API and behaviour may change. + +!!! note + By default the state is not renormalized, so the norm keeps useful information (the + accumulated truncation error in real time, or the decaying weight in imaginary time). + Pass `normalize = true` to `timestep`/`time_evolve` to renormalize after every step + instead. This is independent of `imaginary_evolution`. + +## Fields + +$(TYPEDFIELDS) + +## References + +* Ceruti, Kusch & Lubich, *A parallel rank-adaptive integrator for dynamical low-rank + approximation*, SIAM J. Sci. Comput. **46** (2024). +* Ceruti, Kusch, Lubich & Sulz, *A parallel Basis Update and Galerkin integrator for tree tensor + networks*, arXiv:2412.00858 (2024). +""" +struct ParallelBUG{A, O, T, S, F} <: AbstractParallelBUG + "algorithm used in the exponential solvers" + integrator::A + + "tolerance for gauging algorithm" + tolgauge::Float64 + + "maximal amount of iterations for gauging algorithm" + gaugemaxiter::Int + + "algorithm used to re-orthonormalize the basis after each local update" + alg_orth::O + + "truncation scheme used to cut the augmented bonds back down" + trscheme::T + + "algorithm used for the singular value decomposition" + alg_svd::S + + "safety constant `c` in the step-rejection threshold `h·η > c·ϑ` (paper value ≈ 10)" + c::Float64 + + "maximum number of rejection recomputes per step (0 disables step rejection)" + maxiter_rejection::Int + + "callback function applied after each iteration, of signature `finalize(iter, ψ, H, envs) -> ψ, envs`" + finalize::F +end +function ParallelBUG(; + integrator = Defaults.alg_expsolve(), tolgauge = Defaults.tolgauge, + gaugemaxiter = Defaults.maxiter, alg_orth = Defaults.alg_orth(), + trscheme = notrunc(), alg_svd = Defaults.alg_svd(), + c = 10.0, maxiter_rejection = 0, + finalize = Defaults._finalize + ) + return ParallelBUG( + integrator, tolgauge, gaugemaxiter, alg_orth, trscheme, alg_svd, + c, maxiter_rejection, finalize + ) +end + +""" +$(TYPEDEF) + +Single site MPS time-evolution algorithm based on the *second-order* variant of the *parallel* +Basis-Update & Galerkin (BUG) integrator for tree tensor networks (Kusch 2024, "Variant 2"), +specialized to the linear (`FiniteMPS`) tree. It is to [`ParallelBUG`](@ref) what [`TDVP2`](@ref) is +to [`TDVP`](@ref): a separate, more accurate integrator sharing the same interface. + +The genuine second order comes from **pre-augmenting every bond basis with one `H·ψ₀` application +before evolving**: each bond isometry is enlarged to rank `2r` with the directions opened by a single +effective-Hamiltonian application to the frozen center, the Galerkin K-steps are evolved on those +enriched environments, and the `O(dt²)` content the first-order scheme discards is transported to the +root, keeping the "new–new" corner zero (so the local error is `O(dt³)` rather than `O(dt²)`). +Everything else matches [`ParallelBUG`](@ref): one frozen `t₀` snapshot, mutually independent local +solves, amplitude carried once at the root, `notrunc()` for the fixed-rank variant, and forward-only +evolution (no backward substep). + +!!! warning "Experimental" + This integrator is **work in progress**; the API and behaviour may change. Its single-step error + scales as `O(dt³)` (log–log slope `≈ 3`) versus the first-order `O(dt²)`. Only plain Hamiltonians + are supported (no `LazySum` / time-dependent operators). + +!!! note + By default the state is not renormalized, so the norm keeps useful information (the + accumulated truncation error in real time, or the decaying weight in imaginary time). + Pass `normalize = true` to `timestep`/`time_evolve` to renormalize after every step + instead. This is independent of `imaginary_evolution`. + +## Fields + +$(TYPEDFIELDS) + +## References + +* Kusch, *Second-order robust parallel integrators for dynamical low-rank approximation*, + arXiv:2403.02834 (2024). +* Ceruti, Kusch, Lubich & Sulz, *A parallel Basis Update and Galerkin integrator for tree tensor + networks*, arXiv:2412.00858 (2024). +""" +struct ParallelBUG2{A, O, T, S, F} <: AbstractParallelBUG + "algorithm used in the exponential solvers" + integrator::A + + "tolerance for gauging algorithm" + tolgauge::Float64 + + "maximal amount of iterations for gauging algorithm" + gaugemaxiter::Int + + "algorithm used to re-orthonormalize the basis after each local update" + alg_orth::O + + "truncation scheme used to cut the augmented bonds back down" + trscheme::T + + "algorithm used for the singular value decomposition" + alg_svd::S + + "safety constant `c` in the step-rejection threshold `h·η > c·ϑ` (paper value ≈ 10)" + c::Float64 + + "maximum number of rejection recomputes per step (0 disables step rejection)" + maxiter_rejection::Int + + "callback function applied after each iteration, of signature `finalize(iter, ψ, H, envs) -> ψ, envs`" + finalize::F +end +function ParallelBUG2(; + integrator = Defaults.alg_expsolve(), tolgauge = Defaults.tolgauge, + gaugemaxiter = Defaults.maxiter, alg_orth = Defaults.alg_orth(), + trscheme = notrunc(), alg_svd = Defaults.alg_svd(), + c = 10.0, maxiter_rejection = 0, + finalize = Defaults._finalize + ) + return ParallelBUG2( + integrator, tolgauge, gaugemaxiter, alg_orth, trscheme, alg_svd, + c, maxiter_rejection, finalize + ) +end + +function timestep!( + ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::AbstractParallelBUG, + envs::AbstractMPSEnvironments = environments(ψ, H, ψ); + imaginary_evolution::Bool = false, normalize::Bool = false + ) + L = length(ψ) + if L == 1 # single site: a plain forward center step + AC = integrate( + AC_hamiltonian(1, ψ, H, ψ, envs), ψ.AC[1], t, dt, alg.integrator; imaginary_evolution + ) + normalize && normalize!(AC) + ψ.AC[1] = AC + return ψ, envs + end + + truncates = !(alg.trscheme isa MatrixAlgebraKit.NoTruncation) + Vs = [right_virtualspace(ψ, b) for b in 1:(L - 1)] # pre-step spaces, for the fixed-rank restore + + ϕ, ηh = _pbug_assemble(ψ, H, t, dt, alg; imaginary_evolution, normalize) + augVs = [right_virtualspace(ϕ, b) for b in 1:(L - 1)] # the (doubled) augmented spaces + _pbug_truncate!(ϕ, alg, Vs; normalize) + + # step rejection (opt-in): a bond that kept its full augmented space was under-resolved by a + # single doubling; recompute as two half-steps so one doubling per sub-step suffices. + if truncates && alg.maxiter_rejection > 0 + saturated = any(1:(L - 1)) do b + return right_virtualspace(ϕ, b) == augVs[b] && augVs[b] != Vs[b] + end + @debug "ParallelBUG step" ηh saturated + if saturated + alg′ = _pbug_with_rejections(alg, alg.maxiter_rejection - 1) + timestep!(ψ, H, t, dt / 2, alg′; imaginary_evolution, normalize) + timestep!(ψ, H, t + dt / 2, dt / 2, alg′; imaginary_evolution, normalize) + return ψ, environments(ψ, H, ψ) + end + end + + # overwrite `ψ` in place with the assembled state (adopt `ϕ`'s tensors; identity-keyed envs self-heal) + for f in (:ALs, :ARs, :ACs, :Cs) + copyto!(getfield(ψ, f), getfield(ϕ, f)) + end + return ψ, environments(ψ, H, ψ) +end + +# rebuild a parallel-BUG algorithm with a reduced rejection budget (the structs are immutable) +function _pbug_with_rejections(alg::AbstractParallelBUG, n::Int) + return typeof(alg)( + alg.integrator, alg.tolgauge, alg.gaugemaxiter, alg.alg_orth, + alg.trscheme, alg.alg_svd, alg.c, n, alg.finalize + ) +end + +# ---- second-order (Variant 2) assembly --------------------------------------------------------- + +# old-first LEFT enrichment `Û0[1..L-1]` (rank `2r`) plus the enriched left-environment chain `GLhat` +# (`⟨Û0|H|Û0⟩`): enlarge each `AL⁰[i]` with the range of the frozen derivative image `W[i]=(H·ψ₀)ᵢ`, +# stacked leaves→root with the mixed `⟨new|H|old⟩` coupling so directions opened deep in the chain +# reach the root. Envs are folded by explicit transfer (no `FiniteMPS` round-trip, so the zero-weight +# enriched directions do not collapse under canonicalization). +function _pbug2_left_enrich(ψ₀, H, envs₀, W, alg_orth) + L = length(ψ₀) + Û0 = Vector{typeof(ψ₀.AL[1])}(undef, L - 1) + GLhat = Vector{Any}(undef, L) + GLhat[1] = leftenv(envs₀, 1, ψ₀) + GLmix = leftenv(envs₀, 1, ψ₀) # mixed ⟨new|H|old⟩ chain + local GLnew + for i in 1:(L - 1) + if i == 1 + C⁰, Ĉ = ψ₀.AL[1], W[1] + else + C̃ = MPO_AC_Hamiltonian(GLnew, H[i], rightenv(envs₀, i, ψ₀))(ψ₀.AC[i]) + C⁰ = _pbug_stack_child(ψ₀.AL[i], zerovector!(similar(C̃))) + Ĉ = _pbug_stack_child(W[i], C̃) + end + Ũ, = _pbug_newdirs(C⁰, Ĉ, alg_orth) + Û0[i] = catdomain(C⁰, Ũ) + GLhat[i + 1] = GLhat[i] * TransferMatrix(Û0[i], H[i], Û0[i]) + GLnew = GLmix * TransferMatrix(ψ₀.AL[i], H[i], Ũ) # ket=old, bra=new + i == L - 1 && break + GLmix = GLmix * TransferMatrix(ψ₀.AL[i], H[i], Û0[i]) # ket=old, bra=enriched + end + return Û0, GLhat +end + +# old-first RIGHT enrichment `V̂0[2..L]` (rank `2r`) plus the enriched right-environment chain `GRhat`, +# the mirror of `_pbug2_left_enrich`. The interior K-step freezes this enriched right basis (freezing +# the old right basis instead only yields local slope 2). +function _pbug2_right_enrich(ψ₀, H, envs₀, W, alg_orth) + L = length(ψ₀) + V̂0 = Vector{typeof(ψ₀.AR[1])}(undef, L) + GRhat = Vector{Any}(undef, L) + GRhat[L] = rightenv(envs₀, L, ψ₀) + GRmix = rightenv(envs₀, L, ψ₀) + local GRnew + for i in L:-1:2 + if i == L + C⁰, Ĉ = ψ₀.AR[L], W[L] + else + C̃ = MPO_AC_Hamiltonian(leftenv(envs₀, i, ψ₀), H[i], GRnew)(ψ₀.AC[i]) + C⁰ = catdomain(ψ₀.AR[i], zerovector!(similar(C̃))) + Ĉ = catdomain(W[i], C̃) + end + # new left-bond directions `Ĉ` opens beyond the right-isometry `C⁰` (in `_transpose_tail` form) + N = right_null!(_transpose_tail(C⁰; copy = true)) + _, Q = right_orth(_transpose_tail(Ĉ) * N'; alg = alg_orth) + V̂ = Q * N + V̂0[i] = _transpose_front(catcodomain(_transpose_tail(C⁰), V̂)) + GRhat[i - 1] = TransferMatrix(V̂0[i], H[i], V̂0[i]) * GRhat[i] + GRnew = TransferMatrix(ψ₀.AR[i], H[i], _transpose_front(V̂)) * GRmix # ket=old, bra=new + i == 2 && break + GRmix = TransferMatrix(ψ₀.AR[i], H[i], V̂0[i]) * GRmix # ket=old, bra=enriched + end + return V̂0, GRhat +end + +# Genuine second-order parallel-BUG assembly (Kusch 2024, Variant 2 / `4r`) rooted at site `L`, from a +# single frozen `t₀` snapshot: (1) pre-augment the bond bases to rank `2r` with one `H·ψ₀` application +# (`Û0`/`V̂0`) and fold the enriched envs by explicit transfer; (2) K-step every center on the enriched +# envs, freezing `V̂0`; (3) a leaves→root sweep builds the `4r` isometries `[Û0 | Ũ2]` and transports +# the evolved-amplitude coupling `R = Ũ2ᵀĈ` through the frozen `V̂0`; (4) the root stacks the `2r` +# Galerkin `Kevo[L]` with the transported coupling, keeping the "new–new" corner zero (local `O(dt³)`). +function _pbug_assemble(ψ, H, t, dt, alg::ParallelBUG2; imaginary_evolution::Bool = false, normalize::Bool = false) + L = length(ψ) + ψ.AC[L] # gauge to the root + ψ₀ = copy(ψ) + envs₀ = environments(ψ₀, H, ψ₀) + _pbug_warmup_envs!(envs₀, L, ψ₀) + scheduler = Defaults.scheduler[] + + # frozen `H·ψ₀` images (independent, threaded) + W = Vector{typeof(ψ₀.AC[1])}(undef, L) + tmap!(W, 1:L; scheduler) do i + return AC_hamiltonian(i, ψ₀, H, ψ₀, envs₀)(ψ₀.AC[i]) + end + + Û0, GLhat = _pbug2_left_enrich(ψ₀, H, envs₀, W, alg.alg_orth) + V̂0, GRhat = _pbug2_right_enrich(ψ₀, H, envs₀, W, alg.alg_orth) + + GRof(i) = i < L ? GRhat[i] : rightenv(envs₀, L, ψ₀) + enrL(i) = i == 1 ? left_virtualspace(ψ₀, 1) : + (i <= L - 1 ? space(Û0[i], 1) : only(domain(Û0[L - 1]))) + enrR(i) = i == L ? right_virtualspace(ψ₀, L) : space(V̂0[i + 1], 1) + # embed the old center into the enriched bond spaces with zero weight in the new directions + embed(i) = absorb!( + zerovector!(similar(ψ₀.AC[i], (enrL(i) ⊗ physicalspace(ψ₀, i)) ← enrR(i))), ψ₀.AC[i] + ) + + # K/S steps on the enriched environments (freeze `V̂0`); independent ⇒ threaded + Kevo = Vector{typeof(ψ₀.AC[1])}(undef, L) + tmap!(Kevo, 1:L; scheduler) do i + return integrate( + MPO_AC_Hamiltonian(GLhat[i], H[i], GRof(i)), embed(i), t, dt, alg.integrator; + imaginary_evolution + ) + end + + # leaves→root assembly with the transported evolved-amplitude coupling + As = Vector{typeof(ψ₀.AL[1])}(undef, L) + ηh = zero(real(scalartype(ψ₀))) + local R + for i in 1:(L - 1) + if i == 1 + C⁰, Ĉ¹ = Û0[1], Kevo[1] + else + C̃ = _transpose_front(R * _transpose_tail(V̂0[i])) # transport R through frozen V̂0[i] + ηh = max(ηh, norm(C̃)) + zc0 = zerovector!(similar(C̃, codomain(C̃) ← domain(Û0[i]))) + C⁰ = _pbug_stack_child(Û0[i], zc0) + Ĉ¹ = _pbug_stack_child(Kevo[i], C̃) + end + Ũ2, = _pbug_newdirs(C⁰, Ĉ¹, alg.alg_orth) + As[i] = catdomain(C⁰, Ũ2) + R = Ũ2' * Ĉ¹ # evolved amplitude in the new dirs + end + # root: `2r` Galerkin `Kevo[L]` (Û0 rows) stacked with the transported coupling (Ũ2 rows); the + # "new–new" corner stays zero (right bond trivial at the root), so the local error is `O(dt³)`. + C̃L = _transpose_front(R * _transpose_tail(V̂0[L])) + ηh = max(ηh, norm(C̃L)) + As[L] = _pbug_stack_child(Kevo[L], C̃L) + + return FiniteMPS(As; overwrite = true, normalize), ηh +end + +# First-order parallel-BUG assembly (Ceruti et al. 2024, arXiv:2412.00858, Alg. 1-4) rooted at site +# `L`. Phase 1: Galerkin-evolve every center `AC[i]` from the frozen snapshot (independent ⇒ parallel). +# Phase 2: a leaves→root sweep orthonormalizes each evolved center *stacked with a first-order coupling +# block* `C̃ᵢ` on the previous bond's new rows (`[old │ Ũᵢ]`), propagating deep new directions to the +# root. Interior tensors are pure isometries (amplitude/phase discarded); all first-order + amplitude +# content enters once, at the root `[C̄_L; C̃_L]`. The zero "new–new" corners make it first order. +function _pbug_assemble(ψ, H, t, dt, alg::ParallelBUG; imaginary_evolution::Bool = false, normalize::Bool = false) + L = length(ψ) + ψ.AC[L] # gauge to the root + ψ₀ = copy(ψ) + envs₀ = environments(ψ₀, H, ψ₀) + dt′ = imaginary_evolution ? -dt : -im * dt + + # warm the lazily-cached envs serially so the threaded phase-1 solves below only read them + _pbug_warmup_envs!(envs₀, L, ψ₀) + + # phase 1: frozen-snapshot Galerkin evolutions (independent local solves, threaded) + scheduler = Defaults.scheduler[] + Cevo = Vector{typeof(ψ₀.AC[1])}(undef, L) + tmap!(Cevo, 1:L; scheduler) do i + return integrate( + AC_hamiltonian(i, ψ₀, H, ψ₀, envs₀), ψ₀.AC[i], t, dt, alg.integrator; + imaginary_evolution + ) + end + C̄L = Cevo[L] + + # phase 2: leaves→root augmentation, threading the mixed ⟨augmented|H|old⟩ envs. `ηh = max‖C̃ᵢ‖` is + # the local error estimator (Ceruti-Kusch-Lubich 2024, eq. 6): the frozen derivative projected onto + # the new directions the old basis misses. + As = Vector{typeof(ψ₀.AL[1])}(undef, L) + ηh = zero(real(scalartype(ψ₀))) + GLmix = _pbug_mixedenv_init(H, envs₀, ψ₀) + local GLnew + for i in 1:(L - 1) + if i == 1 + C⁰, Ĉ¹ = ψ₀.AL[1], Cevo[1] + else + # first-order coupling block on the new rows of bond i-1 (two-arg apply at midpoint + # `t + dt/2`, matching what `integrate` freezes; a `TimedOperator` has no one-arg apply) + C̃ = scale(_pbug_coupling_hamiltonian(GLnew, H, i, envs₀, ψ₀)(ψ₀.AC[i], t + dt / 2), dt′) + ηh = max(ηh, norm(C̃)) + C⁰ = _pbug_stack_child(ψ₀.AL[i], zerovector!(similar(C̃))) + Ĉ¹ = _pbug_stack_child(Cevo[i], C̃) + end + Ũ, = _pbug_newdirs(C⁰, Ĉ¹, alg.alg_orth) + As[i] = catdomain(C⁰, Ũ) + GLnew = _pbug_mixedenv_step(GLmix, H, i, ψ₀.AL[i], Ũ) + i == L - 1 && break # the full mixed environment is no longer needed + GLmix = _pbug_mixedenv_step(GLmix, H, i, ψ₀.AL[i], As[i]) + end + # root: the amplitude is carried once, by the evolved center and its coupling row + C̃L = scale(_pbug_coupling_hamiltonian(GLnew, H, L, envs₀, ψ₀)(ψ₀.AC[L], t + dt / 2), dt′) + ηh = max(ηh, norm(C̃L)) + As[L] = _pbug_stack_child(C̄L, C̃L) + + return FiniteMPS(As; overwrite = true, normalize), ηh +end + +# warm the lazily-cached frozen environments serially (leftenv/rightenv on `FiniteEnvironments` mutate +# their cache), so the parallel phase-1 solves only read them. A `LazySum` yields a +# `MultipleEnvironments`, so recurse into its per-summand `FiniteEnvironments`. +function _pbug_warmup_envs!(envs::FiniteEnvironments, L, ψ₀) + for i in 1:L + leftenv(envs, i, ψ₀) + rightenv(envs, i, ψ₀) + end + return envs +end +function _pbug_warmup_envs!(envs::MultipleEnvironments, L, ψ₀) + foreach(e -> _pbug_warmup_envs!(e, L, ψ₀), envs.envs) + return envs +end + +# mixed ⟨augmented|H|old⟩ left environments: initial (trivial-bond) env and one-site transfer step, +# dispatching through `MultipliedOperator`/`LazySum` like the effective-Hamiltonian constructors. +_pbug_mixedenv_init(H, envs, ψ₀) = leftenv(envs, 1, ψ₀) +_pbug_mixedenv_init(H::MultipliedOperator, envs, ψ₀) = _pbug_mixedenv_init(H.op, envs, ψ₀) +function _pbug_mixedenv_init(H::LazySum, envs::MultipleEnvironments, ψ₀) + return map((o, e) -> _pbug_mixedenv_init(o, e, ψ₀), H.ops, envs.envs) +end + +_pbug_mixedenv_step(GL, H, i, above, below) = GL * TransferMatrix(above, H[i], below) +_pbug_mixedenv_step(GL, H::MultipliedOperator, i, above, below) = + _pbug_mixedenv_step(GL, H.op, i, above, below) +function _pbug_mixedenv_step(GLs::Vector, H::LazySum, i, above, below) + return map((gl, o) -> _pbug_mixedenv_step(gl, o, i, above, below), GLs, H.ops) +end + +# one-site effective derivative with the mixed new-direction rows as left env and the frozen old right +# env: applied to the center `AC[i]` this yields the first-order coupling block `C̃ᵢ`. +_pbug_coupling_hamiltonian(GL, H, i, envs, ψ₀) = MPO_AC_Hamiltonian(GL, H[i], rightenv(envs, i, ψ₀)) +function _pbug_coupling_hamiltonian(GL, H::MultipliedOperator, i, envs, ψ₀) + return MultipliedOperator(_pbug_coupling_hamiltonian(GL, H.op, i, envs, ψ₀), H.f) +end +function _pbug_coupling_hamiltonian(GLs::Vector, H::LazySum, i, envs::MultipleEnvironments, ψ₀) + Hs = map((gl, o, e) -> _pbug_coupling_hamiltonian(gl, o, i, e, ψ₀), GLs, H.ops, envs.envs) + elT = Union{D, MultipliedOperator{D}} where {D <: DerivativeOperator} + return LazySum{elT}(Hs) +end + +# new bond directions: the component of the evolved (stacked) candidate orthogonal to the old basis, +# re-orthonormalized. NOTE: keep `left_null` here — swapping it for `project_complement!`+QR adds +# completion columns outside the old-basis complement and breaks two-site exactness / first order. +function _pbug_newdirs(AL, Cevo, alg_orth = Defaults.alg_orth()) + N = left_null(AL) + g = N' * Cevo + Q, _ = left_orth(g; alg = alg_orth) + return N * Q, domain(Q) +end + +# stack two MPS tensors along the child (left-virtual) bond: `[top; bot]`, doubling that bond +_pbug_stack_child(top, bot) = + _transpose_front(catcodomain(_transpose_tail(top), _transpose_tail(bot))) + +# cut the augmented bonds back down: a truncating `trscheme` selects rank-adaptivity, `notrunc()` +# restores the pre-step virtual space of every bond (fixed-rank). Environments self-heal lazily. +function _pbug_truncate!(ϕ, alg::AbstractParallelBUG, Vs; normalize::Bool = false) + if !(alg.trscheme isa MatrixAlgebraKit.NoTruncation) + changebonds!(ϕ, SvdCut(; trscheme = alg.trscheme, alg_svd = alg.alg_svd); normalize) + else + for i in (length(ϕ) - 1):-1:1 + U, S, Vᴴ = svd_trunc(ϕ.C[i]; trunc = truncspace(Vs[i]), alg = alg.alg_svd) + ϕ.AC[i] = (ϕ.AL[i] * U, S) + ϕ.AC[i + 1] = (S, _transpose_front(Vᴴ * _transpose_tail(ϕ.AR[i + 1]))) + end + normalize && normalize!(ϕ) + end + return ϕ +end + +# copying version, shared by both BUG integrators +function timestep( + ψ::AbstractFiniteMPS, H, time::Number, timestep::Number, + alg::Union{BUG, ParallelBUG, ParallelBUG2}, envs::AbstractMPSEnvironments...; + imaginary_evolution::Bool = false, normalize::Bool = false, kwargs... + ) + isreal = (scalartype(ψ) <: Real && !imaginary_evolution) + ψ′ = isreal ? complex(ψ) : copy(ψ) + if length(envs) != 0 && isreal + @warn "Currently cannot reuse real environments for complex evolution" + envs′ = environments(ψ′, H, ψ′) + elseif length(envs) == 1 + envs′ = only(envs) + else + @assert length(envs) == 0 "Invalid signature" + envs′ = environments(ψ′, H, ψ′) + end + return timestep!(ψ′, H, time, timestep, alg, envs′; imaginary_evolution, normalize, kwargs...) +end diff --git a/src/algorithms/timestep/tdvp.jl b/src/algorithms/timestep/tdvp.jl index 993782e54..afaf3ba4e 100644 --- a/src/algorithms/timestep/tdvp.jl +++ b/src/algorithms/timestep/tdvp.jl @@ -11,10 +11,11 @@ the enlarged bond back down (selecting the truncated-SVD gauge). The expansion i state-preserving, as required for a consistent time evolution. !!! note - Real-time evolution preserves the norm: neither the bond expansion nor the truncation - renormalizes, so the state norm reflects the accumulated truncation error. Imaginary-time - evolution instead renormalizes at every step, like a ground-state search. CBE is only - available for finite MPS. + By default the norm is preserved: neither the bond expansion nor the truncation + renormalizes, so the state norm keeps useful information (the accumulated truncation + error in real time, or the decaying weight in imaginary time). Pass `normalize = true` + to `timestep`/`time_evolve` to renormalize at every step instead, like a ground-state + search. This is independent of `imaginary_evolution`. CBE is only available for finite MPS. ## Fields @@ -61,11 +62,14 @@ end function timestep( ψ::InfiniteMPS, H, t::Number, dt::Number, alg::TDVP, envs::AbstractMPSEnvironments = environments(ψ, H, ψ); - leftorthflag = true, imaginary_evolution::Bool = false + leftorthflag = true, imaginary_evolution::Bool = false, normalize::Bool = false ) + # `normalize` is accepted for signature uniformity with the finite integrators, but an + # `InfiniteMPS` is always normalized to norm-1-per-site by the gauge/reconstruction below + # (a structural gauge requirement, not information erasure), so the flag has no effect here. # convert state to complex if necessary if scalartype(ψ) <: Real && (!imaginary_evolution || !isreal(dt)) - return timestep(complex(ψ), H, t, dt, alg, envs; leftorthflag, imaginary_evolution) + return timestep(complex(ψ), H, t, dt, alg, envs; leftorthflag, imaginary_evolution, normalize) end temp_ACs = similar(ψ.AC) @@ -120,23 +124,23 @@ end function timestep!( ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP, envs::AbstractMPSEnvironments = environments(ψ, H, ψ); - imaginary_evolution::Bool = false + imaginary_evolution::Bool = false, normalize::Bool = false ) # sweep left to right for i in 1:(length(ψ) - 1) # 1. optionally expand the bond ahead of the local update (CBE) isnothing(alg.alg_expand) || - changebond!(i, Val(:right), ψ, H, alg.alg_expand, envs; normalize = imaginary_evolution) + changebond!(i, Val(:right), ψ, H, alg.alg_expand, envs; normalize) # 2. evolve the (possibly expanded) center tensor forward Hac = AC_hamiltonian(i, ψ, H, ψ, envs) AC = integrate(Hac, ψ.AC[i], t, dt / 2, alg.integrator; imaginary_evolution) # 3. gauge: split AC -> AL[i], C[i] (QR center-move, or truncated SVD cutting the - # enlarged bond back down) and move the center to i+1. Real-time evolution preserves - # the norm; imaginary-time evolution renormalizes. - left_gauge!(ψ, i, AC, alg.alg_gauge; normalize = imaginary_evolution) + # enlarged bond back down) and move the center to i+1. By default the norm is + # preserved; `normalize` renormalizes. + left_gauge!(ψ, i, AC, alg.alg_gauge; normalize) # 4. evolve the bond tensor backward Hc = C_hamiltonian(i, ψ, H, ψ, envs) @@ -154,7 +158,7 @@ function timestep!( for i in length(ψ):-1:2 # 1. optionally expand the bond ahead of the local update (CBE) isnothing(alg.alg_expand) || - changebond!(i, Val(:left), ψ, H, alg.alg_expand, envs; normalize = imaginary_evolution) + changebond!(i, Val(:left), ψ, H, alg.alg_expand, envs; normalize) # 2. evolve the (possibly expanded) center tensor forward Hac = AC_hamiltonian(i, ψ, H, ψ, envs) @@ -163,9 +167,9 @@ function timestep!( imaginary_evolution ) - # 3. gauge: split AC -> C[i-1], AR[i] and move the center to i-1 (real-time preserves the - # norm; imaginary-time renormalizes) - right_gauge!(ψ, i, AC, alg.alg_gauge; normalize = imaginary_evolution) + # 3. gauge: split AC -> C[i-1], AR[i] and move the center to i-1 (norm preserved by + # default; `normalize` renormalizes) + right_gauge!(ψ, i, AC, alg.alg_gauge; normalize) # 4. evolve the bond tensor backward Hc = C_hamiltonian(i - 1, ψ, H, ψ, envs) @@ -221,7 +225,7 @@ end function timestep!( ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP2, envs::AbstractMPSEnvironments = environments(ψ, H, ψ); - imaginary_evolution::Bool = false + imaginary_evolution::Bool = false, normalize::Bool = false ) # sweep left to right @@ -231,6 +235,9 @@ function timestep!( ac2′ = integrate(Hac2, ac2, t, dt / 2, alg.integrator; imaginary_evolution) nal, nc, nar = svd_trunc!(ac2′; trunc = alg.trscheme, alg = alg.alg_svd) + # `nc` is the norm-carrying bond tensor (`nal`/`nar` are isometries), so normalizing it + # normalizes the whole state, mirroring single-site TDVP's per-gauge renormalization + normalize && normalize!(nc) ψ.AC[i] = (nal, complex(nc)) ψ.AC[i + 1] = (complex(nc), _transpose_front(nar)) @@ -250,6 +257,7 @@ function timestep!( ac2′ = integrate(Hac2, ac2, t + dt / 2, dt / 2, alg.integrator; imaginary_evolution) nal, nc, nar = svd_trunc!(ac2′; trunc = alg.trscheme, alg = alg.alg_svd) + normalize && normalize!(nc) ψ.AC[i - 1] = (nal, complex(nc)) ψ.AC[i] = (complex(nc), _transpose_front(nar)) @@ -269,7 +277,7 @@ end function timestep( ψ::AbstractFiniteMPS, H, time::Number, timestep::Number, alg::Union{TDVP, TDVP2}, envs::AbstractMPSEnvironments...; - imaginary_evolution::Bool = false, kwargs... + imaginary_evolution::Bool = false, normalize::Bool = false, kwargs... ) isreal = (scalartype(ψ) <: Real && !imaginary_evolution) ψ′ = isreal ? complex(ψ) : copy(ψ) @@ -282,5 +290,5 @@ function timestep( @assert length(envs) == 0 "Invalid signature" envs′ = environments(ψ′, H, ψ′) end - return timestep!(ψ′, H, time, timestep, alg, envs′; imaginary_evolution, kwargs...) + return timestep!(ψ′, H, time, timestep, alg, envs′; imaginary_evolution, normalize, kwargs...) end diff --git a/src/algorithms/timestep/time_evolve.jl b/src/algorithms/timestep/time_evolve.jl index 51288abd7..a321520c0 100644 --- a/src/algorithms/timestep/time_evolve.jl +++ b/src/algorithms/timestep/time_evolve.jl @@ -20,6 +20,11 @@ through each of the time points obtained by iterating t_span. instead, (i.e. ``\\exp(-Hdt)`` instead of ``\\exp(-iHdt)``). This can be useful for using this function to compute the ground state of a Hamiltonian, or to compute finite-temperature properties of a system. +- `normalize::Bool=false`: if true, the state is renormalized after every step. This is + independent of `imaginary_evolution`: by default the norm is preserved, so it retains + useful information (the accumulated truncation error in real time, or the decaying weight + in imaginary time). Pass `true` to renormalize each step, e.g. when computing a ground + state or thermal state via imaginary-time evolution. """ function time_evolve end, function time_evolve! end @@ -27,7 +32,7 @@ for (timestep, time_evolve) in zip((:timestep, :timestep!), (:time_evolve, :time @eval function $time_evolve( ψ, H, t_span::AbstractVector{<:Number}, alg, envs = environments(ψ, H, ψ); - verbosity::Int = 0, imaginary_evolution::Bool = false + verbosity::Int = 0, imaginary_evolution::Bool = false, normalize::Bool = false ) log = IterLog("TDVP") LoggingExtras.withlevel(; verbosity) do @@ -36,7 +41,7 @@ for (timestep, time_evolve) in zip((:timestep, :timestep!), (:time_evolve, :time t = t_span[iter] dt = t_span[iter + 1] - t - ψ, envs = $timestep(ψ, H, t, dt, alg, envs; imaginary_evolution) + ψ, envs = $timestep(ψ, H, t, dt, alg, envs; imaginary_evolution, normalize) ψ, envs = alg.finalize(t, ψ, H, envs)::Tuple{typeof(ψ), typeof(envs)} @infov 3 logiter!(log, iter, 0, t) @@ -69,6 +74,11 @@ solving the Schroedinger equation: ``i ∂ψ/∂t = H ψ``. instead, (i.e. ``\\exp(-Hdt)`` instead of ``\\exp(-iHdt)``). This can be useful for using this function to compute the ground state of a Hamiltonian, or to compute finite-temperature properties of a system. +- `normalize::Bool=false`: if true, the state is renormalized after the step. This is + independent of `imaginary_evolution`: by default the norm is preserved, so it retains + useful information (the accumulated truncation error in real time, or the decaying weight + in imaginary time). Pass `true` to renormalize, e.g. when computing a ground state or + thermal state via imaginary-time evolution. """ function timestep end, function timestep! end diff --git a/test/algorithms/bug.jl b/test/algorithms/bug.jl new file mode 100644 index 000000000..d635ebd66 --- /dev/null +++ b/test/algorithms/bug.jl @@ -0,0 +1,527 @@ +println(" +----------------------------- +| BUG time-stepping tests | +----------------------------- +") + +using .TestSetup +using Test, TestExtras +using MPSKit +using TensorKit +using TensorKit: ℙ +using LinearAlgebra: dot, norm +using Random + +@testset "BUG time evolution" verbose = true begin + dt = 0.1 + L = 10 + + H = force_planar(heisenberg_XXX(Float64, Trivial; spin = 1 // 2, L)) + ψ = FiniteMPS(rand, Float64, L, ℙ^2, ℙ^4) + ψ₀, = find_groundstate(ψ, H; verbosity = 0) + E₀ = expectation_value(ψ₀, H) + + # 1. energy conservation + eigenstate phase + @testset "energy conservation" begin + ψ1, envs = timestep(ψ₀, H, 0.0, dt, BUG()) + E1 = expectation_value(ψ1, H, envs) + @test E₀ ≈ E1 atol = 1.0e-2 + @test dot(ψ1, ψ₀) ≈ exp(im * dt * E₀) atol = 1.0e-4 + end + + # 2. agreement with TDVP over a few real-time steps of a random MPS + @testset "agreement with TDVP" begin + Random.seed!(1234) + ψr = complex(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^4)) + δt = 0.01 + ψ_bug, ψ_tdvp = ψr, ψr + for k in 0:4 + ψ_bug, = timestep(ψ_bug, H, k * δt, δt, BUG()) + ψ_tdvp, = timestep(ψ_tdvp, H, k * δt, δt, TDVP()) + end + @test expectation_value(ψ_bug, H) ≈ expectation_value(ψ_tdvp, H) atol = 1.0e-3 + @test abs(dot(ψ_bug, ψ_tdvp)) ≈ 1 atol = 1.0e-3 + end + + # 3. second-order convergence on a small full-rank system (isolates the temporal order) + @testset "second-order convergence" begin + Random.seed!(2) + Lc = 4 + Hc = force_planar(transverse_field_ising(ComplexF64, Trivial; L = Lc)) + ψ_full = FiniteMPS(rand, ComplexF64, Lc, ℙ^2, ℙ^4) # full-rank: 1,2,4,2,1 + + Hmat = convert(TensorMap, Hc) + ψvec = convert(TensorMap, ψ_full) + ψvec /= norm(ψvec) + + T = 0.5 + dts = [0.1, 0.05, 0.025] + errs = map(dts) do δt + n = round(Int, T / δt) + ref = exp(-im * Hmat * (n * δt)) * ψvec + ψ = copy(ψ_full) + envs = environments(ψ, Hc, ψ) + for k in 0:(n - 1) + timestep!(ψ, Hc, k * δt, δt, BUG(), envs) + end + ψout = convert(TensorMap, ψ) + ψout /= norm(ψout) + return 1 - abs(dot(ψout, ref)) + end + + slopes = [ + (log(errs[i + 1]) - log(errs[i])) / (log(dts[i + 1]) - log(dts[i])) + for i in 1:(length(dts) - 1) + ] + @info "BUG convergence" errs slopes + for s in slopes + @test s ≈ 2 atol = 0.3 + end + end + + # 4. imaginary-time evolution lowers the energy toward the ground state (and, having no + # backward substep, stays norm-preserving/stable) + @testset "imaginary-time lowers energy" begin + Random.seed!(5) + ψi = complex(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^4)) + E_start = real(expectation_value(ψi, H)) + E_prev = E_start + for _ in 1:20 + ψi, = timestep(ψi, H, 0.0, 0.1, BUG(); imaginary_evolution = true, normalize = true) + E_now = real(expectation_value(ψi, H)) + @test E_now ≤ E_prev + 1.0e-6 # monotone (non-increasing) energy + E_prev = E_now + end + @test E_prev < E_start - 1.0 # substantial lowering toward the ground state + @test norm(ψi) ≈ 1 atol = 1.0e-6 # imaginary-time BUG renormalizes each step + end + + # 5. LazySum / MultipliedOperator smoke tests + @testset "LazySum" begin + Hlazy = LazySum([3 * H, 1.55 * H, -0.1 * H]) + ψl, envs = timestep(ψ₀, Hlazy, 0.0, dt, BUG()) + E = expectation_value(ψl, Hlazy, envs) + @test (3 + 1.55 - 0.1) * E₀ ≈ E atol = 1.0e-2 + end + + @testset "TimeDependent LazySum" begin + Ht = MultipliedOperator(H, t -> 4) + MultipliedOperator(H, 1.45) + ψa, envsa = timestep(ψ₀, Ht(1.0), 0.0, dt, BUG()) + Ea = expectation_value(ψa, Ht(1.0), envsa) + + ψt, envst = timestep(ψ₀, Ht, 1.0, dt, BUG()) + Et = expectation_value(ψt, Ht(1.0), envst) + @test Ea ≈ Et atol = 1.0e-8 + end +end + +# Rank-adaptive BUG (Stage 2): a truncating `trscheme` enables basis augmentation + truncation, so +# the bond dimension grows and shrinks to track the entanglement of the evolving state. Trivial +# tensors only (symmetry stress is Chunk 2.3). The default `notrunc()` regression is covered by the +# fixed-rank testsets above (they must all still pass unchanged). +@testset "BUG rank-adaptive" verbose = true begin + # 1. bond growth: a tight tolerance grows a low-bond-dim state, a looser tolerance keeps it + # smaller. + @testset "bond growth" begin + Random.seed!(101) + L = 6 + H = force_planar(transverse_field_ising(ComplexF64, Trivial; L)) + ψ₀ = complex(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^2)) # low bond dim (2) + normalize!(ψ₀) + Dstart = maximum(dim(left_virtualspace(ψ₀, k)) for k in 1:L) + + ψtight = ψ₀ + for k in 0:2 + ψtight, = timestep(ψtight, H, k * 0.05, 0.05, BUG(; trscheme = truncerror(; atol = 1.0e-10))) + end + Dtight = maximum(dim(left_virtualspace(ψtight, k)) for k in 1:L) + + ψloose = ψ₀ + for k in 0:2 + ψloose, = timestep(ψloose, H, k * 0.05, 0.05, BUG(; trscheme = truncerror(; atol = 1.0e-2))) + end + Dloose = maximum(dim(left_virtualspace(ψloose, k)) for k in 1:L) + + @info "BUG rank-adaptive bond growth" Dstart Dloose Dtight + @test Dtight > Dstart # rank-adaptivity grows the bond + @test Dloose < Dtight # a looser tolerance keeps a smaller bond + end + + # 2. THE HARD GATE: overlap-error vs the dense exp(-iHT) reference decreases (monotonically, up + # to the plateau at the fixed-dt floor) as the truncation tolerance ϑ shrinks. This proves the + # augmentation actually captures the true dynamics. + @testset "accuracy improves as ϑ decreases" begin + Random.seed!(202) + L = 6 + H = force_planar(transverse_field_ising(ComplexF64, Trivial; L)) + ψ₀ = complex(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^2)) # low-rank start the dynamics grows + normalize!(ψ₀) + + Hmat = convert(TensorMap, H) + ψvec = convert(TensorMap, ψ₀) + ψvec /= norm(ψvec) + + T = 0.2 + dt = 0.05 + n = round(Int, T / dt) + ref = exp(-im * Hmat * (n * dt)) * ψvec + + ϑs = [1.0e-2, 1.0e-4, 1.0e-6] + errs = map(ϑs) do ϑ + alg = BUG(; trscheme = truncerror(; atol = ϑ)) + ψ = copy(ψ₀) + envs = environments(ψ, H, ψ) + for k in 0:(n - 1) + timestep!(ψ, H, k * dt, dt, alg, envs) + end + ψout = convert(TensorMap, ψ) + ψout /= norm(ψout) + return 1 - abs(dot(ψout, ref)) + end + + @info "BUG rank-adaptive accuracy vs ϑ" ϑs errs + for i in 1:(length(ϑs) - 1) + @test errs[i + 1] ≤ 1.5 * errs[i] # monotone within plateau noise near the dt-floor + end + @test errs[end] < errs[1] / 10 # clear net improvement toward the dt-floor + end + + # 3. CBE-style comparison: from a low-rank state, rank-adaptive BUG tracks a bond-adaptive TDVP2 + # reference better than fixed-rank `BUG()` does (mirrors the CBE-TDVP test). + @testset "tracks TDVP2 better than fixed-rank BUG" begin + Random.seed!(303) + L = 8 + H = force_planar(heisenberg_XXX(Float64, Trivial; spin = 1 // 2, L)) + Dstart, Dcap, dt = 2, 16, 0.05 + ψ₀ = complex(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^Dstart)) + + ref, adaptive, fixed = ψ₀, ψ₀, ψ₀ + for _ in 1:6 + ref, = timestep(ref, H, 0.0, dt, TDVP2(; trscheme = truncrank(Dcap))) + adaptive, = timestep(adaptive, H, 0.0, dt, BUG(; trscheme = truncerror(; atol = 1.0e-8))) + fixed, = timestep(fixed, H, 0.0, dt, BUG()) + end + + @test dim(left_virtualspace(adaptive, L ÷ 2)) > Dstart # adaptive grew the bond + @test dim(left_virtualspace(fixed, L ÷ 2)) == Dstart # fixed-rank stuck at Dstart + @test abs(dot(ref, adaptive)) > abs(dot(ref, fixed)) # and tracks the reference better + end + + # 4. imaginary-time ground-state search: from a low bond dim, rank-adaptive imaginary-time BUG + # grows the bond and lowers the energy toward the true ground state. + @testset "imaginary-time grows bond and lowers energy" begin + Random.seed!(404) + L = 8 + H = force_planar(heisenberg_XXX(Float64, Trivial; spin = 1 // 2, L)) + ψgs, = find_groundstate(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^16), H; verbosity = 0) + Egs = real(expectation_value(ψgs, H)) + + ψ = complex(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^2)) # low bond dim + Dstart = maximum(dim(left_virtualspace(ψ, k)) for k in 1:L) + E_start = real(expectation_value(ψ, H)) + for _ in 1:30 + ψ, = timestep(ψ, H, 0.0, 0.1, BUG(; trscheme = truncerror(; atol = 1.0e-8)); imaginary_evolution = true, normalize = true) + end + Dend = maximum(dim(left_virtualspace(ψ, k)) for k in 1:L) + E_end = real(expectation_value(ψ, H)) + + @info "BUG rank-adaptive imaginary-time" Dstart Dend E_start E_end Egs + @test Dend > Dstart # the bond grew as entanglement built up + @test E_end < E_start - 1.0 # substantial lowering + @test E_end ≈ Egs atol = 0.6 # toward the true ground state (loose) + @test norm(ψ) ≈ 1 atol = 1.0e-6 # imaginary-time renormalizes each step + end + + # 5. real-time energy conservation with truncation (loose atol) and norm behaviour: for a small + # tolerance the truncation is negligible, so the norm stays ≈ 1. + @testset "real-time energy conservation and norm" begin + Random.seed!(505) + L = 6 + H = force_planar(transverse_field_ising(ComplexF64, Trivial; L)) + ψ₀ = complex(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^4)) + normalize!(ψ₀) + E₀ = real(expectation_value(ψ₀, H)) + + ψ = ψ₀ + for k in 0:4 + ψ, = timestep(ψ, H, k * 0.05, 0.05, BUG(; trscheme = truncerror(; atol = 1.0e-8))) + end + @test real(expectation_value(ψ, H)) ≈ E₀ atol = 1.0e-2 # energy conserved (loose) + @test norm(ψ) ≈ 1 atol = 1.0e-6 # tiny ϑ ⇒ norm preserved + end +end + +# Charge-sector (symmetric-tensor) coverage for the fixed-rank BUG. These use *genuine* +# symmetric tensors (no `force_planar`), exercising the graded-bond paths flagged in the design +# doc's hsector risk register (H1/H6/H7): the transport-tensor seed `isomorphism(V ← V)`, the +# (co)domain/dual conventions of the inlined transport contractions, and the adjoints carrying +# sector duals. A fixed-rank step must preserve the total charge and the graded structure of every bond. +@testset "BUG symmetric tensors" verbose = true begin + dt = 0.1 + L = 6 + + # 1. U(1)-symmetric Heisenberg, both in the natural total-Sz = 0 sector and in a fixed nonzero + # total-charge (Sz = 1) sector: energy conservation + eigenstate phase + sector preservation. + @testset "U(1) Heisenberg (total Sz = $label)" for (label, right) in + (("0", U1Space(0 => 1)), ("1", U1Space(1 => 1))) + Random.seed!(2718) + H = heisenberg_XXX(ComplexF64, U1Irrep; spin = 1 // 2, L) + maxV = MPSKit.max_virtualspaces(physicalspace(H)) + ψ = FiniteMPS(physicalspace(H), maxV[2:(end - 1)]; right) + ψ₀, = find_groundstate(ψ, H; verbosity = 0) + E₀ = expectation_value(ψ₀, H) + + Vl₀ = left_virtualspace.(Ref(ψ₀), 1:L) + Vr₀ = right_virtualspace.(Ref(ψ₀), 1:L) + + ψ1, envs = timestep(ψ₀, H, 0.0, dt, BUG()) + E1 = expectation_value(ψ1, H, envs) + + @test E₀ ≈ E1 atol = 1.0e-2 + @test imag(E1) ≈ 0 atol = 1.0e-8 + @test dot(ψ1, ψ₀) ≈ exp(im * dt * E₀) atol = 1.0e-4 + # the fixed-rank step preserves the graded structure (sector content) of every bond + @test left_virtualspace.(Ref(ψ1), 1:L) == Vl₀ + @test right_virtualspace.(Ref(ψ1), 1:L) == Vr₀ + end + + # 2. A second symmetry group. Z2 (transverse-field Ising) and SU2 (Heisenberg) both stress the + # graded transport tensor; same assertions (energy conservation + eigenstate phase). + @testset "Z2 transverse-field Ising" begin + Random.seed!(161803) + H = transverse_field_ising(ComplexF64, Z2Irrep; g = 1.0, L) + ψ = FiniteMPS(physicalspace(H), Z2Space(0 => 4, 1 => 4)) + ψ₀, = find_groundstate(ψ, H; verbosity = 0) + E₀ = expectation_value(ψ₀, H) + Vl₀ = left_virtualspace.(Ref(ψ₀), 1:L) + + ψ1, envs = timestep(ψ₀, H, 0.0, dt, BUG()) + E1 = expectation_value(ψ1, H, envs) + @test E₀ ≈ E1 atol = 1.0e-2 + @test dot(ψ1, ψ₀) ≈ exp(im * dt * E₀) atol = 1.0e-4 + @test left_virtualspace.(Ref(ψ1), 1:L) == Vl₀ + end + + @testset "SU(2) Heisenberg" begin + Random.seed!(577215) + H = heisenberg_XXX(ComplexF64, SU2Irrep; spin = 1 // 2, L) + # SU(2) spin-1/2 bonds alternate between integer / half-integer spins, so use the + # model's own full-rank virtual spaces rather than a hand-picked (integer-only) space. + maxV = MPSKit.max_virtualspaces(physicalspace(H)) + ψ = FiniteMPS(physicalspace(H), maxV[2:(end - 1)]) + ψ₀, = find_groundstate(ψ, H; verbosity = 0) + E₀ = expectation_value(ψ₀, H) + Vl₀ = left_virtualspace.(Ref(ψ₀), 1:L) + + ψ1, envs = timestep(ψ₀, H, 0.0, dt, BUG()) + E1 = expectation_value(ψ1, H, envs) + @test E₀ ≈ E1 atol = 1.0e-2 + @test dot(ψ1, ψ₀) ≈ exp(im * dt * E₀) atol = 1.0e-4 + @test left_virtualspace.(Ref(ψ1), 1:L) == Vl₀ + end + + # 3. Imaginary-time symmetric evolution lowers the energy while preserving the sector + norm. + @testset "imaginary-time (U(1))" begin + Random.seed!(141421) + H = heisenberg_XXX(ComplexF64, U1Irrep; spin = 1 // 2, L) + maxV = MPSKit.max_virtualspaces(physicalspace(H)) + ψi = FiniteMPS(physicalspace(H), maxV[2:(end - 1)]) + Vl₀ = left_virtualspace.(Ref(ψi), 1:L) + Vr₀ = right_virtualspace.(Ref(ψi), 1:L) + E_start = real(expectation_value(ψi, H)) + E_prev = E_start + for _ in 1:15 + ψi, = timestep(ψi, H, 0.0, 0.1, BUG(); imaginary_evolution = true, normalize = true) + E_now = real(expectation_value(ψi, H)) + @test E_now ≤ E_prev + 1.0e-6 # monotone (non-increasing) energy + E_prev = E_now + end + @test E_prev < E_start - 0.5 # substantial lowering toward the ground state + @test norm(ψi) ≈ 1 atol = 1.0e-6 # imaginary-time BUG renormalizes each step + # the sector content of every bond is preserved throughout the imaginary-time sweep + @test left_virtualspace.(Ref(ψi), 1:L) == Vl₀ + @test right_virtualspace.(Ref(ψi), 1:L) == Vr₀ + end +end + +# Charge-sector RANK-ADAPTIVE coverage (Chunk 2.3): the "sector-adaptivity" path — genuine +# symmetric tensors under a *truncating* `trscheme`, so bonds grow and shrink per sector. This is +# the hardest part of the design doc; it stresses the H2–H5/H10 pitfalls of the risk register: +# augmentation as a per-sector direct sum (H2/H10), the "old basis first" invariant per sector +# (H3), a global-ϑ truncation dropping a sector to dimension 0 ⇒ dynamic bond grading (H4), and +# total-boundary-charge conservation through the zero-block embeddings (H5). All states start +# genuinely low-rank so the dynamics *must* grow the bonds; no `force_planar`. +@testset "BUG rank-adaptive symmetric" verbose = true begin + # A genuinely low-rank U(1) start. For spin-1/2 the virtual bonds alternate integer / + # half-integer parity, so a single-sector cap collapses the state (no fusion channels); use a + # mixed-parity cap with multiplicity 1 per sector. The full-rank Sz=0 profile is [1,2,4,8,4,2]; + # this cap gives [1,2,3,2,3,2], leaving room (esp. on the middle bond) for per-sector growth. + u1cap = U1Space(-1 // 2 => 1, 1 // 2 => 1, 0 => 1, 1 => 1, -1 => 1) + function low_rank_u1(L; seed = 2718, right = U1Space(0 => 1)) + Random.seed!(seed) + H = heisenberg_XXX(ComplexF64, U1Irrep; spin = 1 // 2, L) + ψ = FiniteMPS(physicalspace(H), u1cap; right) + normalize!(ψ) + return H, ψ + end + maxbond(ψ) = maximum(dim(left_virtualspace(ψ, k)) for k in 1:length(ψ)) + secmults(V) = [c => dim(V, c) for c in sectors(V)] + + # 1. THE CORE GATE: rank-adaptivity under symmetry grows the bond per sector for a tight ϑ and + # keeps it small for a loose ϑ, without any `SpaceMismatch`, while the total boundary charge + # (the fixed `right` virtual space) is preserved. Checked in the natural total-Sz=0 sector and + # in a fixed nonzero-charge (Sz=1) sector. + @testset "rank-adaptivity + total-charge preservation (Sz = $label)" for (label, right) in + (("0", U1Space(0 => 1)), ("1", U1Space(1 => 1))) + L = 6 + H, ψ₀ = low_rank_u1(L; right) + Rtot = right_virtualspace(ψ₀, L) + Dstart = maxbond(ψ₀) + + ψtight = ψ₀ + for k in 0:2 + ψtight, = timestep(ψtight, H, k * 0.05, 0.05, BUG(; trscheme = truncerror(; atol = 1.0e-10))) + end + Dtight = maxbond(ψtight) + + ψloose = ψ₀ + for k in 0:2 + ψloose, = timestep(ψloose, H, k * 0.05, 0.05, BUG(; trscheme = truncerror(; atol = 1.0e-2))) + end + Dloose = maxbond(ψloose) + + @info "BUG rank-adaptive symmetric bond growth (Sz=$label)" Dstart Dloose Dtight + @test Dtight > Dstart # per-sector augmentation grows the bond + @test Dloose < Dtight # a looser tolerance keeps a smaller bond + @test right_virtualspace(ψtight, L) == Rtot # H5: total boundary charge preserved + @test right_virtualspace(ψloose, L) == Rtot + @test norm(ψtight) ≈ 1 atol = 1.0e-6 # tiny ϑ ⇒ negligible truncation + end + + # 2. THE HARD GATE: overlap-error vs a *dense* `exp(-iH·T)` reference decreases as ϑ shrinks + # (monotonically up to the fixed-dt plateau). This proves per-sector augmentation+truncation + # actually captures the true dynamics under symmetry, not just that some bond grows. + @testset "accuracy improves as ϑ decreases (U(1))" begin + L = 6 + H, ψ₀ = low_rank_u1(L; seed = 202) + + Hmat = convert(TensorMap, H) + ψvec = convert(TensorMap, ψ₀) + ψvec /= norm(ψvec) + + T = 0.2 + dt = 0.05 + n = round(Int, T / dt) + ref = exp(-im * Hmat * (n * dt)) * ψvec + + ϑs = [1.0e-2, 1.0e-4, 1.0e-6] + errs = map(ϑs) do ϑ + ψ = copy(ψ₀) + for k in 0:(n - 1) + ψ, = timestep(ψ, H, k * dt, dt, BUG(; trscheme = truncerror(; atol = ϑ))) + end + ψout = convert(TensorMap, ψ) + ψout /= norm(ψout) + return 1 - abs(dot(ψout, ref)) + end + + @info "BUG rank-adaptive symmetric accuracy vs ϑ" ϑs errs + for i in 1:(length(ϑs) - 1) + @test errs[i + 1] ≤ 1.5 * errs[i] # monotone within plateau noise near the dt-floor + end + # Clear net improvement toward the dt-floor. This charge-capped low-rank system saturates + # its (small) per-sector bonds quickly: at ϑ ≲ 1e-4 no directions are truncated, so the two + # tightest tolerances land on the *identical* pure-dt-discretization floor (the ϑ=1e-4 and + # 1e-6 errors coincide to the digit). The net truncation-limited improvement is therefore + # modest (≈2.5×) rather than the orders of magnitude a full-rank trivial system shows. + @test errs[end] < errs[1] / 2 + end + + # 3. DYNAMIC BOND GRADING: an interior bond's per-sector multiplicities + # `[dim(V, c) for c in sectors(V)]` are time-dependent — they change across a single + # rank-adaptive step (the graded structure is not fixed), while the total charge is. + @testset "dynamic bond grading (per-sector multiplicities are time-dependent)" begin + L = 6 + H, ψ₀ = low_rank_u1(L) + mid = 4 + before = secmults(left_virtualspace(ψ₀, mid)) + ψ, = timestep(ψ₀, H, 0.0, 0.05, BUG(; trscheme = truncerror(; atol = 1.0e-10))) + after = secmults(left_virtualspace(ψ, mid)) + @info "BUG rank-adaptive symmetric dynamic grading" before after + @test before != after # graded structure evolved + @test right_virtualspace(ψ, L) == right_virtualspace(ψ₀, L) # ... at fixed total charge + end + + # 4. H4 — a global-ϑ cut can truncate a whole sector to dimension 0 (removing it from the bond), + # and downstream rank-adaptive steps must tolerate the changed / asymmetric grading. + # + # NOTE on the sharpest H4 test (deterministic *drop-and-re-add* in one run): this is NOT + # achievable with single-site BUG, and asserting it would be wrong. The augmentation candidate + # `ACᵢ` in each half-sweep always carries the (already-truncated) bond as its *domain*, so + # `_bug_augment_left`/`_bug_augment_right` can only append directions whose charge sectors are + # already present on that bond — they grow multiplicity within existing sectors but can never + # re-introduce a sector that was truncated away (the graded analog of single-site TDVP's + # inability to change bond quantum numbers; re-adding a sector needs a two-site / CBE-style + # candidate, cf. `OptimalExpand`). We verified empirically that a sector dropped from an + # interior bond does not reappear over many subsequent tight steps. We therefore assert the + # deterministic *drop* and the H4 tolerance requirement (subsequent steps run without + # `SpaceMismatch`, yield a valid normalizable state, and preserve the total charge). + @testset "sector drop-to-zero + dynamic-grading tolerance (H4)" begin + L = 6 + H, ψ₀ = low_rank_u1(L) + # grow a rich interior bond that carries the subdominant ±1 sectors + ψrich = ψ₀ + for k in 0:5 + ψrich, = timestep(ψrich, H, k * 0.05, 0.05, BUG(; trscheme = truncerror(; atol = 1.0e-10))) + end + mid = 3 + Rtot = right_virtualspace(ψrich, L) + @test U1Irrep(-1) in sectors(left_virtualspace(ψrich, mid)) # subdominant sector present + + # a global rank cut pools singular values across all sectors under one threshold, so the + # subdominant sector is truncated to dimension 0 and removed from the bond (H4). + ψdrop, = timestep(ψrich, H, 0.0, 0.05, BUG(; trscheme = truncrank(2))) + @test !(U1Irrep(-1) in sectors(left_virtualspace(ψdrop, mid))) # dropped to dim 0 + @test right_virtualspace(ψdrop, L) == Rtot # ... charge still preserved + + # subsequent rank-adaptive steps must tolerate the reduced / asymmetric grading: no + # SpaceMismatch, a valid normalizable state, and a conserved total charge. + ψcont = ψdrop + for k in 0:3 + ψcont, = timestep(ψcont, H, k * 0.05, 0.05, BUG(; trscheme = truncerror(; atol = 1.0e-10))) + end + @info "BUG rank-adaptive symmetric H4" dropped = sectors(left_virtualspace(ψdrop, mid)) continued = sectors(left_virtualspace(ψcont, mid)) + @test isfinite(real(expectation_value(ψcont, H))) + @test norm(ψcont) > 0 + @test right_virtualspace(ψcont, L) == Rtot + end + + # 5. IMAGINARY-TIME symmetric rank-adaptive ground-state search: from a low bond dim it grows the + # per-sector bonds, lowers the energy toward `find_groundstate`, and preserves both the total + # charge and the (renormalized) norm. + @testset "imaginary-time symmetric grows bond and lowers energy" begin + L = 6 + H, ψ₀ = low_rank_u1(L) + maxV = MPSKit.max_virtualspaces(physicalspace(H)) + ψgs, = find_groundstate(FiniteMPS(physicalspace(H), maxV[2:(end - 1)]), H; verbosity = 0) + Egs = real(expectation_value(ψgs, H)) + + Rtot = right_virtualspace(ψ₀, L) + Dstart = maxbond(ψ₀) + E_start = real(expectation_value(ψ₀, H)) + + ψ = ψ₀ + for _ in 1:30 + ψ, = timestep(ψ, H, 0.0, 0.1, BUG(; trscheme = truncerror(; atol = 1.0e-8)); imaginary_evolution = true, normalize = true) + end + Dend = maxbond(ψ) + E_end = real(expectation_value(ψ, H)) + + @info "BUG rank-adaptive symmetric imaginary-time" Dstart Dend E_start E_end Egs + @test Dend > Dstart # per-sector bonds grew as entanglement built up + @test E_end < E_start - 1.0 # substantial lowering + @test E_end ≈ Egs atol = 0.1 # toward the true ground state + @test norm(ψ) ≈ 1 atol = 1.0e-6 # imaginary-time renormalizes each step + @test right_virtualspace(ψ, L) == Rtot # total charge preserved throughout + end +end diff --git a/test/algorithms/bug_augment.jl b/test/algorithms/bug_augment.jl new file mode 100644 index 000000000..947d2c378 --- /dev/null +++ b/test/algorithms/bug_augment.jl @@ -0,0 +1,136 @@ +println(" +------------------------------------ +| BUG basis-augmentation tests | +------------------------------------ +") + +using .TestSetup +using Test, TestExtras +using MPSKit +using MPSKit: _bug_augment_left, _bug_augment_right, _transpose_tail, _transpose_front, + left_orth, right_gauge +using TensorKit +using TensorKit: ℙ +using LinearAlgebra: I, norm +using Random + +# The augment helpers keep the OLD isometry `U₀` as the leading per-sector block and append the +# component of the evolved candidate `K₁` that is orthogonal to it (no truncation). The four core +# properties (checked below for both sweep directions, trivial + U(1)): +# 1. isometry `Û' Û ≈ 𝟙` +# 2. old-first, per sector `M = Û' U₀ ≈ [𝟙; 0]` block-by-block +# 3. range ⊇ old, candidate `Û (Û' U₀) ≈ U₀` and `Û (Û' K₁) ≈ K₁` +# 4. rank growth `dim(Vr₀) < dim(V̂) ≤ 2·dim(Vr₀)` + +# Assert the per-sector "old-first" invariant on the overlap `M` (codomain = augmented bond, +# domain = old bond `V_old`): every sector block must be `[𝟙; 0]`, i.e. the leading +# `dim(V_old, c)` rows are the identity and the rest vanish. Iterating `blocks(M)` visits exactly +# the sectors common to `M`'s (co)domain, so purely-new sectors (absent from `V_old`) are correctly +# skipped here. +function _check_old_first(M, V_old; tol = 1.0e-10) + for (c, b) in blocks(M) + r0 = dim(V_old, c) + @test b[1:r0, :] ≈ I + if r0 < size(b, 1) + @test norm(b[(r0 + 1):end, :]) < tol + end + end + return nothing +end + +function check_augment_left(U₀, K₁; tol = 1.0e-10) + Û = _bug_augment_left(U₀, K₁) + M = Û' * U₀ # old bond's coordinates in the augmented basis (V̂ ← Vr₀) + # 1. isometry + @test Û' * Û ≈ one(Û' * Û) + # 2. old-first, per sector (M : V̂ ← Vr₀) + _check_old_first(M, domain(U₀)[1]; tol) + # 3. range captures both the old basis and the evolved candidate + @test Û * (Û' * U₀) ≈ U₀ + @test Û * (Û' * K₁) ≈ K₁ + @test Û * M ≈ U₀ # consistency of the returned overlap + # 4. rank growth: strictly bigger than the old bond, at most doubled + r = dim(domain(U₀)) + @test r < dim(domain(Û)) ≤ 2r + return Û, M +end + +function check_augment_right(U₀, K₁; tol = 1.0e-10) + Û = _bug_augment_right(U₀, K₁) + ût = _transpose_tail(Û) # V̂ ← P ⊗ Vr, right-isometric (row space) + u0t = _transpose_tail(U₀) + k1t = _transpose_tail(K₁) + M = ût * u0t' # old bond's coordinates in the augmented basis (V̂ ← Vl₀) + # 1. isometry (right-canonical ⇒ tail has orthonormal rows) + @test ût * ût' ≈ one(ût * ût') + # 2. old-first, per sector (M : V̂ ← Vl₀) + _check_old_first(M, codomain(u0t)[1]; tol) + # 3. row space captures both the old basis and the evolved candidate + @test (u0t * ût') * ût ≈ u0t + @test (k1t * ût') * ût ≈ k1t + @test M' * ût ≈ u0t # consistency of the returned overlap + # 4. rank growth on the left bond + r = dim(codomain(u0t)) + @test r < dim(codomain(ût)) ≤ 2r + return Û, M +end + +@testset "BUG basis augmentation" verbose = true begin + # ------------------------------------------------------------------------------------------- + # Trivial (dense) tensors, ComplexF64 to exercise the adjoints. + # ------------------------------------------------------------------------------------------- + @testset "trivial tensors" begin + Random.seed!(20260707) + Vl = ℂ^2 + P = ℂ^3 + Vr = ℂ^2 + + # left→right: augment the RIGHT bond (domain of a left-isometry) + U₀_L, _ = left_orth(randn(ComplexF64, Vl ⊗ P ← Vr)) # Vl⊗P ← Vr, left-isometric + K₁_L = randn(ComplexF64, Vl ⊗ P ← Vr) # generic evolved candidate + Û_L, _ = check_augment_left(U₀_L, K₁_L) + + # right→left: augment the LEFT bond (codomain of a right-isometry) + _, U₀_R = right_gauge(randn(ComplexF64, Vl ⊗ P ← Vr)) # V ⊗ P ← Vr, right-isometric + K₁_R = randn(ComplexF64, Vl ⊗ P ← Vr) + Û_R, _ = check_augment_right(U₀_R, K₁_R) + end + + # ------------------------------------------------------------------------------------------- + # U(1)-symmetric tensors: per-sector direct sums, old-first sector-by-sector, and a genuinely + # new sector introduced by the candidate on the left augment. + # ------------------------------------------------------------------------------------------- + @testset "U(1) symmetric tensors" begin + Random.seed!(31415926) + + # left→right, WITH a new sector: `Vr₀` deliberately omits sector 0, which `Vl ⊗ P` + # contains and the candidate `K₁` populates ⇒ augmentation must add it to `V̂`. + Vl = U1Space(0 => 1, 1 => 1) + P = U1Space(0 => 1, 1 => 1) # Vl⊗P : 0=>1, 1=>2, 2=>1 + Vr0 = U1Space(1 => 1, 2 => 1) # omits sector 0 + Vr0K = U1Space(0 => 1, 1 => 1, 2 => 1) # candidate populates sector 0 + + U₀_L, _ = left_orth(randn(ComplexF64, Vl ⊗ P ← Vr0)) + @test domain(U₀_L)[1] == Vr0 # left_orth kept the full old bond + K₁_L = randn(ComplexF64, Vl ⊗ P ← Vr0K) + Û_L, _ = check_augment_left(U₀_L, K₁_L) + + # the augmented bond is a per-sector direct sum that gained sector 0 + V̂_L = domain(Û_L)[1] + @test !(U1Irrep(0) in sectors(Vr0)) + @test U1Irrep(0) in sectors(V̂_L) + # every old sector survives with at least its old multiplicity (old-first ⇒ nothing dropped) + for c in sectors(Vr0) + @test dim(V̂_L, c) ≥ dim(Vr0, c) + end + + # right→left augment on a compatible symmetric configuration + Plr = U1Space(0 => 1, 1 => 1) + Vrr = U1Space(0 => 1, 1 => 1) + Vl0 = U1Space(0 => 1, 1 => 1) + _, U₀_R = right_gauge(randn(ComplexF64, Vl0 ⊗ Plr ← Vrr)) + Vl_K = U1Space(0 => 1, 1 => 1) + K₁_R = randn(ComplexF64, Vl_K ⊗ Plr ← Vrr) + Û_R, _ = check_augment_right(U₀_R, K₁_R) + end +end diff --git a/test/algorithms/parallelbug.jl b/test/algorithms/parallelbug.jl new file mode 100644 index 000000000..91a5aa1f2 --- /dev/null +++ b/test/algorithms/parallelbug.jl @@ -0,0 +1,446 @@ +println(" +------------------------------------- +| ParallelBUG time-stepping tests | +------------------------------------- +") + +using .TestSetup +using Test, TestExtras +using MPSKit +using TensorKit +using TensorKit: ℙ +using LinearAlgebra: dot, norm +using Random + +# NOTE (experimental). `ParallelBUG` is the parallel Basis-Update & Galerkin integrator +# (Ceruti et al. 2024) specialized to the caterpillar `FiniteMPS`. It reproduces the exact matrix +# parallel-BUG step for two sites, conserves energy / the eigenstate phase exactly (amplitude carried +# once, at the root), grows bonds adaptively, agrees with `TDVP` over short times, and converges at +# (at least) the documented first order in `dt` toward the dense reference. +@testset "ParallelBUG time evolution" verbose = true begin + dt = 0.1 + L = 6 + + H = force_planar(heisenberg_XXX(Float64, Trivial; spin = 1 // 2, L)) + ψ = FiniteMPS(rand, Float64, L, ℙ^2, ℙ^4) + ψ₀, = find_groundstate(ψ, H; verbosity = 0) + E₀ = expectation_value(ψ₀, H) + + # 1. energy conservation + eigenstate phase (amplitude is carried exactly once, at the root) + @testset "energy conservation + eigenstate phase" begin + ψ1, envs = timestep(ψ₀, H, 0.0, dt, ParallelBUG(; trscheme = truncerror(; atol = 1.0e-12))) + E1 = expectation_value(ψ1, H, envs) + @test E₀ ≈ E1 atol = 1.0e-2 + @test dot(ψ1, ψ₀) ≈ exp(im * dt * E₀) atol = 1.0e-4 + end + + # 2. two sites: exact reproduction of the dense exp(-iH·dt) step (locks the block conventions) + @testset "two-site exactness vs dense reference" begin + Random.seed!(2) + Lc = 2 + Hc = force_planar(transverse_field_ising(ComplexF64, Trivial; L = Lc)) + ψc = FiniteMPS(rand, ComplexF64, Lc, ℙ^2, ℙ^4) + Hmat = convert(TensorMap, Hc) + ψvec = convert(TensorMap, ψc); ψvec /= norm(ψvec) + ref = exp(-im * Hmat * 0.05) * ψvec + ψ1, = timestep(ψc, Hc, 0.0, 0.05, ParallelBUG(; trscheme = truncerror(; atol = 1.0e-12))) + out = convert(TensorMap, ψ1); out /= norm(out) + @test 1 - abs(dot(out, ref)) < 1.0e-10 + end + + # 3. agreement with TDVP over a few short real-time steps of a random MPS + @testset "agreement with TDVP" begin + Random.seed!(1234) + ψr = complex(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^4)) + δt = 0.01 + ψ_p, ψ_tdvp = ψr, ψr + for k in 0:4 + ψ_p, = timestep(ψ_p, H, k * δt, δt, ParallelBUG(; trscheme = truncerror(; atol = 1.0e-12))) + ψ_tdvp, = timestep(ψ_tdvp, H, k * δt, δt, TDVP()) + end + @test expectation_value(ψ_p, H) ≈ expectation_value(ψ_tdvp, H) atol = 1.0e-3 + @test abs(dot(ψ_p, ψ_tdvp)) ≈ 1 atol = 1.0e-3 + end + + # 4. imaginary-time evolution lowers the energy monotonically and stays normalized + @testset "imaginary-time lowers energy" begin + Random.seed!(5) + ψi = complex(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^4)) + E_start = real(expectation_value(ψi, H)) + E_prev = E_start + for _ in 1:20 + ψi, = timestep( + ψi, H, 0.0, 0.1, ParallelBUG(; trscheme = truncerror(; atol = 1.0e-10)); + imaginary_evolution = true, normalize = true + ) + E_now = real(expectation_value(ψi, H)) + @test E_now ≤ E_prev + 1.0e-6 # monotone (non-increasing) energy + E_prev = E_now + end + @test E_prev < E_start - 1.0 # substantial lowering toward the ground state + @test norm(ψi) ≈ 1 atol = 1.0e-6 # imaginary-time renormalizes each step + end + + # 5. rank adaptivity: a low-bond-dim start grows under a tight tolerance and stays small under a + # loose one (the augment-then-SvdCut mechanism injects the new directions the evolution finds). + @testset "bond growth" begin + Random.seed!(3) + Lg = 6 + Hg = force_planar(transverse_field_ising(ComplexF64, Trivial; L = Lg)) + ψg = normalize!(complex(FiniteMPS(rand, Float64, Lg, ℙ^2, ℙ^2))) # bond dim 2 + Dstart = maximum(dim(right_virtualspace(ψg, b)) for b in 1:(Lg - 1)) + + ψtight = ψg + for k in 0:4 + ψtight, = timestep(ψtight, Hg, k * 0.05, 0.05, ParallelBUG(; trscheme = truncerror(; atol = 1.0e-8))) + end + Dtight = maximum(dim(right_virtualspace(ψtight, b)) for b in 1:(Lg - 1)) + + ψloose = ψg + for k in 0:4 + ψloose, = timestep(ψloose, Hg, k * 0.05, 0.05, ParallelBUG(; trscheme = truncerror(; atol = 1.0e-2))) + end + Dloose = maximum(dim(right_virtualspace(ψloose, b)) for b in 1:(Lg - 1)) + + @info "ParallelBUG bond growth" Dstart Dloose Dtight + @test Dtight > Dstart # rank-adaptivity grows the bond + @test Dloose < Dtight # a looser tolerance keeps a smaller bond + end + + # 6. LazySum smoke test + @testset "LazySum" begin + Hlazy = LazySum([3 * H, 1.55 * H, -0.1 * H]) + ψl, envs = timestep(ψ₀, Hlazy, 0.0, dt, ParallelBUG(; trscheme = truncerror(; atol = 1.0e-12))) + E = expectation_value(ψl, Hlazy, envs) + @test (3 + 1.55 - 0.1) * E₀ ≈ E atol = 1.0e-2 + end + + # 6b. TimeDependent operator smoke test: exercises the `TimedOperator` coupling path (which + # has no one-arg apply). A constant-in-`t` coefficient makes midpoint-freezing exact, so + # evolving the time-dependent `Ht` at `t=1.0` must match evolving the pre-evaluated `Ht(1.0)`. + @testset "TimeDependent LazySum" begin + Ht = MultipliedOperator(H, t -> 4) + MultipliedOperator(H, 1.45) + alg = ParallelBUG(; trscheme = truncerror(; atol = 1.0e-12)) + ψa, envsa = timestep(ψ₀, Ht(1.0), 0.0, dt, alg) + Ea = expectation_value(ψa, Ht(1.0), envsa) + + ψt, envst = timestep(ψ₀, Ht, 1.0, dt, alg) + Et = expectation_value(ψt, Ht(1.0), envst) + @test Ea ≈ Et atol = 1.0e-8 + end + + # 7. convergence order: the integrator is documented as (globally) first order in dt; assert at + # least that. (At these bond dimensions the augmented spans are near-exact and the measured + # slope is ≈ 2, so only a lower bound is imposed.) + @testset "first-order accuracy" begin + Random.seed!(2) + Lc = 4 + Hc = force_planar(transverse_field_ising(ComplexF64, Trivial; L = Lc)) + ψf = normalize!(FiniteMPS(rand, ComplexF64, Lc, ℙ^2, ℙ^4)) + Hmat = convert(TensorMap, Hc) + ψvec = convert(TensorMap, ψf); ψvec /= norm(ψvec) + + Tfin = 0.2 + dts = [0.05, 0.025, 0.0125] + errs = map(dts) do δt + n = round(Int, Tfin / δt) + ref = exp(-im * Hmat * (n * δt)) * ψvec + ψc = copy(ψf) + for k in 0:(n - 1) + ψc, = timestep(ψc, Hc, k * δt, δt, ParallelBUG(; trscheme = truncerror(; atol = 1.0e-12))) + end + out = convert(TensorMap, ψc); out /= norm(out) + return 1 - abs(dot(out, ref)) + end + slopes = [ + (log(errs[i + 1]) - log(errs[i])) / (log(dts[i + 1]) - log(dts[i])) + for i in 1:(length(dts) - 1) + ] + @info "ParallelBUG convergence" errs slopes + # documented rate: first order (slope ≈ 1); assert at least that + @test all(s -> s > 0.8, slopes) + end + + # 7b. THE decisive local-order gate (Ceruti et al. 2024; MPSKit design §8): the single-step local + # error ‖step(dt) − exp(−iH·dt)ψ₀‖ against the dense propagator, log–log slope vs dt, at + # L = 3 and L = 4 (2-site is exact and does NOT exercise the dropped "new–new" corner, so + # L ≥ 3 is required). A first-order method has local error O(dt²) ⇒ slope ≈ 2; this is the gate + # that separates genuine second order (slope ≈ 3, see `ParallelBUG2`) from a better-constant + # first-order scheme. Full bonds (ℙ^8) isolate the time-discretization error. + @testset "local order slope L=$Lc" for Lc in (3, 4) + Random.seed!(2) + Hc = force_planar(transverse_field_ising(ComplexF64, Trivial; L = Lc)) + ψf = normalize!(FiniteMPS(rand, ComplexF64, Lc, ℙ^2, ℙ^8)) + Hmat = convert(TensorMap, Hc) + ψvec = convert(TensorMap, ψf); ψvec /= norm(ψvec) + dts = [0.04, 0.02, 0.01, 0.005] + errs = map(dts) do δt + ref = exp(-im * Hmat * δt) * ψvec + ϕ, = timestep(ψf, Hc, 0.0, δt, ParallelBUG(; trscheme = truncerror(; atol = 1.0e-12))) + out = convert(TensorMap, ϕ) + return norm(out - ref) + end + slopes = [ + (log(errs[i + 1]) - log(errs[i])) / (log(dts[i + 1]) - log(dts[i])) + for i in 1:(length(dts) - 1) + ] + @info "ParallelBUG local order" Lc errs slopes + for s in slopes + @test s ≈ 2 atol = 0.3 + end + end + + # 8. the truncation tolerance ϑ maps onto the accumulating `c·n·ϑ` error term: tightening ϑ + # improves the overlap with a ϑ → 0 run of the same integrator (which shares the time + # discretization error, so the comparison isolates the truncation term). + @testset "accuracy improves as ϑ decreases" begin + Random.seed!(202) + Lc = 6 + Hc = force_planar(transverse_field_ising(ComplexF64, Trivial; L = Lc)) + ψ0 = normalize!(complex(FiniteMPS(rand, Float64, Lc, ℙ^2, ℙ^2))) + Tfin = 0.2; δt = 0.05; n = round(Int, Tfin / δt) + evolve(ϑ) = foldl(0:(n - 1); init = ψ0) do ψc, k + first(timestep(ψc, Hc, k * δt, δt, ParallelBUG(; trscheme = truncerror(; atol = ϑ)))) + end + ref = evolve(1.0e-12) + ϑs = [1.0e-2, 1.0e-4, 1.0e-6] + errs = map(ϑs) do ϑ + ψc = evolve(ϑ) + return 1 - abs(dot(ψc, ref)) / (norm(ψc) * norm(ref)) + end + @info "ParallelBUG accuracy vs ϑ" ϑs errs + @test issorted(errs; rev = true) + @test errs[end] < errs[1] / 10 + end + + # 9. step rejection: a small-bond start under a tight tolerance + a large `dt` can saturate the + # doubling cap; `maxiter_rejection > 0` then recomputes the step as half-steps. The recompute + # path must run cleanly, keep a normalized state, and never worsen the overlap with the dense + # reference (sub-stepping only refines it). + @testset "step rejection" begin + Random.seed!(7) + Lr = 4 + Hr = force_planar(transverse_field_ising(ComplexF64, Trivial; L = Lr)) + Hmat = convert(TensorMap, Hr) + ψr = normalize!(complex(FiniteMPS(rand, Float64, Lr, ℙ^2, ℙ^2))) # bond dim 2 (small) + ψvec = convert(TensorMap, ψr); ψvec /= norm(ψvec) + dtbig = 0.3 + ref = exp(-im * Hmat * dtbig) * ψvec + + tol = truncerror(; atol = 1.0e-6) + ψno, = timestep(ψr, Hr, 0.0, dtbig, ParallelBUG(; trscheme = tol, maxiter_rejection = 0)) + ψrej, = timestep(ψr, Hr, 0.0, dtbig, ParallelBUG(; trscheme = tol, maxiter_rejection = 4)) + outno = convert(TensorMap, ψno); outno /= norm(outno) + outrej = convert(TensorMap, ψrej); outrej /= norm(outrej) + errno = 1 - abs(dot(outno, ref)) + errrej = 1 - abs(dot(outrej, ref)) + @info "ParallelBUG step rejection" errno errrej + @test errrej ≤ errno + 1.0e-12 # sub-stepping never worsens the accuracy + @test isfinite(errrej) + end + + # 10. first-order error accumulation (Ceruti et al. 2024, Thm 4.5: `‖error‖ ≲ c·n·ϑ`). At fixed + # `dt` and `ϑ` the global error grows ~linearly in the number of steps `n`. (A `dt`-refinement + # slope is ill-posed here: at fixed `ϑ`, halving `dt` doubles `n` and thus the accumulated + # truncation error, so refining `dt` need not converge — hence this per-`n` test instead.) + @testset "linear error accumulation" begin + Random.seed!(13) + Le = 6 + He = force_planar(transverse_field_ising(ComplexF64, Trivial; L = Le)) + Hmat = convert(TensorMap, He) + ψe = normalize!(complex(FiniteMPS(rand, Float64, Le, ℙ^2, ℙ^2))) + ψvec = convert(TensorMap, ψe); ψvec /= norm(ψvec) + δt = 0.02 + alg = ParallelBUG(; trscheme = truncerror(; atol = 1.0e-4)) + ns = [5, 10, 20] + # state 2-norm error (phase-aligned), the quantity the `c·n·ϑ` bound controls — the overlap + # *infidelity* `1-|⟨·⟩|` is its square, so it would grow ~n² and must not be used here. + errs = map(ns) do n + ref = exp(-im * Hmat * (n * δt)) * ψvec + ψc = copy(ψe) + for k in 0:(n - 1) + ψc, = timestep(ψc, He, k * δt, δt, alg) + end + out = convert(TensorMap, ψc); out /= norm(out) + return sqrt(2 * (1 - abs(dot(out, ref)))) + end + @info "ParallelBUG error accumulation" ns errs + @test issorted(errs) # error grows monotonically with the step count + @test 2.5 * errs[1] < errs[3] < 6 * errs[1] # ~linear in n (5→20 = 4×), not saturating/quadratic + end +end + +# Charge-sector coverage: a fixed-rank symmetric step must preserve the total charge and the graded +# structure of every bond (energy conservation + eigenstate phase carry over from the trivial case). +@testset "ParallelBUG symmetric tensors" verbose = true begin + dt = 0.1 + L = 6 + + @testset "U(1) Heisenberg" begin + Random.seed!(2718) + H = heisenberg_XXX(ComplexF64, U1Irrep; spin = 1 // 2, L) + maxV = MPSKit.max_virtualspaces(physicalspace(H)) + ψ = FiniteMPS(physicalspace(H), maxV[2:(end - 1)]; right = U1Space(0 => 1)) + ψ₀, = find_groundstate(ψ, H; verbosity = 0) + E₀ = expectation_value(ψ₀, H) + Rtot = right_virtualspace(ψ₀, L) + Vl₀ = left_virtualspace.(Ref(ψ₀), 1:L) + + ψ1, envs = timestep(ψ₀, H, 0.0, dt, ParallelBUG()) # fixed-rank (notrunc) + E1 = expectation_value(ψ1, H, envs) + @test E₀ ≈ E1 atol = 1.0e-2 + @test imag(E1) ≈ 0 atol = 1.0e-8 + @test dot(ψ1, ψ₀) ≈ exp(im * dt * E₀) atol = 1.0e-4 + @test right_virtualspace(ψ1, L) == Rtot # total boundary charge preserved + @test left_virtualspace.(Ref(ψ1), 1:L) == Vl₀ # graded structure preserved (fixed rank) + end + + @testset "Z2 transverse-field Ising" begin + Random.seed!(161803) + H = transverse_field_ising(ComplexF64, Z2Irrep; g = 1.0, L) + ψ = FiniteMPS(physicalspace(H), Z2Space(0 => 4, 1 => 4)) + ψ₀, = find_groundstate(ψ, H; verbosity = 0) + E₀ = expectation_value(ψ₀, H) + + ψ1, envs = timestep(ψ₀, H, 0.0, dt, ParallelBUG()) + E1 = expectation_value(ψ1, H, envs) + @test E₀ ≈ E1 atol = 1.0e-2 + @test dot(ψ1, ψ₀) ≈ exp(im * dt * E₀) atol = 1.0e-4 + end +end + +# `ParallelBUG2` (second-order parallel BUG, Kusch 2024 Variant 2 / `4r`, EXPERIMENTAL). The genuine +# second-order assembly enriches every bond basis with one `H·ψ₀` application (`Û0`/`V̂0`), runs the +# `2r` Galerkin K-steps freezing the enriched right basis, and transports the EVOLVED-amplitude +# coupling to the root through the frozen right isometries, keeping the "new–new" corner zero (O(dt³)). +# It reproduces the exact two-site step, conserves energy / the eigenstate phase, agrees with `TDVP`, +# lowers the energy monotonically in imaginary time, preserves the charge sectors, and — the decisive +# gate — has a local error one order higher than the first-order `ParallelBUG`: slope ≥ 3 (vs ≈ 2). +@testset "ParallelBUG2 (second order)" verbose = true begin + dt = 0.1 + L = 6 + H = force_planar(heisenberg_XXX(Float64, Trivial; spin = 1 // 2, L)) + ψ = FiniteMPS(rand, Float64, L, ℙ^2, ℙ^4) + ψ₀, = find_groundstate(ψ, H; verbosity = 0) + E₀ = expectation_value(ψ₀, H) + + # 1. two sites: exact reproduction of the dense exp(-iH·dt) step + @testset "two-site exactness vs dense reference" begin + Random.seed!(2) + Hc = force_planar(transverse_field_ising(ComplexF64, Trivial; L = 2)) + ψc = FiniteMPS(rand, ComplexF64, 2, ℙ^2, ℙ^4) + Hmat = convert(TensorMap, Hc) + ψvec = convert(TensorMap, ψc); ψvec /= norm(ψvec) + ref = exp(-im * Hmat * 0.05) * ψvec + ψ1, = timestep(ψc, Hc, 0.0, 0.05, ParallelBUG2(; trscheme = truncerror(; atol = 1.0e-12))) + out = convert(TensorMap, ψ1); out /= norm(out) + @test 1 - abs(dot(out, ref)) < 1.0e-10 + end + + # 2. energy conservation + eigenstate phase (amplitude carried exactly once, at the root) + @testset "energy conservation + eigenstate phase" begin + ψ1, envs = timestep(ψ₀, H, 0.0, dt, ParallelBUG2(; trscheme = truncerror(; atol = 1.0e-12))) + E1 = expectation_value(ψ1, H, envs) + @test E₀ ≈ E1 atol = 1.0e-2 + @test dot(ψ1, ψ₀) ≈ exp(im * dt * E₀) atol = 1.0e-4 + end + + # 3. THE decisive local-order gate (mirror of the first-order testset above). The single-step + # local error ‖step(dt) − exp(−iH·dt)ψ₀‖ has slope ≈ 2 for first order and slope ≈ 3 for genuine + # second order (local `O(dt³)`). Full bonds (ℙ^8) isolate the time-discretization error. + @testset "local order slope ≈ 3 (full, L=$Lc)" for Lc in (4, 5) + Random.seed!(2) + Hc = force_planar(transverse_field_ising(ComplexF64, Trivial; L = Lc)) + ψf = normalize!(FiniteMPS(rand, ComplexF64, Lc, ℙ^2, ℙ^8)) + Hmat = convert(TensorMap, Hc) + ψvec = convert(TensorMap, ψf); ψvec /= norm(ψvec) + dts = [0.04, 0.02, 0.01, 0.005] + errs = map(dts) do δt + ref = exp(-im * Hmat * δt) * ψvec + ϕ, = timestep(ψf, Hc, 0.0, δt, ParallelBUG2(; trscheme = truncerror(; atol = 1.0e-12))) + return norm(convert(TensorMap, ϕ) - ref) + end + slopes = [ + (log(errs[i + 1]) - log(errs[i])) / (log(dts[i + 1]) - log(dts[i])) + for i in 1:(length(dts) - 1) + ] + @info "ParallelBUG2 local order (full)" Lc errs slopes + for s in slopes + @test s ≈ 3 atol = 0.4 # genuine second order (vs the first-order value 2) + end + end + + # the second-order slope must persist as the chain grows past a single interior bond; L = 6 with + # a reduced bond dimension ℙ^4 has an unsaturated interior bond (a genuine dropped corner). + @testset "local order slope ≈ 3 (reduced bond, L=6)" begin + Random.seed!(2) + Lc = 6 + Hc = force_planar(transverse_field_ising(ComplexF64, Trivial; L = Lc)) + ψf = normalize!(FiniteMPS(rand, ComplexF64, Lc, ℙ^2, ℙ^4)) + Hmat = convert(TensorMap, Hc) + ψvec = convert(TensorMap, ψf); ψvec /= norm(ψvec) + dts = [0.04, 0.02, 0.01, 0.005] + errs = map(dts) do δt + ref = exp(-im * Hmat * δt) * ψvec + ϕ, = timestep(ψf, Hc, 0.0, δt, ParallelBUG2(; trscheme = truncerror(; atol = 1.0e-12))) + return norm(convert(TensorMap, ϕ) - ref) + end + slopes = [ + (log(errs[i + 1]) - log(errs[i])) / (log(dts[i + 1]) - log(dts[i])) + for i in 1:(length(dts) - 1) + ] + @info "ParallelBUG2 local order (reduced)" errs slopes + for s in slopes + @test s ≈ 3 atol = 0.4 + end + end + + # 4. agreement with TDVP over a few short real-time steps of a random MPS + @testset "agreement with TDVP" begin + Random.seed!(1234) + ψr = complex(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^4)) + δt = 0.01 + ψ_p, ψ_tdvp = ψr, ψr + for k in 0:4 + ψ_p, = timestep(ψ_p, H, k * δt, δt, ParallelBUG2(; trscheme = truncerror(; atol = 1.0e-12))) + ψ_tdvp, = timestep(ψ_tdvp, H, k * δt, δt, TDVP()) + end + @test expectation_value(ψ_p, H) ≈ expectation_value(ψ_tdvp, H) atol = 1.0e-3 + @test abs(dot(ψ_p, ψ_tdvp)) ≈ 1 atol = 1.0e-3 + end + + # 5. imaginary-time evolution lowers the energy monotonically and stays normalized + @testset "imaginary-time lowers energy" begin + Random.seed!(5) + ψi = complex(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^4)) + E_start = real(expectation_value(ψi, H)) + E_prev = E_start + for _ in 1:20 + ψi, = timestep( + ψi, H, 0.0, 0.1, ParallelBUG2(; trscheme = truncerror(; atol = 1.0e-10)); + imaginary_evolution = true, normalize = true + ) + E_now = real(expectation_value(ψi, H)) + @test E_now ≤ E_prev + 1.0e-6 # monotone (non-increasing) energy + E_prev = E_now + end + @test E_prev < E_start - 1.0 # substantial lowering toward the ground state + @test norm(ψi) ≈ 1 atol = 1.0e-6 # imaginary-time renormalizes each step + end + + # 6. charge-sector coverage: a fixed-rank (notrunc) Z2 step preserves the total boundary charge + @testset "Z2 symmetry preservation" begin + Random.seed!(161803) + Hz = transverse_field_ising(ComplexF64, Z2Irrep; g = 1.0, L) + ψz = FiniteMPS(physicalspace(Hz), Z2Space(0 => 4, 1 => 4)) + ψz0, = find_groundstate(ψz, Hz; verbosity = 0) + Ez = expectation_value(ψz0, Hz) + Rtot = right_virtualspace(ψz0, L) + + ψ1, envs = timestep(ψz0, Hz, 0.0, dt, ParallelBUG2()) # fixed-rank (notrunc) + E1 = expectation_value(ψ1, Hz, envs) + @test Ez ≈ E1 atol = 1.0e-2 + @test dot(ψ1, ψz0) ≈ exp(im * dt * Ez) atol = 1.0e-4 + @test right_virtualspace(ψ1, L) == Rtot # total boundary charge preserved + end +end diff --git a/test/algorithms/timestep.jl b/test/algorithms/timestep.jl index a846b2c5c..8f0fe730f 100644 --- a/test/algorithms/timestep.jl +++ b/test/algorithms/timestep.jl @@ -9,7 +9,7 @@ using Test, TestExtras using MPSKit using TensorKit using TensorKit: ℙ -using LinearAlgebra: norm +using LinearAlgebra: norm, Diagonal using Random verbosity_full = 5 @@ -122,8 +122,9 @@ end @test abs(dot(ref, cbe)) > abs(dot(ref, plain)) end - # the bond truncation must preserve the norm for real-time evolution (the norm reflects the - # discarded weight) and only renormalize for imaginary-time evolution + # by default (`normalize = false`) the bond truncation preserves the norm, so it reflects the + # discarded weight; `normalize = true` renormalizes each step. This is independent of + # `imaginary_evolution`. @testset "norm handling" begin Random.seed!(6) ψ₀ = complex(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^Dstart)) @@ -132,15 +133,24 @@ end ψrt = ψ₀ for _ in 1:12 - ψrt, = timestep(ψrt, H, 0.0, 0.5, lossy) # real time + ψrt, = timestep(ψrt, H, 0.0, 0.5, lossy) # real time, norm preserved by default end @test norm(ψrt) < 1 - 1.0e-3 # truncation loss is not renormalized away + # imaginary-time, norm preserved by default: the weight is *not* pinned to unit norm + # (imaginary-time evolution rescales the state, so the norm drifts away from 1) ψit = ψ₀ for _ in 1:12 - ψit, = timestep(ψit, H, 0.0, 0.5, lossy; imaginary_evolution = true) # no external normalize! + ψit, = timestep(ψit, H, 0.0, 0.5, lossy; imaginary_evolution = true) end - @test norm(ψit) ≈ 1 atol = 1.0e-6 # imaginary-time renormalizes each step + @test abs(norm(ψit) - 1) > 1.0e-3 + + # imaginary-time with `normalize = true`: renormalized to unit norm each step + ψn = ψ₀ + for _ in 1:12 + ψn, = timestep(ψn, H, 0.0, 0.5, lossy; imaginary_evolution = true, normalize = true) + end + @test norm(ψn) ≈ 1 atol = 1.0e-6 end @testset "imaginary-time lowers energy" begin @@ -150,7 +160,7 @@ end E₀ = real(expectation_value(ψ₀, H)) ψ = ψ₀ for _ in 1:8 - ψ, = timestep(ψ, H, 0.0, 0.1, alg; imaginary_evolution = true) # gauge renormalizes + ψ, = timestep(ψ, H, 0.0, 0.1, alg; imaginary_evolution = true, normalize = true) end @test real(expectation_value(ψ, H)) < E₀ @test dim(left_virtualspace(ψ, L ÷ 2)) > Dstart @@ -182,3 +192,105 @@ end @test E₀ ≈ E atol = 1.0e-2 end end + +# A single spin flip on the fully polarized vacuum of the XX chain (`XY_model` with `g = 0`) +# propagates as a free particle hopping on the open L-site chain, so the real-time correlator +# has a closed form: +# C_nm(t) = ⟨vac| S^x_n e^{-iHt} S^x_m |vac⟩ = ¼ [exp(-i h t)]_{nm}, +# where the single-particle hopping matrix `h` is tridiagonal (`h_{i,i±1} = -1`) and diagonalizes +# analytically on the open chain: eigenvalues `ε_k = -2 cos(kπ/(L+1))` with eigenvectors +# `φ_k(n) = √(2/(L+1)) sin(nkπ/(L+1))`. Being a genuine two-time correlator, this also exercises +# the norm-preserving `normalize = false` default: the evolving state must NOT be renormalized. +@testset "dynamical correlator vs free-particle analytics" begin + L = 16 + m = 8 # source site (chain centre) + ns = 5:11 # measured sink sites + dt = 0.05 + ts = 0.0:dt:1.5 + + H = force_planar(XY_model(ComplexF64, Trivial; g = 0, L = L)) + + # fully polarized product state |↑…↑⟩ (basis index 1 on every site): a zero-energy eigenstate + Vp = ℙ^2 + Vl = oneunit(Vp) + ψvac = FiniteMPS([isometry(ComplexF64, Vl ⊗ Vp, Vl) for _ in 1:L]) + @test norm(H * ψvac) < 1.0e-10 # H annihilates the polarized vacuum ⇒ E₀ = 0, no bra phase + + # apply S^x at a single site (S^x = ½σ^x flips one spin ⇒ a single free particle) + Sx = force_planar(S_x()) + apply_Sx(i) = FiniteMPOHamiltonian(physicalspace(ψvac), i => Sx) * ψvac + ϕ_bra = Dict(n => apply_Sx(n) for n in ns) # static ½|n⟩ + ϕ = apply_Sx(m) # ½|m⟩, evolved below + + # analytic single-particle propagator exp(-i h t) on the open chain + ks = 1:L + ε = [-2 * cos(k * π / (L + 1)) for k in ks] + φ = [sqrt(2 / (L + 1)) * sin(n * k * π / (L + 1)) for n in 1:L, k in ks] + Uref(t) = φ * Diagonal(cis.(-ε .* t)) * transpose(φ) + Cref(n, t) = 0.25 * Uref(t)[n, m] + + # single-particle physics stays rank ≤ 2, so TDVP2 with a generous cap is essentially exact + alg = TDVP2(; trscheme = truncrank(6)) + t_prev = 0.0 + maxerr = 0.0 + for t in ts + t > 0 && ((ϕ, _) = timestep(ϕ, H, t_prev, t - t_prev, alg); t_prev = t) + for n in ns + maxerr = max(maxerr, abs(dot(ϕ_bra[n], ϕ) - Cref(n, t))) + end + end + @test maxerr < 1.0e-3 + + # the correlator retains the un-normalized amplitude: real-time evolution preserves the norm + @test norm(ϕ) ≈ norm(ϕ_bra[m]) atol = 1.0e-6 +end + +# The same free-particle benchmark with an explicitly fermionic (fℤ₂-graded) model: a single +# fermion on the tight-binding chain (`kitaev_model` with `mu = Delta = 0`) is a free particle, so +# the single-particle Green's function is the closed-form propagator +# G_nm(t) = ⟨0| c_n e^{-iHt} c_m† |0⟩ = [exp(-i h t)]_{nm}, h_{i,i±1} = -t/2, +# with `ε_k = -t cos(kπ/(L+1))` and `φ_k(n) = √(2/(L+1)) sin(nkπ/(L+1))`. This exercises the graded +# fermionic tensors (Jordan-Wigner strings handled by the fℤ₂ braiding): from the empty vacuum, +# `c_m†|0⟩` is just the localised single-particle basis state (the string acts trivially on the +# vacuum), built here as a bond-1 MPS with a fermion-parity charge kink at site m. It relies on the +# hopping being Hermitian (`f_hopping`, not `f_plus_f_min + f_min_f_plus`), else the norm drifts. +@testset "fermionic dynamical correlator vs free-particle analytics" begin + L = 12 + m = 6 + ns = 4:8 + dt = 0.05 + ts = 0.0:dt:1.2 + + H = kitaev_model(ComplexF64, Trivial; t = 1.0, mu = 0, Delta = 0, L = L) + P = physicalspace(H)[1] # fermion_space = Vect[fℤ₂](0 => 1, 1 => 1) + S = typeof(P) + Vtriv = oneunit(P) + Vodd = S(c => dim(P, c) for c in sectors(P) if !isone(c)) + + # single-particle basis state |j⟩ = c_j†|0⟩: empty except occupied at j, with the fermion-parity + # kink even→odd at site j (even bond to the left, odd bond to the right, odd right boundary) + onepart(j) = FiniteMPS([ + isometry(ComplexF64, (k <= j ? Vtriv : Vodd) ⊗ P, (k < j ? Vtriv : Vodd)) for k in 1:L + ]) + @test norm(H * onepart(m)) ≈ sqrt(2) * 0.5 atol = 1.0e-8 # hops to both neighbours, amp t/2 + + ϕ_bra = Dict(n => onepart(n) for n in ns) + ϕ = onepart(m) + + ks = 1:L + ε = [-cos(k * π / (L + 1)) for k in ks] # -t cos(kπ/(L+1)) with t = 1 + φ = [sqrt(2 / (L + 1)) * sin(n * k * π / (L + 1)) for n in 1:L, k in ks] + Uref(t) = φ * Diagonal(cis.(-ε .* t)) * transpose(φ) + + alg = TDVP2(; trscheme = truncrank(6)) + t_prev = 0.0 + maxerr = 0.0 + for t in ts + t > 0 && ((ϕ, _) = timestep(ϕ, H, t_prev, t - t_prev, alg); t_prev = t) + for n in ns + maxerr = max(maxerr, abs(dot(ϕ_bra[n], ϕ) - Uref(t)[n, m])) + end + end + @test maxerr < 1.0e-3 + @test norm(ϕ) ≈ 1 atol = 1.0e-8 # real-time evolution preserves the norm +end diff --git a/test/setup/testsetup.jl b/test/setup/testsetup.jl index 053d272f3..8c2bfbd3f 100644 --- a/test/setup/testsetup.jl +++ b/test/setup/testsetup.jl @@ -11,7 +11,7 @@ using BlockTensorKit using LinearAlgebra: Diagonal using Combinatorics: permutations using TensorKitTensors.SpinOperators: S_x, S_y, S_z, S_x_S_x, S_y_S_y, S_z_S_z, S_exchange, S_plus_S_min, S_min_S_plus -using TensorKitTensors.FermionOperators: f_plus_f_min, f_min_f_plus, f_plus_f_plus, f_min_f_min, f_num +using TensorKitTensors.FermionOperators: f_plus_f_min, f_min_f_plus, f_plus_f_plus, f_min_f_min, f_num, f_hopping # exports export S_x, S_y, S_z @@ -154,8 +154,8 @@ function kitaev_model( T::Type{<:Number} = ComplexF64, sym::Type{<:Sector} = Trivial; t = 1.0, mu = 1.0, Delta = 1.0, L = Inf ) - TB = scale!(f_plus_f_min(T, sym) + f_min_f_plus(T, sym), -t / 2) # tight-binding term - SC = scale!(f_plus_f_plus(T, sym) + f_min_f_min(T, sym), Delta / 2) # superconducting term + TB = scale!(f_hopping(T, sym), -t / 2) # tight-binding term + SC = scale!(f_min_f_min(T, sym) - f_plus_f_plus(T, sym), Delta / 2) # superconducting term CP = scale!(f_num(T, sym), -mu) # chemical potential term if L == Inf @@ -165,7 +165,7 @@ function kitaev_model( lattice = fill(space(TB, 1), L) onsite_terms = ((i,) => CP for i in 1:L) twosite_terms = ((i, i + 1) => TB + SC for i in 1:(L - 1)) - terms = Iterators.flatten(twosite_terms, onsite_terms) + terms = Iterators.flatten((twosite_terms, onsite_terms)) return FiniteMPOHamiltonian(lattice, terms) end end