diff --git a/docs/api/bond-update-bug/index.md b/docs/api/bond-update-bug/index.md
new file mode 100644
index 0000000..2879b15
--- /dev/null
+++ b/docs/api/bond-update-bug/index.md
@@ -0,0 +1,40 @@
+# bond_update_bug
+
+Alice's bond_update_bug (Basis-Update & Galerkin) integrator evolves an MPS in real or imaginary time under a nearest-neighbour Hamiltonian. It is the rank-adaptive BUG of Ceruti, Kusch & Lubich ([arXiv:2304.05660](https://arxiv.org/abs/2304.05660)): commuting even/odd Trotter sweeps of *local* K/L/S bond updates. Each update augments the left frame from the evolved **K** factor, augments the right frame from the evolved **L** factor, evolves the small core **S** in the augmented bases (Galerkin), then truncates with an SVD — so the bond dimension adapts to the growing entanglement (the basis augmentation). The local substeps exponentiate the *projected* effective Hamiltonian internally (Krylov `expv`); no pre-formed propagator gate is applied, and the update is exact at full rank.
+
+The bond Hamiltonians are reused directly from the [AutoMPO](../interaction/build-interaction.md) interaction list, so any nearest-neighbour model and symmetry that `build_interaction` supports works unchanged.
+
+## API
+
+| Symbol | Description |
+|--------|-------------|
+| [Options](options.md) | Run options: time step, steps, Trotter order, bond dimension |
+| [Summary](summary.md) | Output: evolved MPS, time/norm history, kept and augmented bond dims |
+| [run](run.md) | Top-level entry point |
+
+## Usage Pattern
+
+```python
+from alice import build_interaction, init_mps
+from alice.algorithm import bond_update_bug
+
+interactions, spc, geo = build_interaction("config.toml")
+mps = init_mps(geo.L, spc, Op, config=[0, 1] * (geo.L // 2), target_qn=0)
+opts = bond_update_bug.Options(dt=0.05, n_steps=40, order='strang', max_bond=128)
+
+summary = bond_update_bug.run(mps, interactions, opts)
+print(summary.max_bond_dims) # kept bond dimension per step
+print(summary.aug_dims) # proposed (pre-truncation) augmentation per step
+```
+
+## Trotter Orders
+
+| Name | Alias | Description |
+|------|-------|-------------|
+| `'strang'` | `'second'`, `'2'` | Symmetric second-order step `U_even(dt/2) U_odd(dt) U_even(dt/2)` |
+| `'lie'` | `'first'`, `'1'` | First-order step `U_even(dt) U_odd(dt)` |
+
+## See Also
+
+- [bond_update_bug.run](run.md) — full parameter reference.
+- [build_interaction](../interaction/build-interaction.md) — build the `interactions` argument.
diff --git a/docs/api/bond-update-bug/options.md b/docs/api/bond-update-bug/options.md
new file mode 100644
index 0000000..95bcd64
--- /dev/null
+++ b/docs/api/bond-update-bug/options.md
@@ -0,0 +1,38 @@
+# Options
+
+bond_update_bug run options.
+
+::: alice.algorithm.bond_update_bug.Options
+ options:
+ heading_level: 2
+
+## TOML Loading
+
+`Options` can be loaded directly from an `[algorithm]` TOML section:
+
+```python
+import tomllib
+from alice.algorithm import bond_update_bug
+
+with open("config.toml", "rb") as f:
+ cfg = tomllib.load(f)
+
+opts = bond_update_bug.Options.from_toml(cfg["heisenberg"]["algorithm"])
+```
+
+Example TOML block:
+
+```toml
+[heisenberg.algorithm]
+dt = 0.05
+n_steps = 40
+order = "strang"
+max_bond = 128
+trunc_thresh = 1e-12
+imaginary_time = false
+```
+
+## See Also
+
+- [Summary](summary.md) — output dataclass.
+- [run](run.md) — pass `Options` here.
diff --git a/docs/api/bond-update-bug/run.md b/docs/api/bond-update-bug/run.md
new file mode 100644
index 0000000..8fc69da
--- /dev/null
+++ b/docs/api/bond-update-bug/run.md
@@ -0,0 +1,13 @@
+# Launch
+
+Evolve an MPS under a nearest-neighbour Hamiltonian with the bond_update_bug integrator.
+
+::: alice.algorithm.bond_update_bug.run
+ options:
+ heading_level: 2
+
+## See Also
+
+- [Options](options.md) — configure the run.
+- [Summary](summary.md) — interpret the output.
+- [build_interaction](../interaction/build-interaction.md) — create the `interactions` argument.
diff --git a/docs/api/bond-update-bug/summary.md b/docs/api/bond-update-bug/summary.md
new file mode 100644
index 0000000..505b134
--- /dev/null
+++ b/docs/api/bond-update-bug/summary.md
@@ -0,0 +1,12 @@
+# Summary
+
+bond_update_bug output.
+
+::: alice.algorithm.bond_update_bug.Summary
+ options:
+ heading_level: 2
+
+## See Also
+
+- [Options](options.md) — configure the run.
+- [run](run.md) — produces this dataclass.
diff --git a/docs/api/index.md b/docs/api/index.md
index 07587d0..f20d3c2 100644
--- a/docs/api/index.md
+++ b/docs/api/index.md
@@ -67,6 +67,27 @@ Ground-state DMRG algorithm.
| [Summary](dmrg/summary.md) | DMRG output dataclass |
| [run](dmrg/run.md) | Top-level DMRG entry point |
+## bond_update_bug
+
+Rank-adaptive Basis-Update & Galerkin time integrator (the discarded-projector
+K/L/S sweep; real and imaginary time).
+
+| Symbol | Description |
+|--------|-------------|
+| [Options](bond-update-bug/options.md) | bond_update_bug run options |
+| [Summary](bond-update-bug/summary.md) | bond_update_bug output dataclass |
+| [run](bond-update-bug/run.md) | Top-level bond_update_bug entry point |
+
+## Two-Site TDVP
+
+Rank-adaptive two-site TDVP time integrator (real and imaginary time).
+
+| Symbol | Description |
+|--------|-------------|
+| [Options](tdvp2/options.md) | TDVP run options |
+| [Summary](tdvp2/summary.md) | TDVP output dataclass |
+| [run](tdvp2/run.md) | Top-level TDVP entry point |
+
## Logging
| Symbol | Description |
diff --git a/docs/api/tdvp2/index.md b/docs/api/tdvp2/index.md
new file mode 100644
index 0000000..488ed24
--- /dev/null
+++ b/docs/api/tdvp2/index.md
@@ -0,0 +1,48 @@
+# Two-Site TDVP
+
+Alice's two-site TDVP (Time-Dependent Variational Principle) integrator evolves an MPS in real or imaginary time under a Hamiltonian MPO. It is the projector-splitting scheme of Haegeman et al. ([arXiv:1408.5056](https://arxiv.org/abs/1408.5056)) with a two-site update so the bond dimension adapts. A forward half-sweep evolves each two-site block forward in time and the carried one-site tensor backward in time (the inverse-free backward correction that removes the double counting of the shared bond); a reverse half-sweep mirrors it; a symmetric step composes the two halves for second-order accuracy. The local substeps exponentiate the *effective Hamiltonian* — the MPS tensor bracketed by the left/right MPO environments — with a Hermitian Krylov `expv`.
+
+TDVP needs the full effective Hamiltonian, so it takes a Hamiltonian MPO built by [`build_hamiltonian`](../hamiltonian/build-hamiltonian.md) — exactly like [DMRG](../dmrg/index.md). It reuses the DMRG environment machinery and effective-Hamiltonian contractions.
+
+## API
+
+| Symbol | Description |
+|--------|-------------|
+| [Options](options.md) | Run options: time step, steps, bond dimension, real/imaginary time |
+| [Summary](summary.md) | Output: evolved MPS, time/norm history, kept bond dims |
+| [run](run.md) | Top-level entry point |
+
+## Usage Pattern
+
+```python
+from alice import build_interaction, build_hamiltonian, init_mps
+from alice.algorithm import tdvp2
+
+interactions, spc, geo = build_interaction("config.toml")
+mpo = build_hamiltonian(interactions, geo.L, spc)
+mps = init_mps(geo.L, spc, Op, config=[0, 1] * (geo.L // 2), target_qn=0)
+opts = tdvp2.Options(dt=0.05, n_steps=40, max_bond=128)
+
+summary = tdvp2.run(mps, mpo, opts)
+print(summary.max_bond_dims) # kept bond dimension per step
+print(summary.norms) # norm per step (≈1 for real time; decays for imaginary)
+```
+
+## Real vs. Imaginary Time
+
+| `imaginary_time` | Propagator | Use |
+|------------------|------------|-----|
+| `False` (default) | `exp(-i dt H)` | unitary real-time dynamics; the norm is conserved |
+| `True` | `exp(-dt H)` | imaginary-time cooling toward the ground state (pair with `normalize=True`) |
+
+!!! note "Convergence at fixed bond dimension"
+ At fixed or adaptively-capped bond dimension, two-site TDVP's error is a
+ *manifold-projection* error that does not vanish as `dt → 0` — it plateaus —
+ rather than the `O(dt²)` state error of a full-rank propagator. Refine the bond
+ dimension (`max_bond`) to reduce the plateau.
+
+## See Also
+
+- [tdvp2.run](run.md) — full parameter reference.
+- [build_hamiltonian](../hamiltonian/build-hamiltonian.md) — build the `mpo` argument.
+- [DMRG](../dmrg/index.md) — ground-state search sharing the same MPO/environment core.
diff --git a/docs/api/tdvp2/options.md b/docs/api/tdvp2/options.md
new file mode 100644
index 0000000..3677825
--- /dev/null
+++ b/docs/api/tdvp2/options.md
@@ -0,0 +1,37 @@
+# Options
+
+Two-site TDVP run options.
+
+::: alice.algorithm.tdvp2.Options
+ options:
+ heading_level: 2
+
+## TOML Loading
+
+`Options` can be loaded directly from an `[algorithm]` TOML section:
+
+```python
+import tomllib
+from alice.algorithm import tdvp2
+
+with open("config.toml", "rb") as f:
+ cfg = tomllib.load(f)
+
+opts = tdvp2.Options.from_toml(cfg["heisenberg"]["algorithm"])
+```
+
+Example TOML block:
+
+```toml
+[heisenberg.algorithm]
+dt = 0.05
+n_steps = 40
+max_bond = 128
+cutoff = 1e-12
+imaginary_time = false
+```
+
+## See Also
+
+- [Summary](summary.md) — output dataclass.
+- [run](run.md) — pass `Options` here.
diff --git a/docs/api/tdvp2/run.md b/docs/api/tdvp2/run.md
new file mode 100644
index 0000000..59a039b
--- /dev/null
+++ b/docs/api/tdvp2/run.md
@@ -0,0 +1,13 @@
+# Launch
+
+Evolve an MPS under a Hamiltonian MPO with the two-site TDVP integrator.
+
+::: alice.algorithm.tdvp2.run
+ options:
+ heading_level: 2
+
+## See Also
+
+- [Options](options.md) — configure the run.
+- [Summary](summary.md) — interpret the output.
+- [build_hamiltonian](../hamiltonian/build-hamiltonian.md) — create the `mpo` argument.
diff --git a/docs/api/tdvp2/summary.md b/docs/api/tdvp2/summary.md
new file mode 100644
index 0000000..d6b0e8c
--- /dev/null
+++ b/docs/api/tdvp2/summary.md
@@ -0,0 +1,12 @@
+# Summary
+
+Two-site TDVP output.
+
+::: alice.algorithm.tdvp2.Summary
+ options:
+ heading_level: 2
+
+## See Also
+
+- [Options](options.md) — configure the run.
+- [run](run.md) — produces this dataclass.
diff --git a/docs/getting-started/changelog.md b/docs/getting-started/changelog.md
index aa301d2..e4971bf 100644
--- a/docs/getting-started/changelog.md
+++ b/docs/getting-started/changelog.md
@@ -1,5 +1,105 @@
# Changelog
+## [Unreleased]
+
+**Two-Site BUG Time Integrator**
+
+Adds `alice.algorithm.bond_update_bug`, the rank-adaptive bond_update_bug
+(Basis-Update & Galerkin) integrator of Ceruti, Kusch & Lubich
+([arXiv:2304.05660](https://arxiv.org/abs/2304.05660)) for real- and
+imaginary-time evolution of an MPS under a nearest-neighbour Hamiltonian. The
+Alice-facing driver is built on the existing Alice/Nicole stack — `MPS`, the
+AutoMPO interaction list, and the PyTorch backend; the symmetry-aware KLS
+local kernel is vendored, Nicole-native, in a private `_kernel` subpackage.
+
+### `alice.algorithm.bond_update_bug`
+
+- **`run(mps, interactions, opts)`** evolves the state with commuting even/odd
+ Trotter sweeps of *local* K/L/S bond updates: each update augments the left and
+ right frames from the evolved K and L factors, evolves the small core in the
+ augmented bases (Galerkin), and truncates with an SVD so the bond dimension
+ adapts (the basis augmentation). The local substeps exponentiate the projected
+ effective Hamiltonian internally (Krylov `expv`) — exact at full rank. Supports
+ first-order (`'lie'`) and symmetric second-order (`'strang'`) steps and
+ imaginary-time cooling.
+- **Bond Hamiltonians** are reused from the AutoMPO interaction list: the leading
+ and terminal MPO tensors of each nearest-neighbour `Interaction2Site` are
+ contracted over their operator channel to form the bare two-site term fed to the
+ KLS kernel. The kernel is symmetry-aware (works with the U(1) charge sectors of
+ the MPS).
+- **`Options`** (TOML-loadable) and **`Summary`** mirror the DMRG interface. The
+ summary records, per step, the kept bond dimension and the *proposed* augmented
+ dimension, so the rank growth and the discarded augmentation are both visible.
+- Validated against exact diagonalization (state fidelity, exact norm
+ conservation, U(1) charge conservation, and second-order Trotter convergence).
+
+**Discarded-Projector BUG Variant**
+
+Adds `alice.algorithm.discarded_bug`, the MPS specialisation of the Lubich
+tree-tensor-network BUG (Ceruti–Lubich–Walach,
+[arXiv:2304.05660](https://arxiv.org/abs/2304.05660)). Like two-site TDVP and DMRG
+it takes a Hamiltonian **MPO** and exponentiates the two-site effective Hamiltonian
+with the left/right MPO environments, reusing the DMRG environment machinery; it is
+inverse-free (no backward substep, no overlap-matrix inverse).
+
+### `alice.algorithm.discarded_bug`
+
+- **`run(mps, mpo, opts)`** evolves the state by **recursive bisection** of the
+ chain — the MPS realisation of the reference's balanced-binary-tree `Step` (whose
+ tree is built by recursive bisection of the 1D modes). Each step updates the
+ central bisection bond, then recurses into the two half-chains, until every bond —
+ every tree node — has had its two-site node update. Because every bond is a node,
+ the bond dimension grows along the whole chain (the full ballistic light cone) as
+ a domain wall melts, matching the bond growth of forward two-site TDVP.
+- **Node update.** At each bisection bond the two-site block is evolved once,
+ `Θ1 = exp(τ H₂) Θ0` (Hermitian → tensor Lanczos); the K-step and L-step grow the
+ left/right frames with the **discarded** projector — `qr([Θ1_left | U0])` /
+ `qr([Θ1_right ; V0])`, the direct sum of the old frame with the evolved block's
+ column/row space — with **no** augmented overlap matrices; the Galerkin core is
+ the projection `Û† Θ1 V̂†` of the already-evolved block, SVD-truncated to set the
+ rank. The frames are read off the *evolved* block so a product-state interface
+ grows its genuine rank-2 entanglement (a frozen-neighbour generator would project
+ it out). Everything stays in the U(1) block-sparse Nicole representation, so the
+ kept bond dimension respects the charge sectors.
+- The step is first order in `dt` (no backward substep); the validated property is
+ the rank growth / light-cone spread. A second-order symmetric composition is left
+ to future work.
+- **`Options`** (TOML-loadable) and **`Summary`** mirror the DMRG interface; the
+ summary records the kept bond dimension per step and the final bond dimensions
+ (the light cone).
+- Validated against exact diagonalization (full light-cone growth tracking forward
+ two-site TDVP, first-order single-step convergence, exact norm and U(1) charge
+ conservation, imaginary-time energy descent).
+
+**Two-Site TDVP Time Integrator**
+
+Adds `alice.algorithm.tdvp2`, a rank-adaptive two-site Time-Dependent Variational
+Principle integrator (Haegeman et al., [arXiv:1408.5056](https://arxiv.org/abs/1408.5056))
+for real- and imaginary-time evolution of an MPS under a Hamiltonian MPO. It is
+built on the shared MPS/MPO core and reuses the DMRG environment machinery and
+effective-Hamiltonian contractions, so it depends only on `alice.network` and
+`alice.algorithm.dmrg`.
+
+### `alice.algorithm.tdvp2`
+
+- **`run(mps, mpo, opts)`** evolves the state with symmetric (Strang) steps: a
+ forward half-sweep evolves each two-site block forward in time and the carried
+ one-site tensor backward in time (the inverse-free backward correction that
+ removes the shared-bond double counting), and a reverse half-sweep mirrors it.
+ The local substeps exponentiate the *effective Hamiltonian* — the MPS tensor
+ bracketed by the left/right MPO environments — with a Hermitian Krylov `expv`.
+ The per-bond SVD truncation makes the bond dimension adapt; real and imaginary
+ time (ground-state cooling) are both supported.
+- Takes a Hamiltonian **MPO** (from `build_hamiltonian`), like DMRG, and reuses
+ the DMRG `Environment` blocks, transfer-matrix steps, and the 1-/2-site
+ effective-Hamiltonian contractions. The local Krylov exponential and the
+ evolution-prefactor handling are self-contained in the package.
+- **`Options`** (TOML-loadable) and **`Summary`** mirror the DMRG interface; the
+ summary records the kept bond dimension and the norm per step.
+- Validated against exact diagonalization on the Heisenberg chain (state fidelity,
+ exact norm conservation, U(1) total-Sz conservation, imaginary-time cooling
+ toward the ground state, and bond-dimension growth as a domain wall melts).
+
## [0.1.6] - 2026-06-10
**MPS Initialization for Odd Chains**
diff --git a/mkdocs.yml b/mkdocs.yml
index a3fa473..f8ba31c 100644
--- a/mkdocs.yml
+++ b/mkdocs.yml
@@ -103,6 +103,16 @@ nav:
- Options: api/dmrg/options.md
- Summary: api/dmrg/summary.md
- Launch: api/dmrg/run.md
+ - bond_update_bug:
+ - Overview: api/bond-update-bug/index.md
+ - Options: api/bond-update-bug/options.md
+ - Summary: api/bond-update-bug/summary.md
+ - Launch: api/bond-update-bug/run.md
+ - Two-Site TDVP:
+ - Overview: api/tdvp2/index.md
+ - Options: api/tdvp2/options.md
+ - Summary: api/tdvp2/summary.md
+ - Launch: api/tdvp2/run.md
- Examples:
- Overview: examples/index.md
- DMRG:
diff --git a/pyproject.toml b/pyproject.toml
index 76a0b7e..816c5e2 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -80,6 +80,12 @@ markers = [
line-length = 100
target-version = "py311"
+[tool.ruff.lint.per-file-ignores]
+# Vendored, Nicole-native bond_update_bug kernel — kept close to its upstream form.
+# Several modules re-export helpers consumed by sibling kernel modules, so the
+# unused-import rule would force churn that breaks those re-exports.
+"src/alice/algorithm/bond_update_bug/_kernel/**" = ["F401"]
+
[tool.mypy]
python_version = "3.11"
check_untyped_defs = true
diff --git a/src/alice/__init__.py b/src/alice/__init__.py
index c7a7ed6..4abe249 100644
--- a/src/alice/__init__.py
+++ b/src/alice/__init__.py
@@ -28,7 +28,7 @@
init_mps,
observe,
)
-from .algorithm import dmrg
+from .algorithm import dmrg, tdvp2, bond_update_bug
from .logging import configure_logging
__version__ = version('alice-net')
@@ -51,6 +51,8 @@
'observe',
# algorithms (as submodules)
'dmrg',
+ 'tdvp2',
+ 'bond_update_bug',
# logging
'configure_logging',
]
diff --git a/src/alice/algorithm/__init__.py b/src/alice/algorithm/__init__.py
index a17a83b..7774205 100644
--- a/src/alice/algorithm/__init__.py
+++ b/src/alice/algorithm/__init__.py
@@ -16,10 +16,19 @@
# along with Alice. If not, see .
-"""Algorithm module: tensor network algorithms built on the network layer."""
+"""Algorithm module: tensor network algorithms built on the network layer.
+`bond_update_bug` is the single Basis-Update & Galerkin time integrator (the
+discarded-projector K/L/S sweep, mirrored by `bond_update_bug!` in BUG-Julia);
+`dmrg` and `tdvp2` are the ground-state and TDVP algorithms.
+"""
+
+from . import bond_update_bug
from . import dmrg
+from . import tdvp2
__all__ = [
+ 'bond_update_bug',
'dmrg',
+ 'tdvp2',
]
diff --git a/src/alice/algorithm/bond_update_bug/__init__.py b/src/alice/algorithm/bond_update_bug/__init__.py
new file mode 100644
index 0000000..32ffd16
--- /dev/null
+++ b/src/alice/algorithm/bond_update_bug/__init__.py
@@ -0,0 +1,43 @@
+# Copyright (C) 2025-2026 Changkai Zhang.
+#
+# This file is part of Alice project.
+#
+# Alice is free software: you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published
+# by the Free Software Foundation, either version 3 of the License,
+# or (at your option) any later version.
+#
+# Alice is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with Alice. If not, see .
+
+
+"""bond_update_bug algorithm package.
+
+Implements the bond_update_bug (Basis-Update & Galerkin) time integrator
+of Ceruti, Kusch & Lubich (arXiv:2304.05660): a nearest-neighbour Hamiltonian is
+evolved by odd/even Trotter sweeps of local K/L/S bond updates. Each update
+augments the left/right frames from the evolved K/L factors, evolves the small
+core in the augmented bases (Galerkin), and truncates with an SVD — exact at
+full rank, rank-adaptive otherwise. The public API includes:
+
+- `Options` — run options (loadable from TOML).
+- `Summary` — output dataclass.
+- `run` — top-level entry point.
+
+The KLS local kernel lives in the vendored, Nicole-native `_kernel`
+subpackage; this package wires it to Alice's `MPS` and AutoMPO bond terms.
+"""
+
+from .bond_update_bug import Options, Summary
+from .bond_update_bug import run
+
+__all__ = [
+ 'Options',
+ 'Summary',
+ 'run',
+]
diff --git a/src/alice/algorithm/bond_update_bug/_kernel/__init__.py b/src/alice/algorithm/bond_update_bug/_kernel/__init__.py
new file mode 100644
index 0000000..3e6c6e7
--- /dev/null
+++ b/src/alice/algorithm/bond_update_bug/_kernel/__init__.py
@@ -0,0 +1,57 @@
+# Copyright (C) 2025-2026 Changkai Zhang.
+#
+# This file is part of Alice project.
+#
+# Alice is free software: you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published
+# by the Free Software Foundation, either version 3 of the License,
+# or (at your option) any later version.
+#
+# Alice is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with Alice. If not, see .
+
+
+"""Vendored bond_update_bug local-bond kernel.
+
+This subpackage is the Nicole-native Basis-Update & Galerkin (BUG) local kernel —
+the discarded-projector K/L/S two-site update (Ceruti–Kusch–Lubich,
+arXiv:2304.05660), ported from the reference Julia `bond_update_bug!`. It is
+symmetry-aware (works with the U(1) charge sectors of an Alice `MPS`) and depends
+only on `nicole` + torch:
+
+- `_kls_local_bond_candidate` — one K/L/S local bond update.
+- `Ix` / `fresh_itag` — lightweight Nicole-index handles used by the kernel.
+- `qr` / `lq` — Nicole-backed decompositions returning `Ix` metadata.
+- `dag` / `tcontract` / `make_tensor` / `to_dense` — Nicole tensor helpers.
+- `with_time_prefactor` / `with_expv_backend` — evolution-prefactor and Krylov
+ backend context managers used to drive the local `expv` substeps.
+
+It is private to `alice.algorithm.bond_update_bug`; the Alice-facing driver in
+`bond_update_bug.py` builds the bond Hamiltonians from AutoMPO and runs the
+odd/even Strang sweep on an Alice `MPS` through this kernel.
+"""
+
+from .indices import Ix, fresh_itag
+from .krylov import with_expv_backend, with_time_prefactor
+from .kls import _kls_local_bond_candidate
+from .linalg import lq, qr
+from .nicole_helpers import dag, make_tensor, tcontract, to_dense
+
+__all__ = [
+ 'Ix',
+ 'fresh_itag',
+ 'with_expv_backend',
+ 'with_time_prefactor',
+ '_kls_local_bond_candidate',
+ 'lq',
+ 'qr',
+ 'dag',
+ 'make_tensor',
+ 'tcontract',
+ 'to_dense',
+]
diff --git a/src/alice/algorithm/bond_update_bug/_kernel/indices.py b/src/alice/algorithm/bond_update_bug/_kernel/indices.py
new file mode 100644
index 0000000..5579295
--- /dev/null
+++ b/src/alice/algorithm/bond_update_bug/_kernel/indices.py
@@ -0,0 +1,274 @@
+# Copyright (C) 2025-2026 Changkai Zhang.
+#
+# This file is part of Alice project.
+#
+# Alice is free software: you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published
+# by the Free Software Foundation, either version 3 of the License,
+# or (at your option) any later version.
+#
+# Alice is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with Alice. If not, see .
+
+"""Index builders and symmetry helpers used throughout ``bug_nicole``.
+
+The Julia code this package was ported from leans heavily on lightweight index
+wrappers and symmetry-aware site constructors. This module keeps that role, but
+spells the ideas out in plain Python so the rest of the code can use readable
+helpers instead of manipulating Nicole indices directly at every call site.
+"""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+from itertools import count
+from typing import Iterable, Literal
+
+from nicole import Direction, Index, Sector, U1Group
+
+SymmetryName = Literal["trivial", "u1"]
+
+_GROUP = U1Group()
+_FRESH = count()
+_SPIN_HALF_U1_SECTORS = (Sector(1, 1), Sector(-1, 1))
+
+__all__ = [
+ "Ix",
+ "SymmetryName",
+ "bond_index",
+ "fresh_itag",
+ "has_nontrivial_symmetry",
+ "idx",
+ "is_trivial_ix",
+ "normalize_symmetry",
+ "resolved_sectors",
+ "siteinds",
+ "spin_half_site_sectors",
+]
+
+
+def normalize_symmetry(symmetry: str) -> SymmetryName:
+ """Normalize public symmetry spellings to the package-internal name.
+
+ Args:
+ symmetry: User-facing symmetry label such as ``"u1"``, ``"sz"``,
+ ``"trivial"`` or ``"dense"``.
+
+ Returns:
+ ``"trivial"`` or ``"u1"``.
+ """
+ key = symmetry.strip().lower()
+ if key in {"trivial", "none", "dense"}:
+ return "trivial"
+ if key in {"u1", "sz", "u(1)"}:
+ return "u1"
+ raise ValueError(f"Unknown symmetry specification: {symmetry!r}")
+
+
+def spin_half_site_sectors(symmetry: str = "u1") -> tuple[Sector, ...]:
+ """Return the canonical spin-1/2 site sectors for one symmetry choice.
+
+ Args:
+ symmetry: Symmetry label understood by :func:`normalize_symmetry`.
+
+ Returns:
+ The sector tuple used for one spin-1/2 physical site.
+ """
+ key = normalize_symmetry(symmetry)
+ if key == "trivial":
+ return (Sector(0, 2),)
+ return _SPIN_HALF_U1_SECTORS
+
+
+@dataclass(frozen=True)
+class Ix:
+ """Lightweight Nicole-index handle used throughout the package.
+
+ Args:
+ itag: Nicole tag.
+ dim: Total index dimension.
+ direction: Nicole direction carried by the leg.
+ sectors: Optional explicit sector tuple. ``None`` means one trivial
+ dense sector of size ``dim``.
+ group: Nicole symmetry group object. The default is the package U(1)
+ group handle.
+ """
+
+ itag: str
+ dim: int
+ direction: Direction
+ sectors: tuple[Sector, ...] | None = None
+ group: object = _GROUP
+
+ def nicole(self) -> Index:
+ """Materialize this wrapper as a Nicole :class:`Index`.
+
+ Returns:
+ A Nicole index with the same tag metadata and sectors.
+ """
+ return Index(self.direction, self.group, resolved_sectors(self))
+
+ def resolved_sectors(self) -> tuple[Sector, ...]:
+ """Return the explicit sector tuple for this index.
+
+ Returns:
+ The stored sectors, or a single trivial sector when the index is
+ dense.
+ """
+ return resolved_sectors(self)
+
+ def reversed(self) -> "Ix":
+ """Return a copy whose Nicole direction is reversed.
+
+ Returns:
+ A new :class:`Ix` with the same metadata and opposite direction.
+ """
+ return Ix(
+ self.itag,
+ self.dim,
+ self.direction.reverse(),
+ self.sectors,
+ self.group,
+ )
+
+ def retag(self, itag: str) -> "Ix":
+ """Return a copy with a different Nicole tag.
+
+ Args:
+ itag: Replacement tag.
+
+ Returns:
+ A new :class:`Ix` with the requested tag.
+ """
+ return Ix(itag, self.dim, self.direction, self.sectors, self.group)
+
+ def is_trivial(self) -> bool:
+ """Return whether this index is one dense neutral sector.
+
+ Returns:
+ ``True`` when the index has only the neutral dense sector.
+ """
+ return is_trivial_ix(self)
+
+
+def resolved_sectors(ix: Ix | Index) -> tuple[Sector, ...]:
+ """Return the explicit sector tuple for an ``Ix`` or Nicole ``Index``.
+
+ Args:
+ ix: Wrapped or native Nicole index.
+
+ Returns:
+ An explicit tuple of Nicole sectors.
+ """
+ if isinstance(ix, Ix):
+ return ix.sectors if ix.sectors is not None else (Sector(0, ix.dim),)
+ return ix.sectors
+
+
+def is_trivial_ix(ix: Ix | Index) -> bool:
+ """Return whether an index carries only the neutral dense sector.
+
+ Args:
+ ix: Wrapped or native Nicole index.
+
+ Returns:
+ ``True`` when the sector structure is trivial.
+ """
+ sectors = resolved_sectors(ix)
+ return len(sectors) == 1 and sectors[0].charge == 0
+
+
+def has_nontrivial_symmetry(ixs: Iterable[Ix | Index]) -> bool:
+ """Return whether any index in a collection carries charge structure.
+
+ Args:
+ ixs: Iterable of wrapped or native Nicole indices.
+
+ Returns:
+ ``True`` when at least one index is not dense-trivial.
+ """
+ return any(not is_trivial_ix(ix) for ix in ixs)
+
+
+def idx(direction: Direction, dim: int, itag: str, **kwargs: object) -> Ix:
+ """Construct an :class:`Ix` using the field order most call sites prefer.
+
+ Args:
+ direction: Nicole direction for the index.
+ dim: Total index dimension.
+ itag: Nicole tag string.
+ **kwargs: Optional ``sectors=...`` and ``group=...`` overrides.
+
+ Returns:
+ A new :class:`Ix` wrapper.
+ """
+
+ sectors = kwargs.pop("sectors", None)
+ group = kwargs.pop("group", _GROUP)
+ if kwargs:
+ unknown = ", ".join(sorted(kwargs))
+ raise TypeError(f"Unknown idx option(s): {unknown}")
+ return Ix(itag, dim, direction, sectors, group)
+
+
+def bond_index(
+ itag: str,
+ direction: Direction,
+ charge_dims: Iterable[tuple[int, int]],
+ *,
+ group: object = _GROUP,
+) -> Ix:
+ """Build a bond index from explicit ``(charge, multiplicity)`` data.
+
+ Args:
+ itag: Nicole tag string.
+ direction: Nicole direction for the bond.
+ charge_dims: Iterable of ``(charge, dim)`` pairs.
+ group: Nicole symmetry group handle.
+
+ Returns:
+ A symmetry-aware :class:`Ix` wrapper for the bond.
+ """
+ sectors = tuple(Sector(int(charge), int(dim)) for charge, dim in charge_dims)
+ return Ix(itag, sum(int(sector.dim) for sector in sectors), direction, sectors, group)
+
+
+def fresh_itag(base: str) -> str:
+ """Generate a unique Nicole tag with a monotone suffix.
+
+ Args:
+ base: Prefix that should remain recognizable in debug output.
+
+ Returns:
+ A fresh tag such as ``"b3#17"``.
+ """
+ return f"{base}#{next(_FRESH)}"
+
+
+def siteinds(n: int, d: int = 2, *, symmetry: str = "trivial") -> list[Ix]:
+ """Build the canonical physical site indices ``s1, s2, ..., sN``.
+
+ Args:
+ n: Number of sites.
+ d: On-site Hilbert-space dimension.
+ symmetry: Symmetry label such as ``"trivial"`` or ``"u1"``.
+
+ Returns:
+ A list of OUT-directed physical site indices.
+ """
+ key = normalize_symmetry(symmetry)
+ if key == "u1":
+ if d != 2:
+ raise NotImplementedError("U(1) site indices currently support only spin-1/2 sites.")
+ sectors = spin_half_site_sectors("u1")
+ else:
+ sectors = (Sector(0, d),)
+
+ # Physical site legs always point outward in the MPS/MPO conventions used
+ # throughout this port.
+ return [Ix(f"s{k}", d, Direction.OUT, sectors) for k in range(1, n + 1)]
diff --git a/src/alice/algorithm/bond_update_bug/_kernel/kls/__init__.py b/src/alice/algorithm/bond_update_bug/_kernel/kls/__init__.py
new file mode 100644
index 0000000..a525147
--- /dev/null
+++ b/src/alice/algorithm/bond_update_bug/_kernel/kls/__init__.py
@@ -0,0 +1,55 @@
+# Copyright (C) 2025-2026 Changkai Zhang.
+#
+# This file is part of Alice project.
+#
+# Alice is free software: you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published
+# by the Free Software Foundation, either version 3 of the License,
+# or (at your option) any later version.
+#
+# Alice is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with Alice. If not, see .
+
+"""Local bond_update_bug K/L/S bond updates for dense and U(1)-symmetric tensors.
+
+This module contains the Python port of the local Lubich-style K/L/S update used
+by the bond_update_bug sweep. The user-facing helper is
+``_kls_local_bond_candidate``. Internally, the code is organized around one
+explicit concept:
+
+- ``LocalBondFrame`` gives names to the tensors and indices on the active bond
+ so the update logic reads like the algorithm rather than a raw dictionary walk.
+"""
+
+from .augment import (
+ _augmented_left_isometry_from_k,
+ _augmented_right_isometry_from_l,
+ _pick_left_update,
+ _pick_right_update,
+ _truncate_quantum_s_step,
+ _truncate_quantum_s_step_reverse,
+)
+from .candidate import _kls_local_bond_candidate
+from .symmetric_completion import (
+ _symmetric_augmented_left_isometry_from_k,
+ _symmetric_augmented_right_isometry_from_l,
+)
+from .frame import LocalBondFrame
+
+__all__ = [
+ "LocalBondFrame",
+ "_augmented_left_isometry_from_k",
+ "_augmented_right_isometry_from_l",
+ "_kls_local_bond_candidate",
+ "_pick_left_update",
+ "_pick_right_update",
+ "_symmetric_augmented_left_isometry_from_k",
+ "_symmetric_augmented_right_isometry_from_l",
+ "_truncate_quantum_s_step",
+ "_truncate_quantum_s_step_reverse",
+]
diff --git a/src/alice/algorithm/bond_update_bug/_kernel/kls/augment.py b/src/alice/algorithm/bond_update_bug/_kernel/kls/augment.py
new file mode 100644
index 0000000..3ed04bc
--- /dev/null
+++ b/src/alice/algorithm/bond_update_bug/_kernel/kls/augment.py
@@ -0,0 +1,459 @@
+# Copyright (C) 2025-2026 Changkai Zhang.
+#
+# This file is part of Alice project.
+#
+# Alice is free software: you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published
+# by the Free Software Foundation, either version 3 of the License,
+# or (at your option) any later version.
+#
+# Alice is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with Alice. If not, see .
+
+"""Augmentation and basis building logic for KLS updates."""
+
+from __future__ import annotations
+
+import math
+from typing import Any
+
+import torch
+from nicole import Sector, Tensor, decomp, einsum
+
+from ..indices import Ix, fresh_itag, resolved_sectors
+from ..krylov import active_time_prefactor, linear_substep, tensor_lanczos_expv
+from ..linalg import (
+ identity_overlap_matrix,
+ qr_column_basis,
+ qr_row_basis,
+)
+from ..nicole_helpers import dag, flatten_fortran, make_tensor, reshape_fortran, to_dense, tcontract
+from .frame import (
+ LocalBondFrame,
+ _apply_gate_named,
+ _clone_tensor_with_ixs,
+ _dense_from_tensor_with_ixs,
+ _left_row_indices_by_flux,
+ _right_col_indices_by_flux,
+ _sector_offsets,
+ _tensor_ix,
+)
+
+
+def _tensor_expv(
+ apply,
+ dt: complex,
+ tensor: Tensor,
+ lanczos_maxiter: int = 30,
+ lanczos_tol: float = 1e-15,
+) -> Tensor:
+ """Apply a Lanczos ``expv`` step to a Nicole tensor.
+
+ Args:
+ apply: Matrix-free tensor action representing the local Hamiltonian.
+ dt: Local timestep for this substep.
+ tensor: Input state tensor.
+ lanczos_maxiter: Maximum Lanczos iterations per local substep.
+ lanczos_tol: Lanczos termination tolerance.
+
+ Returns:
+ The evolved tensor after ``exp(prefactor * dt * H)``.
+ """
+ return tensor_lanczos_expv(
+ apply,
+ active_time_prefactor() * dt,
+ tensor,
+ maxiter=lanczos_maxiter,
+ tol=lanczos_tol,
+ )
+
+
+def _collect_tensor_krylov_directions(
+ seed: Tensor,
+ apply,
+ dt: complex,
+ aug_krylov_depth: int = 1,
+ lanczos_maxiter: int = 30,
+ lanczos_tol: float = 1e-15,
+) -> list[Tensor]:
+ """Collect the K/L Krylov directions used for local basis growth.
+
+ Args:
+ seed: Input tensor for the K or L substep.
+ apply: Matrix-free tensor action for the corresponding projected local
+ Hamiltonian.
+ dt: Local timestep used for the first Krylov direction.
+ aug_krylov_depth: Number of K/L Krylov directions stacked before basis extraction.
+ lanczos_maxiter: Maximum Lanczos iterations per local substep.
+ lanczos_tol: Lanczos termination tolerance.
+
+ Returns:
+ A list containing the evolved first direction followed by repeated
+ projected-Hamiltonian applications when ``aug_krylov_depth > 1``.
+ """
+ first_direction = _tensor_expv(apply, dt, seed, lanczos_maxiter=lanczos_maxiter, lanczos_tol=lanczos_tol)
+ directions = [first_direction]
+ next_direction = first_direction
+ for _ in range(2, aug_krylov_depth + 1):
+ next_direction = apply(next_direction)
+ directions.append(next_direction)
+ return directions
+
+
+def _pick_left_update(
+ U0_mat: torch.Tensor,
+ K1_mat: torch.Tensor,
+ augment: bool = True,
+ max_rank: int | float = math.inf,
+ aug_tol: float = 1e-12,
+):
+ """Choose an augmented left basis and its overlap with ``U0_mat``.
+
+ Args:
+ U0_mat: Current left isometry as a dense matrix.
+ K1_mat: Candidate K-step directions as dense columns.
+ augment: Whether new Krylov directions may enlarge the basis.
+ max_rank: Hard cap on the returned basis rank.
+ aug_tol: Threshold used to discard nearly dependent directions.
+
+ Returns:
+ ``(basis, overlap, n_new)`` for the chosen left basis.
+ """
+
+ if not augment or K1_mat.numel() == 0 or K1_mat.shape[1] == 0:
+ overlap = identity_overlap_matrix(U0_mat.dtype, U0_mat.shape[1], device=U0_mat.device)
+ return U0_mat, overlap, 0
+
+ # Sulz augmented BUG: no pre-filter against U0. Orthonormalise K1 (rank-revealed
+ # at machine eps by qr_column_basis), stack with U0, and take the rank-revealing
+ # QR range basis of [U0 | K1] (rank <= 2r). Redundant/near-dependent directions are
+ # removed by the QR's own non-zero-R-norm rank count; final rank control happens at
+ # the post-S-step SVD truncation (no aug_tol heuristic discard).
+ Qk, _ = qr_column_basis(K1_mat)
+ cand = torch.cat([U0_mat, Qk], dim=1) if Qk.numel() else U0_mat
+ Q, _ = qr_column_basis(cand)
+ if max_rank is not math.inf:
+ Q = Q[:, : min(Q.shape[1], int(max_rank))]
+ overlap = Q.conj().transpose(0, 1) @ U0_mat
+ n_new = max(0, Q.shape[1] - U0_mat.shape[1])
+ return Q, overlap, n_new
+
+
+def _pick_right_update(
+ V0_mat: torch.Tensor,
+ L1_mat: torch.Tensor,
+ augment: bool = True,
+ max_rank: int | float = math.inf,
+ aug_tol: float = 1e-12,
+):
+ """Choose an augmented right basis and its overlap with ``V0_mat``.
+
+ Args:
+ V0_mat: Current right isometry as a dense matrix.
+ L1_mat: Candidate L-step directions as dense rows.
+ augment: Whether new Krylov directions may enlarge the basis.
+ max_rank: Hard cap on the returned basis rank.
+ aug_tol: Threshold used to discard nearly dependent directions.
+
+ Returns:
+ ``(basis, overlap, n_new)`` for the chosen right basis.
+ """
+
+ if not augment or L1_mat.numel() == 0 or L1_mat.shape[0] == 0:
+ overlap = identity_overlap_matrix(V0_mat.dtype, V0_mat.shape[0], device=V0_mat.device)
+ return V0_mat, overlap, 0
+
+ # Sulz augmented BUG (row mirror of _pick_left_update): no pre-filter against V0.
+ Ql, _ = qr_row_basis(L1_mat)
+ cand = torch.cat([V0_mat, Ql], dim=0) if Ql.numel() else V0_mat
+ Q, _ = qr_row_basis(cand)
+ if max_rank is not math.inf:
+ Q = Q[: min(Q.shape[0], int(max_rank)), :]
+ overlap = V0_mat @ Q.conj().transpose(0, 1)
+ n_new = max(0, Q.shape[0] - V0_mat.shape[0])
+ return Q, overlap, n_new
+
+
+def _left_tensor_matrix(U_tens, link_l: Ix, site_l: Ix, mid: Ix):
+ """Reshape a left tensor `(link_l, site_l, mid)` into matrix form."""
+ block = _dense_from_tensor_with_ixs(U_tens, [link_l, site_l, mid]).to(torch.complex128)
+ return reshape_fortran(block, (link_l.dim * site_l.dim, mid.dim))
+
+
+def _right_tensor_matrix(V_tens, mid: Ix, site_r: Ix, link_r: Ix):
+ """Reshape a right tensor `(mid, site_r, link_r)` into matrix form."""
+ block = _dense_from_tensor_with_ixs(V_tens, [mid, site_r, link_r]).to(torch.complex128)
+ return reshape_fortran(block, (mid.dim, site_r.dim * link_r.dim))
+
+
+def _augmented_left_isometry_from_k(
+ U0_tens,
+ K1_tens,
+ *args,
+ link_l=None,
+ site_l=None,
+ old_mid=None,
+ augment: bool = True,
+ max_rank: int | float = math.inf,
+ aug_tol: float = 1e-12,
+ **kwargs: Any,
+):
+ """Build the augmented left isometry tensor from K-step directions.
+
+ Args:
+ U0_tens: Current left canonical factor.
+ K1_tens: Stacked K-step directions.
+ *args: Legacy positional tail ``(link_l, site_l, old_mid)``.
+ link_l: Left bond index.
+ site_l: Left physical site index.
+ old_mid: Current middle bond index.
+ augment: Whether new Krylov directions may enlarge the basis.
+ max_rank: Hard cap on the returned basis rank.
+ aug_tol: Threshold used to discard nearly dependent directions.
+ **kwargs: Keyword overrides for index parameters.
+
+ Returns:
+ ``(U_aug_tens, overlap_tens, n_new)``.
+ """
+
+ if args:
+ if len(args) != 3:
+ raise TypeError("_augmented_left_isometry_from_k expects (link_l, site_l, old_mid) after the tensors.")
+ if any(name in ("link_l", "site_l", "old_mid") for name in kwargs):
+ raise TypeError("Provide left-augmentation indices either positionally or by keyword, not both.")
+ link_l, site_l, old_mid = args
+ if link_l is None or site_l is None or old_mid is None:
+ try:
+ link_l = kwargs.pop("link_l") if link_l is None else link_l
+ site_l = kwargs.pop("site_l") if site_l is None else site_l
+ old_mid = kwargs.pop("old_mid") if old_mid is None else old_mid
+ except KeyError as exc:
+ raise TypeError("Missing left-augmentation index input.") from exc
+ if kwargs:
+ unknown = ", ".join(sorted(kwargs))
+ raise TypeError(f"Unknown left-augmentation option(s): {unknown}")
+
+ U0_mat = _left_tensor_matrix(U0_tens, link_l, site_l, old_mid)
+ K1_mat = _left_tensor_matrix(K1_tens, link_l, site_l, old_mid)
+ U1_mat, overlap_mat, n_new = _pick_left_update(U0_mat, K1_mat, augment=augment, max_rank=max_rank, aug_tol=aug_tol)
+ new_mid = Ix(old_mid.itag, U1_mat.shape[1], old_mid.direction)
+ U1_tens = make_tensor(
+ reshape_fortran(U1_mat, (link_l.dim, site_l.dim, new_mid.dim)),
+ [link_l, site_l, new_mid],
+ dtype=torch.complex128,
+ )
+ overlap_tens = make_tensor(
+ overlap_mat,
+ [Ix(new_mid.itag, new_mid.dim, old_mid.direction), old_mid],
+ dtype=torch.complex128,
+ )
+ return U1_tens, overlap_tens, n_new
+
+
+def _augmented_right_isometry_from_l(
+ V0_tens,
+ L1_tens,
+ *args,
+ old_mid=None,
+ site_r=None,
+ link_r=None,
+ augment: bool = True,
+ max_rank: int | float = math.inf,
+ aug_tol: float = 1e-12,
+ **kwargs: Any,
+):
+ """Build the augmented right isometry tensor from L-step directions.
+
+ Args:
+ V0_tens: Current right canonical factor.
+ L1_tens: Stacked L-step directions.
+ *args: Legacy positional tail ``(old_mid, site_r, link_r)``.
+ old_mid: Current middle bond index.
+ site_r: Right physical site index.
+ link_r: Right bond index.
+ augment: Whether new Krylov directions may enlarge the basis.
+ max_rank: Hard cap on the returned basis rank.
+ aug_tol: Threshold used to discard nearly dependent directions.
+ **kwargs: Keyword overrides for index parameters.
+
+ Returns:
+ ``(V_aug_tens, overlap_tens, n_new)``.
+ """
+
+ if args:
+ if len(args) != 3:
+ raise TypeError("_augmented_right_isometry_from_l expects (old_mid, site_r, link_r) after the tensors.")
+ if any(name in ("old_mid", "site_r", "link_r") for name in kwargs):
+ raise TypeError("Provide right-augmentation indices either positionally or by keyword, not both.")
+ old_mid, site_r, link_r = args
+ if old_mid is None or site_r is None or link_r is None:
+ try:
+ old_mid = kwargs.pop("old_mid") if old_mid is None else old_mid
+ site_r = kwargs.pop("site_r") if site_r is None else site_r
+ link_r = kwargs.pop("link_r") if link_r is None else link_r
+ except KeyError as exc:
+ raise TypeError("Missing right-augmentation index input.") from exc
+ if kwargs:
+ unknown = ", ".join(sorted(kwargs))
+ raise TypeError(f"Unknown right-augmentation option(s): {unknown}")
+
+ V0_mat = _right_tensor_matrix(V0_tens, old_mid, site_r, link_r)
+ L1_mat = _right_tensor_matrix(L1_tens, old_mid, site_r, link_r)
+ V1_mat, overlap_mat, n_new = _pick_right_update(V0_mat, L1_mat, augment=augment, max_rank=max_rank, aug_tol=aug_tol)
+ new_mid = Ix(old_mid.itag, V1_mat.shape[0], old_mid.direction)
+ V1_tens = make_tensor(
+ reshape_fortran(V1_mat, (new_mid.dim, site_r.dim, link_r.dim)),
+ [new_mid, site_r, link_r],
+ dtype=torch.complex128,
+ )
+ overlap_tens = make_tensor(
+ overlap_mat,
+ [old_mid, Ix(new_mid.itag, new_mid.dim, old_mid.direction)],
+ dtype=torch.complex128,
+ )
+ return V1_tens, overlap_tens, n_new
+
+
+def _transported_s_start_from_augmented_bases(U_basis, V_basis, theta0_tens, *args, **kwargs):
+ """Project ``theta0_tens`` into augmented bases to form an initial S tensor.
+
+ Args:
+ U_basis: Left augmented basis matrix.
+ V_basis: Right augmented basis matrix.
+ theta0_tens: Two-site tensor to project.
+ *args: Legacy positional tail ``(link_l, site_l, site_r, link_r)``.
+ **kwargs: Keyword form of the same four indices.
+
+ Returns:
+ Rank-2 Nicole tensor containing the projected S data.
+ """
+
+ if args:
+ if len(args) != 4:
+ raise TypeError("_transported_s_start_from_augmented_bases expects four index arguments.")
+ if any(name in kwargs for name in ("link_l", "site_l", "site_r", "link_r")):
+ raise TypeError("Provide transport indices either positionally or by keyword, not both.")
+ kwargs.update({"link_l": args[0], "site_l": args[1], "site_r": args[2], "link_r": args[3]})
+ try:
+ link_l = kwargs.pop("link_l")
+ site_l = kwargs.pop("site_l")
+ site_r = kwargs.pop("site_r")
+ link_r = kwargs.pop("link_r")
+ except KeyError as exc:
+ raise TypeError("Missing transported-S basis index input.") from exc
+ if kwargs:
+ unknown = ", ".join(sorted(kwargs))
+ raise TypeError(f"Unknown transported-S option(s): {unknown}")
+
+ theta = to_dense(theta0_tens, [link_l.itag, site_l.itag, site_r.itag, link_r.itag]).to(torch.complex128)
+ theta_mat = reshape_fortran(theta, (link_l.dim * site_l.dim, site_r.dim * link_r.dim))
+ S = U_basis.conj().transpose(0, 1) @ theta_mat @ V_basis.conj().transpose(0, 1)
+ return make_tensor(
+ S,
+ [Ix("s_mid_l", U_basis.shape[1], link_l.direction), Ix("s_mid_r", V_basis.shape[0], link_r.direction)],
+ dtype=torch.complex128,
+ )
+
+
+def _advance_s_tensor_in_bases(H_eff_mat, dt: complex, S_old_tens):
+ """Evolve the S tensor with `linear_substep(..., method='expv')`."""
+ block = to_dense(S_old_tens, list(S_old_tens.itags)).to(torch.complex128)
+ s_old = flatten_fortran(block)
+ s_new, numops = linear_substep(
+ H_eff_mat,
+ active_time_prefactor() * dt,
+ s_old,
+ method="expv",
+ lanczos_tol=1e-14,
+ lanczos_maxiter=max(4, len(s_old)),
+ )
+ shape = block.shape
+ out = make_tensor(
+ reshape_fortran(s_new, shape),
+ [
+ Ix(S_old_tens.itags[0], shape[0], S_old_tens.indices[0].direction),
+ Ix(S_old_tens.itags[1], shape[1], S_old_tens.indices[1].direction),
+ ],
+ dtype=torch.complex128,
+ )
+ return out, numops
+
+
+def _truncate_quantum_s_step(S_new_tens, maxdim: int):
+ """Truncate `S_new_tens` by SVD and return split factors for write-back."""
+ block = to_dense(S_new_tens, list(S_new_tens.itags)).to(torch.complex128)
+ U, s, Vh = torch.linalg.svd(block, full_matrices=False)
+ keep = min(int(s.numel()), int(maxdim))
+ U_s = U[:, :keep]
+ SV = torch.diag(s[:keep]) @ Vh[:keep, :]
+ U_tens = make_tensor(
+ U_s,
+ [Ix(S_new_tens.itags[0], U_s.shape[0], S_new_tens.indices[0].direction), Ix("keep", keep, S_new_tens.indices[0].direction.reverse())],
+ dtype=torch.complex128,
+ )
+ SV_tens = make_tensor(
+ SV,
+ [Ix("keep", keep, S_new_tens.indices[1].direction), Ix(S_new_tens.itags[1], SV.shape[1], S_new_tens.indices[1].direction)],
+ dtype=torch.complex128,
+ )
+ return U_tens, SV_tens, keep, s
+
+
+def _truncate_quantum_s_step_reverse(S_new_tens, maxdim: int):
+ """Reverse-sweep alias of `_truncate_quantum_s_step`."""
+ return _truncate_quantum_s_step(S_new_tens, maxdim)
+
+
+def _stack_left_krylov_directions(directions, link_l: Ix, site_l: Ix, mid_k: Ix):
+ if len(directions) == 1:
+ return directions[0], mid_k
+
+ mats = [_left_tensor_matrix(direction, link_l, site_l, mid_k) for direction in directions]
+ sectors = tuple(Sector(int(sec.charge), int(sec.dim * len(directions))) for sec in resolved_sectors(mid_k))
+ ext_mid = Ix(fresh_itag(mid_k.itag), sum(sec.dim for sec in sectors), mid_k.direction, sectors, mid_k.group)
+ old_offsets = _sector_offsets(mid_k)
+ new_offsets = _sector_offsets(ext_mid)
+ stacked = torch.zeros((link_l.dim * site_l.dim, ext_mid.dim), dtype=torch.complex128, device=mats[0].device)
+ for sec in resolved_sectors(mid_k):
+ old_start, old_dim = old_offsets[sec.charge]
+ new_start, _ = new_offsets[sec.charge]
+ old_sl = slice(old_start, old_start + old_dim)
+ for depth, mat in enumerate(mats):
+ new_sl = slice(new_start + depth * old_dim, new_start + (depth + 1) * old_dim)
+ stacked[:, new_sl] = mat[:, old_sl]
+ tensor = make_tensor(
+ reshape_fortran(stacked, (link_l.dim, site_l.dim, ext_mid.dim)),
+ [link_l, site_l, ext_mid],
+ dtype=torch.complex128,
+ )
+ return tensor, ext_mid
+
+
+def _stack_right_krylov_directions(directions, mid_l: Ix, site_r: Ix, link_r: Ix):
+ if len(directions) == 1:
+ return directions[0], mid_l
+
+ mats = [_right_tensor_matrix(direction, mid_l, site_r, link_r) for direction in directions]
+ sectors = tuple(Sector(int(sec.charge), int(sec.dim * len(directions))) for sec in resolved_sectors(mid_l))
+ ext_mid = Ix(fresh_itag(mid_l.itag), sum(sec.dim for sec in sectors), mid_l.direction, sectors, mid_l.group)
+ old_offsets = _sector_offsets(mid_l)
+ new_offsets = _sector_offsets(ext_mid)
+ stacked = torch.zeros((ext_mid.dim, site_r.dim * link_r.dim), dtype=torch.complex128, device=mats[0].device)
+ for sec in resolved_sectors(mid_l):
+ old_start, old_dim = old_offsets[sec.charge]
+ new_start, _ = new_offsets[sec.charge]
+ old_sl = slice(old_start, old_start + old_dim)
+ for depth, mat in enumerate(mats):
+ new_sl = slice(new_start + depth * old_dim, new_start + (depth + 1) * old_dim)
+ stacked[new_sl, :] = mat[old_sl, :]
+ tensor = make_tensor(
+ reshape_fortran(stacked, (ext_mid.dim, site_r.dim, link_r.dim)),
+ [ext_mid, site_r, link_r],
+ dtype=torch.complex128,
+ )
+ return tensor, ext_mid
diff --git a/src/alice/algorithm/bond_update_bug/_kernel/kls/candidate.py b/src/alice/algorithm/bond_update_bug/_kernel/kls/candidate.py
new file mode 100644
index 0000000..3fc2a71
--- /dev/null
+++ b/src/alice/algorithm/bond_update_bug/_kernel/kls/candidate.py
@@ -0,0 +1,243 @@
+# Copyright (C) 2025-2026 Changkai Zhang.
+#
+# This file is part of Alice project.
+#
+# Alice is free software: you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published
+# by the Free Software Foundation, either version 3 of the License,
+# or (at your option) any later version.
+#
+# Alice is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with Alice. If not, see .
+# Author of code: Madhav Menon.
+
+
+"""The `bond_update_bug` local K/L/S bond candidate.
+
+The discarded-projector Basis-Update & Galerkin update for one bond, on the state
+``Θ0 = U0 · S0 · V0``. Two features define it:
+
+1. **Project-before.** The discarded (orthogonal-complement) projector is applied
+ to the K/L *generator* before the exponential, not to the integrated factor.
+ The K generator becomes ``G_K = P⊥_U0 · H_K`` with ``P⊥_U0 = I − U0 U0†`` and
+ the L generator ``G_L = H_L · P⊥_V0`` with ``P⊥_V0 = I − V0† V0``. Because the
+ projected generator is non-Hermitian, the K/L substep uses the general
+ (``issymmetric=False``) Krylov path — a symmetry-preserving tensor Arnoldi
+ exponential — rather than the Hermitian Lanczos.
+
+2. **Act the augmented isometries, no overlap matrices.** The new directions are
+ isolated by the discarded projector and stacked onto the old isometry to form
+ ``Û = [U0 | Qk]`` / ``V̂ = [V0 ; Ql]``. The S-step then projects the *current*
+ two-site tensor directly onto the augmented bases, ``Ŝ0 = Û† Θ0 V̂†``, evolves
+ it in the augmented basis (the Hermitian Galerkin generator), and truncates
+ with an SVD.
+
+This is the Alice realisation of the reference Julia ``bond_update_bug!`` per-bond
+candidate. The Nicole tensor helpers, the Krylov ``expv`` substeps, the QR/SVD
+linear algebra, and the augmented-isometry construction are all shared with the
+rest of the ``_kernel`` subpackage.
+"""
+
+from __future__ import annotations
+
+import math
+from typing import Any
+
+from nicole import Tensor, decomp
+
+from ..indices import Ix, fresh_itag
+from ..krylov import active_time_prefactor
+from ..local_solvers import local_expv
+from ..nicole_helpers import dag, tcontract
+from .frame import (
+ LocalBondFrame,
+ _apply_gate_named,
+ _clone_tensor_with_ixs,
+ _singular_values_from_diag_tensor,
+ _tensor_ix,
+)
+from .symmetric_completion import (
+ _symmetric_augmented_left_isometry_from_k,
+ _symmetric_augmented_right_isometry_from_l,
+)
+
+
+def _discarded_local_bond_candidate(
+ frame: LocalBondFrame,
+ gate: Tensor,
+ dt: complex,
+ maxdim: int = 200,
+ s_dt: complex | None = None,
+ augment: bool = True,
+ aug_krylov_depth: int = 1,
+ aug_tol: float = 1e-12,
+ trunc_thresh: float | None = None,
+ lanczos_tol: float = 1e-15,
+ lanczos_maxiter: int = 30,
+ solver: str = 'krylov',
+ solver_substeps: int = 1,
+ kl_cutoff: float | None = None,
+):
+ """Run one discarded-projector K/L/S local update (see module docstring).
+
+ The K/L/S local exponentials are computed by the selected ``solver`` (see
+ :mod:`alice.algorithm.bond_update_bug._kernel.local_solvers`): ``'krylov'`` is the
+ exact reference, ``'midpoint'``/``'rk4'`` are explicit RK with ``solver_substeps``
+ internal steps, and ``'trapezoid'`` is the A-stable Crank–Nicolson rule. In
+ imaginary time the evolution is non-unitary so any stable integrator is valid.
+ """
+ s_dt_eff = dt if s_dt is None else s_dt
+ augment_left_here = augment and frame.old_rank < frame.left_capacity
+ augment_right_here = augment and frame.old_rank < frame.right_capacity
+ prefactor = active_time_prefactor()
+
+ # ---- K-step: project-before, then integrate K0 = U0·S0 ----
+ # H_K x = V0†-projected gate action; G_K x = P⊥_U0 (H_K x), P⊥_U0 = I − U0 U0†.
+ # The projected generator is NON-Hermitian, so we use a symmetry-preserving
+ # tensor Arnoldi exponential (never densifying to the standard basis, which
+ # would break the U(1) block structure of the Nicole tensor).
+ K0_tens = tcontract(frame.U0_tens, frame.S0_tens) # (link_l, site_l, mid_k)
+ mid_k = _tensor_ix(K0_tens, 2)
+
+ def apply_gk(x_tens: Tensor) -> Tensor:
+ #get discarded projector
+ theta = tcontract(x_tens, frame.V0_tens)
+ evolved = _apply_gate_named(gate, theta, frame.site_l.itag, frame.site_r.itag)
+ HK = tcontract(evolved, dag(frame.V0_tens)) # H_K x on (link_l, site_l, mid_k)
+ # P⊥_U0 on (link_l, site_l): HK − U0 (U0† HK).
+ return HK - tcontract(frame.U0_tens, tcontract(dag(frame.U0_tens), HK))
+
+ K1_tens = local_expv(apply_gk, prefactor * dt, K0_tens,
+ solver=solver, substeps=solver_substeps, hermitian=False,
+ krylov_maxiter=lanczos_maxiter, krylov_tol=lanczos_tol)
+ # Direct sum Û = [U0 | Qk], built per U(1) charge sector so the Nicole block
+ # structure stays valid (a symmetry-blind dense QR would mix sectors and be
+ # rejected). No overlap matrix M̂ is formed — the discarded variant projects
+ # Θ0 onto the augmented bases directly in the S-step below.
+ U_aug_tens, _M_hat, n_new_k = _symmetric_augmented_left_isometry_from_k(
+ frame.U0_tens, K1_tens, frame.link_l, frame.site_l, frame.canon_u0, mid_k,
+ augment=augment_left_here, max_rank=math.inf, aug_tol=aug_tol, kl_cutoff=kl_cutoff)
+
+ # ---- L-step: project-before, then integrate L0 = S0·V0 ----
+ L0_tens = tcontract(frame.S0_tens, frame.V0_tens) # (mid_l, site_r, link_r)
+ mid_l = _tensor_ix(L0_tens, 0)
+
+ def apply_gl(x_tens: Tensor) -> Tensor:
+ theta = tcontract(frame.U0_tens, x_tens)
+ evolved = _apply_gate_named(gate, theta, frame.site_l.itag, frame.site_r.itag)
+ HL = tcontract(dag(frame.U0_tens), evolved) # H_L x on (mid_l, site_r, link_r)
+ # P⊥_V0 on (site_r, link_r): HL − (HL V0†) V0.
+ return HL - tcontract(tcontract(HL, dag(frame.V0_tens)), frame.V0_tens)
+
+ L1_tens = local_expv(apply_gl, prefactor * dt, L0_tens,
+ solver=solver, substeps=solver_substeps, hermitian=False,
+ krylov_maxiter=lanczos_maxiter, krylov_tol=lanczos_tol)
+ V_aug_tens, _N_hat, n_new_l = _symmetric_augmented_right_isometry_from_l(
+ frame.V0_tens, L1_tens, frame.canon_v0, mid_l, frame.site_r, frame.link_r,
+ augment=augment_right_here, max_rank=math.inf, aug_tol=aug_tol, kl_cutoff=kl_cutoff)
+
+ # ---- S-step: project Θ0 directly onto the augmented bases (no M̂/N̂), evolve ----
+ # Ŝ0 = Û† Θ0 V̂† as a tensor contraction. dag(U_aug) exposes the augmented left
+ # mid-leg, dag(V_aug) the augmented right mid-leg, so Ŝ0 is automatically tagged
+ # to contract back with U_aug_tens / V_aug_tens in apply_s_tensor below.
+ theta0_tens = tcontract(tcontract(frame.U0_tens, frame.S0_tens), frame.V0_tens)
+ S_start_tens = tcontract(tcontract(dag(U_aug_tens), theta0_tens), dag(V_aug_tens))
+
+ def apply_s_tensor(x_tens: Tensor) -> Tensor:
+ theta = tcontract(tcontract(U_aug_tens, x_tens), V_aug_tens)
+ evolved = _apply_gate_named(gate, theta, frame.site_l.itag, frame.site_r.itag)
+ projected = tcontract(dag(U_aug_tens), evolved)
+ return tcontract(projected, dag(V_aug_tens))
+
+ # S-step generator is Hermitian (the Galerkin generator on the
+ # augmented bases); imaginary time makes the flow a contraction either way.
+ S_new_tens = local_expv(apply_s_tensor, prefactor * s_dt_eff, S_start_tens,
+ solver=solver, substeps=solver_substeps, hermitian=True,
+ krylov_maxiter=lanczos_maxiter, krylov_tol=lanczos_tol)
+
+ # ---- truncate: SVD sets the new (rank-adaptive) bond dimension ----
+ # Done in the symmetry-blocked Nicole representation (mirrors the
+ # kernel's S-step split), so the kept rank respects the U(1) sectors.
+ final_left_tag = fresh_itag(frame.link_mid.itag)
+ final_right_tag = fresh_itag(frame.link_mid.itag)
+ U_s, Sdiag, Vh = decomp(
+ S_new_tens, 0, mode="SVD",
+ itag=(final_left_tag, final_right_tag),
+ trunc={
+ "nkeep": int(maxdim),
+ "thresh": max(float(aug_tol if trunc_thresh is None else trunc_thresh), 1e-14),
+ },
+ )
+ left_tmp = tcontract(U_aug_tens, U_s)
+ right_tmp = tcontract(tcontract(Sdiag, Vh, axes=([1], [0])), V_aug_tens)
+ left_tmp.retag({final_left_tag: frame.link_mid.itag})
+ right_tmp.retag({final_left_tag: frame.link_mid.itag})
+
+ new_bond = Ix(frame.link_mid.itag, int(left_tmp.indices[2].dim), left_tmp.indices[2].direction,
+ left_tmp.indices[2].sectors, left_tmp.indices[2].group)
+ right_bond = Ix(frame.link_mid.itag, int(right_tmp.indices[0].dim), right_tmp.indices[0].direction,
+ right_tmp.indices[0].sectors, right_tmp.indices[0].group)
+ left_core = _clone_tensor_with_ixs(left_tmp, [frame.link_l, frame.site_l, new_bond])
+ right_core = _clone_tensor_with_ixs(right_tmp, [right_bond, frame.site_r, frame.link_r])
+ svals = _singular_values_from_diag_tensor(Sdiag)
+
+ return {
+ "left_core": left_core,
+ "right_core": right_core,
+ "U_aug_tens": U_aug_tens,
+ "V_aug_tens": V_aug_tens,
+ "S_new": S_new_tens,
+ "n_new_k": int(n_new_k),
+ "n_new_l": int(n_new_l),
+ "keep": int(left_core.indices[2].dim),
+ "svals": svals,
+ }
+
+
+def _kls_local_bond_candidate(
+ bond_data: dict[str, Any],
+ *,
+ gate,
+ dt: complex,
+ maxdim: int = 200,
+ s_dt: complex | None = None,
+ augment: bool = True,
+ aug_krylov_depth: int = 1,
+ aug_tol: float = 1e-12,
+ trunc_thresh: float | None = None,
+ lanczos_tol: float = 1e-15,
+ lanczos_maxiter: int = 30,
+ solver: str = 'krylov',
+ solver_substeps: int = 1,
+ kl_cutoff: float | None = None,
+ **kwargs: Any,
+):
+ """Return the discarded-projector BUG candidate on one bond.
+
+ Mirrors the call surface of
+ :func:`alice.algorithm.bond_update_bug._kernel.kls.candidate._kls_local_bond_candidate`
+ so the odd/even sweep can swap kernels without any other change. ``solver`` and
+ ``solver_substeps`` select the local (imaginary-time) integrator for the K/L/S
+ substeps (see :mod:`alice.algorithm.bond_update_bug._kernel.local_solvers`).
+ """
+ if aug_krylov_depth != 1:
+ raise ValueError("discarded variant currently supports aug_krylov_depth == 1 only.")
+ kwargs.pop("substep_method", None)
+ kwargs.pop("matrixfree_sstep", None)
+ if kwargs:
+ unknown = ", ".join(sorted(kwargs))
+ raise TypeError(f"Unknown discarded variant option(s): {unknown}")
+
+ frame = LocalBondFrame.from_mapping(bond_data)
+ return _discarded_local_bond_candidate(
+ frame, gate, dt,
+ maxdim=maxdim, s_dt=s_dt, augment=augment, aug_krylov_depth=aug_krylov_depth,
+ aug_tol=aug_tol, trunc_thresh=trunc_thresh,
+ lanczos_tol=lanczos_tol, lanczos_maxiter=lanczos_maxiter,
+ solver=solver, solver_substeps=solver_substeps, kl_cutoff=kl_cutoff,
+ )
diff --git a/src/alice/algorithm/bond_update_bug/_kernel/kls/frame.py b/src/alice/algorithm/bond_update_bug/_kernel/kls/frame.py
new file mode 100644
index 0000000..bbaef0d
--- /dev/null
+++ b/src/alice/algorithm/bond_update_bug/_kernel/kls/frame.py
@@ -0,0 +1,192 @@
+# Copyright (C) 2025-2026 Changkai Zhang.
+#
+# This file is part of Alice project.
+#
+# Alice is free software: you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published
+# by the Free Software Foundation, either version 3 of the License,
+# or (at your option) any later version.
+#
+# Alice is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with Alice. If not, see .
+
+"""LocalBondFrame and basic tensor/index utilities for KLS updates."""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+from typing import Any, Mapping
+
+import torch
+from nicole import Tensor, einsum
+
+from ..indices import Ix, has_nontrivial_symmetry, resolved_sectors
+
+
+@dataclass(frozen=True)
+class LocalBondFrame:
+ """Canonical two-site data needed by one local KLS update.
+
+ Args:
+ link_l: Left bond index entering the active two-site block.
+ link_mid: Bond index between the active left and right sites.
+ link_r: Right bond index exiting the active two-site block.
+ site_l: Left physical site index.
+ site_r: Right physical site index.
+ U0_tens: Left canonical isometry.
+ V0_tens: Right canonical isometry.
+ S0_tens: Bond-center tensor between the canonical frames.
+ canon_u0: Middle index carried by ``U0_tens`` and ``S0_tens``.
+ canon_v0: Middle index carried by ``S0_tens`` and ``V0_tens``.
+ theta0_tens: Optional assembled two-site tensor.
+ """
+
+ link_l: Ix
+ link_mid: Ix
+ link_r: Ix
+ site_l: Ix
+ site_r: Ix
+ U0_tens: Tensor
+ V0_tens: Tensor
+ S0_tens: Tensor
+ canon_u0: Ix
+ canon_v0: Ix
+ theta0_tens: Tensor | None = None
+
+ @classmethod
+ def from_mapping(cls, data: Mapping[str, Any]) -> "LocalBondFrame":
+ """Create a frame object from the historical snapshot dictionary.
+
+ Args:
+ data: Snapshot dictionary produced by the bug or environment code.
+
+ Returns:
+ A :class:`LocalBondFrame` instance with the expected fields.
+ """
+ return cls(
+ link_l=data["link_l"],
+ link_mid=data["link_mid"],
+ link_r=data["link_r"],
+ site_l=data["site_l"],
+ site_r=data["site_r"],
+ U0_tens=data["U0_tens"],
+ V0_tens=data["V0_tens"],
+ S0_tens=data["S0_tens"],
+ canon_u0=data["canon_u0"],
+ canon_v0=data["canon_v0"],
+ theta0_tens=data.get("theta0_tens"),
+ )
+
+ @property
+ def old_rank(self) -> int:
+ """Return the current middle-bond rank."""
+ return int(self.link_mid.dim)
+
+ @property
+ def left_capacity(self) -> int:
+ """Return the maximum admissible left-frame rank ``dim(link_l)*dim(site_l)``."""
+ return int(self.link_l.dim * self.site_l.dim)
+
+ @property
+ def right_capacity(self) -> int:
+ """Return the maximum admissible right-frame rank ``dim(site_r)*dim(link_r)``."""
+ return int(self.site_r.dim * self.link_r.dim)
+
+ def has_symmetry(self) -> bool:
+ """Return ``True`` when any leg on the active bond carries nontrivial symmetry."""
+ return has_nontrivial_symmetry([self.link_l, self.link_mid, self.link_r, self.site_l, self.site_r])
+
+
+def _clone_tensor_with_ixs(tensor: Tensor, ixs: list[Ix]) -> Tensor:
+ indices = tuple(ix.nicole() for ix in ixs)
+ itags = tuple(ix.itag for ix in ixs)
+ data = {tuple(key): block.clone() for key, block in tensor.data.items()}
+ intw = None if tensor.intw is None else {tuple(key): bridge.clone() for key, bridge in tensor.intw.items()}
+ return Tensor(indices=indices, itags=itags, data=data, intw=intw, dtype=tensor.dtype, label=tensor.label)
+
+
+def _tensor_ix(tensor: Tensor, axis: int) -> Ix:
+ idx = tensor.indices[axis]
+ return Ix(tensor.itags[axis], int(idx.dim), idx.direction, idx.sectors, idx.group)
+
+
+def _sector_offsets(ix: Ix) -> dict[object, tuple[int, int]]:
+ offsets: dict[object, tuple[int, int]] = {}
+ cursor = 0
+ for sector in resolved_sectors(ix):
+ offsets[sector.charge] = (cursor, sector.dim)
+ cursor += sector.dim
+ return offsets
+
+
+def _dense_from_tensor_with_ixs(tensor: Tensor, ixs: list[Ix]) -> torch.Tensor:
+ shape = tuple(ix.dim for ix in ixs)
+ device = next(iter(tensor.data.values())).device if tensor.data else torch.device("cpu")
+ dense = torch.zeros(shape, dtype=tensor.dtype, device=device)
+ offsets = [_sector_offsets(ix) for ix in ixs]
+ for key, block in tensor.data.items():
+ slices = tuple(slice(offsets[axis][key[axis]][0], offsets[axis][key[axis]][0] + offsets[axis][key[axis]][1]) for axis in range(len(key)))
+ dense[slices] = block
+ return dense
+
+
+def _left_row_indices_by_flux(link_l: Ix, site_l: Ix) -> dict[object, list[int]]:
+ rows: dict[object, list[int]] = {}
+ link_offsets = _sector_offsets(link_l)
+ site_offsets = _sector_offsets(site_l)
+ dl = int(link_l.dim)
+ d_link = int(link_l.direction)
+ d_site = int(site_l.direction)
+
+ for q_link, (link_start, link_dim) in link_offsets.items():
+ for q_site, (site_start, site_dim) in site_offsets.items():
+ flux = -(d_link * q_link + d_site * q_site)
+ block_rows = rows.setdefault(flux, [])
+ for site_local in range(site_dim):
+ for link_local in range(link_dim):
+ block_rows.append((link_start + link_local) + dl * (site_start + site_local))
+ return rows
+
+
+def _right_col_indices_by_flux(site_r: Ix, link_r: Ix) -> dict[object, list[int]]:
+ cols: dict[object, list[int]] = {}
+ site_offsets = _sector_offsets(site_r)
+ link_offsets = _sector_offsets(link_r)
+ ds = int(site_r.dim)
+ d_site = int(site_r.direction)
+ d_link = int(link_r.direction)
+
+ for q_site, (site_start, site_dim) in site_offsets.items():
+ for q_link, (link_start, link_dim) in link_offsets.items():
+ flux = -(d_site * q_site + d_link * q_link)
+ block_cols = cols.setdefault(flux, [])
+ for link_local in range(link_dim):
+ for site_local in range(site_dim):
+ block_cols.append((site_start + site_local) + ds * (link_start + link_local))
+ return cols
+
+
+def _apply_gate_named(gate, theta, site_l_tag: str, site_r_tag: str):
+ out = einsum("LRlr,aLRb->alrb", gate, theta)
+ out.retag({f"{site_l_tag}*": site_l_tag, f"{site_r_tag}*": site_r_tag})
+ return out
+
+
+def _singular_values_from_diag_tensor(S) -> torch.Tensor:
+ """Extract singular values from a diagonal Nicole tensor.
+
+ Args:
+ S: Diagonal tensor returned by Nicole's SVD.
+
+ Returns:
+ A flat torch tensor containing every block-diagonal entry.
+ """
+ vals = [torch.diagonal(block) for block in S.data.values()]
+ if not vals:
+ return torch.empty((0,), dtype=torch.complex128)
+ return torch.cat(vals)
diff --git a/src/alice/algorithm/bond_update_bug/_kernel/kls/symmetric_completion.py b/src/alice/algorithm/bond_update_bug/_kernel/kls/symmetric_completion.py
new file mode 100644
index 0000000..59d4337
--- /dev/null
+++ b/src/alice/algorithm/bond_update_bug/_kernel/kls/symmetric_completion.py
@@ -0,0 +1,357 @@
+# Copyright (C) 2025-2026 Changkai Zhang.
+#
+# This file is part of Alice project.
+#
+# Alice is free software: you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published
+# by the Free Software Foundation, either version 3 of the License,
+# or (at your option) any later version.
+#
+# Alice is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with Alice. If not, see .
+
+"""Symmetric U(1)-aware augmented isometry functions for BUG/KLS updates."""
+
+from __future__ import annotations
+
+import math
+from typing import Any
+
+import torch
+from nicole import Sector
+
+from ..indices import Ix, fresh_itag, resolved_sectors
+from ..nicole_helpers import make_tensor, reshape_fortran
+from .augment import (
+ _left_row_indices_by_flux,
+ _left_tensor_matrix,
+ _pick_left_update,
+ _pick_right_update,
+ _right_col_indices_by_flux,
+ _right_tensor_matrix,
+ _sector_offsets,
+)
+
+
+def _random_orthonormal_columns(m: int, n: int, dtype, device, generator) -> torch.Tensor:
+ """``m x n`` random complex orthonormal block (``n`` capped at ``m``).
+
+ Used by the missing-quantum-number fill: a reachable charge sector that neither
+ ``U0`` nor ``K1`` populate is opened with a minimal random orthonormal seed so the
+ Galerkin S-step can rotate physical weight into it. The seed is drawn from a
+ fixed-seed generator so a run is reproducible; the S-step and the post-S SVD make
+ the result independent of the seed's orientation (an unpopulated sector is pruned).
+ """
+ n = min(int(n), int(m))
+ if n <= 0:
+ return torch.zeros((m, 0), dtype=dtype, device=device)
+ real = torch.randn((m, n), dtype=torch.float64, device=device, generator=generator)
+ imag = torch.randn((m, n), dtype=torch.float64, device=device, generator=generator)
+ q, _ = torch.linalg.qr(torch.complex(real, imag).to(dtype), mode="reduced")
+ return q[:, :n]
+
+
+def _kl_truncated_left(U0_sub, K1_sub, kl_cutoff: float, max_rank):
+ """Augmented left block ``[U0 | SVD-weight-truncated discarded complement of K1]``.
+
+ Instead of completing ``U0`` to the full local capacity (``d*r``), keep ``U0``
+ EXACTLY and admit only the discarded directions ``(I - U0 U0+) K1`` whose singular
+ value clears ``kl_cutoff`` (relative to the top one) — the per-bond analogue of the
+ global discarded_bug's SVD-truncated complement. This caps the augmented rank
+ between ``r`` and ``d*r`` instead of always ``d*r``. Returns ``(Q, overlap, n_new)``
+ with ``Q = [U0 | Q_new]`` orthonormal (``Q_new`` is orthogonal to ``U0`` by
+ construction). ``overlap`` (the M-hat block) is computed for signature parity but is
+ unused by the discarded variant, which seeds the S-step from ``Û† Θ0 V̂†`` directly.
+ """
+ r = U0_sub.shape[1]
+ eye = U0_sub.conj().transpose(0, 1) @ U0_sub
+ if K1_sub.numel() == 0 or K1_sub.shape[1] == 0:
+ return U0_sub, eye, 0
+ k_perp = K1_sub - U0_sub @ (U0_sub.conj().transpose(0, 1) @ K1_sub)
+ u_k, s_k, _ = torch.linalg.svd(k_perp, full_matrices=False)
+ if s_k.numel() == 0 or float(s_k[0]) == 0.0:
+ return U0_sub, eye, 0
+ keep = int((s_k > kl_cutoff * float(s_k[0])).sum().item())
+ if max_rank is not math.inf:
+ keep = min(keep, max(0, int(max_rank) - r))
+ if keep <= 0:
+ return U0_sub, eye, 0
+ Q = torch.cat([U0_sub, u_k[:, :keep]], dim=1)
+ overlap = Q.conj().transpose(0, 1) @ U0_sub
+ return Q, overlap, keep
+
+
+def _kl_truncated_right(V0_sub, L1_sub, kl_cutoff: float, max_rank):
+ """Augmented right block ``[V0 ; SVD-weight-truncated discarded complement of L1]`` (row-wise mirror)."""
+ r = V0_sub.shape[0]
+ eye = V0_sub @ V0_sub.conj().transpose(0, 1)
+ if L1_sub.numel() == 0 or L1_sub.shape[0] == 0:
+ return V0_sub, eye, 0
+ l_perp = L1_sub - (L1_sub @ V0_sub.conj().transpose(0, 1)) @ V0_sub
+ _, s_l, vh_l = torch.linalg.svd(l_perp, full_matrices=False)
+ if s_l.numel() == 0 or float(s_l[0]) == 0.0:
+ return V0_sub, eye, 0
+ keep = int((s_l > kl_cutoff * float(s_l[0])).sum().item())
+ if max_rank is not math.inf:
+ keep = min(keep, max(0, int(max_rank) - r))
+ if keep <= 0:
+ return V0_sub, eye, 0
+ B = torch.cat([V0_sub, vh_l[:keep, :]], dim=0)
+ overlap = V0_sub @ B.conj().transpose(0, 1)
+ return B, overlap, keep
+
+
+def _symmetric_augmented_left_isometry_from_k(
+ U0_tens,
+ K1_tens,
+ *args,
+ link_l=None,
+ site_l=None,
+ old_mid_u=None,
+ k_mid=None,
+ augment: bool = True,
+ max_rank: int | float = math.inf,
+ aug_tol: float = 1e-12,
+ aug_missing_fill: int = 1,
+ kl_cutoff: float | None = None,
+ **kwargs: Any,
+):
+ if args:
+ if len(args) != 4:
+ raise TypeError(
+ "_symmetric_augmented_left_isometry_from_k expects (link_l, site_l, old_mid_u, k_mid) after the tensors."
+ )
+ if any(name in ("link_l", "site_l", "old_mid_u", "k_mid") for name in kwargs):
+ raise TypeError("Provide symmetric left-augmentation indices either positionally or by keyword, not both.")
+ link_l, site_l, old_mid_u, k_mid = args
+ if link_l is None or site_l is None or old_mid_u is None or k_mid is None:
+ try:
+ link_l = kwargs.pop("link_l") if link_l is None else link_l
+ site_l = kwargs.pop("site_l") if site_l is None else site_l
+ old_mid_u = kwargs.pop("old_mid_u") if old_mid_u is None else old_mid_u
+ k_mid = kwargs.pop("k_mid") if k_mid is None else k_mid
+ except KeyError as exc:
+ raise TypeError("Missing symmetric left-augmentation index input.") from exc
+ if kwargs:
+ unknown = ", ".join(sorted(kwargs))
+ raise TypeError(f"Unknown symmetric left-augmentation option(s): {unknown}")
+
+ dtype = torch.complex128
+ device = next(iter(U0_tens.data.values())).device if U0_tens.data else torch.device("cpu")
+ U0_mat = _left_tensor_matrix(U0_tens, link_l, site_l, old_mid_u)
+ K1_mat = _left_tensor_matrix(K1_tens, link_l, site_l, k_mid)
+ row_blocks = _left_row_indices_by_flux(link_l, site_l)
+ u_offsets = _sector_offsets(old_mid_u)
+ k_offsets = _sector_offsets(k_mid)
+ d_old = int(old_mid_u.direction)
+ d_k = int(k_mid.direction)
+
+ fluxes: list[object] = []
+ for sec in resolved_sectors(old_mid_u):
+ flux = d_old * sec.charge
+ if flux not in fluxes:
+ fluxes.append(flux)
+ for sec in resolved_sectors(k_mid):
+ flux = d_k * sec.charge
+ if flux not in fluxes:
+ fluxes.append(flux)
+ for flux in row_blocks:
+ if flux not in fluxes:
+ fluxes.append(flux)
+
+ rng = torch.Generator(device=device)
+ rng.manual_seed(0x5EED)
+ pieces: list[tuple[object, list[int], torch.Tensor, torch.Tensor]] = []
+ total_dim = 0
+ n_new_total = 0
+ for flux in fluxes:
+ rows = row_blocks.get(flux, [])
+ if not rows:
+ continue
+
+ old_charge = flux // d_old
+ k_charge = flux // d_k
+ old_slice = u_offsets.get(old_charge)
+ k_slice = k_offsets.get(k_charge)
+ U0_sub = U0_mat[rows, old_slice[0] : old_slice[0] + old_slice[1]] if old_slice else torch.zeros((len(rows), 0), dtype=dtype, device=device)
+ K1_sub = K1_mat[rows, k_slice[0] : k_slice[0] + k_slice[1]] if k_slice else torch.zeros((len(rows), 0), dtype=dtype, device=device)
+ if kl_cutoff is not None and augment:
+ # Efficient path: keep U0 exact, admit only the SVD-weight-significant
+ # discarded directions of K1 (no full d*r completion).
+ Q_block, overlap_block, n_new = _kl_truncated_left(U0_sub, K1_sub, kl_cutoff, max_rank)
+ else:
+ # Sulz augmented basis of this sector: orth([U0 | K1]) at rank <= 2r (no
+ # complete_column_basis padding to the full d*r local dimension).
+ Q_block, overlap_block, n_new = _pick_left_update(U0_sub, K1_sub, augment=augment, max_rank=max_rank, aug_tol=aug_tol)
+ if augment and Q_block.shape[1] == 0 and len(rows) > 0:
+ # Missing-quantum-number fill (replaces complete_column_basis). This
+ # locally reachable charge sector is populated by neither U0 nor K1 --
+ # U(1) charge conservation with the frozen right frame keeps K1 in U0's
+ # sector, so orth([U0|K1]) can never OPEN a new sector the way the dense
+ # Sulz BUG's K1 does. Seed a minimal random orthonormal block; the S-step
+ # rotates physical weight into it and the post-S SVD prunes it if the
+ # dynamics leaves it empty. (complete_column_basis instead filled the
+ # sector to its FULL local dim -> the d*r augmented-rank blow-up.)
+ n_seed = min(len(rows), max(1, int(aug_missing_fill)))
+ Q_block = _random_orthonormal_columns(len(rows), n_seed, dtype, device, rng)
+ overlap_block = Q_block.conj().transpose(0, 1) @ U0_sub
+ n_new = int(Q_block.shape[1])
+ if Q_block.shape[1] == 0:
+ continue
+ pieces.append((old_charge, rows, Q_block, overlap_block))
+ total_dim += int(Q_block.shape[1])
+ n_new_total += int(n_new)
+
+ if total_dim == 0:
+ raise ValueError("Left symmetric K-step augmentation produced zero rank.")
+
+ sectors = tuple(Sector(int(charge), int(block.shape[1])) for charge, _, block, _ in pieces)
+ new_mid = Ix(fresh_itag(old_mid_u.itag), total_dim, old_mid_u.direction, sectors, old_mid_u.group)
+ U_aug_mat = torch.zeros((link_l.dim * site_l.dim, total_dim), dtype=dtype, device=device)
+ M_hat_mat = torch.zeros((total_dim, old_mid_u.dim), dtype=dtype, device=device)
+
+ cursor = 0
+ for charge, rows, block, overlap_block in pieces:
+ width = int(block.shape[1])
+ sl = slice(cursor, cursor + width)
+ U_aug_mat[rows, sl] = block
+ old_slice = u_offsets.get(charge)
+ if old_slice is not None and overlap_block.numel():
+ old_sl = slice(old_slice[0], old_slice[0] + old_slice[1])
+ M_hat_mat[sl, old_sl] = overlap_block
+ cursor += width
+
+ U_aug_tens = make_tensor(
+ reshape_fortran(U_aug_mat, (link_l.dim, site_l.dim, new_mid.dim)),
+ [link_l, site_l, new_mid],
+ dtype=dtype,
+ )
+ M_hat_tens = make_tensor(M_hat_mat, [new_mid.reversed(), old_mid_u], dtype=dtype)
+ return U_aug_tens, M_hat_tens, n_new_total
+
+
+def _symmetric_augmented_right_isometry_from_l(
+ V0_tens,
+ L1_tens,
+ *args,
+ old_mid_v=None,
+ l_mid=None,
+ site_r=None,
+ link_r=None,
+ augment: bool = True,
+ max_rank: int | float = math.inf,
+ aug_tol: float = 1e-12,
+ aug_missing_fill: int = 1,
+ kl_cutoff: float | None = None,
+ **kwargs: Any,
+):
+ if args:
+ if len(args) != 4:
+ raise TypeError(
+ "_symmetric_augmented_right_isometry_from_l expects (old_mid_v, l_mid, site_r, link_r) after the tensors."
+ )
+ if any(name in ("old_mid_v", "l_mid", "site_r", "link_r") for name in kwargs):
+ raise TypeError("Provide symmetric right-augmentation indices either positionally or by keyword, not both.")
+ old_mid_v, l_mid, site_r, link_r = args
+ if old_mid_v is None or l_mid is None or site_r is None or link_r is None:
+ try:
+ old_mid_v = kwargs.pop("old_mid_v") if old_mid_v is None else old_mid_v
+ l_mid = kwargs.pop("l_mid") if l_mid is None else l_mid
+ site_r = kwargs.pop("site_r") if site_r is None else site_r
+ link_r = kwargs.pop("link_r") if link_r is None else link_r
+ except KeyError as exc:
+ raise TypeError("Missing symmetric right-augmentation index input.") from exc
+ if kwargs:
+ unknown = ", ".join(sorted(kwargs))
+ raise TypeError(f"Unknown symmetric right-augmentation option(s): {unknown}")
+
+ dtype = torch.complex128
+ device = next(iter(V0_tens.data.values())).device if V0_tens.data else torch.device("cpu")
+ V0_mat = _right_tensor_matrix(V0_tens, old_mid_v, site_r, link_r)
+ L1_mat = _right_tensor_matrix(L1_tens, l_mid, site_r, link_r)
+ col_blocks = _right_col_indices_by_flux(site_r, link_r)
+ v_offsets = _sector_offsets(old_mid_v)
+ l_offsets = _sector_offsets(l_mid)
+ d_old = int(old_mid_v.direction)
+ d_l = int(l_mid.direction)
+
+ fluxes: list[object] = []
+ for sec in resolved_sectors(old_mid_v):
+ flux = d_old * sec.charge
+ if flux not in fluxes:
+ fluxes.append(flux)
+ for sec in resolved_sectors(l_mid):
+ flux = d_l * sec.charge
+ if flux not in fluxes:
+ fluxes.append(flux)
+ for flux in col_blocks:
+ if flux not in fluxes:
+ fluxes.append(flux)
+
+ rng = torch.Generator(device=device)
+ rng.manual_seed(0x5EED)
+ pieces: list[tuple[object, list[int], torch.Tensor, torch.Tensor]] = []
+ total_dim = 0
+ n_new_total = 0
+ for flux in fluxes:
+ cols = col_blocks.get(flux, [])
+ if not cols:
+ continue
+
+ old_charge = flux // d_old
+ l_charge = flux // d_l
+ old_slice = v_offsets.get(old_charge)
+ l_slice = l_offsets.get(l_charge)
+ V0_sub = V0_mat[old_slice[0] : old_slice[0] + old_slice[1], cols] if old_slice else torch.zeros((0, len(cols)), dtype=dtype, device=device)
+ L1_sub = L1_mat[l_slice[0] : l_slice[0] + l_slice[1], cols] if l_slice else torch.zeros((0, len(cols)), dtype=dtype, device=device)
+ if kl_cutoff is not None and augment:
+ B_block, overlap_block, n_new = _kl_truncated_right(V0_sub, L1_sub, kl_cutoff, max_rank)
+ else:
+ # Sulz augmented basis of this sector: orth([V0 ; L1]) at rank <= 2r (no
+ # complete_row_basis padding to the full d*r local dimension).
+ B_block, overlap_block, n_new = _pick_right_update(V0_sub, L1_sub, augment=augment, max_rank=max_rank, aug_tol=aug_tol)
+ if augment and B_block.shape[0] == 0 and len(cols) > 0:
+ # Missing-quantum-number fill (row mirror of the K-step; replaces
+ # complete_row_basis): seed a minimal random row-orthonormal block so the
+ # S-step can open this reachable-but-empty charge sector.
+ n_seed = min(len(cols), max(1, int(aug_missing_fill)))
+ B_block = _random_orthonormal_columns(len(cols), n_seed, dtype, device, rng).transpose(0, 1)
+ overlap_block = V0_sub @ B_block.conj().transpose(0, 1)
+ n_new = int(B_block.shape[0])
+ if B_block.shape[0] == 0:
+ continue
+ pieces.append((old_charge, cols, B_block, overlap_block))
+ total_dim += int(B_block.shape[0])
+ n_new_total += int(n_new)
+
+ if total_dim == 0:
+ raise ValueError("Right symmetric L-step augmentation produced zero rank.")
+
+ sectors = tuple(Sector(int(charge), int(block.shape[0])) for charge, _, block, _ in pieces)
+ new_mid = Ix(fresh_itag(old_mid_v.itag), total_dim, old_mid_v.direction, sectors, old_mid_v.group)
+ V_aug_mat = torch.zeros((total_dim, site_r.dim * link_r.dim), dtype=dtype, device=device)
+ N_hat_mat = torch.zeros((old_mid_v.dim, total_dim), dtype=dtype, device=device)
+
+ cursor = 0
+ for charge, cols, block, overlap_block in pieces:
+ height = int(block.shape[0])
+ sl = slice(cursor, cursor + height)
+ V_aug_mat[sl, cols] = block
+ old_slice = v_offsets.get(charge)
+ if old_slice is not None and overlap_block.numel():
+ old_sl = slice(old_slice[0], old_slice[0] + old_slice[1])
+ N_hat_mat[old_sl, sl] = overlap_block
+ cursor += height
+
+ V_aug_tens = make_tensor(
+ reshape_fortran(V_aug_mat, (new_mid.dim, site_r.dim, link_r.dim)),
+ [new_mid, site_r, link_r],
+ dtype=dtype,
+ )
+ N_hat_tens = make_tensor(N_hat_mat, [old_mid_v, new_mid.reversed()], dtype=dtype)
+ return V_aug_tens, N_hat_tens, n_new_total
diff --git a/src/alice/algorithm/bond_update_bug/_kernel/krylov.py b/src/alice/algorithm/bond_update_bug/_kernel/krylov.py
new file mode 100644
index 0000000..55e69e5
--- /dev/null
+++ b/src/alice/algorithm/bond_update_bug/_kernel/krylov.py
@@ -0,0 +1,579 @@
+# Copyright (C) 2025-2026 Changkai Zhang.
+#
+# This file is part of Alice project.
+#
+# Alice is free software: you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published
+# by the Free Software Foundation, either version 3 of the License,
+# or (at your option) any later version.
+#
+# Alice is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with Alice. If not, see .
+
+"""Krylov and micro-step helpers for dense vectors and Nicole tensors.
+
+The original Julia code threaded method names, tolerances, iteration caps, and
+backend choices through many low-level calls. This module keeps the numerical
+behavior but packages the runtime controls into small option objects so the
+higher-level BUG code can call it in a more readable way.
+"""
+
+from __future__ import annotations
+
+import contextlib
+from collections.abc import Callable
+from dataclasses import dataclass, replace
+from typing import Any, Iterable
+
+import torch
+from nicole import Tensor, conj as _nconj, einsum as _neinsum
+
+from .nicole_helpers import flatten_fortran, to_dense
+
+BUG_DEFAULT_EXPV_BACKEND = "krylovkit"
+BUG_ALLOWED_EXPV_BACKENDS = ("krylovkit", "native_hermitian_lanczos")
+_ACTIVE_BACKEND = [BUG_DEFAULT_EXPV_BACKEND]
+_ACTIVE_PREFACTOR = [complex(0.0, -1.0)]
+
+
+@dataclass(frozen=True)
+class LanczosOptions:
+ """Options for Hermitian Lanczos exponentiation.
+
+ Args:
+ tol: Termination tolerance for the Lanczos recurrence.
+ krylovdim: Maximum Krylov basis size.
+ """
+
+ tol: float = 1e-13
+ krylovdim: int = 60
+
+
+@dataclass(frozen=True)
+class LinearSubstepOptions:
+ """Options for one dense or matrix-free linear micro-step.
+
+ Args:
+ method: Local integrator name. Supported values are ``"expv"``,
+ ``"euler"``, and ``"rk4"``.
+ lanczos_tol: Tolerance passed to Lanczos-based ``expv`` routines.
+ lanczos_maxiter: Maximum Lanczos basis size.
+ restart: Reserved compatibility flag retained from the Julia API.
+ issymmetric: Optional hint used by :func:`general_linear_substep`.
+ """
+
+ method: str = "expv"
+ lanczos_tol: float = 1e-13
+ lanczos_maxiter: int = 60
+ restart: int = 1
+ issymmetric: bool = True
+
+ def lanczos_options(self) -> LanczosOptions:
+ """Return the matching :class:`LanczosOptions` view.
+
+ Returns:
+ A :class:`LanczosOptions` instance derived from the substep fields.
+ """
+ return LanczosOptions(tol=self.lanczos_tol, krylovdim=self.lanczos_maxiter)
+
+
+_LANCZOS_OPTION_FIELDS = {field.name for field in LanczosOptions.__dataclass_fields__.values()}
+_SUBSTEP_OPTION_FIELDS = {field.name for field in LinearSubstepOptions.__dataclass_fields__.values()}
+_LANCZOS_ALIASES = {"maxiter": "krylovdim"}
+_SUBSTEP_ALIASES = {"tol": "lanczos_tol", "maxiter": "lanczos_maxiter", "krylovdim": "lanczos_maxiter"}
+
+__all__ = [
+ "BUG_ALLOWED_EXPV_BACKENDS",
+ "BUG_DEFAULT_EXPV_BACKEND",
+ "LanczosOptions",
+ "LinearSubstepOptions",
+ "active_expv_backend",
+ "active_time_prefactor",
+ "complex_tensor_array",
+ "complex_tensor_vec",
+ "general_linear_substep",
+ "hermitian_tridiagonal_exp_coeffs",
+ "linear_substep",
+ "native_hermitian_lanczos_exponentiate",
+ "tensor_inner",
+ "tensor_lanczos_expv",
+ "with_expv_backend",
+ "with_time_prefactor",
+]
+
+
+def _coerce_lanczos_options(options: LanczosOptions | None = None, **kwargs: Any) -> LanczosOptions:
+ """Normalize Lanczos options from an object or legacy kwargs.
+
+ Args:
+ options: Existing :class:`LanczosOptions` instance.
+ **kwargs: Field overrides or legacy aliases such as ``maxiter``.
+
+ Returns:
+ A normalized :class:`LanczosOptions` instance.
+ """
+ overrides: dict[str, Any] = {}
+ for name, value in kwargs.items():
+ normalized = _LANCZOS_ALIASES.get(name, name)
+ if normalized not in _LANCZOS_OPTION_FIELDS:
+ raise TypeError(f"Unknown Lanczos option: {name}")
+ overrides[normalized] = value
+ if options is None:
+ return LanczosOptions(**overrides)
+ return replace(options, **overrides)
+
+
+def _coerce_substep_options(options: LinearSubstepOptions | None = None, **kwargs: Any) -> LinearSubstepOptions:
+ """Normalize linear-substep options from an object or legacy kwargs.
+
+ Args:
+ options: Existing :class:`LinearSubstepOptions` instance.
+ **kwargs: Field overrides or legacy aliases such as ``tol``.
+
+ Returns:
+ A normalized :class:`LinearSubstepOptions` instance.
+ """
+ overrides: dict[str, Any] = {}
+ for name, value in kwargs.items():
+ normalized = _SUBSTEP_ALIASES.get(name, name)
+ if normalized not in _SUBSTEP_OPTION_FIELDS:
+ raise TypeError(f"Unknown linear-substep option: {name}")
+ overrides[normalized] = value
+ if options is None:
+ return LinearSubstepOptions(**overrides)
+ return replace(options, **overrides)
+
+
+def _as_complex_tensor(x) -> torch.Tensor:
+ """Convert input data to a complex128 torch tensor.
+
+ Args:
+ x: Tensor-like object.
+
+ Returns:
+ A ``torch.complex128`` tensor.
+ """
+ if isinstance(x, torch.Tensor):
+ return x.to(dtype=torch.complex128)
+ return torch.as_tensor(x, dtype=torch.complex128)
+
+
+def hermitian_tridiagonal_exp_coeffs(
+ alpha: Iterable[float] | torch.Tensor,
+ beta: Iterable[float] | torch.Tensor,
+ dt: complex,
+) -> torch.Tensor:
+ """Compute the Krylov coefficients for ``exp(dt*T)e1``.
+
+ Args:
+ alpha: Diagonal entries of the Hermitian tridiagonal matrix.
+ beta: Off-diagonal entries.
+ dt: Scalar prefactor used in the exponential.
+
+ Returns:
+ The coefficient vector in the Lanczos basis.
+ """
+ alpha_t = torch.as_tensor(tuple(alpha) if not isinstance(alpha, torch.Tensor) else alpha, dtype=torch.float64)
+ beta_t = torch.as_tensor(tuple(beta) if not isinstance(beta, torch.Tensor) else beta, dtype=torch.float64)
+ if alpha_t.numel() == 0:
+ return torch.empty((0,), dtype=torch.complex128)
+
+ tridiagonal = torch.diag(alpha_t)
+ if beta_t.numel() > 0:
+ tridiagonal = tridiagonal + torch.diag(beta_t, diagonal=1) + torch.diag(beta_t, diagonal=-1)
+ evals, evecs = torch.linalg.eigh(tridiagonal)
+ evecs_c = evecs.to(torch.complex128)
+ weights = torch.exp(dt * evals.to(torch.complex128)) * evecs_c[0, :]
+ return evecs_c @ weights
+
+
+def native_hermitian_lanczos_exponentiate(
+ matvec: Callable[[torch.Tensor], torch.Tensor | object],
+ dt: complex,
+ x,
+ *,
+ options: LanczosOptions | None = None,
+ **kwargs: Any,
+) -> tuple[torch.Tensor, int]:
+ """Apply ``exp(dt * H)`` to ``x`` using a native Hermitian Lanczos solve.
+
+ Args:
+ matvec: Matrix-free Hermitian action on dense vectors.
+ dt: Scalar prefactor used in the exponential.
+ x: Input vector.
+ options: Optional :class:`LanczosOptions` instance.
+ **kwargs: Legacy option overrides such as ``tol=...``.
+
+ Returns:
+ A pair ``(y, numops)`` containing the evolved vector and the number of
+ matrix-vector products performed.
+ """
+ options = _coerce_lanczos_options(options, **kwargs)
+ x_work = _as_complex_tensor(x).reshape(-1).clone()
+ n = int(x_work.numel())
+ if n == 0:
+ return torch.empty((0,), dtype=torch.complex128), 0
+
+ norm_x = torch.linalg.norm(x_work)
+ if norm_x == 0:
+ return torch.zeros_like(x_work), 0
+
+ mmax = min(max(int(options.krylovdim), 1), n)
+ basis = torch.empty((n, mmax), dtype=torch.complex128)
+ alpha = torch.empty((mmax,), dtype=torch.float64)
+ beta = torch.empty((max(mmax - 1, 0),), dtype=torch.float64)
+
+ basis[:, 0] = x_work / norm_x
+ numops = 0
+ final_dim = 1
+
+ # Standard Hermitian Lanczos recurrence on dense vectors.
+ for j in range(mmax):
+ vj = basis[:, j]
+ work = _as_complex_tensor(matvec(vj)).reshape(-1)
+ numops += 1
+
+ if j > 0:
+ work = work - beta[j - 1] * basis[:, j - 1]
+
+ alpha[j] = torch.real(torch.vdot(vj, work))
+ work = work - alpha[j] * vj
+
+ if j == mmax - 1:
+ final_dim = j + 1
+ break
+
+ beta_j = torch.linalg.norm(work)
+ if float(beta_j) <= options.tol:
+ final_dim = j + 1
+ break
+
+ beta[j] = beta_j.real
+ basis[:, j + 1] = work / beta_j
+ final_dim = j + 2
+
+ coeff = hermitian_tridiagonal_exp_coeffs(alpha[:final_dim], beta[: max(final_dim - 1, 0)], dt)
+ y = norm_x.to(torch.complex128) * (basis[:, :final_dim] @ coeff)
+ return y, numops
+
+
+@contextlib.contextmanager
+def with_expv_backend(backend: str):
+ """Temporarily set the active expv backend name.
+
+ Args:
+ backend: Backend label from :data:`BUG_ALLOWED_EXPV_BACKENDS`.
+
+ Returns:
+ A context manager that restores the previous backend on exit.
+ """
+ if backend not in BUG_ALLOWED_EXPV_BACKENDS:
+ raise ValueError(f"Unknown expv backend: {backend}")
+ previous = _ACTIVE_BACKEND[0]
+ _ACTIVE_BACKEND[0] = backend
+ try:
+ yield
+ finally:
+ _ACTIVE_BACKEND[0] = previous
+
+
+def active_expv_backend() -> str:
+ """Return the currently active expv backend label.
+
+ Returns:
+ The active backend name.
+ """
+ return _ACTIVE_BACKEND[0]
+
+
+@contextlib.contextmanager
+def with_time_prefactor(c: complex):
+ """Temporarily override the global evolution prefactor.
+
+ Args:
+ c: New complex prefactor.
+
+ Returns:
+ A context manager that restores the previous prefactor on exit.
+ """
+ previous = _ACTIVE_PREFACTOR[0]
+ _ACTIVE_PREFACTOR[0] = complex(c)
+ try:
+ yield
+ finally:
+ _ACTIVE_PREFACTOR[0] = previous
+
+
+def active_time_prefactor() -> complex:
+ """Return the currently active evolution prefactor.
+
+ Returns:
+ The active complex prefactor.
+ """
+ return _ACTIVE_PREFACTOR[0]
+
+
+def _matrix_linear_substep(
+ H,
+ dt: complex,
+ x,
+ *,
+ options: LinearSubstepOptions,
+) -> tuple[torch.Tensor, int]:
+ """Apply one micro-step when the operator is available as a dense matrix.
+
+ Args:
+ H: Dense matrix.
+ dt: Step size or exponential prefactor.
+ x: Input vector.
+ options: Linear-substep configuration.
+
+ Returns:
+ A pair ``(y, numops)`` describing the updated vector and the number of
+ explicit matvecs counted for Krylov methods.
+ """
+ H_dense = _as_complex_tensor(H)
+ x_dense = _as_complex_tensor(x).reshape(-1)
+ if options.method == "expv":
+ if active_expv_backend() == "native_hermitian_lanczos":
+ return native_hermitian_lanczos_exponentiate(
+ lambda v: H_dense @ v,
+ dt,
+ x_dense,
+ options=options.lanczos_options(),
+ )
+ return torch.linalg.matrix_exp(dt * H_dense) @ x_dense, 0
+
+ return linear_substep(
+ lambda v: H_dense @ v,
+ dt,
+ x_dense,
+ options=options,
+ )
+
+
+def linear_substep(
+ H_or_matvec,
+ dt: complex,
+ x,
+ *,
+ options: LinearSubstepOptions | None = None,
+ **kwargs: Any,
+) -> tuple[torch.Tensor, int]:
+ """Advance one dense or matrix-free micro-step.
+
+ Args:
+ H_or_matvec: Dense matrix or matrix-free action.
+ dt: Step size or exponential prefactor.
+ x: Input vector.
+ options: Optional :class:`LinearSubstepOptions` instance.
+ **kwargs: Legacy option overrides such as ``method="expv"``.
+
+ Returns:
+ A pair ``(y, numops)`` containing the updated vector and the counted
+ operator applications.
+ """
+ options = _coerce_substep_options(options, **kwargs)
+ if callable(H_or_matvec):
+ x_vec = _as_complex_tensor(x).reshape(-1)
+ matvec = H_or_matvec
+
+ if options.method == "expv":
+ return native_hermitian_lanczos_exponentiate(matvec, dt, x_vec, options=options.lanczos_options())
+
+ if options.method == "euler":
+ return x_vec + dt * _as_complex_tensor(matvec(x_vec)).reshape(-1), 1
+
+ if options.method == "rk4":
+ # Keep the explicit stages readable; the dimensions are tiny compared
+ # with the conceptual cost of understanding hidden helper machinery.
+ k1 = _as_complex_tensor(matvec(x_vec)).reshape(-1)
+ k2 = _as_complex_tensor(matvec(x_vec + (dt / 2) * k1)).reshape(-1)
+ k3 = _as_complex_tensor(matvec(x_vec + (dt / 2) * k2)).reshape(-1)
+ k4 = _as_complex_tensor(matvec(x_vec + dt * k3)).reshape(-1)
+ return x_vec + (dt / 6) * (k1 + 2 * k2 + 2 * k3 + k4), 4
+
+ raise ValueError(f"Unknown substep method: {options.method}. Supported: expv, euler, rk4.")
+
+ return _matrix_linear_substep(H_or_matvec, dt, x, options=options)
+
+
+def general_linear_substep(
+ matvec,
+ dt: complex,
+ x,
+ *,
+ options: LinearSubstepOptions | None = None,
+ **kwargs: Any,
+) -> tuple[torch.Tensor, int]:
+ """Variant of :func:`linear_substep` with explicit symmetry dispatch.
+
+ Args:
+ matvec: Matrix-free operator action.
+ dt: Step size or exponential prefactor.
+ x: Input vector.
+ options: Optional :class:`LinearSubstepOptions` instance.
+ **kwargs: Legacy option overrides such as ``issymmetric=False``.
+
+ Returns:
+ A pair ``(y, numops)`` containing the updated vector and the counted
+ operator applications.
+ """
+ options = _coerce_substep_options(options, **kwargs)
+ x_vec = _as_complex_tensor(x).reshape(-1)
+ if options.method != "expv":
+ return linear_substep(matvec, dt, x_vec, options=options)
+
+ if options.issymmetric:
+ return native_hermitian_lanczos_exponentiate(
+ matvec,
+ dt,
+ x_vec,
+ options=LanczosOptions(
+ tol=options.lanczos_tol,
+ krylovdim=min(len(x_vec), max(options.lanczos_maxiter, 4)),
+ ),
+ )
+
+ # The nonsymmetric fallback is intentionally dense because the current port
+ # only needs it for very small diagnostic problems.
+ n = len(x_vec)
+ eye = torch.eye(n, dtype=torch.complex128)
+ dense = torch.empty((n, n), dtype=torch.complex128)
+ for col in range(n):
+ dense[:, col] = _as_complex_tensor(matvec(eye[:, col])).reshape(-1)
+ return torch.linalg.matrix_exp(dt * dense) @ x_vec, n
+
+
+def tensor_inner(a: Tensor, b: Tensor) -> complex:
+ """Return the canonical inner product ```` for same-shape tensors.
+
+ Args:
+ a: Left tensor.
+ b: Right tensor.
+
+ Returns:
+ The complex scalar inner product.
+ """
+ equation = "".join(chr(97 + axis) for axis in range(len(a.itags)))
+ result = _neinsum(f"{equation},{equation}->", _nconj(a), b)
+ # A structurally-zero inner product (no matching charge blocks) is a scalar
+ # tensor with no block; nicole's .item() rejects that. It IS exactly zero --
+ # e.g. for a purely off-diagonal (pure XX flip-flop) H on an Sz-basis
+ # product state. Return 0 rather than requiring a diagonal regularizer.
+ if not getattr(result, "data", None):
+ return 0.0 + 0.0j
+ return result.item()
+
+
+# Opt-in Krylov-depth instrumentation (off by default => zero overhead). When
+# enabled, every tensor_lanczos_expv call appends its Krylov dimension (number of
+# matrix-free H applications) to KRYLOV_LOG, for the N_Krylov diagnostic.
+KRYLOV_LOG: list[int] = []
+_KRYLOV_RECORD = False
+
+
+def enable_krylov_log() -> None:
+ global _KRYLOV_RECORD
+ _KRYLOV_RECORD = True
+ KRYLOV_LOG.clear()
+
+
+def disable_krylov_log() -> None:
+ global _KRYLOV_RECORD
+ _KRYLOV_RECORD = False
+
+
+def get_krylov_log() -> list[int]:
+ return list(KRYLOV_LOG)
+
+
+def tensor_lanczos_expv(
+ apply: Callable[[Tensor], Tensor],
+ dt: complex,
+ x: Tensor,
+ *,
+ options: LanczosOptions | None = None,
+ **kwargs: Any,
+) -> Tensor:
+ """Return ``exp(dt * H) @ x`` for Hermitian Nicole tensor actions.
+
+ Args:
+ apply: Matrix-free Hermitian action on Nicole tensors.
+ dt: Scalar prefactor used in the exponential.
+ x: Input Nicole tensor.
+ options: Optional :class:`LanczosOptions` instance.
+ **kwargs: Legacy option overrides such as ``maxiter=60``.
+
+ Returns:
+ The evolved Nicole tensor.
+ """
+ options = _coerce_lanczos_options(options, **kwargs)
+ beta0 = x.norm()
+ if beta0 == 0:
+ if _KRYLOV_RECORD:
+ KRYLOV_LOG.append(0)
+ return x
+
+ v = (1.0 / beta0) * x
+ basis = [v]
+ alpha: list[float] = []
+ betas: list[float] = []
+
+ w = apply(v)
+ a = tensor_inner(v, w).real
+ alpha.append(a)
+ w = w + (-a) * v
+
+ # This is the same Hermitian recurrence as the dense version, but each basis
+ # vector is now a Nicole tensor instead of a flat torch vector.
+ for _ in range(1, options.krylovdim):
+ b = w.norm()
+ if float(b) < options.tol:
+ break
+ betas.append(float(b))
+ v = (1.0 / b) * w
+ basis.append(v)
+ w = apply(v)
+ a = tensor_inner(v, w).real
+ alpha.append(a)
+ w = w + (-a) * v + (-b) * basis[-2]
+
+ if _KRYLOV_RECORD:
+ KRYLOV_LOG.append(len(alpha))
+ coeff = hermitian_tridiagonal_exp_coeffs(alpha, betas, dt) * beta0
+ out = coeff[0] * basis[0]
+ for idx in range(1, len(alpha)):
+ out = out + coeff[idx] * basis[idx]
+ return out
+
+
+def complex_tensor_array(tensor: Tensor, itag_order):
+ """Convert a Nicole tensor to a dense complex128 torch array.
+
+ Args:
+ tensor: Nicole tensor to densify.
+ itag_order: Tag order passed to :func:`bug_nicole.nicole_helpers.to_dense`.
+
+ Returns:
+ A dense ``torch.complex128`` tensor.
+ """
+ return to_dense(tensor, itag_order).to(dtype=torch.complex128)
+
+
+def complex_tensor_vec(tensor: Tensor, itag_order):
+ """Flatten :func:`complex_tensor_array` in Fortran/column-major order.
+
+ Args:
+ tensor: Nicole tensor to flatten.
+ itag_order: Tag order used for densification.
+
+ Returns:
+ A one-dimensional dense vector.
+ """
+ return flatten_fortran(complex_tensor_array(tensor, itag_order))
diff --git a/src/alice/algorithm/bond_update_bug/_kernel/linalg.py b/src/alice/algorithm/bond_update_bug/_kernel/linalg.py
new file mode 100644
index 0000000..81717b3
--- /dev/null
+++ b/src/alice/algorithm/bond_update_bug/_kernel/linalg.py
@@ -0,0 +1,405 @@
+# Copyright (C) 2025-2026 Changkai Zhang.
+#
+# This file is part of Alice project.
+#
+# Alice is free software: you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published
+# by the Free Software Foundation, either version 3 of the License,
+# or (at your option) any later version.
+#
+# Alice is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with Alice. If not, see .
+
+"""Linear-algebra helpers for Nicole tensors and their dense projections.
+
+The higher-level BUG code talks in terms of QR, LQ, SVD, and orthonormal basis
+completion. This module wraps Nicole's decompositions with small Python helpers
+that preserve the richer :class:`bug_nicole.indices.Ix` metadata and expose a
+friendlier, more explicit surface to the rest of the package.
+"""
+
+from __future__ import annotations
+
+import math
+from dataclasses import dataclass, replace
+from typing import Any, Sequence
+
+import torch
+from nicole import Tensor
+from nicole import decomp as _ndecomp
+
+from .indices import Ix, fresh_itag
+from .nicole_helpers import tcontract
+
+__all__ = [
+ "SVDOptions",
+ "identity_overlap_matrix",
+ "lq",
+ "qr",
+ "qr_column_basis",
+ "qr_nonzero_diagonal_rank",
+ "qr_row_basis",
+ "random_unitary",
+ "reconstruct_from_svd",
+ "svd",
+ "truncate",
+]
+
+
+@dataclass(frozen=True)
+class SVDOptions:
+ """Options controlling the tensor SVD wrapper.
+
+ Args:
+ maxdim: Maximum kept bond dimension. ``math.inf`` means no explicit cap.
+ cutoff: Relative singular-value cutoff.
+ tag: Base tag used for the newly created bond indices.
+ """
+
+ maxdim: int | float = math.inf
+ cutoff: float = 0.0
+ tag: str = "b"
+
+
+_SVD_OPTION_FIELDS = {field.name for field in SVDOptions.__dataclass_fields__.values()}
+
+
+def _coerce_svd_options(options: SVDOptions | None = None, **kwargs: Any) -> SVDOptions:
+ """Normalize SVD options from an object or legacy keyword arguments.
+
+ Args:
+ options: Existing :class:`SVDOptions` instance.
+ **kwargs: Field overrides such as ``maxdim=64``.
+
+ Returns:
+ A normalized :class:`SVDOptions` instance.
+ """
+ overrides = {name: kwargs.pop(name) for name in list(kwargs) if name in _SVD_OPTION_FIELDS}
+ if kwargs:
+ unknown = ", ".join(sorted(kwargs))
+ raise TypeError(f"Unknown SVD option(s): {unknown}")
+ if options is None:
+ return SVDOptions(**overrides)
+ return replace(options, **overrides)
+
+
+def _axis_positions(tensor: Tensor, ixs: Sequence[Ix]) -> list[int]:
+ """Map a sequence of ``Ix`` handles to their axis positions in a tensor.
+
+ Args:
+ tensor: Nicole tensor whose itags should be searched.
+ ixs: Sequence of :class:`Ix` handles.
+
+ Returns:
+ The matching axis positions in ``tensor``.
+ """
+ positions: list[int] = []
+ for ix in ixs:
+ try:
+ positions.append(tensor.itags.index(ix.itag))
+ except ValueError as exc:
+ raise ValueError(f"itag {ix.itag!r} is missing from tensor tags {tensor.itags}.") from exc
+ return positions
+
+
+def _axes_arg(positions: Sequence[int]) -> int | list[int]:
+ """Convert one or many positions into Nicole's decomp ``axes`` argument.
+
+ Args:
+ positions: Axis positions selected for one side of a decomposition.
+
+ Returns:
+ Either a single integer or a list of integers.
+ """
+ return positions[0] if len(positions) == 1 else list(positions)
+
+
+def _remaining_positions(tensor: Tensor, selected: Sequence[int]) -> list[int]:
+ """Return every axis position not present in ``selected``.
+
+ Args:
+ tensor: Nicole tensor whose axes are being partitioned.
+ selected: Selected axis positions.
+
+ Returns:
+ The complementary axis positions.
+ """
+ selected_set = set(selected)
+ return [axis for axis in range(len(tensor.indices)) if axis not in selected_set]
+
+
+def _bond_ix(tensor: Tensor, axis: int) -> Ix:
+ """Wrap one Nicole tensor leg as an :class:`Ix`.
+
+ Args:
+ tensor: Nicole tensor.
+ axis: Axis to wrap.
+
+ Returns:
+ An :class:`Ix` view of that tensor leg.
+ """
+ index = tensor.indices[axis]
+ return Ix(tensor.itags[axis], int(index.dim), index.direction, index.sectors, index.group)
+
+
+def _clone_tensor_with_ixs(tensor: Tensor, ixs: Sequence[Ix]) -> Tensor:
+ """Clone a Nicole tensor and replace its index metadata with ``ixs``.
+
+ Args:
+ tensor: Tensor whose data should be preserved.
+ ixs: Replacement wrapped indices.
+
+ Returns:
+ A cloned Nicole tensor with the requested tags and indices.
+ """
+ if len(ixs) != len(tensor.indices):
+ raise ValueError(f"Cannot reattach {len(ixs)} indices to rank-{len(tensor.indices)} tensor.")
+ indices = tuple(ix.nicole() for ix in ixs)
+ itags = tuple(ix.itag for ix in ixs)
+ data = {tuple(key): block.clone() for key, block in tensor.data.items()}
+ intw = None if tensor.intw is None else {tuple(key): bridge.clone() for key, bridge in tensor.intw.items()}
+ return Tensor(indices=indices, itags=itags, data=data, intw=intw, dtype=tensor.dtype, label=tensor.label)
+
+
+def qr(A: Tensor, Qixs: Sequence[Ix], tag: str = "b", positive: bool = False):
+ """Split a tensor into ``(Q, R, new_bond)`` using Nicole's QR.
+
+ Args:
+ A: Tensor to split.
+ Qixs: Legs that should remain on the ``Q`` side.
+ tag: Base tag for the new bond.
+ positive: Preserved for API compatibility. Nicole's native phase
+ convention is used unchanged.
+
+ Returns:
+ ``(Q, R, bond)`` where ``bond`` is the new :class:`Ix` wrapper.
+ """
+ positions = _axis_positions(A, Qixs)
+ bond_tag = fresh_itag(tag)
+ Q, R = _ndecomp(A, _axes_arg(positions), mode="QR", itag=bond_tag)
+
+ # The port never requests positive QR phases, but we keep the parameter so
+ # callers can remain close to the Julia API.
+ if positive:
+ pass
+
+ bond = _bond_ix(Q, len(Q.indices) - 1)
+ remaining = [_bond_ix(A, axis) for axis in _remaining_positions(A, positions)]
+ q_ixs = [*Qixs, bond]
+ r_bond = Ix(bond.itag, bond.dim, R.indices[0].direction, bond.sectors, bond.group)
+ r_ixs = [r_bond, *remaining]
+ return _clone_tensor_with_ixs(Q, q_ixs), _clone_tensor_with_ixs(R, r_ixs), bond
+
+
+def lq(A: Tensor, Qixs: Sequence[Ix], tag: str = "b"):
+ """Split a tensor into a left factor and right isometry via Nicole ``LV``.
+
+ Args:
+ A: Tensor to split.
+ Qixs: Legs that should remain on the right-isometric factor.
+ tag: Base tag for the new bond.
+
+ Returns:
+ ``(L, Q, bond)`` where ``Q`` is right-isometric and ``bond`` is the new
+ :class:`Ix` wrapper.
+ """
+ q_positions = _axis_positions(A, Qixs)
+ left_positions = _remaining_positions(A, q_positions)
+ bond_tag = fresh_itag(tag)
+ L, Q = _ndecomp(A, _axes_arg(left_positions), mode="LV", itag=bond_tag)
+
+ left_ixs = [_bond_ix(A, axis) for axis in left_positions]
+ bond = _bond_ix(Q, 0)
+ l_bond = Ix(bond.itag, bond.dim, L.indices[-1].direction, bond.sectors, bond.group)
+ return _clone_tensor_with_ixs(L, [*left_ixs, l_bond]), _clone_tensor_with_ixs(Q, [bond, *Qixs]), bond
+
+
+def truncate(s: torch.Tensor, maxdim: int | float, cutoff: float) -> int:
+ """Compute how many singular values should be kept.
+
+ Args:
+ s: Singular values sorted in descending order.
+ maxdim: Explicit cap on the kept rank.
+ cutoff: Relative cutoff measured against ``abs(s[0])``.
+
+ Returns:
+ The kept rank after applying the cutoff and cap.
+ """
+ if s.numel() == 0:
+ return 0
+ if maxdim is None or maxdim == math.inf:
+ maxdim_int = int(s.numel())
+ else:
+ maxdim_int = max(1, int(maxdim))
+
+ thresh = float(cutoff) * float(torch.abs(s[0]))
+ keep = int((torch.abs(s) > thresh).sum().item())
+ if keep == 0:
+ keep = 1
+ return min(keep, maxdim_int, int(s.numel()))
+
+
+def svd(
+ A: Tensor,
+ Uixs: Sequence[Ix],
+ *,
+ options: SVDOptions | None = None,
+ **kwargs: Any,
+):
+ """Split a tensor into ``(U, S, V, bond_u, bond_v)`` using Nicole's SVD.
+
+ Args:
+ A: Tensor to split.
+ Uixs: Legs that should remain on the left factor ``U``.
+ options: Optional :class:`SVDOptions` instance.
+ **kwargs: Legacy overrides such as ``maxdim=64`` or ``cutoff=1e-12``.
+
+ Returns:
+ ``(U, S, V, bond_u, bond_v)`` with the richer :class:`Ix` metadata
+ restored on the tensor factors.
+ """
+ options = _coerce_svd_options(options, **kwargs)
+ positions = _axis_positions(A, Uixs)
+ left_tag = fresh_itag(options.tag)
+ right_tag = fresh_itag(f"{options.tag}r")
+
+ trunc_spec: dict[str, int | float] = {}
+ if options.maxdim is not None and options.maxdim != math.inf:
+ trunc_spec["nkeep"] = int(options.maxdim)
+ if options.cutoff > 0:
+ trunc_spec["thresh"] = float(options.cutoff)
+
+ U, S, V = _ndecomp(
+ A,
+ _axes_arg(positions),
+ mode="SVD",
+ itag=(left_tag, right_tag),
+ trunc=trunc_spec or None,
+ )
+
+ bond_u = _bond_ix(U, len(U.indices) - 1)
+ bond_v = _bond_ix(V, 0)
+ remaining = [_bond_ix(A, axis) for axis in _remaining_positions(A, positions)]
+ return (
+ _clone_tensor_with_ixs(U, [*Uixs, bond_u]),
+ S,
+ _clone_tensor_with_ixs(V, [bond_v, *remaining]),
+ bond_u,
+ bond_v,
+ )
+
+
+def random_unitary(m: int, n: int | None = None, dtype: torch.dtype = torch.complex128) -> torch.Tensor:
+ """Sample a matrix with orthonormal columns.
+
+ Args:
+ m: Ambient row dimension.
+ n: Number of orthonormal columns. Defaults to ``m``.
+ dtype: Output dtype.
+
+ Returns:
+ An ``m x n`` matrix whose columns are orthonormal.
+ """
+ if n is None:
+ n = m
+ if n > m:
+ raise ValueError(f"n must satisfy n<=m; got n={n}, m={m}")
+
+ real = torch.randn((m, m), dtype=torch.float64)
+ if dtype.is_complex:
+ imag = torch.randn((m, m), dtype=torch.float64)
+ mat = (real + 1j * imag).to(dtype=dtype)
+ else:
+ mat = real.to(dtype=dtype)
+
+ q, r = torch.linalg.qr(mat, mode="reduced")
+
+ # Normalize the QR phases so the result is invariant under the arbitrary QR
+ # sign/phase convention returned by torch.
+ diag = torch.diagonal(r)
+ phases = torch.ones_like(diag)
+ nonzero = diag != 0
+ phases[nonzero] = diag[nonzero] / torch.abs(diag[nonzero])
+ q = q * phases.conj().unsqueeze(0)
+ return q[:, :n]
+
+
+def qr_nonzero_diagonal_rank(rmat: torch.Tensor, tol: float | None = None) -> int:
+ """Estimate the numerical rank of a QR ``R`` factor.
+
+ Args:
+ rmat: Upper-triangular QR factor.
+ tol: Optional magnitude threshold.
+
+ Returns:
+ The number of diagonal entries above ``tol``.
+ """
+ diag = torch.abs(torch.diagonal(rmat))
+ if diag.numel() == 0:
+ return 0
+ if tol is None:
+ tol = max(rmat.shape) * torch.finfo(diag.dtype).eps * float(diag.max())
+ return int(torch.count_nonzero(diag > tol).item())
+
+
+def qr_column_basis(a: torch.Tensor, tol: float | None = None) -> tuple[torch.Tensor, int]:
+ """Return an orthonormal basis for the column space of ``a``.
+
+ Args:
+ a: Dense matrix.
+ tol: Optional QR rank tolerance.
+
+ Returns:
+ ``(basis, rank)``.
+ """
+ q, r = torch.linalg.qr(a, mode="reduced")
+ rank = qr_nonzero_diagonal_rank(r, tol)
+ return q[:, :rank], rank
+
+
+def qr_row_basis(a: torch.Tensor, tol: float | None = None) -> tuple[torch.Tensor, int]:
+ """Return a row-orthonormal basis for the row space of ``a``.
+
+ Args:
+ a: Dense matrix.
+ tol: Optional QR rank tolerance.
+
+ Returns:
+ ``(basis, rank)`` where the basis rows span the row space of ``a``.
+ """
+ q, r = torch.linalg.qr(a.transpose(-2, -1), mode="reduced")
+ rank = qr_nonzero_diagonal_rank(r, tol)
+ return q[:, :rank].transpose(-2, -1), rank
+
+
+def identity_overlap_matrix(dtype: torch.dtype, n: int, *, device: torch.device | None = None) -> torch.Tensor:
+ """Return an ``n x n`` identity matrix for no-augmentation overlaps.
+
+ Args:
+ dtype: Matrix dtype.
+ n: Matrix size.
+ device: Optional torch device.
+
+ Returns:
+ An identity matrix.
+ """
+ return torch.eye(n, dtype=dtype, device=device)
+
+
+def reconstruct_from_svd(U: Tensor, S: Tensor, V: Tensor) -> Tensor:
+ """Reconstruct ``U * S * V`` in tensor form.
+
+ Args:
+ U: Left SVD tensor.
+ S: Diagonal singular-value tensor.
+ V: Right SVD tensor.
+
+ Returns:
+ The contracted reconstruction.
+ """
+ return tcontract(tcontract(U, S), V)
diff --git a/src/alice/algorithm/bond_update_bug/_kernel/local_solvers.py b/src/alice/algorithm/bond_update_bug/_kernel/local_solvers.py
new file mode 100644
index 0000000..ffc7a5e
--- /dev/null
+++ b/src/alice/algorithm/bond_update_bug/_kernel/local_solvers.py
@@ -0,0 +1,243 @@
+# Copyright (C) 2025-2026 Changkai Zhang.
+#
+# This file is part of Alice project.
+#
+# Alice is free software: you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published
+# by the Free Software Foundation, either version 3 of the License,
+# or (at your option) any later version.
+#
+# Alice is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with Alice. If not, see .
+# Author of code: Madhav Menon.
+
+
+"""Pluggable local time-integrators for the (imaginary-time) BUG substeps.
+
+Every BUG local substep computes ``y = exp(tau * A) x`` for a matrix-free tensor
+action ``A`` (``apply``) and a (generally complex) local timestep ``tau``. In
+**unitary** real-time evolution this must be an *exact* exponential, so the
+kernel uses a Krylov ``expv``. In **imaginary time** (cooling toward the
+ground state) the evolution is no longer unitary, and ``y = exp(tau A) x`` is just
+the exact flow of the linear ODE ``x'(s) = A x(s)`` over ``s in [0, tau]`` — *any*
+stable integrator of that ODE may be used. This module provides a family of them
+behind one uniform ``(apply, tau, x)`` call surface so the discarded-projector BUG
+and the bond_update_bug (``variant='discarded'``) can swap the local solver:
+
+ * ``'krylov'`` : Lanczos (Hermitian) / Arnoldi (general) exponential — ``≈`` exact.
+ * ``'midpoint'`` : explicit midpoint (RK2), ``substeps`` internal steps — 2nd order,
+ explicit (fast), conditionally stable.
+ * ``'rk4'`` : classical Runge–Kutta 4, ``substeps`` internal steps — 4th order,
+ explicit, conditionally stable.
+ * ``'trapezoid'`` : implicit trapezoidal / Crank–Nicolson, ``substeps`` internal steps —
+ 2nd order, A-stable (unconditionally stable), one matrix-free
+ linear solve per substep.
+
+For the substepped integrators the local exponential error is ``O((tau/substeps)^p)``
+(``p = 2`` for midpoint/trapezoid, ``p = 4`` for rk4); raising ``substeps`` converges
+monotonically to the exact action (and, for the two explicit schemes, also restores
+stability when ``|tau| * ||A||`` is large). Everything stays in the symmetry-blocked
+Nicole tensor representation (only ``apply``, tensor add/scale, and inner products are
+used), so no admissible-block structure is ever broken.
+"""
+
+from __future__ import annotations
+
+import torch
+from nicole import Tensor
+
+from .krylov import tensor_inner, tensor_lanczos_expv
+
+# Solver names accepted by :func:`local_expv`.
+LOCAL_SOLVERS = ('krylov', 'midpoint', 'rk4', 'trapezoid')
+
+
+def _norm(x: Tensor) -> float:
+ n = x.norm()
+ return float(n.real if hasattr(n, 'real') else n)
+
+
+# ---------------------------------------------------------------------------
+# Krylov (general / non-Hermitian) exponential — the Arnoldi counterpart of the
+# Hermitian tensor Lanczos in :mod:`.krylov`.
+# ---------------------------------------------------------------------------
+
+def tensor_arnoldi_expv(apply, tau: complex, x: Tensor, *, maxiter: int = 30, tol: float = 1e-15) -> Tensor:
+ """Return ``exp(tau * A) @ x`` for a NON-Hermitian Nicole-tensor action ``apply``.
+
+ A tensor-native Arnoldi (modified Gram–Schmidt) exponential: builds an
+ orthonormal Krylov basis of Nicole tensors and a small dense upper-Hessenberg
+ matrix ``H``, then forms ``y = beta * V * exp(tau H) e1``. Stays in the
+ symmetry-blocked representation throughout (the non-Hermitian counterpart of
+ :func:`alice.algorithm.bond_update_bug._kernel.krylov.tensor_lanczos_expv`).
+ """
+ beta0 = _norm(x)
+ if beta0 == 0.0:
+ return x
+ m = max(int(maxiter), 1)
+ basis = [(1.0 / beta0) * x]
+ H = torch.zeros((m, m), dtype=torch.complex128)
+ used = 1
+ for j in range(m):
+ w = apply(basis[j])
+ for i in range(j + 1):
+ hij = tensor_inner(basis[i], w)
+ H[i, j] = hij
+ w = w + (-hij) * basis[i]
+ used = j + 1
+ nrm = _norm(w)
+ if nrm <= tol or j == m - 1:
+ break
+ H[j + 1, j] = nrm
+ basis.append((1.0 / nrm) * w)
+
+ Hk = H[:used, :used]
+ coeff = torch.linalg.matrix_exp(tau * Hk)[:, 0] * beta0
+ out = coeff[0] * basis[0]
+ for idx in range(1, used):
+ out = out + coeff[idx] * basis[idx]
+ return out
+
+
+# ---------------------------------------------------------------------------
+# Explicit Runge–Kutta exponential actions (midpoint / RK4)
+# ---------------------------------------------------------------------------
+
+def tensor_rk_expv(apply, tau: complex, x: Tensor, *, order: int, substeps: int) -> Tensor:
+ """Approximate ``exp(tau A) x`` by explicit RK integration of ``x' = A x``.
+
+ Integrates the linear ODE over ``s in [0, tau]`` with ``substeps`` equal steps
+ ``h = tau / substeps``. ``order=2`` is the explicit midpoint rule (RK2),
+ ``order=4`` the classical RK4. Explicit and matrix-free (only ``apply`` and
+ tensor arithmetic), so it is fast but conditionally stable: the local error is
+ ``O(h^order)`` and stability needs ``|h| * ||A||`` inside the method's stability
+ region — raise ``substeps`` if either is violated.
+ """
+ n = max(int(substeps), 1)
+ h = tau / n
+ y = x
+ if order == 2:
+ for _ in range(n):
+ k1 = apply(y)
+ k2 = apply(y + (0.5 * h) * k1)
+ y = y + h * k2
+ elif order == 4:
+ for _ in range(n):
+ k1 = apply(y)
+ k2 = apply(y + (0.5 * h) * k1)
+ k3 = apply(y + (0.5 * h) * k2)
+ k4 = apply(y + h * k3)
+ y = y + (h / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)
+ else:
+ raise ValueError(f"tensor_rk_expv: unsupported order {order!r} (use 2 or 4)")
+ return y
+
+
+# ---------------------------------------------------------------------------
+# Implicit trapezoidal (Crank–Nicolson) exponential action + matrix-free GMRES
+# ---------------------------------------------------------------------------
+
+def tensor_gmres(linop, b: Tensor, *, tol: float = 1e-12, maxiter: int = 60) -> Tensor:
+ """Solve ``linop(y) = b`` for a matrix-free Nicole-tensor linear operator.
+
+ Full (non-restarted) GMRES from the zero initial guess, working entirely in the
+ symmetry-blocked tensor representation (Arnoldi + a small dense least-squares on
+ the Hessenberg matrix). The local systems here are tiny, so a handful of
+ iterations reach ``tol``; ``maxiter`` is the hard cap.
+ """
+ bnorm = _norm(b)
+ if bnorm == 0.0:
+ return b
+ m = max(int(maxiter), 1)
+ V = [(1.0 / bnorm) * b]
+ H = torch.zeros((m + 1, m), dtype=torch.complex128)
+ e1 = torch.zeros(m + 1, dtype=torch.complex128)
+ e1[0] = bnorm
+ for j in range(m):
+ w = linop(V[j])
+ for i in range(j + 1):
+ H[i, j] = tensor_inner(V[i], w)
+ w = w + (-H[i, j]) * V[i]
+ hjj = _norm(w)
+ H[j + 1, j] = hjj
+ # Least-squares solve of the (j+2, j+1) Hessenberg system for the residual.
+ y, *_ = torch.linalg.lstsq(H[:j + 2, :j + 1], e1[:j + 2].unsqueeze(1))
+ y = y.squeeze(1)
+ resid = float(torch.linalg.norm(e1[:j + 2] - H[:j + 2, :j + 1] @ y).real)
+ if hjj <= tol * bnorm or resid <= tol * bnorm or j == m - 1:
+ sol = y[0] * V[0]
+ for idx in range(1, j + 1):
+ sol = sol + y[idx] * V[idx]
+ return sol
+ V.append((1.0 / hjj) * w)
+ # Unreachable: the loop always returns at j == m - 1.
+ raise RuntimeError("tensor_gmres did not return")
+
+
+def tensor_trapezoid_expv(apply, tau: complex, x: Tensor, *, substeps: int,
+ tol: float = 1e-12, maxiter: int = 60) -> Tensor:
+ """Approximate ``exp(tau A) x`` by the implicit trapezoidal rule (Crank–Nicolson).
+
+ Each of the ``substeps`` steps ``h = tau / substeps`` advances
+ ``(I - (h/2) A) y_{k+1} = (I + (h/2) A) y_k`` — the (1,1)-Padé approximant of
+ ``exp(h A)``. It is 2nd order and A-stable (the stability function maps the left
+ half-plane into the unit disc), so it never blows up however large ``|h| * ||A||``
+ is; raising ``substeps`` drives the ``O(h^2)`` error down. The implicit solve is
+ a matrix-free GMRES (:func:`tensor_gmres`).
+ """
+ n = max(int(substeps), 1)
+ h = tau / n
+ c = 0.5 * h
+ y = x
+ for _ in range(n):
+ rhs = y + c * apply(y) # (I + (h/2) A) y_k
+ y = tensor_gmres(lambda z: z + (-c) * apply(z), rhs, tol=tol, maxiter=maxiter)
+ return y
+
+
+# ---------------------------------------------------------------------------
+# Dispatcher
+# ---------------------------------------------------------------------------
+
+def local_expv(apply, tau: complex, x: Tensor, *, solver: str = 'krylov', substeps: int = 1,
+ hermitian: bool = False, krylov_maxiter: int = 30, krylov_tol: float = 1e-15) -> Tensor:
+ """Compute ``exp(tau A) x`` with the requested local integrator.
+
+ Parameters
+ ----------
+ apply:
+ Matrix-free tensor action ``x -> A x``.
+ tau:
+ Local timestep (already including the evolution prefactor, e.g. ``-dt`` for
+ imaginary time, ``-1j*dt`` for real time).
+ x:
+ Input tensor.
+ solver:
+ One of :data:`LOCAL_SOLVERS`.
+ substeps:
+ Number of internal steps for the substepped integrators (ignored by
+ ``'krylov'``).
+ hermitian:
+ Whether ``A`` is Hermitian — selects Lanczos vs Arnoldi for ``'krylov'``
+ (ignored by the other solvers).
+ krylov_maxiter, krylov_tol:
+ Krylov dimension cap and tolerance for ``'krylov'`` (also used as the GMRES
+ tolerance / cap for ``'trapezoid'``).
+ """
+ if solver == 'krylov':
+ if hermitian:
+ return tensor_lanczos_expv(apply, tau, x, maxiter=krylov_maxiter, tol=krylov_tol)
+ return tensor_arnoldi_expv(apply, tau, x, maxiter=krylov_maxiter, tol=krylov_tol)
+ if solver == 'midpoint':
+ return tensor_rk_expv(apply, tau, x, order=2, substeps=substeps)
+ if solver == 'rk4':
+ return tensor_rk_expv(apply, tau, x, order=4, substeps=substeps)
+ if solver == 'trapezoid':
+ return tensor_trapezoid_expv(apply, tau, x, substeps=substeps,
+ tol=krylov_tol, maxiter=max(krylov_maxiter, 60))
+ raise ValueError(f"unknown local solver {solver!r}; recognised values are: {', '.join(LOCAL_SOLVERS)}")
diff --git a/src/alice/algorithm/bond_update_bug/_kernel/nicole_helpers.py b/src/alice/algorithm/bond_update_bug/_kernel/nicole_helpers.py
new file mode 100644
index 0000000..1a50706
--- /dev/null
+++ b/src/alice/algorithm/bond_update_bug/_kernel/nicole_helpers.py
@@ -0,0 +1,508 @@
+# Copyright (C) 2025-2026 Changkai Zhang.
+#
+# This file is part of Alice project.
+#
+# Alice is free software: you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published
+# by the Free Software Foundation, either version 3 of the License,
+# or (at your option) any later version.
+#
+# Alice is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with Alice. If not, see .
+
+"""Nicole tensor convenience helpers used across the BUG stack.
+
+These helpers keep the tensor-manipulation code readable by centralizing a few
+repeated chores: Fortran-order reshaping, dense/block conversion, explicit
+identity construction, and a small collection of Nicole-flavored conjugation and
+contraction utilities.
+"""
+
+from __future__ import annotations
+
+import math
+from dataclasses import dataclass, replace
+from typing import Any, Iterable, Mapping, Sequence
+
+import torch
+from nicole import Direction, Tensor
+from nicole import conj as _nconj
+from nicole import contract as _ncontract
+from nicole import identity as _nidentity
+from nicole import inv as _ninv
+from nicole.blocks import BlockSchema
+
+from .indices import Ix, resolved_sectors
+
+__all__ = [
+ "IdentityOptions",
+ "conj",
+ "dag",
+ "delta",
+ "diag_tensor",
+ "flatten_fortran",
+ "identity_tensor",
+ "inv_diag",
+ "make_tensor",
+ "norm",
+ "prime_bra",
+ "reshape_fortran",
+ "scalar",
+ "star_itags",
+ "tcontract",
+ "to_dense",
+]
+
+
+@dataclass(frozen=True)
+class IdentityOptions:
+ """Options for raw-tag identity and relabeling tensors.
+
+ Args:
+ dim: Dense index dimension used when the source is given only as a tag.
+ sectors: Optional explicit sector tuple for raw-tag construction.
+ group: Optional Nicole symmetry group handle for raw-tag construction.
+ direction: Nicole direction of the source leg when only raw tag metadata
+ is provided.
+ """
+
+ dim: int | None = None
+ sectors: Sequence[object] | None = None
+ group: object | None = None
+ direction: Direction = Direction.IN
+
+
+_IDENTITY_OPTION_FIELDS = {field.name for field in IdentityOptions.__dataclass_fields__.values()}
+
+
+def _coerce_identity_options(options: IdentityOptions | None = None, **kwargs: Any) -> IdentityOptions:
+ """Normalize identity-construction options from an object or kwargs.
+
+ Args:
+ options: Existing options object to start from.
+ **kwargs: Field overrides for :class:`IdentityOptions`.
+
+ Returns:
+ A normalized :class:`IdentityOptions` instance.
+ """
+ overrides = {name: kwargs.pop(name) for name in list(kwargs) if name in _IDENTITY_OPTION_FIELDS}
+ if kwargs:
+ unknown = ", ".join(sorted(kwargs))
+ raise TypeError(f"Unknown identity option(s): {unknown}")
+ if options is None:
+ return IdentityOptions(**overrides)
+ return replace(options, **overrides)
+
+
+def _materialize_indices(ixs: Sequence[Ix]) -> tuple:
+ """Convert wrapped indices into Nicole indices.
+
+ Args:
+ ixs: Sequence of :class:`Ix` wrappers.
+
+ Returns:
+ A tuple of Nicole :class:`Index` objects.
+ """
+ return tuple(ix.nicole() for ix in ixs)
+
+
+def _index_offsets(ix: Ix | object) -> tuple[dict[object, tuple[int, int]], int]:
+ """Map each sector charge to its dense offset along one axis.
+
+ Args:
+ ix: Wrapped or native Nicole index.
+
+ Returns:
+ A pair ``(offsets, total_dim)``.
+ """
+ offsets: dict[object, tuple[int, int]] = {}
+ cursor = 0
+ for sector in resolved_sectors(ix):
+ offsets[sector.charge] = (cursor, sector.dim)
+ cursor += sector.dim
+ return offsets, cursor
+
+
+def reshape_fortran(tensor: torch.Tensor, shape: Sequence[int]) -> torch.Tensor:
+ """Return the torch equivalent of ``reshape(..., order='F')``.
+
+ Args:
+ tensor: Input torch tensor.
+ shape: Target shape interpreted in Fortran/column-major order.
+
+ Returns:
+ A reshaped tensor with the requested shape.
+ """
+ target = tuple(int(dim) for dim in shape)
+ if math.prod(target) != int(tensor.numel()):
+ raise ValueError(f"Cannot reshape tensor with {tensor.numel()} entries into {target}.")
+ if len(target) == 0:
+ return tensor.reshape(())
+ if tensor.ndim == 0:
+ return tensor.reshape(target)
+
+ # Reverse, reshape, then reverse back to emulate column-major memory order.
+ rev_in = tuple(reversed(range(tensor.ndim)))
+ rev_out = tuple(reversed(range(len(target))))
+ reshaped = tensor.permute(rev_in).contiguous().reshape(tuple(reversed(target)))
+ return reshaped.permute(rev_out).contiguous()
+
+
+def flatten_fortran(tensor: torch.Tensor) -> torch.Tensor:
+ """Return the torch equivalent of ``reshape(-1, order='F')``.
+
+ Args:
+ tensor: Input torch tensor.
+
+ Returns:
+ A one-dimensional tensor flattened in Fortran/column-major order.
+ """
+ if tensor.ndim <= 1:
+ return tensor.reshape(-1)
+ return tensor.permute(tuple(reversed(range(tensor.ndim)))).contiguous().reshape(-1)
+
+
+def _dense_to_block_data(
+ dense: torch.Tensor,
+ ixs: Sequence[Ix],
+ *,
+ tol: float = 1e-12,
+) -> dict[tuple[object, ...], torch.Tensor]:
+ """Project a dense tensor into the admissible Nicole block dictionary.
+
+ Args:
+ dense: Dense torch tensor with one axis per index in ``ixs``.
+ ixs: Wrapped indices describing the target block structure.
+ tol: Numerical tolerance used when discarding zero blocks and checking
+ for amplitudes outside the admissible symmetry support.
+
+ Returns:
+ A Nicole-style block dictionary keyed by sector charges.
+ """
+ indices = _materialize_indices(ixs)
+ if dense.ndim != len(indices):
+ raise ValueError(f"Dense tensor rank {dense.ndim} does not match index rank {len(indices)}.")
+
+ expected_shape = tuple(int(ix.dim) for ix in indices)
+ if tuple(dense.shape) != expected_shape:
+ raise ValueError(f"Dense tensor shape {tuple(dense.shape)} does not match index dims {expected_shape}.")
+
+ offsets = [_index_offsets(ix) for ix in ixs]
+ reconstructed = torch.zeros(expected_shape, dtype=dense.dtype, device=dense.device)
+ data: dict[tuple[object, ...], torch.Tensor] = {}
+
+ # Nicole only stores symmetry-admissible blocks, so we reconstruct the
+ # admissible support and confirm nothing significant lives outside it.
+ for key in BlockSchema.iter_admissible_keys(indices):
+ if not BlockSchema.charges_conserved(indices, key):
+ continue
+ slices = tuple(
+ slice(offsets[axis][0][key[axis]][0], offsets[axis][0][key[axis]][0] + offsets[axis][0][key[axis]][1])
+ for axis in range(len(key))
+ )
+ block = dense[slices].clone().contiguous()
+ reconstructed[slices] = block
+ if block.numel() == 0:
+ continue
+ if torch.max(torch.abs(block)).item() > tol:
+ data[tuple(key)] = block
+
+ residual = dense - reconstructed
+ if residual.numel() and torch.max(torch.abs(residual)).item() > tol:
+ raise ValueError("Dense tensor contains amplitudes outside the admissible symmetry blocks.")
+
+ return data
+
+
+def make_tensor(
+ array_or_blocks: torch.Tensor | Mapping[tuple[object, ...], object] | object,
+ ixs: Sequence[Ix],
+ *,
+ dtype: torch.dtype | None = None,
+ tol: float = 1e-12,
+) -> Tensor:
+ """Build a Nicole tensor from dense data or explicit block data.
+
+ Args:
+ array_or_blocks: Dense tensor-like data or a Nicole block dictionary.
+ ixs: Wrapped indices describing the target tensor legs.
+ dtype: Optional dtype override.
+ tol: Tolerance forwarded to dense-to-block conversion.
+
+ Returns:
+ A Nicole :class:`Tensor` with the requested indices and tags.
+ """
+ indices = _materialize_indices(ixs)
+ itags = tuple(ix.itag for ix in ixs)
+
+ if isinstance(array_or_blocks, Mapping):
+ data = {tuple(key): torch.as_tensor(value) for key, value in array_or_blocks.items()}
+ if dtype is not None:
+ data = {key: value.to(dtype=dtype) for key, value in data.items()}
+ inferred_dtype = next(iter(data.values())).dtype if data else (dtype or torch.complex128)
+ return Tensor(indices=indices, itags=itags, data=data, dtype=inferred_dtype)
+
+ dense = torch.as_tensor(array_or_blocks)
+ if dtype is not None:
+ dense = dense.to(dtype=dtype)
+
+ if len(ixs) == 0:
+ scalar_value = dense.reshape(())
+ return Tensor(indices=(), itags=(), data={(): scalar_value}, dtype=scalar_value.dtype)
+
+ data = _dense_to_block_data(dense, ixs, tol=tol)
+ return Tensor(indices=indices, itags=itags, data=data, dtype=dense.dtype)
+
+
+def to_dense(tensor: Tensor, itag_order: Sequence[str]) -> torch.Tensor:
+ """Assemble a dense tensor in the requested itag order.
+
+ Args:
+ tensor: Nicole tensor to densify.
+ itag_order: Desired order of the tensor tags in the dense output.
+
+ Returns:
+ A dense torch tensor with axes permuted to ``itag_order``.
+ """
+ if len(tensor.indices) == 0:
+ return next(iter(tensor.data.values())).reshape(())
+ if len(itag_order) != len(tensor.itags):
+ raise ValueError(f"itag_order length {len(itag_order)} does not match tensor rank {len(tensor.itags)}.")
+
+ offsets = [_index_offsets(index) for index in tensor.indices]
+ shape = tuple(total_dim for _, total_dim in offsets)
+ full = torch.zeros(shape, dtype=tensor.dtype, device=tensor.device)
+ for key, block in tensor.data.items():
+ slices = tuple(
+ slice(offsets[axis][0][key[axis]][0], offsets[axis][0][key[axis]][0] + offsets[axis][0][key[axis]][1])
+ for axis in range(len(key))
+ )
+ full[slices] = block
+
+ positions: dict[str, list[int]] = {}
+ for axis, tag in enumerate(tensor.itags):
+ positions.setdefault(tag, []).append(axis)
+
+ used: dict[str, int] = {}
+ permutation: list[int] = []
+ for tag in itag_order:
+ taken = used.get(tag, 0)
+ axes = positions.get(tag)
+ if axes is None or taken >= len(axes):
+ raise ValueError(f"itag '{tag}' is missing from tensor tags {tensor.itags}.")
+ permutation.append(axes[taken])
+ used[tag] = taken + 1
+
+ return full.permute(permutation).contiguous()
+
+
+def dag(tensor: Tensor) -> Tensor:
+ """Return Nicole's conjugated tensor with flipped directions.
+
+ Args:
+ tensor: Input Nicole tensor.
+
+ Returns:
+ The Nicole ``dag``/conjugation result.
+ """
+ return _nconj(tensor)
+
+
+def conj(tensor: Tensor) -> Tensor:
+ """Alias for Nicole's conjugation helper.
+
+ Args:
+ tensor: Input Nicole tensor.
+
+ Returns:
+ The conjugated Nicole tensor.
+ """
+ return _nconj(tensor)
+
+
+def star_itags(tensor: Tensor) -> Tensor:
+ """Clone a tensor and append ``*`` to every Nicole tag.
+
+ Args:
+ tensor: Input Nicole tensor.
+
+ Returns:
+ A cloned tensor with starred tags.
+ """
+ retagged = tensor.clone()
+ retagged.retag({tag: f"{tag}*" for tag in retagged.itags})
+ return retagged
+
+
+def prime_bra(tensor: Tensor) -> Tensor:
+ """Return the Nicole analogue of ``dag(prime(x))``.
+
+ Args:
+ tensor: Input Nicole tensor.
+
+ Returns:
+ A conjugated tensor whose itags have been starred.
+ """
+ bra = _nconj(tensor)
+ bra.retag({tag: f"{tag}*" for tag in bra.itags})
+ return bra
+
+
+def tcontract(
+ left: Tensor,
+ right: Tensor,
+ axes: tuple[int, int] | tuple[Sequence[int], Sequence[int]] | None = None,
+) -> Tensor:
+ """Contract two Nicole tensors with a gentle outer-product fallback.
+
+ Args:
+ left: Left tensor.
+ right: Right tensor.
+ axes: Optional explicit contraction axes.
+
+ Returns:
+ The Nicole contraction result.
+ """
+ if axes is not None:
+ return _ncontract(left, right, axes=axes)
+ try:
+ return _ncontract(left, right)
+ except ValueError as exc:
+ if "No valid contraction pairs found" not in str(exc):
+ raise
+ return _ncontract(left, right, axes=([], []))
+
+
+def scalar(tensor: Tensor):
+ """Extract a Python scalar from a scalar-like Nicole tensor.
+
+ Args:
+ tensor: Scalar tensor or tensor whose open indices all have dimension 1.
+
+ Returns:
+ The scalar value stored in the tensor.
+ """
+ if tensor.is_scalar():
+ return tensor.item()
+ if all(int(index.dim) == 1 for index in tensor.indices):
+ return to_dense(tensor, list(tensor.itags)).reshape(-1)[0].item()
+ raise ValueError("scalar() requires a scalar tensor or all-dimension-1 open indices.")
+
+
+def norm(tensor: Tensor):
+ """Return Nicole's norm for one tensor.
+
+ Args:
+ tensor: Input Nicole tensor.
+
+ Returns:
+ The tensor norm in Nicole's backend dtype.
+ """
+ return tensor.norm()
+
+
+def delta(
+ source: Ix | str,
+ out_itag: str,
+ dim: int | None = None,
+ *,
+ options: IdentityOptions | None = None,
+ **kwargs: Any,
+) -> Tensor:
+ """Construct a relabeling identity tensor.
+
+ Args:
+ source: Source index wrapper or raw itag string.
+ out_itag: Output Nicole tag.
+ dim: Optional raw-tag dimension. This is ignored when ``source`` is an
+ :class:`Ix`.
+ options: Optional :class:`IdentityOptions` instance.
+ **kwargs: Legacy option overrides such as ``sectors=...`` or
+ ``direction=Direction.IN``.
+
+ Returns:
+ A Nicole identity tensor that relabels one leg to ``out_itag``.
+ """
+ options = _coerce_identity_options(options, dim=dim, **kwargs)
+ if isinstance(source, Ix):
+ return _nidentity(source.nicole(), itags=(source.itag, out_itag))
+ if options.dim is None:
+ raise ValueError("dim is required when building a delta tensor from raw itag metadata.")
+ ix = Ix(source, options.dim, options.direction, None if options.sectors is None else tuple(options.sectors), options.group)
+ return _nidentity(ix.nicole(), itags=(source, out_itag))
+
+
+def identity_tensor(
+ left: Ix | str,
+ right: str | None = None,
+ dim: int | None = None,
+ *,
+ options: IdentityOptions | None = None,
+ **kwargs: Any,
+) -> Tensor:
+ """Construct an identity tensor used to extend local operators.
+
+ Args:
+ left: Source index wrapper or raw itag string.
+ right: Optional output tag. When ``left`` is an :class:`Ix`, the default
+ is ``f"{left.itag}*"``.
+ dim: Optional raw-tag dimension.
+ options: Optional :class:`IdentityOptions` instance.
+ **kwargs: Legacy option overrides forwarded to :func:`delta`.
+
+ Returns:
+ A rank-2 Nicole identity tensor.
+ """
+ options = _coerce_identity_options(options, dim=dim, **kwargs)
+ if isinstance(left, Ix):
+ if right is None:
+ right = f"{left.itag}*"
+ return delta(left, right, options=options)
+ if right is None:
+ raise ValueError("identity_tensor requires a right itag when using raw-tag construction.")
+ return delta(left, right, options=options)
+
+
+def diag_tensor(
+ vec: Iterable[complex] | torch.Tensor,
+ left: Ix | str,
+ right: Ix | str,
+ *,
+ dtype: torch.dtype = torch.float64,
+) -> Tensor:
+ """Build a diagonal Nicole tensor from explicit left and right legs.
+
+ Args:
+ vec: Diagonal values.
+ left: Left index wrapper or raw left tag.
+ right: Right index wrapper or raw right tag.
+ dtype: Output tensor dtype.
+
+ Returns:
+ A rank-2 Nicole tensor whose dense form is ``diag(vec)``.
+ """
+ values = torch.as_tensor(tuple(vec) if not isinstance(vec, torch.Tensor) else vec, dtype=dtype)
+ mat = torch.diag(values)
+ if isinstance(left, Ix) and isinstance(right, Ix):
+ return make_tensor(mat, [left, right], dtype=dtype)
+ if isinstance(left, str) and isinstance(right, str):
+ size = int(values.numel())
+ return make_tensor(mat, [Ix(left, size, Direction.OUT), Ix(right, size, Direction.IN)], dtype=dtype)
+ raise TypeError("diag_tensor requires either two Ix handles or two itag strings.")
+
+
+def inv_diag(diag_tensor_obj: Tensor) -> Tensor:
+ """Invert a diagonal Nicole tensor using Nicole's native helper.
+
+ Args:
+ diag_tensor_obj: Diagonal Nicole tensor.
+
+ Returns:
+ The Nicole inverse tensor.
+ """
+ return _ninv(diag_tensor_obj)
diff --git a/src/alice/algorithm/bond_update_bug/bond.py b/src/alice/algorithm/bond_update_bug/bond.py
new file mode 100644
index 0000000..8735267
--- /dev/null
+++ b/src/alice/algorithm/bond_update_bug/bond.py
@@ -0,0 +1,216 @@
+# Copyright (C) 2025-2026 Changkai Zhang.
+#
+# This file is part of Alice project.
+#
+# Alice is free software: you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published
+# by the Free Software Foundation, either version 3 of the License,
+# or (at your option) any later version.
+#
+# Alice is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with Alice. If not, see .
+# Author of code: Madhav Menon.
+
+
+"""Nearest-neighbour bond Hamiltonians for the bond_update_bug integrator.
+
+The Basis-Update & Galerkin (BUG) integrator (Ceruti, Kusch & Lubich,
+*BIT* 2022; arXiv:2304.05660) evolves an `MPS` under a nearest-neighbour
+Hamiltonian split into commuting odd/even bond groups. Each bond carries the
+*bare* two-site Hamiltonian term `h_{i,i+1}` — not a pre-exponentiated gate. The
+KLS local update (see :mod:`alice.algorithm.bond_update_bug.kls`) exponentiates the
+*projected* effective Hamiltonian internally; this module only supplies the bond
+terms.
+
+The bond Hamiltonian for bond *(i, i+1)* is reused directly from the AutoMPO
+interaction list (`build_interaction`): the leading and terminal MPO tensors of
+an `Interaction2Site` term are contracted over their shared operator channel,
+exactly as `build_hamiltonian` would, so no new operator algebra is introduced.
+
+Index convention (shared with the rest of Alice): a bond Hamiltonian `h` is a
+4-index tensor with axes `(bra_i, ket_i, bra_{i+1}, ket_{i+1})`, the physical
+`bra`/`ket` directions matching the MPS physical index and its dual.
+"""
+
+from __future__ import annotations
+
+from typing import List, Optional
+
+import torch
+from nicole import Tensor, contract, permute
+
+from alice.network.interaction import Interaction, Interaction1Site, Interaction2Site
+
+
+def to_complex(tensor: Tensor) -> Tensor:
+ """Return a copy of `tensor` with every block cast to `complex128`.
+
+ The KLS update exponentiates Hamiltonian terms, so the state must share the
+ `complex128` dtype of the PyTorch backend.
+
+ Parameters
+ ----------
+ tensor:
+ Nicole tensor with real or complex blocks.
+
+ Returns
+ -------
+ Tensor
+ Tensor with identical indices and itags but `complex128` block data.
+ """
+ new_intw = None
+ if tensor.intw is not None:
+ new_intw = {
+ key: bridge.to(tensor.device, dtype=torch.complex128)
+ for key, bridge in tensor.intw.items()
+ }
+ return Tensor(
+ indices=tensor.indices,
+ itags=tensor.itags,
+ data={key: block.to(torch.complex128) for key, block in tensor.data.items()},
+ intw=new_intw,
+ dtype=torch.complex128,
+ )
+
+
+def bond_hamiltonian(intr: Interaction2Site) -> Tensor:
+ """Build the bare two-site bond Hamiltonian of one nearest-neighbour term.
+
+ Contracts the leading and terminal MPO tensors of `intr` over their shared
+ operator channel and drops the two trivial boundary bonds, returning the
+ physical two-site operator scaled by the coupling `intr.cpl`. This mirrors
+ the contraction `build_hamiltonian` performs, so the bond Hamiltonian is
+ exactly the term that enters the AutoMPO Hamiltonian.
+
+ Parameters
+ ----------
+ intr:
+ Nearest-neighbour two-site interaction with populated `leading_tnsr`
+ and `terminal_tnsr` and `terminal_site == leading_site + 1`.
+
+ Returns
+ -------
+ Tensor
+ 4-index bond Hamiltonian with axes
+ `(bra_i, ket_i, bra_{i+1}, ket_{i+1})`.
+
+ Raises
+ ------
+ ValueError
+ If `intr` is not nearest-neighbour, or its tensors are not populated.
+ """
+ if intr.terminal_site != intr.leading_site + 1:
+ raise ValueError(
+ "bond_hamiltonian requires a nearest-neighbour term "
+ f"(terminal_site == leading_site + 1), got leading_site="
+ f"{intr.leading_site}, terminal_site={intr.terminal_site}"
+ )
+ if intr.leading_tnsr is None or intr.terminal_tnsr is None:
+ raise ValueError(
+ "bond_hamiltonian requires populated leading_tnsr and terminal_tnsr; "
+ "build the interaction list with build_interaction first"
+ )
+
+ # leading_tnsr: (L_trivial_IN, op_OUT, bra_i, ket_i)
+ # terminal_tnsr: (op_IN, R_trivial_OUT, bra_{i+1}, ket_{i+1})
+ # Contract the shared operator channel (leading axis 1, terminal axis 0).
+ h = contract(intr.leading_tnsr, intr.terminal_tnsr, axes=(1, 0))
+ # h axes: (L_trivial, bra_i, ket_i, R_trivial, bra_{i+1}, ket_{i+1}).
+ h.squeeze(0) # drop L_trivial -> (bra_i, ket_i, R_trivial, bra_{i+1}, ket_{i+1})
+ h.squeeze(2) # drop R_trivial -> (bra_i, ket_i, bra_{i+1}, ket_{i+1})
+ return h * intr.cpl
+
+
+def build_bond_generators(interactions: List[Interaction], length: int) -> List[Optional[Tensor]]:
+ """Accumulate per-bond Hamiltonians from an AutoMPO interaction list.
+
+ Sums every nearest-neighbour `Interaction2Site` term onto its bond. Bonds
+ with no term are left as `None`. This yields the bond decomposition
+ `H = Σ_b h_b` used by the odd/even Trotter split.
+
+ Parameters
+ ----------
+ interactions:
+ Interaction list from `build_interaction`. Every active term must be a
+ nearest-neighbour `Interaction2Site`.
+ length:
+ Number of sites `L`; there are `L - 1` bonds.
+
+ Returns
+ -------
+ list of (Tensor or None)
+ Length `L - 1`. Entry *b* is the bond Hamiltonian for bond
+ *(b, b+1)*, or `None` if no term acts on that bond.
+
+ Raises
+ ------
+ NotImplementedError
+ If a non-nearest-neighbour two-site term or a one-site term with a
+ non-zero coupling is present (the bond_update_bug integrator targets
+ nearest-neighbour Hamiltonians).
+ """
+ generators: List[Optional[Tensor]] = [None] * (length - 1)
+ for intr in interactions:
+ if isinstance(intr, Interaction1Site):
+ if intr.cpl != 0.0:
+ raise NotImplementedError(
+ "bond_update_bug currently supports nearest-neighbour two-site "
+ f"Hamiltonians only; found a one-site term on site {intr.site}"
+ )
+ continue
+ if isinstance(intr, Interaction2Site):
+ if intr.cpl == 0.0:
+ continue
+ if intr.terminal_site != intr.leading_site + 1:
+ raise NotImplementedError(
+ "bond_update_bug supports nearest-neighbour terms only; found a "
+ f"term coupling sites {intr.leading_site} and {intr.terminal_site}"
+ )
+ bond = intr.leading_site
+ term = bond_hamiltonian(intr)
+ generators[bond] = term if generators[bond] is None else generators[bond] + term
+ return generators
+
+
+def kernel_gate(h: Tensor, site_l_itag: str, site_r_itag: str) -> Tensor:
+ """Relabel a bond Hamiltonian into the local-KLS kernel's gate convention.
+
+ The KLS kernel applies a bare two-site term `g` to a two-site block
+ `theta` with `einsum('LRlr,aLRb->alrb', g, theta)`, then strips the trailing
+ ``*`` from the output physical itags. It therefore expects `g` with axes
+ `(ket_i, ket_j, bra_i, bra_j)`: the *ket* legs (`L`, `R`) carry the two site
+ itags and contract `theta`'s physical legs, while the *bra* legs (`l`, `r`)
+ carry the starred itags `('{si}*', '{sj}*')` and become the updated legs.
+
+ `bond_hamiltonian` returns the term with axes
+ `(bra_i, ket_i, bra_j, ket_j)`; this permutes to `(ket_i, ket_j, bra_i,
+ bra_j)` and retags the four legs with the two sites' physical itags so the
+ gate contracts the actual MPS physical indices.
+
+ Parameters
+ ----------
+ h:
+ Bond Hamiltonian with axes `(bra_i, ket_i, bra_j, ket_j)` (from
+ :func:`bond_hamiltonian`).
+ site_l_itag:
+ Physical itag of the left site `i` in the MPS.
+ site_r_itag:
+ Physical itag of the right site `i+1` in the MPS.
+
+ Returns
+ -------
+ Tensor
+ Complex gate with axes `(ket_i, ket_j, bra_i, bra_j)` and itags
+ `(site_l, site_r, '{site_l}*', '{site_r}*')`.
+ """
+ gate = permute(to_complex(h), [1, 3, 0, 2])
+ gate.retag(
+ [0, 1, 2, 3],
+ [site_l_itag, site_r_itag, f'{site_l_itag}*', f'{site_r_itag}*'],
+ )
+ return gate
diff --git a/src/alice/algorithm/bond_update_bug/bond_update_bug.py b/src/alice/algorithm/bond_update_bug/bond_update_bug.py
new file mode 100644
index 0000000..5df01a1
--- /dev/null
+++ b/src/alice/algorithm/bond_update_bug/bond_update_bug.py
@@ -0,0 +1,446 @@
+# Copyright (C) 2025-2026 Changkai Zhang.
+#
+# This file is part of Alice project.
+#
+# Alice is free software: you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published
+# by the Free Software Foundation, either version 3 of the License,
+# or (at your option) any later version.
+#
+# Alice is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with Alice. If not, see .
+
+
+"""Top-level bond_update_bug driver: options, summary, and entry-point function.
+
+The Basis-Update & Galerkin (BUG) integrator (Ceruti, Kusch & Lubich,
+arXiv:2304.05660) evolves an `MPS` under a nearest-neighbour Hamiltonian by
+odd/even Trotter sweeps of *local* two-site updates. Each bond update is the
+rank-adaptive K/L/S step: it augments the left frame from the evolved K factor,
+augments the right frame from the evolved L factor, evolves the small core S in
+the augmented bases (Galerkin), and truncates with an SVD. The local substeps
+exponentiate the *projected* effective Hamiltonian internally (Krylov `expv`) —
+no pre-formed gate is applied — so the step is the KLS update, exact at
+full rank. Bond Hamiltonians are reused directly from the AutoMPO interaction
+list, so any nearest-neighbour model and symmetry that `build_interaction`
+supports works unchanged.
+
+Typical usage:
+
+ from alice import build_interaction, init_mps
+ from alice.algorithm import bond_update_bug
+
+ interactions, spc, geo = build_interaction(cfg)
+ mps = init_mps(geo.L, spc, Op, config=[0, 1] * (geo.L // 2), target_qn=0)
+ opts = bond_update_bug.Options(dt=0.05, n_steps=20, order='strang', max_bond=64)
+ summary = bond_update_bug.run(mps, interactions, opts)
+ print(summary.bond_dims)
+"""
+
+from __future__ import annotations
+
+import logging
+from dataclasses import dataclass, field
+from typing import Dict, List, Optional
+
+from alice.network import MPS
+from alice.network.interaction import Interaction
+from alice.network.network import Network
+
+from ..interface import AlgorithmOptions, AlgorithmSummary
+from ._kernel import with_expv_backend, with_time_prefactor
+from ._kernel.local_solvers import LOCAL_SOLVERS
+from .bond import build_bond_generators, kernel_gate, to_complex
+from .scheme import parity_sweep
+
+logger = logging.getLogger(__name__)
+
+# Sentinel bond cap used when `Options.max_bond is None` (keep every singular
+# value at the post-S-step SVD, i.e. unlimited growth up to the local capacity).
+_UNLIMITED_BOND = 1 << 30
+
+
+# ---------------------------------------------------------------------------
+# Order alias resolution
+# ---------------------------------------------------------------------------
+
+_ORDER_ALIASES: Dict[str, str] = {
+ 'lie': 'lie',
+ 'first': 'lie',
+ '1': 'lie',
+ 'strang': 'strang',
+ 'second': 'strang',
+ '2': 'strang',
+}
+
+
+def _resolve_order(alias: str) -> str:
+ """Normalise a Trotter-order alias to its canonical name.
+
+ Parameters
+ ----------
+ alias:
+ User-provided order string.
+
+ Returns
+ -------
+ str
+ Canonical order name (`'lie'` or `'strang'`).
+
+ Raises
+ ------
+ ValueError
+ If `alias` is not a recognised order name.
+ """
+ canonical = _ORDER_ALIASES.get(alias.lower())
+ if canonical is None:
+ known = ', '.join(sorted(_ORDER_ALIASES))
+ raise ValueError(f"unknown Trotter order {alias!r}; recognised values are: {known}")
+ return canonical
+
+
+# ---------------------------------------------------------------------------
+# Options
+# ---------------------------------------------------------------------------
+
+@dataclass
+class Options(AlgorithmOptions):
+ """bond_update_bug run options.
+
+ All fields have sensible defaults so `Options()` is a valid minimal
+ configuration. Use `Options.from_toml` to load from an `[algorithm]` TOML
+ section, or `Options.load_toml` to read directly from a file.
+
+ Parameters
+ ----------
+ dt:
+ Time step. Interpreted as real time (evolution operator `exp(-i dt H)`)
+ unless `imaginary_time` is set.
+ n_steps:
+ Number of time steps to perform.
+ order:
+ Trotter order. Canonical values and their aliases:
+
+ - `'strang'` / `'second'` / `'2'`: symmetric second-order step
+ `U_even(dt/2) · U_odd(dt) · U_even(dt/2)`.
+ - `'lie'` / `'first'` / `'1'`: first-order step `U_even(dt) · U_odd(dt)`.
+ solver:
+ Local integrator for the K/L/S substeps — `'krylov'` (exact, default),
+ `'midpoint'` (explicit RK2), `'rk4'`, or `'trapezoid'` (A-stable
+ Crank–Nicolson). See
+ :mod:`alice.algorithm.bond_update_bug._kernel.local_solvers`.
+ solver_substeps:
+ Number of internal substeps for `'midpoint'`/`'rk4'`/`'trapezoid'` (local
+ error `O((dt/solver_substeps)^p)`; ignored by `'krylov'`).
+ max_bond:
+ Maximum bond dimension kept by the post-S-step SVD truncation. `None`
+ means no explicit cap (rank adapts up to the local capacity).
+ trunc_thresh:
+ Singular-value threshold of the post-S-step SVD. Each bond keeps only the
+ directions whose weight exceeds it, so the rank grows only as far as the
+ state's entanglement requires — the discarded-weight control of the
+ rank adaptation.
+ augment:
+ If `True` (default), the local KLS update may grow the bond basis from
+ the evolved K/L directions. If `False`, the bond dimension is held fixed
+ (parallel basis update without rank adaptation).
+ aug_krylov_depth:
+ Number of K/L Krylov directions stacked before the augmented basis is
+ extracted (`1` is the standard rank-adaptive BUG).
+ lanczos_tol:
+ Termination tolerance of the local Lanczos `expv` solves.
+ lanczos_maxiter:
+ Maximum Lanczos iterations per local substep.
+ imaginary_time:
+ If `True`, evolve with `exp(-dt H)` (imaginary time) instead of
+ `exp(-i dt H)`. Combined with `normalize`, this cools the state toward
+ the ground state.
+ normalize:
+ If `True` (default), renormalise the state after every step. Required
+ for imaginary-time evolution; harmless for real time (it only removes
+ the small norm leakage from truncation).
+ """
+
+ dt: float = 0.05
+ n_steps: int = 10
+ order: str = 'strang'
+ solver: str = 'krylov'
+ solver_substeps: int = 1
+ max_bond: Optional[int] = None
+ trunc_thresh: float = 1e-12
+ augment: bool = True
+ aug_krylov_depth: int = 1
+ lanczos_tol: float = 1e-15
+ lanczos_maxiter: int = 30
+ imaginary_time: bool = False
+ normalize: bool = True
+
+ def __post_init__(self) -> None:
+ self.order = _resolve_order(self.order)
+ # Validate eagerly so a bad solver name fails at construction.
+ if self.solver not in LOCAL_SOLVERS:
+ raise ValueError(
+ f"unknown local solver {self.solver!r}; recognised values are: "
+ f"{', '.join(LOCAL_SOLVERS)}")
+
+
+# ---------------------------------------------------------------------------
+# Summary
+# ---------------------------------------------------------------------------
+
+@dataclass
+class Summary(AlgorithmSummary):
+ """bond_update_bug output.
+
+ Attributes
+ ----------
+ state:
+ Evolved MPS after all steps (orthogonality center at site 0).
+ n_steps:
+ Number of steps performed.
+ times:
+ Cumulative evolution time recorded after each step (length `n_steps`).
+ norms:
+ State norm measured after each step *before* any renormalisation
+ (length `n_steps`). For real time these stay near 1; for imaginary time
+ they decay.
+ bond_dims:
+ Bond dimensions of `state` after the final step (length `L - 1`).
+ max_bond_dims:
+ Maximum *kept* bond dimension after each step (length `n_steps`).
+ aug_dims:
+ Maximum *proposed* (pre-truncation) augmented bond dimension over the
+ bonds of each step (length `n_steps`). This is the rank the K/L
+ augmentation reaches before the truncated S-step split; comparing it
+ with `max_bond_dims` shows how much rank growth the truncation discards.
+ disc_weights:
+ Maximum relative discarded weight over the bonds of each step (length
+ `n_steps`) — the fraction of bond weight the `trunc_thresh` S-step SVD
+ throws away. Near zero means the kept rank captures the state faithfully.
+ """
+
+ state: MPS
+ n_steps: int = 0
+ times: List[float] = field(default_factory=list)
+ norms: List[float] = field(default_factory=list)
+ bond_dims: List[int] = field(default_factory=list)
+ max_bond_dims: List[int] = field(default_factory=list)
+ aug_dims: List[int] = field(default_factory=list)
+ aug_k_dims: List[int] = field(default_factory=list)
+ aug_l_dims: List[int] = field(default_factory=list)
+ disc_weights: List[float] = field(default_factory=list)
+
+ def serialize(self) -> Dict:
+ """Serialize the summary to a plain dict compatible with `torch.save`.
+
+ Returns
+ -------
+ Dict
+ Serialized summary with keys `"version"`, `"n_steps"`, `"times"`,
+ `"norms"`, `"bond_dims"`, `"max_bond_dims"`, `"aug_dims"`,
+ `"disc_weights"`, and `"state"`.
+ """
+ return {
+ 'version': 1,
+ 'n_steps': self.n_steps,
+ 'times': self.times,
+ 'norms': self.norms,
+ 'bond_dims': self.bond_dims,
+ 'max_bond_dims': self.max_bond_dims,
+ 'aug_dims': self.aug_dims,
+ 'aug_k_dims': self.aug_k_dims,
+ 'aug_l_dims': self.aug_l_dims,
+ 'disc_weights': self.disc_weights,
+ 'state': self.state.serialize(),
+ }
+
+ @classmethod
+ def deserialize(cls, data: Dict, device: str = 'cpu') -> Summary:
+ """Reconstruct a `Summary` from a dict produced by `serialize`.
+
+ Parameters
+ ----------
+ data:
+ Dict previously returned by `serialize`.
+ device:
+ Device to place all MPS tensor blocks on. Defaults to `'cpu'`.
+
+ Returns
+ -------
+ Summary
+ Reconstructed summary with the MPS state placed on `device`.
+
+ Raises
+ ------
+ ValueError
+ If `data["version"]` is not `1`.
+ """
+ version = data.get('version', 1)
+ if version != 1:
+ raise ValueError(f"Unsupported Summary serialization version: {version!r}")
+ return cls(
+ state=Network.deserialize(data['state'], device=device),
+ n_steps=data['n_steps'],
+ times=data['times'],
+ norms=data['norms'],
+ bond_dims=data['bond_dims'],
+ max_bond_dims=data['max_bond_dims'],
+ aug_dims=data.get('aug_dims', []),
+ aug_k_dims=data.get('aug_k_dims', []),
+ aug_l_dims=data.get('aug_l_dims', []),
+ disc_weights=data.get('disc_weights', []),
+ )
+
+
+# ---------------------------------------------------------------------------
+# Top-level entry point
+# ---------------------------------------------------------------------------
+
+def run(mps: MPS, interactions: List[Interaction], opts: Optional[Options] = None) -> Summary:
+ """Evolve an MPS under a nearest-neighbour Hamiltonian with the bond_update_bug integrator.
+
+ Builds the per-bond Hamiltonian terms once from the AutoMPO interaction list,
+ then applies `opts.n_steps` odd/even Trotter steps of the K/L/S local
+ update. The state is canonicalised to `center = 0` before the first step and
+ returned with `center = 0`.
+
+ Parameters
+ ----------
+ mps:
+ Initial MPS state. Promoted to `complex128` and canonicalised in-place to
+ `center = 0` first.
+ interactions:
+ Interaction list from `build_interaction`. Every active term must be a
+ nearest-neighbour `Interaction2Site` (see
+ :func:`alice.algorithm.bond_update_bug.bond.build_bond_generators`).
+ opts:
+ Run options. Defaults to `Options()` if `None`.
+
+ Returns
+ -------
+ Summary
+ Evolved state, time/norm/bond-dimension history, and step count.
+
+ Raises
+ ------
+ ValueError
+ If `mps` has fewer than two sites.
+ """
+ if opts is None:
+ opts = Options()
+ if mps.L < 2:
+ raise ValueError(f"bond_update_bug evolution requires at least 2 sites, got L={mps.L}")
+
+ maxdim = opts.max_bond if opts.max_bond is not None else _UNLIMITED_BOND
+ # Real-time evolution uses exp(-i dt H); imaginary time uses exp(-dt H). The
+ # kernel multiplies its local timestep by this prefactor internally.
+ prefactor: complex = -1.0 if opts.imaginary_time else -1j
+
+ # Promote the state to complex128 so every local exponential shares the
+ # PyTorch backend dtype, then bring the center to site 0.
+ for site in range(mps.L):
+ mps[site] = to_complex(mps[site])
+ mps.canonical(0)
+
+ # Bare per-bond Hamiltonian terms, relabelled into the local-KLS kernel's
+ # gate convention against the MPS physical itags. Built once and reused for
+ # every sweep (the kernel exponentiates the projected term per substep).
+ generators = build_bond_generators(interactions, mps.L)
+ gates = [
+ None if h is None else kernel_gate(h, mps[b].itags[2], mps[b + 1].itags[2])
+ for b, h in enumerate(generators)
+ ]
+
+ def sweep(parity: str, tau: float):
+ return parity_sweep(
+ mps, gates, parity, tau, maxdim,
+ opts.augment, opts.aug_krylov_depth, opts.trunc_thresh,
+ opts.lanczos_tol, opts.lanczos_maxiter,
+ solver=opts.solver, solver_substeps=opts.solver_substeps,
+ )
+
+ times: List[float] = []
+ norms: List[float] = []
+ max_bond_dims: List[int] = []
+ aug_dims: List[int] = []
+ aug_k_dims: List[int] = []
+ aug_l_dims: List[int] = []
+ disc_weights: List[float] = []
+
+ n_active = sum(1 for h in generators if h is not None)
+ logger.info("─" * 60)
+ logger.info("Commencing: bond_update_bug Time Evolution".center(60))
+ logger.info("─" * 60)
+ logger.info("")
+ logger.info(" order : %s", opts.order)
+ logger.info(" local solver : %s (substeps %d)", opts.solver, opts.solver_substeps)
+ logger.info(" chain length : %d", mps.L)
+ logger.info(" active bonds : %d / %d", n_active, mps.L - 1)
+ logger.info(" time step : %g", opts.dt)
+ logger.info(" steps : %d", opts.n_steps)
+ logger.info(" evolution : %s", "imaginary" if opts.imaginary_time else "real")
+ logger.info(" max bond dim : %s", opts.max_bond if opts.max_bond is not None else 'unlimited')
+ logger.info(" augment : %s", opts.augment)
+ logger.info("")
+
+ w = len(str(opts.n_steps))
+ with with_time_prefactor(prefactor), with_expv_backend('native_hermitian_lanczos'):
+ for step in range(opts.n_steps):
+ if opts.order == 'strang':
+ # Symmetric Strang step: U_even(dt/2) · U_odd(dt) · U_even(dt/2).
+ results = [
+ sweep('even', 0.5 * opts.dt),
+ sweep('odd', opts.dt),
+ sweep('even', 0.5 * opts.dt),
+ ]
+ else:
+ # First-order Lie step: U_even(dt) · U_odd(dt).
+ results = [
+ sweep('even', opts.dt),
+ sweep('odd', opts.dt),
+ ]
+ aug_k = max(ak for ak, _, _ in results)
+ aug_l = max(al for _, al, _ in results)
+ augmented = max(aug_k, aug_l)
+ discarded = max(disc for _, _, disc in results)
+
+ norm = mps.norm()
+ if opts.normalize:
+ mps.normalize()
+
+ times.append((step + 1) * opts.dt)
+ norms.append(norm)
+ max_bond_dims.append(max(mps.bond_dims) if mps.bond_dims else 1)
+ aug_dims.append(augmented)
+ aug_k_dims.append(aug_k)
+ aug_l_dims.append(aug_l)
+ disc_weights.append(discarded)
+
+ logger.info(
+ "step %*d / %d: t = %g, norm = %.10f, kept bond = %d, aug(K,L) = (%d,%d), disc = %.2e",
+ w, step + 1, opts.n_steps, times[-1], norm, max_bond_dims[-1], aug_k, aug_l, discarded,
+ )
+
+ # Ensure the returned state has the center at site 0 for a well-defined norm.
+ if mps.center != 0:
+ mps.canonical(0)
+
+ logger.info("")
+
+ return Summary(
+ state=mps,
+ n_steps=opts.n_steps,
+ times=times,
+ norms=norms,
+ bond_dims=list(mps.bond_dims),
+ max_bond_dims=max_bond_dims,
+ aug_dims=aug_dims,
+ aug_k_dims=aug_k_dims,
+ aug_l_dims=aug_l_dims,
+ disc_weights=disc_weights,
+ )
diff --git a/src/alice/algorithm/bond_update_bug/scheme.py b/src/alice/algorithm/bond_update_bug/scheme.py
new file mode 100644
index 0000000..d681654
--- /dev/null
+++ b/src/alice/algorithm/bond_update_bug/scheme.py
@@ -0,0 +1,302 @@
+# Copyright (C) 2025-2026 Changkai Zhang.
+#
+# This file is part of Alice project.
+#
+# Alice is free software: you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published
+# by the Free Software Foundation, either version 3 of the License,
+# or (at your option) any later version.
+#
+# Alice is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with Alice. If not, see .
+
+
+"""Odd/even parity sweeps driving the KLS local bond update.
+
+The chain Hamiltonian splits into two commuting groups — bonds with an even
+left-site index (0, 2, 4, …) and bonds with an odd left-site index (1, 3, 5, …).
+Gates within one group act on disjoint site pairs, so a parity sweep applies
+them as an exact factor of the Trotter step; the splitting error lives only
+between the two groups. This is the odd/even BUG sweep used for the domain-wall
+XX chain.
+
+Each bond update is the Ceruti–Kusch–Lubich K/L/S step from
+:mod:`alice.algorithm.bond_update_bug._kernel` (Basis-Update & Galerkin).
+This module is the thin Alice adapter: it brings the orthogonality center onto
+the active bond, takes a canonical two-site snapshot of the Alice `MPS`, calls
+the vendored kernel, and writes the updated cores back. The kernel works in the
+`(link_l, site, link_r)` tensor layout; Alice stores `(left, right, phys)`, so
+the snapshot and writeback transpose between the two.
+"""
+
+from __future__ import annotations
+
+from typing import Callable, Dict, List, Optional, Tuple
+
+import torch
+from nicole import Tensor, permute
+
+from alice.network import MPS
+
+from ._kernel import (
+ Ix,
+ _kls_local_bond_candidate,
+ lq,
+ qr,
+ tcontract,
+ to_dense,
+)
+
+
+def _discarded_weight(s_new: Tensor, keep: int) -> float:
+ """Relative Frobenius weight discarded when the S-step core is cut to `keep`.
+
+ The kernel evolves the (small) two-site core `s_new` and then keeps `keep`
+ singular values. This recomputes the full singular spectrum of `s_new` and
+ returns `sqrt(Σ_{i≥keep} σ_i² / Σ_i σ_i²)` — the fraction of the bond's weight
+ the truncation throws away, the standard MPS discarded-weight diagnostic.
+ """
+ matrix = to_dense(s_new, list(s_new.itags))
+ svals = torch.linalg.svdvals(matrix.reshape(matrix.shape[0], -1))
+ total = float((svals ** 2).sum())
+ if total == 0.0 or keep >= svals.numel():
+ return 0.0
+ tail = float((svals[keep:] ** 2).sum())
+ return (tail / total) ** 0.5
+
+
+def _ix(tensor: Tensor, axis: int) -> Ix:
+ """Wrap one leg of a Nicole tensor as a kernel `Ix` handle."""
+ index = tensor.indices[axis]
+ return Ix(tensor.itags[axis], int(index.dim), index.direction, index.sectors, index.group)
+
+
+def _to_kernel_layout(site: Tensor) -> Tensor:
+ """Transpose an Alice MPS tensor `(left, right, phys)` to `(left, phys, right)`."""
+ return permute(site, [0, 2, 1])
+
+
+def _to_mps_layout(core: Tensor) -> Tensor:
+ """Transpose a kernel core `(left, phys, right)` back to Alice `(left, right, phys)`."""
+ return permute(core, [0, 2, 1])
+
+
+def bond_snapshot(mps: MPS, i: int) -> Dict[str, object]:
+ """Take the canonical two-site snapshot the KLS kernel consumes at bond *(i, i+1)*.
+
+ `mps` must already have its orthogonality center on site *i*. A QR of site
+ *i* gives the left isometry `U0` and an LQ of site *i+1* gives the right
+ isometry `V0`; their inner factors contract to the bond center `S0`. These
+ are exact, truncation-free moves, so the subsequent KLS update — and only it
+ — is responsible for the rank adaptation.
+
+ Parameters
+ ----------
+ mps:
+ State with `center == i`.
+ i:
+ Left site of the bond.
+
+ Returns
+ -------
+ dict
+ Snapshot mapping consumed by `LocalBondFrame.from_mapping`: the five leg
+ handles (`link_l`, `site_l`, `link_mid`, `site_r`, `link_r`), the
+ canonical factors (`U0_tens`, `V0_tens`, `S0_tens`), and the middle bond
+ handles on each side (`canon_u0`, `canon_v0`).
+ """
+ left = _to_kernel_layout(mps[i]) # (link_l, site_l, link_mid)
+ right = _to_kernel_layout(mps[i + 1]) # (link_mid, site_r, link_r)
+
+ link_l = _ix(left, 0)
+ site_l = _ix(left, 1)
+ link_mid = _ix(left, 2)
+ site_r = _ix(right, 1)
+ link_r = _ix(right, 2)
+
+ U0_tens, s_left, canon_u0 = qr(left, [link_l, site_l], tag=link_mid.itag)
+ s_right, V0_tens, canon_v0 = lq(right, [site_r, link_r], tag=link_mid.itag)
+ S0_tens = tcontract(s_left, s_right)
+
+ return {
+ 'link_l': link_l,
+ 'site_l': site_l,
+ 'link_mid': link_mid,
+ 'site_r': site_r,
+ 'link_r': link_r,
+ 'U0_tens': U0_tens,
+ 'V0_tens': V0_tens,
+ 'S0_tens': S0_tens,
+ 'canon_u0': canon_u0,
+ 'canon_v0': canon_v0,
+ }
+
+
+def kls_bond(
+ mps: MPS,
+ i: int,
+ gate: Tensor,
+ tau: float,
+ maxdim: int,
+ augment: bool,
+ aug_krylov_depth: int,
+ trunc_thresh: float,
+ lanczos_tol: float,
+ lanczos_maxiter: int,
+ candidate_fn: Callable = _kls_local_bond_candidate,
+ solver: str = 'krylov',
+ solver_substeps: int = 1,
+) -> Tuple[int, int, float]:
+ """Apply one KLS update to sites *(i, i+1)* of `mps`, in place.
+
+ Moves the orthogonality center onto site *i* (truncation-free), snapshots the
+ bond, runs the vendored K/L/S local update for time `tau` (the active
+ evolution prefactor — `-1j` for real time, `-1` for imaginary — is applied by
+ the kernel), and writes the two updated cores back. After the call
+ `mps.center == i + 1`.
+
+ Parameters
+ ----------
+ mps:
+ State to update in place.
+ i:
+ Left site of the bond.
+ gate:
+ Bare two-site bond Hamiltonian in the kernel convention (see
+ :func:`alice.algorithm.bond_update_bug.bond.kernel_gate`).
+ tau:
+ Real time advanced by this local step.
+ maxdim:
+ Bond-dimension cap kept by the post-S-step SVD truncation.
+ augment, aug_krylov_depth, lanczos_tol, lanczos_maxiter:
+ KLS controls forwarded to the kernel.
+
+ trunc_thresh:
+ Singular-value threshold for the post-S-step SVD: the bond keeps only the
+ directions whose weight exceeds it, so the rank grows only as far as the
+ entanglement of the state requires (the rank-adaptive truncation).
+
+ Returns
+ -------
+ int
+ Proposed augmented **K** bond dimension (old rank + new K directions) at
+ this bond, before the truncated split.
+ int
+ Proposed augmented **L** bond dimension (old rank + new L directions).
+ float
+ Relative weight discarded by this bond's S-step truncation.
+ """
+ mps.canonical(i, trunc=None)
+ bond_data = bond_snapshot(mps, i)
+ old_rank = int(bond_data['link_mid'].dim)
+
+ candidate = candidate_fn(
+ bond_data,
+ gate=gate,
+ dt=tau,
+ maxdim=maxdim,
+ augment=augment,
+ aug_krylov_depth=aug_krylov_depth,
+ trunc_thresh=trunc_thresh,
+ lanczos_tol=lanczos_tol,
+ lanczos_maxiter=lanczos_maxiter,
+ solver=solver,
+ solver_substeps=solver_substeps,
+ )
+
+ mps[i] = _to_mps_layout(candidate['left_core'])
+ mps[i + 1] = _to_mps_layout(candidate['right_core'])
+ mps._center = i + 1
+
+ aug_k = old_rank + int(candidate['n_new_k'])
+ aug_l = old_rank + int(candidate['n_new_l'])
+ discarded = _discarded_weight(candidate['S_new'], int(candidate['keep']))
+ return aug_k, aug_l, discarded
+
+
+def parity_bonds(length: int, parity: str) -> List[int]:
+ """Return the left-site indices of all bonds in one commuting group.
+
+ Parameters
+ ----------
+ length:
+ Number of sites `L`.
+ parity:
+ `'even'` for bonds with an even left-site index (0, 2, 4, …) or `'odd'`
+ for bonds with an odd left-site index (1, 3, 5, …).
+
+ Returns
+ -------
+ list of int
+ Left-site indices of the bonds in the requested group.
+
+ Raises
+ ------
+ ValueError
+ If `parity` is not `'even'` or `'odd'`.
+ """
+ if parity == 'even':
+ return list(range(0, length - 1, 2))
+ if parity == 'odd':
+ return list(range(1, length - 1, 2))
+ raise ValueError(f"parity must be 'even' or 'odd', got {parity!r}")
+
+
+def parity_sweep(
+ mps: MPS,
+ gates: List[Optional[Tensor]],
+ parity: str,
+ tau: float,
+ maxdim: int,
+ augment: bool,
+ aug_krylov_depth: int,
+ trunc_thresh: float,
+ lanczos_tol: float,
+ lanczos_maxiter: int,
+ candidate_fn: Callable = _kls_local_bond_candidate,
+ solver: str = 'krylov',
+ solver_substeps: int = 1,
+) -> Tuple[int, int, float]:
+ """Apply every bond gate of one commuting group to `mps`, in place.
+
+ Bonds of the chosen parity act on disjoint site pairs, so the group is an
+ exact factor of the Trotter step. Bonds whose gate is `None` (no Hamiltonian
+ term) are skipped.
+
+ Parameters
+ ----------
+ mps:
+ State to update in place.
+ gates:
+ Per-bond kernel gates of length `L - 1`; entry *b* acts on bond *(b, b+1)*.
+ parity:
+ `'even'` or `'odd'` — selects the commuting bond group.
+ tau:
+ Real time advanced by each local KLS step in this group.
+ maxdim, augment, aug_krylov_depth, trunc_thresh, lanczos_tol, lanczos_maxiter:
+ KLS controls forwarded to each bond update.
+
+ Returns
+ -------
+ int
+ Largest proposed augmented bond dimension over the bonds of this group
+ (0 if the group has no active bonds).
+ float
+ Largest relative discarded weight over the bonds of this group.
+ """
+ aug_k = aug_l = 0
+ discarded = 0.0
+ for i in parity_bonds(mps.L, parity):
+ if gates[i] is not None:
+ ak, al, disc = kls_bond(mps, i, gates[i], tau, maxdim, augment,
+ aug_krylov_depth, trunc_thresh, lanczos_tol, lanczos_maxiter,
+ candidate_fn, solver, solver_substeps)
+ aug_k = max(aug_k, ak)
+ aug_l = max(aug_l, al)
+ discarded = max(discarded, disc)
+ return aug_k, aug_l, discarded
diff --git a/src/alice/algorithm/tdvp2/__init__.py b/src/alice/algorithm/tdvp2/__init__.py
new file mode 100644
index 0000000..6168890
--- /dev/null
+++ b/src/alice/algorithm/tdvp2/__init__.py
@@ -0,0 +1,44 @@
+# Copyright (C) 2025-2026 Changkai Zhang.
+#
+# This file is part of Alice project.
+#
+# Alice is free software: you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published
+# by the Free Software Foundation, either version 3 of the License,
+# or (at your option) any later version.
+#
+# Alice is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with Alice. If not, see .
+# Author of code: Madhav Menon.
+
+
+"""Two-site TDVP algorithm package.
+
+A 2-site Time-Dependent Variational Principle integrator for an `MPS` evolving
+under a Hamiltonian `MPO` (Haegeman et al., arXiv:1408.5056). A forward half-sweep
+evolves each two-site block forward and the carried one-site tensor backward
+(inverse-free backward correction); a reverse half-sweep mirrors it; a symmetric
+step composes the two halves for second-order accuracy. The bond dimension adapts
+through the per-bond SVD truncation.
+
+Reuses the DMRG environment machinery and effective-Hamiltonian contractions
+(`alice.algorithm.dmrg`) and the Hermitian Krylov exponential vendored with the
+two-site BUG kernel. Public API:
+
+- `Options` — run options (loadable from TOML).
+- `Summary` — output dataclass.
+- `run` — top-level entry point ``run(mps, mpo, opts)``.
+"""
+
+from .tdvp2 import Options, Summary, run
+
+__all__ = [
+ 'Options',
+ 'Summary',
+ 'run',
+]
diff --git a/src/alice/algorithm/tdvp2/_krylov.py b/src/alice/algorithm/tdvp2/_krylov.py
new file mode 100644
index 0000000..67d74cf
--- /dev/null
+++ b/src/alice/algorithm/tdvp2/_krylov.py
@@ -0,0 +1,259 @@
+# Copyright (C) 2025-2026 Changkai Zhang.
+#
+# This file is part of Alice project.
+#
+# Alice is free software: you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published
+# by the Free Software Foundation, either version 3 of the License,
+# or (at your option) any later version.
+#
+# Alice is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with Alice. If not, see .
+# Author of code: Madhav Menon.
+
+
+"""Self-contained Krylov exponential for the 2-site TDVP local substeps.
+
+TDVP advances a site or bond tensor by ``exp(prefactor * dt * H_eff)`` where the
+effective Hamiltonian ``H_eff`` is Hermitian and available only as a matrix-free
+action on a Nicole `Tensor` (the contraction of the MPO environments with the MPO
+site tensors). This module provides a Hermitian tensor Lanczos exponential for
+exactly that situation, plus the evolution-prefactor context manager (``-1j`` for
+real time, ``-1`` for imaginary time).
+
+Keeping these helpers inside the ``tdvp2`` package makes the integrator depend
+only on the shared `network` layer and the DMRG effective-Hamiltonian
+contractions — no other algorithm package — so it can be reviewed and merged on
+its own.
+"""
+
+from __future__ import annotations
+
+import contextlib
+from collections.abc import Callable
+from typing import Iterable
+
+import torch
+from nicole import Tensor, conj as _nconj, einsum as _neinsum
+
+_ACTIVE_PREFACTOR = [complex(0.0, -1.0)]
+
+
+@contextlib.contextmanager
+def with_time_prefactor(c: complex):
+ """Temporarily set the global evolution prefactor.
+
+ Parameters
+ ----------
+ c:
+ New complex prefactor (``-1j`` for real time, ``-1`` for imaginary time).
+
+ Returns
+ -------
+ contextlib._GeneratorContextManager
+ A context manager that restores the previous prefactor on exit.
+ """
+ previous = _ACTIVE_PREFACTOR[0]
+ _ACTIVE_PREFACTOR[0] = complex(c)
+ try:
+ yield
+ finally:
+ _ACTIVE_PREFACTOR[0] = previous
+
+
+def active_time_prefactor() -> complex:
+ """Return the currently active evolution prefactor.
+
+ Returns
+ -------
+ complex
+ The active complex prefactor set by :func:`with_time_prefactor`
+ (default ``-1j``).
+ """
+ return _ACTIVE_PREFACTOR[0]
+
+
+def to_complex(tensor: Tensor) -> Tensor:
+ """Return a copy of ``tensor`` with every block cast to ``complex128``.
+
+ Real-time evolution exponentiates the effective Hamiltonian, so the state, the
+ MPO, and the environment boundary blocks must all share the ``complex128``
+ backend dtype.
+
+ Parameters
+ ----------
+ tensor:
+ Nicole tensor with real or complex blocks.
+
+ Returns
+ -------
+ Tensor
+ Tensor with identical indices and itags but ``complex128`` block data.
+ """
+ new_intw = None
+ if tensor.intw is not None:
+ new_intw = {
+ key: bridge.to(tensor.device, dtype=torch.complex128)
+ for key, bridge in tensor.intw.items()
+ }
+ return Tensor(
+ indices=tensor.indices,
+ itags=tensor.itags,
+ data={key: block.to(torch.complex128) for key, block in tensor.data.items()},
+ intw=new_intw,
+ dtype=torch.complex128,
+ )
+
+
+def _tensor_inner(a: Tensor, b: Tensor) -> complex:
+ """Return the canonical inner product ```` for two same-shape tensors.
+
+ Parameters
+ ----------
+ a:
+ Left (bra) tensor.
+ b:
+ Right (ket) tensor with the same index structure as ``a``.
+
+ Returns
+ -------
+ complex
+ The scalar inner product ``sum(conj(a) * b)``.
+ """
+ equation = "".join(chr(97 + axis) for axis in range(len(a.itags)))
+ return _neinsum(f"{equation},{equation}->", _nconj(a), b).item()
+
+
+def _tridiagonal_exp_first_column(
+ alpha: Iterable[float],
+ beta: Iterable[float],
+ dt: complex,
+) -> torch.Tensor:
+ """Return ``exp(dt * T) e_1`` for the Hermitian tridiagonal Lanczos matrix ``T``.
+
+ Parameters
+ ----------
+ alpha:
+ Diagonal entries of the tridiagonal matrix.
+ beta:
+ Off-diagonal entries (length ``len(alpha) - 1``).
+ dt:
+ Scalar prefactor in the exponential.
+
+ Returns
+ -------
+ torch.Tensor
+ The first column of ``exp(dt * T)``, as a complex vector of length
+ ``len(alpha)``.
+ """
+ alpha_t = torch.as_tensor(tuple(alpha), dtype=torch.float64)
+ beta_t = torch.as_tensor(tuple(beta), dtype=torch.float64)
+ if alpha_t.numel() == 0:
+ return torch.empty((0,), dtype=torch.complex128)
+ tridiagonal = torch.diag(alpha_t)
+ if beta_t.numel() > 0:
+ tridiagonal = tridiagonal + torch.diag(beta_t, 1) + torch.diag(beta_t, -1)
+ evals, evecs = torch.linalg.eigh(tridiagonal)
+ evecs_c = evecs.to(torch.complex128)
+ weights = torch.exp(dt * evals.to(torch.complex128)) * evecs_c[0, :]
+ return evecs_c @ weights
+
+
+# Opt-in Krylov-depth instrumentation (off by default => zero overhead). When
+# enabled, every tensor_lanczos_expv call appends its Krylov dimension (number of
+# matrix-free H applications) to KRYLOV_LOG, for the N_Krylov diagnostic.
+KRYLOV_LOG: list[int] = []
+_KRYLOV_RECORD = False
+
+
+def enable_krylov_log() -> None:
+ global _KRYLOV_RECORD
+ _KRYLOV_RECORD = True
+ KRYLOV_LOG.clear()
+
+
+def disable_krylov_log() -> None:
+ global _KRYLOV_RECORD
+ _KRYLOV_RECORD = False
+
+
+def get_krylov_log() -> list[int]:
+ return list(KRYLOV_LOG)
+
+
+def tensor_lanczos_expv(
+ apply: Callable[[Tensor], Tensor],
+ dt: complex,
+ x: Tensor,
+ *,
+ maxiter: int = 30,
+ tol: float = 1e-13,
+) -> Tensor:
+ """Return ``exp(dt * H) @ x`` for a Hermitian Nicole-tensor action ``apply``.
+
+ Builds an orthonormal Krylov basis of Nicole tensors via the Hermitian Lanczos
+ three-term recurrence, exponentiates the small tridiagonal projection, and
+ recombines the basis. Everything stays in the symmetry-blocked Nicole
+ representation; the operator is never materialised as a dense matrix.
+
+ Parameters
+ ----------
+ apply:
+ Matrix-free Hermitian action ``H`` on a Nicole tensor, returning a tensor
+ with the same index structure as its input.
+ dt:
+ Scalar prefactor in the exponential (already including any evolution
+ prefactor such as ``-1j``).
+ x:
+ Input Nicole tensor.
+ maxiter:
+ Maximum Krylov dimension (number of Lanczos steps).
+ tol:
+ Off-diagonal threshold at which the Lanczos recurrence terminates early.
+
+ Returns
+ -------
+ Tensor
+ The evolved tensor ``exp(dt * H) @ x`` with the same index structure as
+ ``x``.
+ """
+ beta0 = x.norm()
+ if float(abs(beta0)) == 0.0:
+ if _KRYLOV_RECORD:
+ KRYLOV_LOG.append(0)
+ return x
+
+ v = (1.0 / beta0) * x
+ basis = [v]
+ alpha: list[float] = []
+ betas: list[float] = []
+
+ w = apply(v)
+ a = _tensor_inner(v, w).real
+ alpha.append(a)
+ w = w + (-a) * v
+
+ for _ in range(1, maxiter):
+ b = w.norm()
+ if float(b) < tol:
+ break
+ betas.append(float(b))
+ v = (1.0 / b) * w
+ basis.append(v)
+ w = apply(v)
+ a = _tensor_inner(v, w).real
+ alpha.append(a)
+ w = w + (-a) * v + (-b) * basis[-2]
+
+ if _KRYLOV_RECORD:
+ KRYLOV_LOG.append(len(alpha))
+ coeff = _tridiagonal_exp_first_column(alpha, betas, dt) * beta0
+ evolved = coeff[0] * basis[0]
+ for idx in range(1, len(alpha)):
+ evolved = evolved + coeff[idx] * basis[idx]
+ return evolved
diff --git a/src/alice/algorithm/tdvp2/local.py b/src/alice/algorithm/tdvp2/local.py
new file mode 100644
index 0000000..46b3e57
--- /dev/null
+++ b/src/alice/algorithm/tdvp2/local.py
@@ -0,0 +1,122 @@
+# Copyright (C) 2025-2026 Changkai Zhang.
+#
+# This file is part of Alice project.
+#
+# Alice is free software: you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published
+# by the Free Software Foundation, either version 3 of the License,
+# or (at your option) any later version.
+#
+# Alice is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with Alice. If not, see .
+# Author of code: Madhav Menon.
+
+
+"""Local real-time substeps for the 2-site TDVP integrator.
+
+2-site TDVP advances the orthogonality window by exponentiating the *effective
+Hamiltonian* — the MPS bond tensor evolved under ``exp(prefactor·dt·H_eff)`` with
+the left/right MPO environments held fixed — rather than applying a pre-formed
+gate. The two effective Hamiltonians are exactly the DMRG ones, so this module
+reuses the DMRG contractions verbatim:
+
+- the 2-site action ``H_eff^{(2)}`` is :func:`alice.algorithm.dmrg.scheme_2s.matvec_2s`
+ (``E_left · W_i · W_{i+1} · E_right`` applied to Θ), and
+- the 1-site action ``H_eff^{(1)}`` is :func:`alice.algorithm.dmrg.scheme_1s.matvec`,
+
+each fed to the Hermitian tensor Lanczos exponential in this package's
+:mod:`alice.algorithm.tdvp2._krylov` (``tensor_lanczos_expv``). H_eff is Hermitian,
+so the Lanczos path is the right one; the evolution prefactor (``-1j`` real time,
+``-1`` imaginary) is the active prefactor of that module.
+
+The forward sweep evolves each 2-site block forward by ``dt`` then evolves the
+carried 1-site tensor *backward* by ``dt`` (the inverse-free single-site
+correction that prevents double counting the shared bond); the reverse sweep
+mirrors it. A symmetric (Strang) step composes a forward half-sweep and a reverse
+half-sweep.
+"""
+
+from __future__ import annotations
+
+from functools import partial
+
+from nicole import Tensor
+
+from ..dmrg.scheme_1s import matvec as _matvec_1s
+from ..dmrg.scheme_2s import matvec_2s as _matvec_2s
+from ._krylov import active_time_prefactor, tensor_lanczos_expv
+
+
+def evolve_two_site(
+ theta: Tensor,
+ W_i: Tensor,
+ W_i1: Tensor,
+ E_left: Tensor,
+ E_right: Tensor,
+ dt: complex,
+ *,
+ lanczos_tol: float,
+ lanczos_maxiter: int,
+) -> Tensor:
+ """Return ``exp(prefactor·dt·H_eff^{(2)})|Θ⟩`` for the 2-site bond tensor Θ.
+
+ Parameters
+ ----------
+ theta:
+ Bond tensor with axes ``(ket_left, ket_right, phys_ket_i, phys_ket_{i+1})``.
+ W_i, W_i1:
+ MPO tensors at sites ``i`` and ``i+1``.
+ E_left, E_right:
+ Left/right MPO environments bracketing the two-site window.
+ dt:
+ Real time advanced by this substep (multiplied by the active evolution
+ prefactor internally).
+ lanczos_tol, lanczos_maxiter:
+ Lanczos termination tolerance and maximum Krylov dimension.
+ """
+ mv = partial(_matvec_2s, W_i=W_i, W_i1=W_i1, E_left=E_left, E_right=E_right)
+ return tensor_lanczos_expv(
+ mv, active_time_prefactor() * dt, theta,
+ maxiter=lanczos_maxiter, tol=lanczos_tol,
+ )
+
+
+def evolve_one_site(
+ site: Tensor,
+ W: Tensor,
+ E_left: Tensor,
+ E_right: Tensor,
+ dt: complex,
+ *,
+ lanczos_tol: float,
+ lanczos_maxiter: int,
+) -> Tensor:
+ """Return ``exp(prefactor·dt·H_eff^{(1)})|M⟩`` for the 1-site tensor M.
+
+ Used with a *negative* ``dt`` for the TDVP backward correction on the carried
+ bond tensor between two 2-site updates.
+
+ Parameters
+ ----------
+ site:
+ Center site tensor with axes ``(ket_left, ket_right, phys_ket)``.
+ W:
+ MPO tensor at that site.
+ E_left, E_right:
+ Left/right MPO environments bracketing the site.
+ dt:
+ Real time advanced by this substep (multiplied by the active evolution
+ prefactor internally). The caller passes ``-tau`` for the backward step.
+ lanczos_tol, lanczos_maxiter:
+ Lanczos termination tolerance and maximum Krylov dimension.
+ """
+ mv = partial(_matvec_1s, W=W, E_left=E_left, E_right=E_right)
+ return tensor_lanczos_expv(
+ mv, active_time_prefactor() * dt, site,
+ maxiter=lanczos_maxiter, tol=lanczos_tol,
+ )
diff --git a/src/alice/algorithm/tdvp2/sweep.py b/src/alice/algorithm/tdvp2/sweep.py
new file mode 100644
index 0000000..12f2a21
--- /dev/null
+++ b/src/alice/algorithm/tdvp2/sweep.py
@@ -0,0 +1,152 @@
+# Copyright (C) 2025-2026 Changkai Zhang.
+#
+# This file is part of Alice project.
+#
+# Alice is free software: you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published
+# by the Free Software Foundation, either version 3 of the License,
+# or (at your option) any later version.
+#
+# Alice is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with Alice. If not, see .
+# Author of code: Madhav Menon.
+
+
+"""Forward / reverse half-sweeps for the 2-site TDVP integrator.
+
+These sweeps reuse the DMRG 2-site machinery wholesale — the `Environment` blocks
+and their transfer-matrix updates (`step_left_env` / `step_right_env`), the bond
+contraction `build_bulk`, and the truncating SVD splits `split_forward` /
+`split_backward`. The only differences from the DMRG sweep are:
+
+- the local update is a *real-time evolution* of the effective Hamiltonian
+ (`evolve_two_site`) rather than a Davidson eigensolve, and
+- after each forward 2-site step the carried one-site tensor is evolved
+ *backward* in time (`evolve_one_site` with ``-tau``) — the inverse-free TDVP
+ backward correction that removes the double counting of the shared bond.
+
+A forward half-sweep advances every bond left-to-right, leaving the orthogonality
+center at site ``L-1``; the reverse half-sweep mirrors it back to site ``0``. A
+symmetric (Strang) TDVP step is ``forward(dt/2)`` followed by ``reverse(dt/2)``.
+"""
+
+from __future__ import annotations
+
+from typing import Optional
+
+from alice.network import MPS, MPO
+
+from ..dmrg.environ import step_left_env, step_right_env
+from ..dmrg.scheme_2s import build_bulk, split_backward, split_forward
+from .local import evolve_one_site, evolve_two_site
+
+
+def _trunc_dict(maxdim: Optional[int], cutoff: float) -> Optional[dict]:
+ """Assemble the SVD truncation dict consumed by `split_forward`/`split_backward`."""
+ trunc: dict = {'thresh': max(float(cutoff), 0.0)}
+ if maxdim is not None:
+ trunc['nkeep'] = int(maxdim)
+ return trunc
+
+
+def forward_sweep(
+ mps: MPS,
+ mpo: MPO,
+ env_left,
+ env_right,
+ tau: float,
+ *,
+ maxdim: Optional[int],
+ cutoff: float,
+ lanczos_tol: float,
+ lanczos_maxiter: int,
+) -> None:
+ """Left-to-right 2-site TDVP half-sweep advancing the state by time ``tau``.
+
+ Requires ``mps.center == 0`` and every ``env_right`` block populated. After
+ the call ``mps.center == mps.L - 1``.
+
+ For each bond ``(i, i+1)``: contract the two cores, evolve forward by ``tau``,
+ SVD-split (truncating to ``maxdim``/``cutoff``) leaving ``mps[i]`` left-isometric
+ and the singular values carried right, advance the left environment, and — for
+ every bond except the last — evolve the carried one-site tensor backward by
+ ``tau`` before it is absorbed into the next bond.
+ """
+ L = mps.L
+ trunc = _trunc_dict(maxdim, cutoff)
+
+ for i in range(0, L - 1):
+ E_left = env_left.fetch(i)
+ E_right = env_right.fetch(i + 1)
+
+ theta = build_bulk(mps[i], mps[i + 1])
+ theta = evolve_two_site(theta, mpo[i], mpo[i + 1], E_left, E_right, tau,
+ lanczos_tol=lanczos_tol, lanczos_maxiter=lanczos_maxiter)
+
+ itag = mps._bond_itag(i + 1)
+ mps[i], carry = split_forward(theta, itag, trunc) # mps[i] left-iso, carry = S·V
+ mps._center = i + 1
+
+ if i == L - 1 - 1:
+ mps[i + 1] = carry
+ continue
+
+ # Advance the left environment with the freshly fixed left-isometric mps[i],
+ # then evolve the carried bond tensor backward in time on site i+1.
+ env_left[i + 1] = step_left_env(E_left, mps[i], mpo[i])
+ mps[i + 1] = evolve_one_site(
+ carry, mpo[i + 1], env_left[i + 1], E_right, -tau,
+ lanczos_tol=lanczos_tol, lanczos_maxiter=lanczos_maxiter,
+ )
+
+
+def reverse_sweep(
+ mps: MPS,
+ mpo: MPO,
+ env_left,
+ env_right,
+ tau: float,
+ *,
+ maxdim: Optional[int],
+ cutoff: float,
+ lanczos_tol: float,
+ lanczos_maxiter: int,
+) -> None:
+ """Right-to-left 2-site TDVP half-sweep advancing the state by time ``tau``.
+
+ Requires ``mps.center == mps.L - 1`` and every ``env_left`` block populated.
+ After the call ``mps.center == 0``. Mirror of :func:`forward_sweep`: the SVD
+ leaves ``mps[i+1]`` right-isometric and carries the singular values left, and
+ the carried one-site tensor is evolved backward by ``tau`` on site ``i``.
+ """
+ L = mps.L
+ trunc = _trunc_dict(maxdim, cutoff)
+
+ for i in range(L - 2, -1, -1):
+ E_left = env_left.fetch(i)
+ E_right = env_right.fetch(i + 1)
+
+ theta = build_bulk(mps[i], mps[i + 1])
+ theta = evolve_two_site(theta, mpo[i], mpo[i + 1], E_left, E_right, tau,
+ lanczos_tol=lanczos_tol, lanczos_maxiter=lanczos_maxiter)
+
+ itag = mps._bond_itag(i + 1)
+ carry, mps[i + 1] = split_backward(theta, itag, trunc) # mps[i+1] right-iso, carry = U·S
+ mps._center = i
+
+ if i == 0:
+ mps[i] = carry
+ continue
+
+ # Advance the right environment with the freshly fixed right-isometric
+ # mps[i+1], then evolve the carried bond tensor backward on site i.
+ env_right[i] = step_right_env(E_right, mps[i + 1], mpo[i + 1])
+ mps[i] = evolve_one_site(
+ carry, mpo[i], E_left, env_right[i], -tau,
+ lanczos_tol=lanczos_tol, lanczos_maxiter=lanczos_maxiter,
+ )
diff --git a/src/alice/algorithm/tdvp2/tdvp2.py b/src/alice/algorithm/tdvp2/tdvp2.py
new file mode 100644
index 0000000..da088e3
--- /dev/null
+++ b/src/alice/algorithm/tdvp2/tdvp2.py
@@ -0,0 +1,305 @@
+# Copyright (C) 2025-2026 Changkai Zhang.
+#
+# This file is part of Alice project.
+#
+# Alice is free software: you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published
+# by the Free Software Foundation, either version 3 of the License,
+# or (at your option) any later version.
+#
+# Alice is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with Alice. If not, see .
+# Author of code: Madhav Menon.
+
+
+"""Top-level 2-site TDVP driver: options, summary, and entry-point function.
+
+Two-site Time-Dependent Variational Principle (TDVP) integrator on an `MPS`,
+following the Haegeman et al. projector-splitting scheme (arXiv:1408.5056) with a
+2-site update so the bond dimension can adapt. It is the Alice counterpart of the
+reference Julia `tdvp2_step!` (`../../../../src/TDVP/tdvp2_sweep.jl`): a forward
+half-sweep evolves each 2-site block forward by ``dt`` and the carried one-site
+tensor backward by ``dt`` (inverse-free backward correction), a reverse half-sweep
+mirrors it, and a symmetric step composes ``forward(dt/2)`` + ``reverse(dt/2)`` for
+second-order accuracy.
+
+Unlike the BUG integrators (which apply bare two-site gates), TDVP exponentiates
+the full *effective Hamiltonian* with the left/right MPO environments, so it takes
+a Hamiltonian `MPO` (from `build_hamiltonian`) — exactly like `alice.algorithm.dmrg`
+— and reuses the DMRG environment machinery and 2-site/1-site contractions.
+
+Typical usage::
+
+ from alice import build_interaction, build_hamiltonian, init_mps
+ from alice.algorithm import tdvp2
+
+ interactions, spc, geo = build_interaction(cfg)
+ mpo = build_hamiltonian(interactions, geo.L, spc)
+ mps = init_mps(geo.L, spc, Op, config=[0, 1] * (geo.L // 2), target_qn=0)
+ opts = tdvp2.Options(dt=0.05, n_steps=20, max_bond=64)
+ summary = tdvp2.run(mps, mpo, opts)
+"""
+
+from __future__ import annotations
+
+import logging
+from dataclasses import dataclass, field
+from typing import Dict, List, Optional
+
+from alice.network import MPS, MPO
+from alice.network.network import Network
+
+from ..interface import AlgorithmOptions, AlgorithmSummary
+from ..dmrg.environ import (
+ Environment,
+ left_env_boundary,
+ right_env_boundary,
+ step_left_env,
+ step_right_env,
+)
+from ._krylov import to_complex, with_time_prefactor
+from .sweep import forward_sweep, reverse_sweep
+
+logger = logging.getLogger(__name__)
+
+
+# ---------------------------------------------------------------------------
+# Options
+# ---------------------------------------------------------------------------
+
+@dataclass
+class Options(AlgorithmOptions):
+ """2-site TDVP run options.
+
+ Parameters
+ ----------
+ dt:
+ Time step. Real time (``exp(-i dt H)``) unless ``imaginary_time`` is set.
+ n_steps:
+ Number of time steps to perform.
+ max_bond:
+ Maximum bond dimension kept by the per-bond SVD truncation. ``None`` means
+ no explicit cap (the bond grows up to the local capacity).
+ cutoff:
+ Singular-value threshold of the per-bond SVD truncation.
+ lanczos_tol:
+ Termination tolerance of the local Lanczos ``expv`` solves.
+ lanczos_maxiter:
+ Maximum Lanczos iterations per local substep.
+ imaginary_time:
+ If ``True``, evolve with ``exp(-dt H)`` (imaginary time) instead of
+ ``exp(-i dt H)``. Combined with ``normalize`` this cools toward the
+ ground state.
+ normalize:
+ If ``True`` (default), renormalise the state after every step.
+ """
+
+ dt: float = 0.05
+ n_steps: int = 10
+ max_bond: Optional[int] = None
+ cutoff: float = 1e-12
+ lanczos_tol: float = 1e-15
+ lanczos_maxiter: int = 30
+ imaginary_time: bool = False
+ normalize: bool = True
+
+
+# ---------------------------------------------------------------------------
+# Summary
+# ---------------------------------------------------------------------------
+
+@dataclass
+class Summary(AlgorithmSummary):
+ """2-site TDVP output.
+
+ Attributes
+ ----------
+ state:
+ Evolved MPS after all steps (orthogonality center at site 0).
+ n_steps:
+ Number of steps performed.
+ times:
+ Cumulative evolution time after each step (length ``n_steps``).
+ norms:
+ State norm after each step *before* renormalisation (length ``n_steps``).
+ bond_dims:
+ Bond dimensions of ``state`` after the final step (length ``L - 1``).
+ max_bond_dims:
+ Maximum kept bond dimension after each step (length ``n_steps``).
+ """
+
+ state: MPS
+ n_steps: int = 0
+ times: List[float] = field(default_factory=list)
+ norms: List[float] = field(default_factory=list)
+ bond_dims: List[int] = field(default_factory=list)
+ max_bond_dims: List[int] = field(default_factory=list)
+
+ def serialize(self) -> Dict:
+ """Serialize the summary to a plain dict compatible with ``torch.save``."""
+ return {
+ 'version': 1,
+ 'n_steps': self.n_steps,
+ 'times': self.times,
+ 'norms': self.norms,
+ 'bond_dims': self.bond_dims,
+ 'max_bond_dims': self.max_bond_dims,
+ 'state': self.state.serialize(),
+ }
+
+ @classmethod
+ def deserialize(cls, data: Dict, device: str = 'cpu') -> Summary:
+ """Reconstruct a `Summary` from a dict produced by `serialize`."""
+ version = data.get('version', 1)
+ if version != 1:
+ raise ValueError(f"Unsupported Summary serialization version: {version!r}")
+ return cls(
+ state=Network.deserialize(data['state'], device=device),
+ n_steps=data['n_steps'],
+ times=data['times'],
+ norms=data['norms'],
+ bond_dims=data['bond_dims'],
+ max_bond_dims=data.get('max_bond_dims', []),
+ )
+
+
+# ---------------------------------------------------------------------------
+# Half-sweep environment preparation
+# ---------------------------------------------------------------------------
+
+def _do_forward(mps: MPS, mpo: MPO, tau: float, maxdim, cutoff, lanczos_tol, lanczos_maxiter):
+ """Right-canonicalise, build all right environments, run a forward half-sweep.
+
+ The right environments are built here (rather than via the DMRG bulk builder)
+ so the dim-1 boundary block can be promoted to ``complex128`` — for real-time
+ evolution the state and MPO are complex, and the transfer contractions require
+ all three tensors to share a dtype.
+ """
+ L = mps.L
+ mps.canonical(0)
+ env_left = Environment(L, fetch_lo=0, fetch_hi=L - 2)
+ env_right = Environment(L, fetch_lo=1, fetch_hi=L - 1)
+ env_left[0] = to_complex(left_env_boundary(mps, mpo))
+ env_right[L - 1] = to_complex(right_env_boundary(mps, mpo))
+ for i in range(L - 2, 0, -1): # env_right[i] accumulates sites i+1 … L-1
+ env_right[i] = step_right_env(env_right[i + 1], mps[i + 1], mpo[i + 1])
+ forward_sweep(mps, mpo, env_left, env_right, tau,
+ maxdim=maxdim, cutoff=cutoff,
+ lanczos_tol=lanczos_tol, lanczos_maxiter=lanczos_maxiter)
+
+
+def _do_reverse(mps: MPS, mpo: MPO, tau: float, maxdim, cutoff, lanczos_tol, lanczos_maxiter):
+ """Left-canonicalise, build all left environments, run a reverse half-sweep."""
+ L = mps.L
+ mps.canonical(L - 1)
+ env_left = Environment(L, fetch_lo=0, fetch_hi=L - 2)
+ env_right = Environment(L, fetch_lo=1, fetch_hi=L - 1)
+ env_left[0] = to_complex(left_env_boundary(mps, mpo))
+ env_right[L - 1] = to_complex(right_env_boundary(mps, mpo))
+ for i in range(L - 1): # env_left[i+1] accumulates sites 0 … i
+ env_left[i + 1] = step_left_env(env_left[i], mps[i], mpo[i])
+ reverse_sweep(mps, mpo, env_left, env_right, tau,
+ maxdim=maxdim, cutoff=cutoff,
+ lanczos_tol=lanczos_tol, lanczos_maxiter=lanczos_maxiter)
+
+
+# ---------------------------------------------------------------------------
+# Top-level entry point
+# ---------------------------------------------------------------------------
+
+def run(mps: MPS, mpo: MPO, opts: Optional[Options] = None) -> Summary:
+ """Evolve an MPS under a Hamiltonian MPO with the 2-site TDVP integrator.
+
+ Performs ``opts.n_steps`` symmetric (Strang) steps: each step is a forward
+ half-sweep of duration ``dt/2`` followed by a reverse half-sweep of ``dt/2``.
+ The state is returned with ``center == 0``.
+
+ Parameters
+ ----------
+ mps:
+ Initial MPS state. Promoted to ``complex128`` and canonicalised in-place.
+ mpo:
+ Hamiltonian MPO of the same length as ``mps``.
+ opts:
+ Run options. Defaults to ``Options()`` if ``None``.
+
+ Returns
+ -------
+ Summary
+ Evolved state and time/norm/bond-dimension history.
+
+ Raises
+ ------
+ ValueError
+ If ``mps`` has fewer than two sites, or ``mps`` and ``mpo`` differ in length.
+ """
+ if opts is None:
+ opts = Options()
+ if mps.L < 2:
+ raise ValueError(f"2-site TDVP evolution requires at least 2 sites, got L={mps.L}")
+ if mps.L != mpo.L:
+ raise ValueError(f"mps and mpo must have the same length, got {mps.L} and {mpo.L}")
+
+ maxdim = opts.max_bond
+ prefactor: complex = -1.0 if opts.imaginary_time else -1j
+
+ # Promote both state and Hamiltonian to complex128 so every effective-H
+ # contraction and local exponential shares the backend dtype (the DMRG path
+ # keeps these real; real-time TDVP needs the complex exponential).
+ for site in range(mps.L):
+ mps[site] = to_complex(mps[site])
+ mpo = MPO([to_complex(mpo[b]) for b in range(mpo.L)])
+ mps.canonical(0)
+
+ half = 0.5 * opts.dt
+ times: List[float] = []
+ norms: List[float] = []
+ max_bond_dims: List[int] = []
+
+ logger.info("─" * 60)
+ logger.info("Commencing: Two-Site TDVP Time Evolution".center(60))
+ logger.info("─" * 60)
+ logger.info("")
+ logger.info(" chain length : %d", mps.L)
+ logger.info(" time step : %g", opts.dt)
+ logger.info(" steps : %d", opts.n_steps)
+ logger.info(" evolution : %s", "imaginary" if opts.imaginary_time else "real")
+ logger.info(" max bond dim : %s", opts.max_bond if opts.max_bond is not None else 'unlimited')
+ logger.info("")
+
+ w = len(str(opts.n_steps))
+ with with_time_prefactor(prefactor):
+ for step in range(opts.n_steps):
+ # Symmetric Strang step: forward(dt/2) then reverse(dt/2).
+ _do_forward(mps, mpo, half, maxdim, opts.cutoff, opts.lanczos_tol, opts.lanczos_maxiter)
+ _do_reverse(mps, mpo, half, maxdim, opts.cutoff, opts.lanczos_tol, opts.lanczos_maxiter)
+
+ norm = mps.norm()
+ if opts.normalize:
+ mps.normalize()
+
+ times.append((step + 1) * opts.dt)
+ norms.append(norm)
+ max_bond_dims.append(max(mps.bond_dims) if mps.bond_dims else 1)
+
+ logger.info("step %*d / %d: t = %g, norm = %.10f, kept bond = %d",
+ w, step + 1, opts.n_steps, times[-1], norm, max_bond_dims[-1])
+
+ if mps.center != 0:
+ mps.canonical(0)
+
+ logger.info("")
+
+ return Summary(
+ state=mps,
+ n_steps=opts.n_steps,
+ times=times,
+ norms=norms,
+ bond_dims=list(mps.bond_dims),
+ max_bond_dims=max_bond_dims,
+ )
diff --git a/tests/algorithm/bond_update_bug/__init__.py b/tests/algorithm/bond_update_bug/__init__.py
new file mode 100644
index 0000000..e414966
--- /dev/null
+++ b/tests/algorithm/bond_update_bug/__init__.py
@@ -0,0 +1,19 @@
+# Copyright (C) 2025-2026 Changkai Zhang.
+#
+# This file is part of Alice project.
+#
+# Alice is free software: you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published
+# by the Free Software Foundation, either version 3 of the License,
+# or (at your option) any later version.
+#
+# Alice is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with Alice. If not, see .
+
+
+"""Tests for alice.algorithm.bond_update_bug: gate-based two-site BUG integrator."""
diff --git a/tests/algorithm/bond_update_bug/conftest.py b/tests/algorithm/bond_update_bug/conftest.py
new file mode 100644
index 0000000..825989d
--- /dev/null
+++ b/tests/algorithm/bond_update_bug/conftest.py
@@ -0,0 +1,281 @@
+# Copyright (C) 2025-2026 Changkai Zhang.
+#
+# This file is part of Alice project.
+#
+# Alice is free software: you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published
+# by the Free Software Foundation, either version 3 of the License,
+# or (at your option) any later version.
+#
+# Alice is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with Alice. If not, see .
+
+
+"""Pytest fixtures and exact-diagonalization helpers for bond_update_bug tests.
+
+The helpers build a dense Heisenberg Hamiltonian, dense product states, and a
+dense vector from an MPS — all in the same physical basis ordering as Nicole's
+spin-1/2 U(1) space — so the integrator can be checked against exact
+diagonalization.
+"""
+
+from __future__ import annotations
+
+import functools
+from typing import Dict, List, Tuple
+
+import pytest
+import torch
+from nicole import Index, Tensor, load_space
+
+from alice.network import MPS, build_interaction
+
+
+@pytest.fixture(autouse=True)
+def _isolate_cwd(tmp_path, monkeypatch):
+ """Run every test in a fresh working directory."""
+ monkeypatch.chdir(tmp_path)
+
+
+@pytest.fixture(scope='session')
+def spin_space() -> Tuple[Index, Dict[str, Tensor]]:
+ """Spin-1/2 U(1) physical space and operators (shared across the session)."""
+ return load_space('Spin', 'U1', {'J': 0.5})
+
+
+def heisenberg_chain(length: int, coupling: float = 1.0):
+ """Build the Heisenberg interaction list, physical index, and geometry.
+
+ Parameters
+ ----------
+ length:
+ Number of sites.
+ coupling:
+ Isotropic exchange coupling `J`.
+
+ Returns
+ -------
+ tuple
+ `(interactions, spc, geo)` from `build_interaction`.
+ """
+ cfg = {
+ 'geometry': {'lattice': 'chain', 'lx': length, 'bcx': 'OBC', 'n2x': True},
+ 'model': {
+ 'category': 'bosonic', 'label': 'Heisenberg',
+ 'symmetry': 'U1', 'spin': 0.5, 'J': coupling,
+ },
+ }
+ return build_interaction(cfg)
+
+
+def _spin_matrices(charges: List[int]):
+ """Return dense `(Sz, Sp, Sm)` in the sector order given by `charges`."""
+ sz = torch.diag(torch.tensor([c / 2.0 for c in charges], dtype=torch.complex128))
+ up = 0 if charges[0] > charges[1] else 1
+ sp = torch.zeros((2, 2), dtype=torch.complex128)
+ sp[up, 1 - up] = 1.0
+ return sz, sp, sp.conj().T.contiguous()
+
+
+def _embed(op: torch.Tensor, site: int, length: int) -> torch.Tensor:
+ """Embed a single-site operator into the full `2**length` Hilbert space."""
+ eye = torch.eye(2, dtype=torch.complex128)
+ factors = [op if k == site else eye for k in range(length)]
+ return functools.reduce(lambda a, b: torch.kron(a.contiguous(), b.contiguous()), factors)
+
+
+def dense_heisenberg(length: int, charges: List[int], coupling: float = 1.0) -> torch.Tensor:
+ """Build the dense Heisenberg Hamiltonian matching Alice's spin basis.
+
+ Parameters
+ ----------
+ length:
+ Number of sites.
+ charges:
+ Sector charges of the physical index, in dense order (from
+ `Spc.sectors`), used to fix the single-site basis ordering.
+ coupling:
+ Isotropic exchange coupling `J`.
+
+ Returns
+ -------
+ torch.Tensor
+ Dense `(2**length, 2**length)` Hamiltonian.
+ """
+ sz, sp, sm = _spin_matrices(charges)
+ dim = 2 ** length
+ ham = torch.zeros((dim, dim), dtype=torch.complex128)
+ for i in range(length - 1):
+ ham = ham + coupling * (
+ _embed(sz, i, length) @ _embed(sz, i + 1, length)
+ + 0.5 * (_embed(sp, i, length) @ _embed(sm, i + 1, length))
+ + 0.5 * (_embed(sm, i, length) @ _embed(sp, i + 1, length))
+ )
+ return ham
+
+
+def dense_total_sz(length: int, charges: List[int]) -> torch.Tensor:
+ """Build the dense total-`S_z` operator matching Alice's spin basis."""
+ sz, _, _ = _spin_matrices(charges)
+ return sum(_embed(sz, i, length) for i in range(length))
+
+
+def dense_sz_profile(vec: torch.Tensor, length: int, charges: List[int]) -> List[float]:
+ """Per-site `` of a dense state vector, in Alice's spin basis.
+
+ Built from the same `_spin_matrices`/`_embed` machinery as `dense_total_sz`, so
+ it is convention-exact against `mps_to_vector`'s basis ordering rather than
+ relying on a hand-written site embedding.
+ """
+ sz, _, _ = _spin_matrices(charges)
+ return [torch.vdot(vec, _embed(sz, i, length) @ vec).real.item()
+ for i in range(length)]
+
+
+def dense_hamiltonian(interactions, length: int, charges: List[int]) -> torch.Tensor:
+ """Assemble the full `d**L` dense Hamiltonian from Alice's own bond terms.
+
+ Densifies each nearest-neighbour `Interaction2Site` bond Hamiltonian exactly as
+ the integrator consumes it (`build_bond_generators`) and lifts it to the full
+ Hilbert space. This is convention-exact — the dense operator is, by
+ construction, the same Hamiltonian the MPS evolves under — so it avoids any
+ basis/normalisation mismatch a hand-written model matrix could introduce.
+
+ Parameters
+ ----------
+ interactions:
+ Interaction list from `build_interaction`.
+ length:
+ Number of sites `L`.
+ charges:
+ Charges of the physical space in dense order (fixes the local basis).
+
+ Returns
+ -------
+ torch.Tensor
+ Dense `(d**L, d**L)` Hamiltonian, `d = len(charges)`.
+ """
+ from alice.algorithm.bond_update_bug._kernel import to_dense
+ from alice.algorithm.bond_update_bug.bond import build_bond_generators
+
+ generators = build_bond_generators(interactions, length)
+ d = len(charges)
+ dim = d ** length
+ ham = torch.zeros((dim, dim), dtype=torch.complex128)
+ eye = torch.eye(d, dtype=torch.complex128)
+ for bond, h in enumerate(generators):
+ if h is None:
+ continue
+ # h axes: (bra_i, ket_i, bra_j, ket_j). Densify, then reorder to the
+ # operator matrix [(bra_i, bra_j), (ket_i, ket_j)].
+ dense = to_dense(h, [h.itags[0], h.itags[1], h.itags[2], h.itags[3]]).to(torch.complex128)
+ local = dense.permute(0, 2, 1, 3).reshape(d * d, d * d)
+ factors: List[torch.Tensor] = []
+ site = 0
+ while site < length:
+ if site == bond:
+ factors.append(local)
+ site += 2
+ else:
+ factors.append(eye)
+ site += 1
+ lifted = factors[0]
+ for factor in factors[1:]:
+ lifted = torch.kron(lifted.contiguous(), factor.contiguous())
+ ham = ham + lifted
+ return ham
+
+
+def product_vector(config: List[int], charges: List[int]) -> torch.Tensor:
+ """Build the dense product-state vector for a sector-index configuration.
+
+ Parameters
+ ----------
+ config:
+ Per-site sector index (0 or 1) — the same `config` passed to `init_mps`.
+ charges:
+ Sector charges in dense order (unused beyond fixing length-2 basis).
+
+ Returns
+ -------
+ torch.Tensor
+ Dense state vector of length `2**len(config)`.
+ """
+ basis = [
+ torch.tensor([1.0, 0.0], dtype=torch.complex128),
+ torch.tensor([0.0, 1.0], dtype=torch.complex128),
+ ]
+ return functools.reduce(
+ lambda a, b: torch.kron(a.contiguous(), b.contiguous()),
+ [basis[c] for c in config],
+ )
+
+
+def exact_evolve(ham: torch.Tensor, psi0: torch.Tensor, t: float) -> torch.Tensor:
+ """Return `exp(-i t H) |psi0>` via dense eigendecomposition."""
+ evals, evecs = torch.linalg.eigh(ham)
+ return evecs @ (torch.exp(-1j * t * evals) * (evecs.conj().T @ psi0))
+
+
+def _core_dense(core: Tensor, charges: List[int]) -> torch.Tensor:
+ """Densify a 3-index MPS core `(left, right, phys)` to a dense torch tensor.
+
+ The bonds are densified to their own (symmetry-restricted) dimensions, but the
+ physical axis is *embedded into the full local basis* of size `len(charges)`:
+ a symmetric core only stores the physical sectors its charge structure allows
+ (e.g. a boundary site pinned to one charge has a dim-1 physical leg), so each
+ present physical charge `q` is placed at its global basis index
+ `charges.index(q)` and the rest is zero. This makes the contracted state live
+ in the full `len(charges)**L` space the dense ED helpers use.
+ """
+ phys_table = {q: (charges.index(q), 1) for q in charges}
+ offsets = []
+ for axis, index in enumerate(core.indices):
+ if axis == 2: # physical leg → full local basis
+ offsets.append((phys_table, len(charges)))
+ continue
+ table = {}
+ cursor = 0
+ for sector in index.sectors:
+ table[sector.charge] = (cursor, sector.dim)
+ cursor += sector.dim
+ offsets.append((table, cursor))
+ dense = torch.zeros([total for _, total in offsets], dtype=torch.complex128)
+ for key, block in core.data.items():
+ slices = tuple(
+ slice(offsets[axis][0][key[axis]][0],
+ offsets[axis][0][key[axis]][0] + offsets[axis][0][key[axis]][1])
+ for axis in range(3)
+ )
+ dense[slices] = block.to(torch.complex128)
+ return dense
+
+
+def mps_to_vector(mps: MPS, charges: List[int]) -> torch.Tensor:
+ """Contract an OBC MPS into a dense state vector in the full physical basis.
+
+ Parameters
+ ----------
+ mps:
+ MPS with trivial (dimension-1) boundary bonds.
+ charges:
+ Charges of the full physical space, in dense order (from `Spc.sectors`).
+ Each site's physical leg is embedded into this `len(charges)`-dimensional
+ basis (see `_core_dense`), so the result has length `len(charges)**L`
+ regardless of which charges each site's symmetric core actually carries.
+
+ Returns
+ -------
+ torch.Tensor
+ Dense state vector of length `len(charges)**L`.
+ """
+ psi = _core_dense(mps[0], charges)[0] # drop trivial left bond -> (right, phys_0)
+ for site in range(1, mps.L):
+ psi = torch.tensordot(psi, _core_dense(mps[site], charges), dims=([0], [0]))
+ psi = psi.movedim(-2, 0) # keep the open right bond at the front
+ return psi[0].reshape(-1) # drop trivial right bond
diff --git a/tests/algorithm/bond_update_bug/test_bond.py b/tests/algorithm/bond_update_bug/test_bond.py
new file mode 100644
index 0000000..8556847
--- /dev/null
+++ b/tests/algorithm/bond_update_bug/test_bond.py
@@ -0,0 +1,57 @@
+# Copyright (C) 2025-2026 Changkai Zhang.
+#
+# This file is part of Alice project.
+#
+# Alice is free software: you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published
+# by the Free Software Foundation, either version 3 of the License,
+# or (at your option) any later version.
+#
+# Alice is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with Alice. If not, see .
+
+
+"""Tests for nearest-neighbour bond Hamiltonian extraction and gate relabelling."""
+
+from __future__ import annotations
+
+from alice.algorithm.bond_update_bug.bond import bond_hamiltonian, build_bond_generators, kernel_gate
+
+from .conftest import heisenberg_chain
+
+
+def _first_generator(length=4):
+ """Return the first active two-site bond Hamiltonian of a Heisenberg chain."""
+ interactions, _, geo = heisenberg_chain(length)
+ generators = build_bond_generators(interactions, geo.L)
+ bond = next(i for i, g in enumerate(generators) if g is not None)
+ return generators[bond]
+
+
+def test_bond_hamiltonian_is_four_index():
+ """A nearest-neighbour bond term has axes (bra_i, ket_i, bra_{i+1}, ket_{i+1})."""
+ h = _first_generator()
+ assert len(h.indices) == 4
+
+
+def test_kernel_gate_itags_and_axes():
+ """`kernel_gate` relabels the term into the kernel's (ket_i, ket_j, bra_i*, bra_j*) order."""
+ h = _first_generator()
+ gate = kernel_gate(h, 's00', 's01')
+ # Axes are (ket_i, ket_j, bra_i, bra_j) with the bra (output) legs starred.
+ assert list(gate.itags) == ['s00', 's01', 's00*', 's01*']
+ assert len(gate.indices) == 4
+
+
+def test_kernel_gate_is_complex():
+ """The gate must be complex so the local exponentials share the backend dtype."""
+ import torch
+
+ h = _first_generator()
+ gate = kernel_gate(h, 's00', 's01')
+ assert all(block.dtype == torch.complex128 for block in gate.data.values())
diff --git a/tests/algorithm/bond_update_bug/test_bond_update_bug.py b/tests/algorithm/bond_update_bug/test_bond_update_bug.py
new file mode 100644
index 0000000..21edab6
--- /dev/null
+++ b/tests/algorithm/bond_update_bug/test_bond_update_bug.py
@@ -0,0 +1,261 @@
+# Copyright (C) 2025-2026 Changkai Zhang.
+#
+# This file is part of Alice project.
+#
+# Alice is free software: you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published
+# by the Free Software Foundation, either version 3 of the License,
+# or (at your option) any later version.
+#
+# Alice is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with Alice. If not, see .
+
+
+"""Tests for the bond_update_bug integrator (Options, Summary, run)."""
+
+from __future__ import annotations
+
+import dataclasses
+
+import pytest
+import torch
+from nicole import Index, Tensor
+
+from alice import init_mps
+from alice.algorithm import bond_update_bug
+from alice.algorithm.bond_update_bug.bond import build_bond_generators
+from alice.network.interaction import Interaction2Site
+
+from .conftest import (
+ dense_hamiltonian,
+ dense_total_sz,
+ exact_evolve,
+ heisenberg_chain,
+ mps_to_vector,
+)
+
+
+def _domain_wall(length, spin_space):
+ """Return `(mps, interactions, charges, psi0)` for a full-phys Heisenberg domain wall.
+
+ The state is the Sz=0 domain wall `|↓…↓↑…↑⟩`. `init_mps(config=...)` builds it
+ as a product state but pins each physical leg to its single occupied charge
+ (dim-1 phys), which freezes the dynamics and cannot densify to the full `2**L`
+ space. Each physical leg is therefore inflated to the full spin-1/2 index (the
+ occupied-charge block is kept, the empty charge added) so spins can flip and the
+ state densifies to `2**L`. `psi0` is the dense initial vector (one nonzero
+ amplitude), in the same physical basis order as the dense ED helpers.
+ """
+ _, operators = spin_space
+ interactions, spc, _ = heisenberg_chain(length)
+ charges = [sector.charge for sector in spc.sectors]
+ config = [0] * (length // 2) + [1] * (length - length // 2)
+ target = sum(charges[c] for c in config)
+ mps = init_mps(length, spc, operators, config=config, target_qn=target)
+ # Inflate each pinned (dim-1) physical leg to the full local space.
+ for i in range(mps.L):
+ core = mps[i]
+ full_phys = Index(core.indices[2].direction, core.indices[2].group, spc.sectors)
+ mps[i] = Tensor(
+ indices=(core.indices[0], core.indices[1], full_phys),
+ itags=core.itags,
+ data={key: block.clone() for key, block in core.data.items()},
+ dtype=core.dtype,
+ )
+ psi0 = mps_to_vector(mps, charges)
+ return mps, interactions, charges, psi0
+
+
+# ---------------------------------------------------------------------------
+# Options
+# ---------------------------------------------------------------------------
+
+class TestOptions:
+ """Tests for the Options dataclass."""
+
+ def test_default_order(self):
+ assert bond_update_bug.Options().order == 'strang'
+
+ @pytest.mark.parametrize('alias,canonical', [
+ ('strang', 'strang'), ('second', 'strang'), ('2', 'strang'),
+ ('lie', 'lie'), ('first', 'lie'), ('1', 'lie'),
+ ])
+ def test_order_aliases(self, alias, canonical):
+ assert bond_update_bug.Options(order=alias).order == canonical
+
+ def test_unknown_order_raises(self):
+ with pytest.raises(ValueError, match='unknown Trotter order'):
+ bond_update_bug.Options(order='leapfrog')
+
+ def test_from_toml(self):
+ opts = bond_update_bug.Options.from_toml(
+ {'dt': 0.02, 'n_steps': 50, 'order': 'second', 'max_bond': 32}
+ )
+ assert opts.dt == 0.02
+ assert opts.n_steps == 50
+ assert opts.order == 'strang'
+ assert opts.max_bond == 32
+
+ def test_to_toml_round_trip(self, tmp_path):
+ original = bond_update_bug.Options(dt=0.01, n_steps=7, order='lie', max_bond=16)
+ path = tmp_path / 'opts.toml'
+ original.to_toml(path)
+ loaded = bond_update_bug.Options.load_toml(path)
+ assert loaded.dt == 0.01
+ assert loaded.n_steps == 7
+ assert loaded.order == 'lie'
+ assert loaded.max_bond == 16
+
+
+# ---------------------------------------------------------------------------
+# Summary
+# ---------------------------------------------------------------------------
+
+class TestSummary:
+ """Tests for the Summary dataclass."""
+
+ def test_serialize_round_trip(self, spin_space):
+ mps, interactions, _, _ = _domain_wall(6, spin_space)
+ summary = bond_update_bug.run(
+ mps, interactions, bond_update_bug.Options(dt=0.05, n_steps=3, max_bond=16)
+ )
+ restored = bond_update_bug.Summary.deserialize(summary.serialize())
+ assert restored.n_steps == summary.n_steps
+ assert restored.bond_dims == summary.bond_dims
+ assert restored.times == pytest.approx(summary.times)
+ assert restored.state.L == summary.state.L
+
+
+# ---------------------------------------------------------------------------
+# Generators / error handling
+# ---------------------------------------------------------------------------
+
+class TestGenerators:
+ """Tests for bond-generator extraction from the interaction list."""
+
+ def test_all_bonds_populated_for_heisenberg(self):
+ interactions, _, geo = heisenberg_chain(5)
+ generators = build_bond_generators(interactions, geo.L)
+ assert len(generators) == geo.L - 1
+ assert all(g is not None for g in generators)
+
+ def test_long_range_term_raises(self):
+ # A synthetic non-nearest-neighbour term must be rejected.
+ interactions, _, geo = heisenberg_chain(4)
+ far = dataclasses.replace(
+ next(i for i in interactions if isinstance(i, Interaction2Site)),
+ leading_site=0, terminal_site=2,
+ )
+ with pytest.raises(NotImplementedError, match='nearest-neighbour'):
+ build_bond_generators([far], geo.L)
+
+
+# ---------------------------------------------------------------------------
+# Dynamics
+# ---------------------------------------------------------------------------
+
+class TestDynamics:
+ """Physical correctness of the time evolution."""
+
+ def test_norm_conserved_real_time(self, spin_space):
+ mps, interactions, _, _ = _domain_wall(6, spin_space)
+ summary = bond_update_bug.run(
+ mps, interactions,
+ bond_update_bug.Options(dt=0.05, n_steps=10, max_bond=64, normalize=False),
+ )
+ for norm in summary.norms:
+ assert abs(norm - 1.0) < 1e-10
+
+ def test_total_sz_conserved(self, spin_space):
+ mps, interactions, charges, psi0 = _domain_wall(6, spin_space)
+ sz_total = dense_total_sz(6, charges)
+ sz_before = (psi0.conj() @ sz_total @ psi0).real.item() / psi0.norm().item() ** 2
+ summary = bond_update_bug.run(
+ mps, interactions, bond_update_bug.Options(dt=0.05, n_steps=10, max_bond=64)
+ )
+ vec = mps_to_vector(summary.state, charges)
+ sz_after = (vec.conj() @ sz_total @ vec).real.item() / vec.norm().item() ** 2
+ assert abs(sz_after - sz_before) < 1e-10
+
+ def test_fidelity_matches_exact_diagonalization(self, spin_space):
+ length = 6
+ mps, interactions, charges, psi0 = _domain_wall(length, spin_space)
+ ham = dense_hamiltonian(interactions, length, charges)
+ psi0 = psi0 / psi0.norm()
+ dt, n_steps = 0.05, 20
+ summary = bond_update_bug.run(
+ mps, interactions,
+ bond_update_bug.Options(dt=dt, n_steps=n_steps, max_bond=64, normalize=False),
+ )
+ evolved = mps_to_vector(summary.state, charges)
+ evolved = evolved / evolved.norm()
+ exact = exact_evolve(ham, psi0, dt * n_steps)
+ exact = exact / exact.norm()
+ fidelity = abs(torch.vdot(exact, evolved)).item()
+ assert 1.0 - fidelity < 1e-6
+
+ def test_strang_converges_second_order(self, spin_space):
+ length = 6
+ _, interactions, charges, psi0 = _domain_wall(length, spin_space)
+ ham = dense_hamiltonian(interactions, length, charges)
+ psi0 = psi0 / psi0.norm()
+
+ def infidelity(dt, n_steps):
+ mps, _, _, _ = _domain_wall(length, spin_space)
+ summary = bond_update_bug.run(
+ mps, interactions,
+ bond_update_bug.Options(dt=dt, n_steps=n_steps, max_bond=64, normalize=False),
+ )
+ evolved = mps_to_vector(summary.state, charges)
+ evolved = evolved / evolved.norm()
+ exact = exact_evolve(ham, psi0, dt * n_steps)
+ exact = exact / exact.norm()
+ return 1.0 - abs(torch.vdot(exact, evolved)).item()
+
+ coarse = infidelity(0.10, 10)
+ fine = infidelity(0.05, 20)
+ # Strang state error is O(dt^2), so the infidelity is O(dt^4): halving dt
+ # cuts it by ~16. Allow a generous band around the asymptotic ratio.
+ assert coarse / fine > 8.0
+
+ def test_strang_beats_lie(self, spin_space):
+ length = 6
+ _, interactions, charges, _ = _domain_wall(length, spin_space)
+ ham = dense_hamiltonian(interactions, length, charges)
+
+ def infidelity(order):
+ mps, interactions, charges, psi0 = _domain_wall(length, spin_space)
+ psi0 = psi0 / psi0.norm()
+ summary = bond_update_bug.run(
+ mps, interactions,
+ bond_update_bug.Options(dt=0.1, n_steps=10, order=order, max_bond=64, normalize=False),
+ )
+ evolved = mps_to_vector(summary.state, charges)
+ evolved = evolved / evolved.norm()
+ exact = exact_evolve(ham, psi0, 1.0)
+ exact = exact / exact.norm()
+ return 1.0 - abs(torch.vdot(exact, evolved)).item()
+
+ assert infidelity('strang') < infidelity('lie')
+
+ def test_imaginary_time_lowers_energy(self, spin_space):
+ length = 6
+ mps, interactions, charges, psi0 = _domain_wall(length, spin_space)
+ ham = dense_hamiltonian(interactions, length, charges)
+ ground = torch.linalg.eigvalsh(ham)[0].item()
+ psi0 = psi0 / psi0.norm()
+ energy_before = (psi0.conj() @ ham @ psi0).real.item()
+ summary = bond_update_bug.run(
+ mps, interactions,
+ bond_update_bug.Options(dt=0.05, n_steps=40, imaginary_time=True, max_bond=64),
+ )
+ vec = mps_to_vector(summary.state, charges)
+ vec = vec / vec.norm()
+ energy_after = (vec.conj() @ ham @ vec).real.item()
+ assert energy_after < energy_before
+ assert energy_after > ground - 1e-9
diff --git a/tests/algorithm/bond_update_bug/test_kernel.py b/tests/algorithm/bond_update_bug/test_kernel.py
new file mode 100644
index 0000000..f19a271
--- /dev/null
+++ b/tests/algorithm/bond_update_bug/test_kernel.py
@@ -0,0 +1,197 @@
+# Copyright (C) 2025-2026 Changkai Zhang.
+#
+# This file is part of Alice project.
+#
+# Alice is free software: you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published
+# by the Free Software Foundation, either version 3 of the License,
+# or (at your option) any later version.
+#
+# Alice is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with Alice. If not, see .
+# Author of code: Madhav Menon.
+
+
+"""Tests for the `bond_update_bug` local K/L/S bond update.
+
+The update (see :mod:`alice.algorithm.bond_update_bug._kernel.kls.candidate`)
+projects the K/L generators by the discarded (orthogonal-complement) projector
+*before* the exponential and acts the augmented **isometries directly** in the
+S-step (``Ŝ0 = Û† Θ0 V̂†``), forming no overlap matrices. These tests check it
+against exact diagonalization, the conservation laws (norm and total Sz), the
+second-order Strang convergence, and imaginary-time cooling to the ground state.
+"""
+
+from __future__ import annotations
+
+import pytest
+import torch
+from nicole import Index, Tensor
+
+from alice import init_mps
+from alice.algorithm import bond_update_bug
+
+from .conftest import (
+ dense_hamiltonian,
+ dense_sz_profile,
+ dense_total_sz,
+ exact_evolve,
+ heisenberg_chain,
+ mps_to_vector,
+)
+
+
+def _neel(length, spin_space):
+ """Return ``(mps, interactions, charges, psi0)`` for a full-phys Néel state.
+
+ The Néel state ``|↑↓↑↓…⟩`` (alternating ``config=[0,1,0,1,…]``) lives in the
+ Sz=0 sector for even ``length``. Each physical leg is inflated to the full
+ spin-1/2 index so spins can flip and the state densifies to ``2**length``;
+ ``psi0`` is the dense initial vector in the ED helpers' basis order.
+ """
+ _, operators = spin_space
+ interactions, spc, _ = heisenberg_chain(length)
+ charges = [sector.charge for sector in spc.sectors]
+ config = [0, 1] * (length // 2)
+ target = sum(charges[c] for c in config)
+ mps = init_mps(length, spc, operators, config=config, target_qn=target)
+ for i in range(mps.L):
+ core = mps[i]
+ full_phys = Index(core.indices[2].direction, core.indices[2].group, spc.sectors)
+ mps[i] = Tensor(
+ indices=(core.indices[0], core.indices[1], full_phys),
+ itags=core.itags,
+ data={key: block.clone() for key, block in core.data.items()},
+ dtype=core.dtype,
+ )
+ psi0 = mps_to_vector(mps, charges)
+ return mps, interactions, charges, psi0
+
+
+def _opts(**kwargs):
+ """`bond_update_bug` options with sensible test defaults."""
+ return bond_update_bug.Options(**kwargs)
+
+
+# ---------------------------------------------------------------------------
+# Conservation
+# ---------------------------------------------------------------------------
+
+class TestConservation:
+ """Norm (real time) and total Sz are conserved by the discarded S-step."""
+
+ def test_norm_conserved_real_time(self, spin_space):
+ mps, interactions, _, _ = _neel(6, spin_space)
+ summary = bond_update_bug.run(
+ mps, interactions, _opts(dt=0.05, n_steps=10, max_bond=64, normalize=False)
+ )
+ for norm in summary.norms:
+ assert abs(norm - 1.0) < 1e-10
+
+ def test_total_sz_conserved(self, spin_space):
+ mps, interactions, charges, psi0 = _neel(6, spin_space)
+ sz_total = dense_total_sz(6, charges)
+ sz_before = (psi0.conj() @ sz_total @ psi0).real.item() / psi0.norm().item() ** 2
+ summary = bond_update_bug.run(
+ mps, interactions, _opts(dt=0.05, n_steps=10, max_bond=64)
+ )
+ vec = mps_to_vector(summary.state, charges)
+ sz_after = (vec.conj() @ sz_total @ vec).real.item() / vec.norm().item() ** 2
+ assert abs(sz_after - sz_before) < 1e-10
+
+
+# ---------------------------------------------------------------------------
+# Accuracy: the discarded-projector + augmented-isometry S-step is correct
+# ---------------------------------------------------------------------------
+
+class TestAccuracy:
+ """Real-time accuracy of the discarded S-step against ED and the kernel."""
+
+ def test_fidelity_matches_exact_diagonalization(self, spin_space):
+ length = 6
+ mps, interactions, charges, psi0 = _neel(length, spin_space)
+ ham = dense_hamiltonian(interactions, length, charges)
+ psi0 = psi0 / psi0.norm()
+ dt, n_steps = 0.05, 20
+ summary = bond_update_bug.run(
+ mps, interactions, _opts(dt=dt, n_steps=n_steps, max_bond=64, normalize=False)
+ )
+ evolved = mps_to_vector(summary.state, charges)
+ evolved = evolved / evolved.norm()
+ exact = exact_evolve(ham, psi0, dt * n_steps)
+ exact = exact / exact.norm()
+ fidelity = abs(torch.vdot(exact, evolved)).item()
+ # At full bond dimension the only error is the Strang splitting (O(dt^2)).
+ assert 1.0 - fidelity < 1e-6
+
+ def test_strang_converges_second_order(self, spin_space):
+ """Strang state error is O(dt^2) -> infidelity O(dt^4): halving dt cuts ~16x.
+
+ Measured in the ASYMPTOTIC regime. At dt=0.1 the higher-order Trotter terms
+ are still large enough to contaminate the ratio (it reads ~6 at every system
+ size, L=4/6/8 alike), which measures how far dt is from asymptotia rather
+ than the method's order. Halving into dt=0.05/0.025 recovers the expected
+ behaviour. Verified against a dt/L scan: the ratio rises monotonically
+ towards 16 as dt shrinks (L=6: 6.33 -> 10.48 -> 14.63 at T=1.0/0.5/0.2),
+ which is the signature of a genuine second-order method; a rank-projection
+ floor would push the ratio DOWN as the Trotter error vanished, not up.
+ """
+ length = 6
+ _, interactions, charges, psi0 = _neel(length, spin_space)
+ ham = dense_hamiltonian(interactions, length, charges)
+ psi0 = psi0 / psi0.norm()
+
+ def infidelity(dt, n_steps):
+ mps, _, _, _ = _neel(length, spin_space)
+ summary = bond_update_bug.run(
+ mps, interactions, _opts(dt=dt, n_steps=n_steps, max_bond=64, normalize=False)
+ )
+ evolved = mps_to_vector(summary.state, charges)
+ evolved = evolved / evolved.norm()
+ exact = exact_evolve(ham, psi0, dt * n_steps)
+ exact = exact / exact.norm()
+ return 1.0 - abs(torch.vdot(exact, evolved)).item()
+
+ coarse = infidelity(0.05, 10)
+ fine = infidelity(0.025, 20)
+ assert coarse / fine > 8.0
+
+
+# ---------------------------------------------------------------------------
+# Imaginary time
+# ---------------------------------------------------------------------------
+
+class TestImaginaryTime:
+ """Imaginary-time cooling toward the exact ground state."""
+
+ def test_imaginary_time_reaches_ground_state(self, spin_space):
+ length = 6
+ mps, interactions, charges, psi0 = _neel(length, spin_space)
+ ham = dense_hamiltonian(interactions, length, charges)
+ evals, evecs = torch.linalg.eigh(ham)
+ ground_energy = evals[0].item()
+ ground_vec = evecs[:, 0]
+ psi0 = psi0 / psi0.norm()
+ err_before = 1.0 - abs(torch.vdot(ground_vec, psi0)).item()
+
+ summary = bond_update_bug.run(
+ mps, interactions,
+ _opts(dt=0.05, n_steps=160, imaginary_time=True, max_bond=64),
+ )
+ vec = mps_to_vector(summary.state, charges)
+ vec = vec / vec.norm()
+ energy_after = (vec.conj() @ ham @ vec).real.item()
+ err_after = 1.0 - abs(torch.vdot(ground_vec, vec)).item()
+
+ # Variational lower bound, substantial cooling, and tight final overlap.
+ assert energy_after > ground_energy - 1e-9
+ assert energy_after - ground_energy < 1e-2
+ assert err_after < err_before
+ assert err_after < 1e-2
+
+
diff --git a/tests/algorithm/tdvp2/__init__.py b/tests/algorithm/tdvp2/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/tests/algorithm/tdvp2/conftest.py b/tests/algorithm/tdvp2/conftest.py
new file mode 100644
index 0000000..4d53526
--- /dev/null
+++ b/tests/algorithm/tdvp2/conftest.py
@@ -0,0 +1,271 @@
+# Copyright (C) 2025-2026 Changkai Zhang.
+#
+# This file is part of Alice project.
+#
+# Alice is free software: you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published
+# by the Free Software Foundation, either version 3 of the License,
+# or (at your option) any later version.
+#
+# Alice is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with Alice. If not, see .
+# Author of code: Madhav Menon.
+
+
+"""Pytest fixtures and exact-diagonalization helpers for 2-site TDVP tests.
+
+These helpers build a dense Heisenberg Hamiltonian, a dense total-Sz operator,
+and a dense state vector from an MPS — all in the same physical basis ordering as
+Nicole's spin-1/2 U(1) space — so the integrator can be checked against exact
+diagonalization (the analytical reference for the small chains tested here).
+"""
+
+from __future__ import annotations
+
+import functools
+from typing import Dict, List, Tuple
+
+import pytest
+import torch
+from nicole import Index, Tensor, load_space
+
+from alice.network import MPS, build_interaction
+
+
+@pytest.fixture(autouse=True)
+def _isolate_cwd(tmp_path, monkeypatch):
+ """Run every test in a fresh working directory.
+
+ Parameters
+ ----------
+ tmp_path:
+ Pytest per-test temporary directory.
+ monkeypatch:
+ Pytest fixture used to change the working directory for the test.
+ """
+ monkeypatch.chdir(tmp_path)
+
+
+@pytest.fixture(scope='session')
+def spin_space() -> Tuple[Index, Dict[str, Tensor]]:
+ """Spin-1/2 U(1) physical space and operators (shared across the session).
+
+ Returns
+ -------
+ tuple
+ The `(spc, operators)` pair from `load_space('Spin', 'U1', {'J': 0.5})`.
+ """
+ return load_space('Spin', 'U1', {'J': 0.5})
+
+
+def heisenberg_chain(length: int, coupling: float = 1.0):
+ """Build the Heisenberg interaction list, physical index, and geometry.
+
+ Parameters
+ ----------
+ length:
+ Number of sites.
+ coupling:
+ Isotropic exchange coupling `J`.
+
+ Returns
+ -------
+ tuple
+ `(interactions, spc, geo)` from `build_interaction`.
+ """
+ cfg = {
+ 'geometry': {'lattice': 'chain', 'lx': length, 'bcx': 'OBC', 'n2x': True},
+ 'model': {
+ 'category': 'bosonic', 'label': 'Heisenberg',
+ 'symmetry': 'U1', 'spin': 0.5, 'J': coupling,
+ },
+ }
+ return build_interaction(cfg)
+
+
+def _spin_matrices(charges: List[int]):
+ """Return dense `(Sz, Sp, Sm)` spin-1/2 operators in the given sector order.
+
+ Parameters
+ ----------
+ charges:
+ Sector charges of the physical index in dense order, fixing the
+ single-site basis ordering.
+
+ Returns
+ -------
+ tuple
+ Dense `(Sz, Sp, Sm)` matrices.
+ """
+ sz = torch.diag(torch.tensor([c / 2.0 for c in charges], dtype=torch.complex128))
+ up = 0 if charges[0] > charges[1] else 1
+ sp = torch.zeros((2, 2), dtype=torch.complex128)
+ sp[up, 1 - up] = 1.0
+ return sz, sp, sp.conj().T.contiguous()
+
+
+def _embed(op: torch.Tensor, site: int, length: int) -> torch.Tensor:
+ """Embed a single-site operator into the full `2**length` Hilbert space.
+
+ Parameters
+ ----------
+ op:
+ Single-site `2 x 2` operator.
+ site:
+ Site index the operator acts on.
+ length:
+ Number of sites.
+
+ Returns
+ -------
+ torch.Tensor
+ The operator embedded as a `2**length x 2**length` dense matrix.
+ """
+ eye = torch.eye(2, dtype=torch.complex128)
+ factors = [op if k == site else eye for k in range(length)]
+ return functools.reduce(lambda a, b: torch.kron(a.contiguous(), b.contiguous()), factors)
+
+
+def dense_heisenberg(length: int, charges: List[int], coupling: float = 1.0) -> torch.Tensor:
+ """Build the dense Heisenberg Hamiltonian matching Alice's spin basis.
+
+ This is the same nearest-neighbour `J (Sz Sz + ½(S+ S- + S- S+))` chain that
+ `build_hamiltonian` assembles as an MPO from the Heisenberg interaction list,
+ written directly in the dense `2**length` basis so it can serve as the
+ exact-diagonalization reference.
+
+ Parameters
+ ----------
+ length:
+ Number of sites.
+ charges:
+ Sector charges of the physical index, in dense order, fixing the basis.
+ coupling:
+ Isotropic exchange coupling `J`.
+
+ Returns
+ -------
+ torch.Tensor
+ Dense `(2**length, 2**length)` Hamiltonian.
+ """
+ sz, sp, sm = _spin_matrices(charges)
+ dim = 2 ** length
+ ham = torch.zeros((dim, dim), dtype=torch.complex128)
+ for i in range(length - 1):
+ ham = ham + coupling * (
+ _embed(sz, i, length) @ _embed(sz, i + 1, length)
+ + 0.5 * (_embed(sp, i, length) @ _embed(sm, i + 1, length))
+ + 0.5 * (_embed(sm, i, length) @ _embed(sp, i + 1, length))
+ )
+ return ham
+
+
+def dense_total_sz(length: int, charges: List[int]) -> torch.Tensor:
+ """Build the dense total-`S_z` operator matching Alice's spin basis.
+
+ Parameters
+ ----------
+ length:
+ Number of sites.
+ charges:
+ Sector charges of the physical index, in dense order, fixing the basis.
+
+ Returns
+ -------
+ torch.Tensor
+ Dense `(2**length, 2**length)` total-`S_z` operator.
+ """
+ sz, _, _ = _spin_matrices(charges)
+ return sum(_embed(sz, i, length) for i in range(length))
+
+
+def exact_evolve(ham: torch.Tensor, psi0: torch.Tensor, t: float) -> torch.Tensor:
+ """Return `exp(-i t H) |psi0>` via dense eigendecomposition.
+
+ Parameters
+ ----------
+ ham:
+ Dense Hermitian Hamiltonian.
+ psi0:
+ Dense initial state vector.
+ t:
+ Evolution time.
+
+ Returns
+ -------
+ torch.Tensor
+ The exactly evolved dense state vector.
+ """
+ evals, evecs = torch.linalg.eigh(ham)
+ return evecs @ (torch.exp(-1j * t * evals) * (evecs.conj().T @ psi0))
+
+
+def _core_dense(core: Tensor, charges: List[int]) -> torch.Tensor:
+ """Densify a 3-index MPS core `(left, right, phys)` to a dense torch tensor.
+
+ The physical axis is embedded into the full local basis of size
+ `len(charges)`: a symmetric core only stores the physical sectors its charge
+ structure allows, so each present physical charge `q` is placed at its global
+ basis index `charges.index(q)` and the rest is zero.
+
+ Parameters
+ ----------
+ core:
+ MPS site tensor with axes `(left_bond, right_bond, physical)`.
+ charges:
+ Charges of the full physical space, in dense order.
+
+ Returns
+ -------
+ torch.Tensor
+ Dense `(left, right, len(charges))` tensor.
+ """
+ phys_table = {q: (charges.index(q), 1) for q in charges}
+ offsets = []
+ for axis, index in enumerate(core.indices):
+ if axis == 2:
+ offsets.append((phys_table, len(charges)))
+ continue
+ table = {}
+ cursor = 0
+ for sector in index.sectors:
+ table[sector.charge] = (cursor, sector.dim)
+ cursor += sector.dim
+ offsets.append((table, cursor))
+ dense = torch.zeros([total for _, total in offsets], dtype=torch.complex128)
+ for key, block in core.data.items():
+ slices = tuple(
+ slice(offsets[axis][0][key[axis]][0],
+ offsets[axis][0][key[axis]][0] + offsets[axis][0][key[axis]][1])
+ for axis in range(3)
+ )
+ dense[slices] = block.to(torch.complex128)
+ return dense
+
+
+def mps_to_vector(mps: MPS, charges: List[int]) -> torch.Tensor:
+ """Contract an OBC MPS into a dense state vector in the full physical basis.
+
+ Parameters
+ ----------
+ mps:
+ MPS with trivial (dimension-1) boundary bonds.
+ charges:
+ Charges of the full physical space, in dense order. Each site's physical
+ index is embedded into this `len(charges)`-dimensional basis.
+
+ Returns
+ -------
+ torch.Tensor
+ Dense state vector of length `len(charges)**L`.
+ """
+ psi = _core_dense(mps[0], charges)[0]
+ for site in range(1, mps.L):
+ psi = torch.tensordot(psi, _core_dense(mps[site], charges), dims=([0], [0]))
+ psi = psi.movedim(-2, 0)
+ return psi[0].reshape(-1)
diff --git a/tests/algorithm/tdvp2/test_tdvp2.py b/tests/algorithm/tdvp2/test_tdvp2.py
new file mode 100644
index 0000000..c9f93d0
--- /dev/null
+++ b/tests/algorithm/tdvp2/test_tdvp2.py
@@ -0,0 +1,210 @@
+# Copyright (C) 2025-2026 Changkai Zhang.
+#
+# This file is part of Alice project.
+#
+# Alice is free software: you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published
+# by the Free Software Foundation, either version 3 of the License,
+# or (at your option) any later version.
+#
+# Alice is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with Alice. If not, see .
+# Author of code: Madhav Menon.
+
+
+"""Tests for the 2-site TDVP integrator (Options, Summary, run).
+
+The 2-site TDVP integrator (see :mod:`alice.algorithm.tdvp2`) evolves an `MPS`
+under a Hamiltonian `MPO` by symmetric forward/reverse half-sweeps of effective-
+Hamiltonian exponentials, with an inverse-free one-site backward correction. These
+tests check, on the symmetric (isotropic) Heisenberg chain — which conserves total
+Sz and is available by exact diagonalization — that the integrator:
+
+* grows the bond dimension as a domain wall melts (rank adaptivity),
+* tracks the analytical (exact-diagonalization) solution to high accuracy,
+* conserves the state norm to machine precision (real time; TDVP is unitary),
+* conserves total Sz,
+* cools toward the ground state in imaginary time.
+
+A note on convergence: at fixed/adaptive bond dimension, 2-site TDVP's error is a
+projection (manifold) error that does **not** vanish as ``dt → 0`` — it plateaus,
+unlike the Strang state error of a full-rank propagator. The accuracy tests
+therefore assert a bounded, non-increasing error rather than a strict O(dt^2) ratio.
+"""
+
+from __future__ import annotations
+
+import pytest
+import torch
+from nicole import Index, Tensor
+
+from alice import build_hamiltonian, init_mps
+from alice.algorithm import tdvp2
+
+from .conftest import (
+ dense_heisenberg,
+ dense_total_sz,
+ exact_evolve,
+ heisenberg_chain,
+ mps_to_vector,
+)
+
+
+def _domain_wall(length, spin_space):
+ """Return `(mps, mpo, charges, psi0)` for a full-phys Heisenberg domain wall."""
+ _, operators = spin_space
+ interactions, spc, _ = heisenberg_chain(length)
+ charges = [sector.charge for sector in spc.sectors]
+ config = [0] * (length // 2) + [1] * (length - length // 2)
+ target = sum(charges[c] for c in config)
+ mps = init_mps(length, spc, operators, config=config, target_qn=target)
+ for i in range(mps.L):
+ core = mps[i]
+ full_phys = Index(core.indices[2].direction, core.indices[2].group, spc.sectors)
+ mps[i] = Tensor(
+ indices=(core.indices[0], core.indices[1], full_phys),
+ itags=core.itags,
+ data={key: block.clone() for key, block in core.data.items()},
+ dtype=core.dtype,
+ )
+ mpo = build_hamiltonian(interactions, length, spc)
+ psi0 = mps_to_vector(mps, charges)
+ return mps, mpo, charges, psi0
+
+
+# ---------------------------------------------------------------------------
+# Options / Summary
+# ---------------------------------------------------------------------------
+
+class TestOptions:
+ """Options/Summary basics."""
+
+ def test_defaults(self):
+ opts = tdvp2.Options()
+ assert opts.dt == 0.05
+ assert opts.max_bond is None
+ assert opts.normalize is True
+
+ def test_serialize_round_trip(self, spin_space):
+ mps, mpo, _, _ = _domain_wall(6, spin_space)
+ summary = tdvp2.run(mps, mpo, tdvp2.Options(dt=0.05, n_steps=3, max_bond=16))
+ restored = tdvp2.Summary.deserialize(summary.serialize())
+ assert restored.n_steps == summary.n_steps
+ assert restored.bond_dims == summary.bond_dims
+ assert restored.times == pytest.approx(summary.times)
+
+
+# ---------------------------------------------------------------------------
+# Error handling
+# ---------------------------------------------------------------------------
+
+class TestErrors:
+ def test_length_mismatch_raises(self, spin_space):
+ mps6, mpo6, _, _ = _domain_wall(6, spin_space)
+ _, mpo4, _, _ = _domain_wall(4, spin_space)
+ with pytest.raises(ValueError, match='same length'):
+ tdvp2.run(mps6, mpo4, tdvp2.Options(dt=0.05, n_steps=1))
+
+
+# ---------------------------------------------------------------------------
+# Rank adaptivity
+# ---------------------------------------------------------------------------
+
+class TestRankAdaptivity:
+ def test_bond_dimension_grows(self, spin_space):
+ mps, mpo, _, _ = _domain_wall(6, spin_space)
+ assert max(mps.bond_dims) == 1
+ summary = tdvp2.run(mps, mpo, tdvp2.Options(dt=0.05, n_steps=10, max_bond=64))
+ assert max(summary.bond_dims) > 1
+ assert max(summary.max_bond_dims) >= 4
+
+ def test_max_bond_cap_respected(self, spin_space):
+ mps, mpo, _, _ = _domain_wall(6, spin_space)
+ cap = 4
+ summary = tdvp2.run(mps, mpo, tdvp2.Options(dt=0.05, n_steps=10, max_bond=cap))
+ assert max(summary.bond_dims) <= cap
+
+
+# ---------------------------------------------------------------------------
+# Accuracy vs exact diagonalization
+# ---------------------------------------------------------------------------
+
+class TestAccuracy:
+ def test_fidelity_matches_exact_diagonalization(self, spin_space):
+ length = 6
+ mps, mpo, charges, psi0 = _domain_wall(length, spin_space)
+ ham = dense_heisenberg(length, charges)
+ psi0 = psi0 / psi0.norm()
+ dt, n_steps = 0.02, 25
+ summary = tdvp2.run(
+ mps, mpo, tdvp2.Options(dt=dt, n_steps=n_steps, max_bond=64, normalize=False)
+ )
+ evolved = mps_to_vector(summary.state, charges)
+ evolved = evolved / evolved.norm()
+ exact = exact_evolve(ham, psi0, dt * n_steps)
+ exact = exact / exact.norm()
+ assert 1.0 - abs(torch.vdot(exact, evolved)).item() < 1e-5
+
+ def test_error_is_bounded_and_non_increasing(self, spin_space):
+ """2-site TDVP's error is a manifold-projection error: it does NOT vanish as
+ dt→0 (it plateaus), but it must stay small and not GROW as dt shrinks."""
+ length = 6
+ mps0, mpo, charges, psi0 = _domain_wall(length, spin_space)
+ ham = dense_heisenberg(length, charges)
+ psi0 = psi0 / psi0.norm()
+
+ def infidelity(dt, n):
+ mps, _, _, _ = _domain_wall(length, spin_space)
+ tdvp2.run(mps, mpo, tdvp2.Options(dt=dt, n_steps=n, max_bond=64, normalize=False))
+ v = mps_to_vector(mps, charges)
+ v = v / v.norm()
+ ex = exact_evolve(ham, psi0, dt * n)
+ ex = ex / ex.norm()
+ return 1.0 - abs(torch.vdot(ex, v)).item()
+
+ coarse = infidelity(0.1, 5)
+ fine = infidelity(0.05, 10)
+ assert coarse < 1e-4 and fine < 1e-4 # both small
+ assert fine <= coarse * 1.5 # does not grow as dt shrinks
+
+
+# ---------------------------------------------------------------------------
+# Conservation laws
+# ---------------------------------------------------------------------------
+
+class TestConservation:
+ def test_norm_conserved_real_time(self, spin_space):
+ mps, mpo, _, _ = _domain_wall(6, spin_space)
+ summary = tdvp2.run(mps, mpo, tdvp2.Options(dt=0.05, n_steps=10, max_bond=64, normalize=False))
+ for norm in summary.norms:
+ assert abs(norm - 1.0) < 1e-9
+
+ def test_total_sz_conserved(self, spin_space):
+ mps, mpo, charges, psi0 = _domain_wall(6, spin_space)
+ sz_total = dense_total_sz(6, charges)
+ sz_before = (psi0.conj() @ sz_total @ psi0).real.item() / psi0.norm().item() ** 2
+ summary = tdvp2.run(mps, mpo, tdvp2.Options(dt=0.05, n_steps=10, max_bond=64))
+ vec = mps_to_vector(summary.state, charges)
+ sz_after = (vec.conj() @ sz_total @ vec).real.item() / vec.norm().item() ** 2
+ assert abs(sz_after - sz_before) < 1e-9
+
+ def test_imaginary_time_lowers_energy(self, spin_space):
+ length = 6
+ mps, mpo, charges, psi0 = _domain_wall(length, spin_space)
+ ham = dense_heisenberg(length, charges)
+ ground = torch.linalg.eigvalsh(ham)[0].item()
+ psi0 = psi0 / psi0.norm()
+ energy_before = (psi0.conj() @ ham @ psi0).real.item()
+ summary = tdvp2.run(
+ mps, mpo, tdvp2.Options(dt=0.05, n_steps=40, imaginary_time=True, max_bond=64)
+ )
+ vec = mps_to_vector(summary.state, charges)
+ vec = vec / vec.norm()
+ energy_after = (vec.conj() @ ham @ vec).real.item()
+ assert energy_after < energy_before
+ assert energy_after > ground - 1e-9
diff --git a/tests/algorithm/test_imaginary_time_groundstate.py b/tests/algorithm/test_imaginary_time_groundstate.py
new file mode 100644
index 0000000..b161732
--- /dev/null
+++ b/tests/algorithm/test_imaginary_time_groundstate.py
@@ -0,0 +1,139 @@
+# Copyright (C) 2025-2026 Changkai Zhang.
+#
+# This file is part of Alice project.
+#
+# Alice is free software: you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published
+# by the Free Software Foundation, either version 3 of the License,
+# or (at your option) any later version.
+#
+# Alice is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with Alice. If not, see .
+# Author of code: Madhav Menon.
+
+
+"""Cross-method imaginary-time ground-state convergence.
+
+The headline quantity of the BUG-vs-TDVP study is the **overlap error of the
+imaginary-time-cooled state with the exact ground state**. This module checks, on
+a small Heisenberg chain where the ground state is available by exact
+diagonalization, that *every* integrator under comparison cools a Néel product
+state toward that exact ground state:
+
+* bond_update_bug (``bond_update_bug``),
+* two-site TDVP (``tdvp2``).
+
+For each method the final state must have a small overlap error with the exact
+ground state, a near-degenerate energy, and clear cooling relative to the Néel
+start. This is the unit-level guard for the imaginary-time pipeline the full
+``L = 26`` campaign runs.
+"""
+
+from __future__ import annotations
+
+import pytest
+import torch
+from nicole import Index, Tensor, load_space
+
+from alice import build_hamiltonian, init_mps
+from alice.algorithm import tdvp2, bond_update_bug
+
+from tests.algorithm.bond_update_bug.conftest import (
+ dense_hamiltonian,
+ heisenberg_chain,
+ mps_to_vector,
+)
+
+# Imaginary-time schedule: dt matches the production setup; beta = dt * n_steps is
+# made long enough that, at full bond dimension (no truncation on L = 6), the only
+# residual is the O(dt^2) Strang/splitting bias. Kept modest so the test is fast.
+_DT = 0.05
+_N_STEPS = 200
+_LENGTH = 6
+
+# This unit test validates that the inverse-free BUG family converges to the exact
+# ground state in imaginary time. Two-site TDVP is the comparison baseline whose
+# imaginary-time instability (it stalls under truncation and blows up) is the very
+# phenomenon the study figure exhibits — so it is driven by the study harness and
+# its own test module, and is deliberately not asserted as a convergence invariant
+# here.
+_BUG_METHODS = ['bug']
+
+
+@pytest.fixture(autouse=True)
+def _isolate_cwd(tmp_path, monkeypatch):
+ monkeypatch.chdir(tmp_path)
+
+
+@pytest.fixture(scope='module')
+def spin_space():
+ return load_space('Spin', 'U1', {'J': 0.5})
+
+
+def _neel(length, spin_space):
+ """Build the full-phys Néel MPS plus the interaction list, MPO, and ED data."""
+ _, operators = spin_space
+ interactions, spc, _ = heisenberg_chain(length)
+ mpo = build_hamiltonian(interactions, length, spc)
+ charges = [sector.charge for sector in spc.sectors]
+ config = [0, 1] * (length // 2)
+ target = sum(charges[c] for c in config)
+ mps = init_mps(length, spc, operators, config=config, target_qn=target)
+ for i in range(mps.L):
+ core = mps[i]
+ full_phys = Index(core.indices[2].direction, core.indices[2].group, spc.sectors)
+ mps[i] = Tensor(
+ indices=(core.indices[0], core.indices[1], full_phys),
+ itags=core.itags,
+ data={key: block.clone() for key, block in core.data.items()},
+ dtype=core.dtype,
+ )
+ psi0 = mps_to_vector(mps, charges)
+ return mps, interactions, mpo, charges, psi0
+
+
+def _cool(method, mps, interactions, mpo):
+ """Run one method in imaginary time and return its evolved MPS state."""
+ if method == 'bug':
+ return bond_update_bug.run(
+ mps, interactions,
+ bond_update_bug.Options(dt=_DT, n_steps=_N_STEPS,
+ imaginary_time=True, max_bond=64),
+ ).state
+ if method == 'tdvp2':
+ return tdvp2.run(
+ mps, mpo,
+ tdvp2.Options(dt=_DT, n_steps=_N_STEPS, imaginary_time=True, max_bond=64),
+ ).state
+ raise ValueError(f"unknown method {method!r}")
+
+
+@pytest.mark.parametrize('method', _BUG_METHODS)
+def test_cools_neel_to_exact_ground_state(method, spin_space):
+ mps, interactions, mpo, charges, psi0 = _neel(_LENGTH, spin_space)
+ ham = dense_hamiltonian(interactions, _LENGTH, charges)
+ evals, evecs = torch.linalg.eigh(ham)
+ ground_energy = evals[0].item()
+ ground_vec = evecs[:, 0]
+
+ psi0 = psi0 / psi0.norm()
+ energy_before = (psi0.conj() @ ham @ psi0).real.item()
+ err_before = 1.0 - abs(torch.vdot(ground_vec, psi0)).item()
+
+ state = _cool(method, mps, interactions, mpo)
+ vec = mps_to_vector(state, charges)
+ vec = vec / vec.norm()
+ energy_after = (vec.conj() @ ham @ vec).real.item()
+ err_after = 1.0 - abs(torch.vdot(ground_vec, vec)).item()
+
+ # Variational lower bound, genuine cooling, and convergence to the exact GS.
+ assert energy_after > ground_energy - 1e-9, f"{method}: energy below ED ground state"
+ assert energy_after < energy_before - 1e-6, f"{method}: energy did not decrease"
+ assert err_after < err_before, f"{method}: did not cool toward ground state"
+ assert energy_after - ground_energy < 1e-2, f"{method}: energy not converged"
+ assert err_after < 1e-2, f"{method}: final overlap error {err_after:.2e} too large"
diff --git a/tests/algorithm/test_local_solvers.py b/tests/algorithm/test_local_solvers.py
new file mode 100644
index 0000000..d9a9f8a
--- /dev/null
+++ b/tests/algorithm/test_local_solvers.py
@@ -0,0 +1,158 @@
+# Copyright (C) 2025-2026 Changkai Zhang.
+#
+# This file is part of Alice project.
+#
+# Alice is free software: you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published
+# by the Free Software Foundation, either version 3 of the License,
+# or (at your option) any later version.
+#
+# Alice is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with Alice. If not, see .
+# Author of code: Madhav Menon.
+
+
+"""Pluggable local-solver tests for the (imaginary-time) discarded-projector BUGs.
+
+In imaginary time the local update ``y = exp(tau A) x`` is the exact flow of a
+linear ODE, so it may be computed by any stable integrator instead of the exact
+Krylov exponential. Both the bond_update_bug and the global
+``bond_update_bug`` exposes ``solver`` / ``solver_substeps`` for this. These tests
+check, end-to-end through the real symmetry-blocked tensor machinery, that:
+
+* the substepped integrators (``midpoint``/``rk4``/``trapezoid``) reproduce the exact
+ ``krylov`` evolution as ``solver_substeps`` grows (the screenshot's "increase n"),
+ validating the explicit RK actions *and* the implicit Crank–Nicolson GMRES solve;
+* every solver still cools a Néel state toward the exact ground state; and
+* bad solver names are rejected at ``Options`` construction.
+"""
+
+from __future__ import annotations
+
+import pytest
+import torch
+from nicole import Index, Tensor, load_space
+
+from alice import build_hamiltonian, build_interaction, init_mps
+from alice.algorithm import bond_update_bug
+from alice.algorithm.bond_update_bug._kernel.local_solvers import LOCAL_SOLVERS
+
+from tests.algorithm.bond_update_bug.conftest import (
+ dense_hamiltonian,
+ heisenberg_chain,
+ mps_to_vector,
+)
+
+_LENGTH = 6
+
+
+@pytest.fixture(autouse=True)
+def _isolate_cwd(tmp_path, monkeypatch):
+ monkeypatch.chdir(tmp_path)
+
+
+@pytest.fixture(scope='module')
+def spin_space():
+ return load_space('Spin', 'U1', {'J': 0.5})
+
+
+def _neel(length, spin_space):
+ _, operators = spin_space
+ interactions, spc, _ = heisenberg_chain(length)
+ mpo = build_hamiltonian(interactions, length, spc)
+ charges = [sector.charge for sector in spc.sectors]
+ config = [0, 1] * (length // 2)
+ target = sum(charges[c] for c in config)
+ mps = init_mps(length, spc, operators, config=config, target_qn=target)
+ for i in range(mps.L):
+ core = mps[i]
+ full_phys = Index(core.indices[2].direction, core.indices[2].group, spc.sectors)
+ mps[i] = Tensor(
+ indices=(core.indices[0], core.indices[1], full_phys),
+ itags=core.itags,
+ data={key: block.clone() for key, block in core.data.items()},
+ dtype=core.dtype,
+ )
+ return mps, interactions, mpo, charges
+
+
+def _two_site_state(spin_space, *, solver, substeps, n_steps=4, dt=0.05):
+ mps, interactions, _, charges = _neel(_LENGTH, spin_space)
+ state = bond_update_bug.run(
+ mps, interactions,
+ bond_update_bug.Options(solver=solver, solver_substeps=substeps,
+ dt=dt, n_steps=n_steps, imaginary_time=True, max_bond=64),
+ ).state
+ vec = mps_to_vector(state, charges)
+ return vec / vec.norm()
+
+
+def _overlap_err(a, b):
+ # Clamp at 0: when two states agree to machine precision, || can round to
+ # just above 1 and give a tiny negative "error".
+ return max(0.0, 1.0 - abs(torch.vdot(a, b)).item())
+
+
+# ---------------------------------------------------------------------------
+# Option validation
+# ---------------------------------------------------------------------------
+
+class TestSolverOptions:
+
+ def test_known_solvers(self):
+ assert set(LOCAL_SOLVERS) == {'krylov', 'midpoint', 'rk4', 'trapezoid'}
+
+ def test_default_is_krylov(self):
+ assert bond_update_bug.Options().solver == 'krylov'
+
+ @pytest.mark.parametrize('factory', [bond_update_bug.Options])
+ def test_unknown_solver_raises(self, factory):
+ with pytest.raises(ValueError, match='unknown local solver'):
+ factory(solver='euler')
+
+
+# ---------------------------------------------------------------------------
+# bond_update_bug: K/L/S solves
+# ---------------------------------------------------------------------------
+
+class TestTwoSiteSolvers:
+ """Substepped solvers reproduce the exact Krylov evolution as n grows."""
+
+ @pytest.mark.parametrize('solver', ['midpoint', 'rk4', 'trapezoid'])
+ def test_converges_to_krylov_with_substeps(self, solver, spin_space):
+ ref = _two_site_state(spin_space, solver='krylov', substeps=1, n_steps=3)
+ coarse = _two_site_state(spin_space, solver=solver, substeps=2, n_steps=3)
+ fine = _two_site_state(spin_space, solver=solver, substeps=10, n_steps=3)
+ err_coarse = _overlap_err(ref, coarse)
+ err_fine = _overlap_err(ref, fine)
+ # More substeps -> at least as close to the exact Krylov action (rk4 already
+ # hits machine precision at n=2, so allow equality at the FP floor), and tight.
+ assert err_fine <= err_coarse + 1e-12
+ assert err_fine < 1e-4, f"{solver}: err_fine {err_fine:.2e}"
+
+
+# ---------------------------------------------------------------------------
+# Every solver cools toward the ground state
+# ---------------------------------------------------------------------------
+
+class TestCoolsWithEverySolver:
+
+ @pytest.mark.slow
+ @pytest.mark.parametrize('solver', ['krylov', 'midpoint', 'rk4', 'trapezoid'])
+ def test_two_site_discarded_cools(self, solver, spin_space):
+ mps, interactions, _, charges = _neel(_LENGTH, spin_space)
+ ham = dense_hamiltonian(interactions, _LENGTH, charges)
+ evals, evecs = torch.linalg.eigh(ham)
+ ground_vec = evecs[:, 0]
+ psi0 = mps_to_vector(mps, charges)
+ psi0 = psi0 / psi0.norm()
+ err_before = _overlap_err(ground_vec, psi0)
+ vec = _two_site_state(spin_space, solver=solver, substeps=8, n_steps=100, dt=0.05)
+ err_after = _overlap_err(ground_vec, vec)
+ assert err_after < err_before
+ assert err_after < 5e-2, f"{solver}: final overlap error {err_after:.2e}"