diff --git a/src/MPSKit.jl b/src/MPSKit.jl
index a4e1fe605..b29174216 100644
--- a/src/MPSKit.jl
+++ b/src/MPSKit.jl
@@ -37,7 +37,7 @@ export FiniteExcited, QuasiparticleAnsatz, ChepigaAnsatz, ChepigaAnsatz2
 export time_evolve, timestep, timestep!, make_time_mpo
 export TDVP, TDVP2, WI, WII, TaylorCluster
 export changebonds, changebonds!
-export VUMPSSvdCut, OptimalExpand, SvdCut, RandExpand
+export VUMPSSvdCut, OptimalExpand, SvdCut, RandExpand, SketchedExpand
 export propagator
 export DynamicalDMRG, NaiveInvert, Jeckelmann
 export exact_diagonalization, fidelity_susceptibility
@@ -156,6 +156,7 @@ include("algorithms/changebonds/optimalexpand.jl")
 include("algorithms/changebonds/vumpssvd.jl")
 include("algorithms/changebonds/svdcut.jl")
 include("algorithms/changebonds/randexpand.jl")
+include("algorithms/changebonds/sketchedexpand.jl")
 
 include("algorithms/timestep/tdvp.jl")
 include("algorithms/timestep/taylorcluster.jl")
diff --git a/src/algorithms/changebonds/changebonds.jl b/src/algorithms/changebonds/changebonds.jl
index de549ca24..3a5dbeb5e 100644
--- a/src/algorithms/changebonds/changebonds.jl
+++ b/src/algorithms/changebonds/changebonds.jl
@@ -12,27 +12,41 @@ See also: [`SvdCut`](@ref), [`RandExpand`](@ref), [`VUMPSSvdCut`](@ref), [`Optim
 function changebonds end
 function changebonds! end
 
+@doc """
+    changebond(site, dir, ψ, [H], alg, [envs]) -> ψ
+    changebond!(site, dir, ψ, [H], alg, [envs]) -> ψ
+
+Expand a single bond of `ψ` in place by adding directions orthogonal to the current state, keeping the state in mixed-canonical form around the enriched bond.
+The sweep direction `dir` is a `Val(:right)` or `Val(:left)` used for dispatch.
+For `Val(:right)` the bond `(site, site + 1)` is enriched on the right tensor (`ψ.AR[site + 1]`) with zero weight added at `ψ.AC[site]`, so that a subsequent single-site optimization of `site` sees the new directions;
+for `Val(:left)` the mirror is applied to bond `(site - 1, site)`.
+
+See also [`changebonds`](@ref), [`changebonds!`](@ref).
+""" changebond, changebond!
+function changebond end
+function changebond! end
+
 _expand(ψ, AL′, AR′) = _expand!(copy(ψ), AL′, AR′)
 function _expand!(ψ::InfiniteMPS, AL′::PeriodicVector, AR′::PeriodicVector)
     for i in 1:length(ψ)
         # update AL: add vectors, make room for new vectors:
         # AL -> [AL expansion; 0 0]
         al′ = _transpose_tail(catdomain(ψ.AL[i], AL′[i]))
-        lz = zerovector!(similar(al′, _lastspace(AL′[i - 1])' ← domain(al′)))
-        ψ.AL[i] = _transpose_front(catcodomain(al′, lz))
+        al_space = (codomain(al′)[1] ⊕ _lastspace(AL′[i - 1])') ← domain(al′)
+        ψ.AL[i] = _transpose_front(absorb!(zerovector!(similar(al′, al_space)), al′))
 
         # update AR: add vectors, make room for new vectors:
         # AR -> [AR 0; expansion 0]
         ar′ = _transpose_front(catcodomain(_transpose_tail(ψ.AR[i + 1]), AR′[i + 1]))
-        rz = zerovector!(similar(ar′, codomain(ar′) ← _firstspace(AR′[i + 2])))
-        ψ.AR[i + 1] = catdomain(ar′, rz)
+        ar_space = codomain(ar′) ← (domain(ar′)[1] ⊕ _firstspace(AR′[i + 2]))
+        ψ.AR[i + 1] = absorb!(zerovector!(similar(ar′, ar_space)), ar′)
 
         # update C: add vectors, make room for new vectors:
         # C -> [C 0; 0 expansion]
-        l = zerovector!(similar(ψ.C[i], codomain(ψ.C[i]) ← _firstspace(AR′[i + 1])))
-        ψ.C[i] = catdomain(ψ.C[i], l)
-        r = zerovector!(similar(ψ.C[i], _lastspace(AL′[i])' ← domain(ψ.C[i])))
-        ψ.C[i] = catcodomain(ψ.C[i], r)
+        c_dom = codomain(ψ.C[i]) ← (domain(ψ.C[i])[1] ⊕ _firstspace(AR′[i + 1]))
+        ψ.C[i] = absorb!(zerovector!(similar(ψ.C[i], c_dom)), ψ.C[i])
+        c_cod = (codomain(ψ.C[i])[1] ⊕ _lastspace(AL′[i])') ← domain(ψ.C[i])
+        ψ.C[i] = absorb!(zerovector!(similar(ψ.C[i], c_cod)), ψ.C[i])
 
         # update AC: recalculate
         ψ.AC[i] = ψ.AL[i] * ψ.C[i]
@@ -45,3 +59,7 @@ function _expand!(ψ::MultilineMPS, AL′::PeriodicMatrix, AR′::PeriodicMatrix
     end
     return ψ
 end
+
+function changebond(site::Int, dir::Val, ψ::AbstractFiniteMPS, H, alg::Algorithm, envs)
+    return changebond!(site, dir, copy(ψ), H, alg, envs)
+end
diff --git a/src/algorithms/changebonds/optimalexpand.jl b/src/algorithms/changebonds/optimalexpand.jl
index 194808704..73da47899 100644
--- a/src/algorithms/changebonds/optimalexpand.jl
+++ b/src/algorithms/changebonds/optimalexpand.jl
@@ -5,6 +5,9 @@ An algorithm that expands the given mps as described in
 [Zauner-Stauber et al. Phys. Rev. B 97 (2018)](@cite zauner-stauber2018), by selecting the
 dominant contributions of a two-site updated MPS tensor, orthogonal to the original ψ.
 
+The expansion does not alter the state: the added directions are connected through a zero block,
+so that the expanded state is identical to the original one (as required for e.g. TDVP).
+
 !!! note
     [`changebonds!`](@ref) is only defined for `FiniteMPS`, and modifies both the state and its environment.
 
@@ -81,37 +84,71 @@ function changebonds(ψ::MultilineMPS, H, alg::OptimalExpand, envs = environment
     return newψ, envs
 end
 
-function changebonds(ψ::AbstractFiniteMPS, H, alg::OptimalExpand, envs = environments(ψ, H, ψ))
-    return changebonds!(copy(ψ), H, alg, envs)
-end
 
-function changebonds!(ψ::AbstractFiniteMPS, H, alg::OptimalExpand, envs = environments(ψ, H, ψ))
-    #inspired by the infinite mps algorithm, alternative is to use https://arxiv.org/pdf/1501.05504.pdf
+# Finite system
+# -------------
+function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::OptimalExpand, envs; normalize::Bool = true)
+    bond = site
+    left = ψ.AC[site]
+    right = ψ.AR[site + 1]
+    NL = left_null(left)
+    NR = right_null!(_transpose_tail(right; copy = true))
+
+    # two-site update from the projected effective Hamiltonian
+    AC2 = AC2_projection(bond, ψ, H, ψ, envs)
+
+    # select the dominant directions in the complement of the current state
+    g2 = adjoint(NL) * AC2 * adjoint(NR)
+    _, _, Vᴴ = svd_trunc!(normalize!(g2); trunc = alg.trscheme, alg = alg.alg_svd)
+
+    # optimal vectors at site+1
+    ar_re = Vᴴ * NR
+    # embed `left` into the enlarged domain (zero weight in the new directions), leaving the state
+    # unchanged
+    nal_space = codomain(left) ← (only(domain(left)) ⊕ space(Vᴴ, 1))
+    nal, nc = left_gauge(absorb!(zerovector!(similar(left, nal_space)), left))
+    nar = _transpose_front(catcodomain(_transpose_tail(right), ar_re))
+
+    normalize && normalize!(nc)
+    ψ.AC[site] = (nal, nc)
+    ψ.AC[site + 1] = (nc, nar)
+    return ψ
+end
+function changebond!(site::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, alg::OptimalExpand, envs; normalize::Bool = true)
+    bond = site - 1
+    left = ψ.AL[site - 1]
+    right = ψ.AC[site]
+    NL = left_null(left)
+    NR = right_null!(_transpose_tail(right; copy = true))
+
+    # two-site update from the projected effective Hamiltonian
+    AC2 = AC2_projection(bond, ψ, H, ψ, envs)
+
+    # select the dominant directions in the complement of the current state
+    g2 = adjoint(NL) * AC2 * adjoint(NR)
+    U, _, _ = svd_trunc!(normalize!(g2); trunc = alg.trscheme, alg = alg.alg_svd)
+
+    # optimal vectors at site-1
+    Q = NL * U
+    right_tail = _transpose_tail(right)
+    # embed `_transpose_tail(right)` into the enlarged codomain (zero weight in the new
+    # directions), leaving the state unchanged
+    nc_space = (codomain(right_tail)[1] ⊕ space(Q, 3)') ← domain(right_tail)
+    nc, Qr = lq_compact!(absorb!(zerovector!(similar(right_tail, nc_space)), right_tail))
+    AL_exp = catdomain(left, Q)
+
+    normalize && normalize!(nc)
+    ψ.AC[site] = (nc, _transpose_front(Qr))
+    ψ.AC[site - 1] = (AL_exp, nc)
+    return ψ
+end
 
-    #the idea is that we always want to expand the state in such a way that there are zeros at site i
-    #but "optimal vectors" at site i+1
-    #so during optimization of site i, you have access to these optimal vectors :)
+changebonds(ψ::AbstractFiniteMPS, H, alg::OptimalExpand, envs = environments(ψ, H, ψ)) =
+    changebonds!(copy(ψ), H, alg, envs)
 
+function changebonds!(ψ::AbstractFiniteMPS, H, alg::OptimalExpand, envs = environments(ψ, H, ψ))
     for i in 1:(length(ψ) - 1)
-        AC2 = AC2_projection(i, ψ, H, ψ, envs)
-
-        #Calculate nullspaces for left and right
-        NL = left_null(ψ.AC[i])
-        NR = right_null!(_transpose_tail(ψ.AR[i + 1]; copy = true))
-
-        #Use this nullspaces and SVD decomposition to determine the optimal expansion space
-        intermediate = normalize!(adjoint(NL) * AC2 * adjoint(NR))
-        _, _, V, = svd_trunc!(intermediate; trunc = alg.trscheme, alg = alg.alg_svd)
-
-        ar_re = V * NR
-        ar_le = zerovector!(similar(ar_re, codomain(ψ.AC[i]) ← space(V, 1)))
-
-        nal, nc = qr_compact!(catdomain(ψ.AC[i], ar_le))
-        nar = _transpose_front(catcodomain(_transpose_tail(ψ.AR[i + 1]), ar_re))
-
-        ψ.AC[i] = (nal, nc)
-        ψ.AC[i + 1] = (nc, nar)
+        changebond!(i, Val(:right), ψ, H, alg, envs)
     end
-
-    return (ψ, envs)
+    return ψ, envs
 end
diff --git a/src/algorithms/changebonds/randexpand.jl b/src/algorithms/changebonds/randexpand.jl
index 346755e4c..183830e36 100644
--- a/src/algorithms/changebonds/randexpand.jl
+++ b/src/algorithms/changebonds/randexpand.jl
@@ -68,30 +68,69 @@ end
 
 
 function changebonds!(ψ::AbstractFiniteMPS, alg::RandExpand)
+    # the expansion directions are sampled from a randomized two-site update, so no Hamiltonian
+    # or environments are required
     for i in 1:(length(ψ) - 1)
-        AC2 = randomize!(_transpose_front(ψ.AC[i]) * _transpose_tail(ψ.AR[i + 1]))
-
-        #Calculate nullspaces for left and right
-        NL = left_null(ψ.AC[i])
-        NR = right_null!(_transpose_tail(ψ.AR[i + 1]; copy = true))
-
-        #Use this nullspaces and SVD decomposition to determine the optimal expansion space
-        intermediate = normalize!(adjoint(NL) * AC2 * adjoint(NR))
-        _, _, Vᴴ = svd_trunc!(intermediate; trunc = alg.trscheme, alg = alg.alg_svd)
-
-        ar_re = Vᴴ * NR
-        ar_le = zerovector!(similar(ar_re, codomain(ψ.AC[i]) ← space(Vᴴ, 1)))
-
-        nal, nc = qr_compact!(catdomain(ψ.AC[i], ar_le))
-        nar = _transpose_front(catcodomain(_transpose_tail(ψ.AR[i + 1]), ar_re))
-
-        ψ.AC[i] = (nal, nc)
-        ψ.AC[i + 1] = (nc, nar)
+        changebond!(i, Val(:right), ψ, nothing, alg, nothing)
     end
 
     return normalize!(ψ)
 end
 
+function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::RandExpand, envs; normalize::Bool = true)
+    bond = site
+    left = ψ.AC[site]
+    right = ψ.AR[site + 1]
+    NL = left_null(left)
+    NR = right_null!(_transpose_tail(right; copy = true))
+
+    # randomized two-site update; H and envs are unused
+    ac2 = randomize!(AC2(ψ, bond))
+
+    # select the dominant directions in the complement of the current state
+    g2 = adjoint(NL) * ac2 * adjoint(NR)
+    _, _, Vᴴ = svd_trunc!(normalize!(g2); trunc = alg.trscheme, alg = alg.alg_svd)
+
+    # optimal vectors at site+1, zero weight at site
+    ar_re = Vᴴ * NR
+    # embed `left` into the enlarged domain (zero weight in the new directions)
+    nal_space = codomain(left) ← (only(domain(left)) ⊕ space(Vᴴ, 1))
+    nal, nc = left_gauge(absorb!(zerovector!(similar(left, nal_space)), left))
+    nar = _transpose_front(catcodomain(_transpose_tail(right), ar_re))
+
+    normalize && normalize!(nc)
+    ψ.AC[site] = (nal, nc)
+    ψ.AC[site + 1] = (nc, nar)
+    return ψ
+end
+function changebond!(site::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, alg::RandExpand, envs; normalize::Bool = true)
+    bond = site - 1
+    left = ψ.AL[site - 1]
+    right = ψ.AC[site]
+    NL = left_null(left)
+    NR = right_null!(_transpose_tail(right; copy = true))
+
+    # randomized two-site update; H and envs are unused
+    ac2 = randomize!(AC2(ψ, bond))
+
+    # select the dominant directions in the complement of the current state
+    g2 = adjoint(NL) * ac2 * adjoint(NR)
+    U, _, _ = svd_trunc!(normalize!(g2); trunc = alg.trscheme, alg = alg.alg_svd)
+
+    # optimal vectors at site-1, zero weight at site
+    Q = NL * U
+    # embed `_transpose_tail(right)` into the enlarged codomain (zero weight in the new directions)
+    right_tail = _transpose_tail(right)
+    nc_space = (codomain(right_tail)[1] ⊕ space(Q, 3)') ← domain(right_tail)
+    nc, Qr = lq_compact!(absorb!(zerovector!(similar(right_tail, nc_space)), right_tail))
+    AL_exp = catdomain(left, Q)
+
+    normalize && normalize!(nc)
+    ψ.AC[site] = (nc, _transpose_front(Qr))
+    ψ.AC[site - 1] = (AL_exp, nc)
+    return ψ
+end
+
 """
     sample_space(V, strategy)
 
diff --git a/src/algorithms/changebonds/sketchedexpand.jl b/src/algorithms/changebonds/sketchedexpand.jl
new file mode 100644
index 000000000..9e89b3b81
--- /dev/null
+++ b/src/algorithms/changebonds/sketchedexpand.jl
@@ -0,0 +1,141 @@
+"""
+$(TYPEDEF)
+
+An algorithm that expands the bond dimension like [`OptimalExpand`](@ref) — selecting the
+dominant directions of the projected two-site update orthogonal to the current state — but at
+single-site cost using the randomized "shrewd selection" of Controlled Bond Expansion
+(Gleis et al. Phys. Rev. Lett. 130, 246402 (2023)). A random sketch of the orthogonal complement
+is folded into the effective environment, collapsing the large bond before the two-site update is
+ever formed, and the dominant directions are read off a small singular value decomposition.
+
+The state-preserving behaviour matches [`OptimalExpand`](@ref).
+
+!!! note
+    Only defined for `FiniteMPS` (through [`changebond!`](@ref)), so it can be used standalone or
+    as the `alg_expand` strategy of [`DMRG`](@ref). The reported `ϵ_2site` is a randomized
+    estimate, and the folded application does not exploit `JordanMPO` sparsity.
+
+## Fields
+
+$(TYPEDFIELDS)
+"""
+@kwdef struct SketchedExpand{S} <: Algorithm
+    "algorithm used to orthonormalize the sketched complement (passed as the `alg` of `left_orth!`/`right_orth!`); `nothing` selects QR without oversampling and an SVD-based decomposition otherwise"
+    alg_orth::S = nothing
+
+    "algorithm used for truncating the expanded space"
+    trscheme::TruncationStrategy
+
+    "number of extra sketch columns drawn beyond the target rank (range-finder oversampling)"
+    oversampling::Int = 0
+end
+
+"""
+    sketch_space(V, alg::SketchedExpand) -> Vℓ, Vk
+
+The random-sketch space `Vℓ` drawn from the complement `V`, together with its oversampling-free
+target `Vk` (selected by `alg.trscheme`). `Vℓ` enlarges `Vk` by `alg.oversampling` extra
+directions (capped by `V`); the selection is truncated back to `Vk`.
+"""
+function sketch_space(V, alg::SketchedExpand)
+    Vk = sample_space(V, alg.trscheme)
+    alg.oversampling == 0 && return Vk, Vk
+    Vp = typeof(V)(c => min(dim(V, c), alg.oversampling) for c in sectors(V))
+    return infimum(V, Vk ⊕ Vp), Vk
+end
+
+# Finite system
+# -------------
+# Rather than forming the two-site update, a random sketch of the orthogonal complement is folded
+# into the effective environment, and the dominant directions are read off a small SVD of the
+# sketched gradient. The complement projectors act with the isometric MPS tensors directly, which
+# leaves the MPS gauge (and any incrementally-maintained environments) untouched.
+function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::SketchedExpand, envs; normalize::Bool = true)
+    left = ψ.AC[site]
+    right = ψ.AR[site + 1]
+    AL, _ = left_gauge(left)        # local left-isometric form
+    ARtt = _transpose_tail(right)   # AR is already right-isometric
+
+    # nothing to add when either complement is empty (e.g. edge bonds)
+    compL = fuse(codomain(AL)) ⊖ only(domain(AL))
+    compR = fuse(domain(ARtt)) ⊖ only(codomain(ARtt))
+    (dim(compL) == 0 || dim(compR) == 0) && return ψ
+
+    # random sketch of the left complement, folded into the left environment
+    Vℓ, Vk = sketch_space(compL, alg)
+    Ω = randisometry(scalartype(left), codomain(AL) ← Vℓ)
+    Q, _ = qr_compact!(project_complement!(Ω, AL))
+    GL = leftenv(envs, site, ψ) * TransferMatrix(left, H[site], Q)
+    Hac = MPO_AC_Hamiltonian(GL, H[site + 1], rightenv(envs, site + 1, ψ))
+    Y = Hac * right
+
+    # orthonormalize the sketched right complement, truncating away the oversampling padding
+    B = project_complement_right!(_transpose_tail(Y), ARtt)
+    if dim(Vℓ) == dim(Vk)        # no oversampling: a plain QR/LQ suffices
+        _, ar_re = right_orth!(B; alg = alg.alg_orth)
+    else
+        _, ar_re = right_orth!(B; trunc = truncspace(Vk), alg = alg.alg_orth)
+    end
+
+    # optimal vectors at site+1
+    # embed `left` into the enlarged domain (zero weight in the new directions), leaving the state
+    # unchanged
+    nal_space = codomain(left) ← (only(domain(left)) ⊕ space(ar_re, 1))
+    nal, nc = left_gauge(absorb!(zerovector!(similar(left, nal_space)), left))
+    nar = _transpose_front(catcodomain(_transpose_tail(right), ar_re))
+
+    normalize && normalize!(nc)
+    ψ.AC[site] = (nal, nc)
+    ψ.AC[site + 1] = (nc, nar)
+    return ψ
+end
+function changebond!(site::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, alg::SketchedExpand, envs; normalize::Bool = true)
+    left = ψ.AL[site - 1]
+    right = ψ.AC[site]
+    _, ARtt = right_orth!(_transpose_tail(right; copy = true); trunc = notrunc())  # local right-isometric form
+
+    # nothing to add when either complement is empty (e.g. edge bonds)
+    compL = fuse(codomain(left)) ⊖ only(domain(left))
+    compR = fuse(domain(ARtt)) ⊖ only(codomain(ARtt))
+    (dim(compL) == 0 || dim(compR) == 0) && return ψ
+
+    # random sketch of the right complement, folded into the right environment
+    Vℓ, Vk = sketch_space(compR, alg)
+    Ω = adjoint(randisometry(scalartype(right), domain(ARtt) ← Vℓ))
+    _, Qr_o = lq_compact!(project_complement_right!(Ω, ARtt))
+    Qr = _transpose_front(Qr_o)
+    GR = TransferMatrix(right, H[site], Qr) * rightenv(envs, site, ψ)
+    Hac = MPO_AC_Hamiltonian(leftenv(envs, site - 1, ψ), H[site - 1], GR)
+    Y = Hac * left
+
+    # orthonormalize the sketched left complement, truncating away the oversampling padding
+    B = project_complement!(Y, left)
+    if dim(Vℓ) == dim(Vk)        # no oversampling: a plain QR/LQ suffices
+        Q, _ = left_orth!(B; alg = alg.alg_orth)
+    else
+        Q, _ = left_orth!(B; trunc = truncspace(Vk), alg = alg.alg_orth)
+    end
+
+    # optimal vectors at site-1
+    right_tail = _transpose_tail(right)
+    # embed `_transpose_tail(right)` into the enlarged codomain (zero weight in the new
+    # directions), leaving the state unchanged
+    nc_space = (codomain(right_tail)[1] ⊕ space(Q, 3)') ← domain(right_tail)
+    nc, Qr2 = lq_compact!(absorb!(zerovector!(similar(right_tail, nc_space)), right_tail))
+    AL_exp = catdomain(left, Q)
+
+    normalize && normalize!(nc)
+    ψ.AC[site] = (nc, _transpose_front(Qr2))
+    ψ.AC[site - 1] = (AL_exp, nc)
+    return ψ
+end
+
+changebonds(ψ::AbstractFiniteMPS, H, alg::SketchedExpand, envs = environments(ψ, H, ψ)) =
+    changebonds!(copy(ψ), H, alg, envs)
+
+function changebonds!(ψ::AbstractFiniteMPS, H, alg::SketchedExpand, envs = environments(ψ, H, ψ))
+    for i in 1:(length(ψ) - 1)
+        changebond!(i, Val(:right), ψ, H, alg, envs)
+    end
+    return ψ, envs
+end
diff --git a/src/algorithms/groundstate/dmrg.jl b/src/algorithms/groundstate/dmrg.jl
index ca8b09d18..61f3097a0 100644
--- a/src/algorithms/groundstate/dmrg.jl
+++ b/src/algorithms/groundstate/dmrg.jl
@@ -1,13 +1,31 @@
+# whether the gauge algorithm truncates the bond (SVD-based) or preserves it (QR-based, a plain
+# center-move as in textbook single-site DMRG)
+_truncates(::MatrixAlgebraKit.AbstractAlgorithm) = false
+_truncates(::MatrixAlgebraKit.TruncatedAlgorithm) = true
+
 """
 $(TYPEDEF)
 
-Single-site DMRG algorithm for finding the dominant eigenvector.
+Density Matrix Renormalization Group algorithm for finding the dominant eigenvector.
+
+Each site update is, in order: (1) an optional bond expansion (`alg_expand`), (2) a single-site
+eigensolve, and (3) a gauge step (`alg_gauge`). With the defaults (`alg_expand = nothing` and a
+non-truncating QR `alg_gauge`) this is textbook single-site DMRG, which cannot change the bond
+dimension. Setting `alg_expand` to a bond-expansion algorithm (e.g. [`OptimalExpand`](@ref),
+[`RandExpand`](@ref), [`SketchedExpand`](@ref)) enriches the bond with directions orthogonal to
+the current state ahead of each eigensolve, recovering Controlled Bond Expansion (CBE) DMRG; a
+truncating `alg_gauge` is then desirable to cut the enlarged bond back down.
+
+The gauge algorithm is selected in the keyword constructor from the `trscheme` argument: when it
+is `notrunc()` the gauge is a QR decomposition (`alg_orth`, [`Householder`](@extref
+MatrixAlgebraKit.Householder) by default), otherwise it is a truncated SVD (`alg_svd` with the
+given `trscheme`).
 
 ## Fields
 
 $(TYPEDFIELDS)
 """
-struct DMRG{A, F} <: Algorithm
+struct DMRG{A, F, E, G} <: Algorithm
     "tolerance for convergence criterium"
     tol::Float64
 
@@ -22,14 +40,28 @@ struct DMRG{A, F} <: Algorithm
 
     "callback function applied after each iteration, of signature `finalize(iter, ψ, H, envs) -> ψ, envs`"
     finalize::F
+
+    "algorithm used to expand the bond ahead of each local update, or `nothing` for none"
+    alg_expand::E
+
+    "factorization used for the post-update gauge: a QR algorithm (no truncation) or a truncated SVD"
+    alg_gauge::G
 end
 function DMRG(;
         tol = Defaults.tol, maxiter = Defaults.maxiter, alg_eigsolve = (;),
-        verbosity = Defaults.verbosity, finalize = Defaults._finalize
+        verbosity = Defaults.verbosity, finalize = Defaults._finalize,
+        alg_expand = nothing, trscheme = notrunc(),
+        alg_svd = Defaults.alg_svd(), alg_orth = Defaults.alg_orth()
     )
     alg_eigsolve′ = alg_eigsolve isa NamedTuple ? Defaults.alg_eigsolve(; alg_eigsolve...) :
         alg_eigsolve
-    return DMRG(tol, maxiter, verbosity, alg_eigsolve′, finalize)
+    # a no-truncation `trscheme` selects a (bond-preserving) QR gauge, anything else a truncated SVD
+    alg_gauge = trscheme isa MatrixAlgebraKit.NoTruncation ? alg_orth :
+        MatrixAlgebraKit.TruncatedAlgorithm(alg_svd, trscheme)
+    if !isnothing(alg_expand) && !_truncates(alg_gauge)
+        @warn "DMRG with `alg_expand` but no truncation (`trscheme = notrunc()`): the bond dimension will grow unboundedly each sweep."
+    end
+    return DMRG(tol, maxiter, verbosity, alg_eigsolve′, finalize, alg_expand, alg_gauge)
 end
 
 function find_groundstate!(ψ::AbstractFiniteMPS, H, alg::DMRG, envs = environments(ψ, H, ψ))
@@ -46,14 +78,44 @@ function find_groundstate!(ψ::AbstractFiniteMPS, H, alg::DMRG, envs = environme
 
             zerovector!(ϵs)
             @timeit timeroutput "sweep" begin
-                for pos in [1:(length(ψ) - 1); length(ψ):-1:2]
-                    local vec
+                # left-to-right
+                for pos in 1:(length(ψ) - 1)
+                    local AC′
+                    # convergence: pre-expansion single-site Galerkin error
+                    ϵs[pos] = max(ϵs[pos], calc_galerkin(pos, ψ, H, ψ, envs))
+
+                    # 1. expand
+                    isnothing(alg.alg_expand) ||
+                        @timeit timeroutput "expand" changebond!(pos, Val(:right), ψ, H, alg.alg_expand, envs)
+
+                    # 2. local update
                     @timeit timeroutput "AC_eigsolve" begin
-                        h = AC_hamiltonian(pos, ψ, H, ψ, envs)
-                        _, vec = fixedpoint(h, ψ.AC[pos], :SR, alg_eigsolve)
+                        Hac = AC_hamiltonian(pos, ψ, H, ψ, envs)
+                        _, AC′ = fixedpoint(Hac, ψ.AC[pos], :SR, alg_eigsolve)
                     end
+
+                    # 3. gauge (QR center-move or truncated SVD, selected by `alg_gauge`)
+                    @timeit timeroutput "gauge" left_gauge!(ψ, pos, AC′, alg.alg_gauge; normalize = true)
+                end
+
+                # right-to-left
+                for pos in length(ψ):-1:2
+                    local AC′
+                    # convergence: pre-expansion single-site Galerkin error
                     ϵs[pos] = max(ϵs[pos], calc_galerkin(pos, ψ, H, ψ, envs))
-                    @timeit timeroutput "AC_update" ψ.AC[pos] = vec
+
+                    # 1. expand
+                    isnothing(alg.alg_expand) ||
+                        @timeit timeroutput "expand" changebond!(pos, Val(:left), ψ, H, alg.alg_expand, envs)
+
+                    # 2. local update
+                    @timeit timeroutput "AC_eigsolve" begin
+                        Hac = AC_hamiltonian(pos, ψ, H, ψ, envs)
+                        _, AC′ = fixedpoint(Hac, ψ.AC[pos], :SR, alg_eigsolve)
+                    end
+
+                    # 3. gauge (QR center-move or truncated SVD, selected by `alg_gauge`)
+                    @timeit timeroutput "gauge" right_gauge!(ψ, pos, AC′, alg.alg_gauge; normalize = true)
                 end
             end
             ϵ = maximum(ϵs)
@@ -78,6 +140,7 @@ function find_groundstate!(ψ::AbstractFiniteMPS, H, alg::DMRG, envs = environme
     return ψ, envs, ϵ
 end
 
+
 """
 $(TYPEDEF)
 
diff --git a/src/algorithms/timestep/tdvp.jl b/src/algorithms/timestep/tdvp.jl
index a6fc56c2f..993782e54 100644
--- a/src/algorithms/timestep/tdvp.jl
+++ b/src/algorithms/timestep/tdvp.jl
@@ -3,6 +3,19 @@ $(TYPEDEF)
 
 Single site MPS time-evolution algorithm based on the Time-Dependent Variational Principle.
 
+For finite MPS, setting `alg_expand` to a bond-expansion algorithm (e.g. [`OptimalExpand`](@ref),
+[`SketchedExpand`](@ref)) enriches the bond with directions orthogonal to the current state
+ahead of each local integration, recovering Controlled Bond Expansion (CBE) TDVP and lifting the
+fixed-bond limitation of plain single-site TDVP. A truncating `trscheme` is then required to cut
+the enlarged bond back down (selecting the truncated-SVD gauge). The expansion is
+state-preserving, as required for a consistent time evolution.
+
+!!! note
+    Real-time evolution preserves the norm: neither the bond expansion nor the truncation
+    renormalizes, so the state norm reflects the accumulated truncation error. Imaginary-time
+    evolution instead renormalizes at every step, like a ground-state search. CBE is only
+    available for finite MPS.
+
 ## Fields
 
 $(TYPEDFIELDS)
@@ -11,18 +24,38 @@ $(TYPEDFIELDS)
 
 * [Haegeman et al. Phys. Rev. Lett. 107 (2011)](@cite haegeman2011)
 """
-@kwdef struct TDVP{A, F} <: Algorithm
+struct TDVP{A, E, G, F} <: Algorithm
     "algorithm used in the exponential solvers"
-    integrator::A = Defaults.alg_expsolve()
+    integrator::A
 
     "tolerance for gauging algorithm"
-    tolgauge::Float64 = Defaults.tolgauge
+    tolgauge::Float64
 
     "maximal amount of iterations for gauging algorithm"
-    gaugemaxiter::Int = Defaults.maxiter
+    gaugemaxiter::Int
+
+    "algorithm used to expand the bond ahead of each local update, or `nothing` for none (finite CBE-TDVP)"
+    alg_expand::E
+
+    "factorization used for the post-update gauge: a QR algorithm (no truncation) or a truncated SVD"
+    alg_gauge::G
 
     "callback function applied after each iteration, of signature `finalize(iter, ψ, H, envs) -> ψ, envs`"
-    finalize::F = Defaults._finalize
+    finalize::F
+end
+function TDVP(;
+        integrator = Defaults.alg_expsolve(), tolgauge = Defaults.tolgauge,
+        gaugemaxiter = Defaults.maxiter, finalize = Defaults._finalize,
+        alg_expand = nothing, trscheme = notrunc(),
+        alg_svd = Defaults.alg_svd(), alg_orth = Defaults.alg_orth()
+    )
+    # a no-truncation `trscheme` selects a (bond-preserving) QR gauge, anything else a truncated SVD
+    alg_gauge = trscheme isa MatrixAlgebraKit.NoTruncation ? alg_orth :
+        MatrixAlgebraKit.TruncatedAlgorithm(alg_svd, trscheme)
+    if !isnothing(alg_expand) && !_truncates(alg_gauge)
+        @warn "TDVP with `alg_expand` but no truncation (`trscheme = notrunc()`): the bond dimension will grow unboundedly each sweep."
+    end
+    return TDVP(integrator, tolgauge, gaugemaxiter, alg_expand, alg_gauge, finalize)
 end
 
 function timestep(
@@ -92,9 +125,20 @@ function timestep!(
 
     # sweep left to right
     for i in 1:(length(ψ) - 1)
+        # 1. optionally expand the bond ahead of the local update (CBE)
+        isnothing(alg.alg_expand) ||
+            changebond!(i, Val(:right), ψ, H, alg.alg_expand, envs; normalize = imaginary_evolution)
+
+        # 2. evolve the (possibly expanded) center tensor forward
         Hac = AC_hamiltonian(i, ψ, H, ψ, envs)
-        ψ.AC[i] = integrate(Hac, ψ.AC[i], t, dt / 2, alg.integrator; imaginary_evolution)
+        AC = integrate(Hac, ψ.AC[i], t, dt / 2, alg.integrator; imaginary_evolution)
 
+        # 3. gauge: split AC -> AL[i], C[i] (QR center-move, or truncated SVD cutting the
+        #    enlarged bond back down) and move the center to i+1. Real-time evolution preserves
+        #    the norm; imaginary-time evolution renormalizes.
+        left_gauge!(ψ, i, AC, alg.alg_gauge; normalize = imaginary_evolution)
+
+        # 4. evolve the bond tensor backward
         Hc = C_hamiltonian(i, ψ, H, ψ, envs)
         ψ.C[i] = integrate(
             Hc, ψ.C[i], t + dt / 2, -dt / 2, alg.integrator;
@@ -108,12 +152,22 @@ function timestep!(
 
     # sweep right to left
     for i in length(ψ):-1:2
+        # 1. optionally expand the bond ahead of the local update (CBE)
+        isnothing(alg.alg_expand) ||
+            changebond!(i, Val(:left), ψ, H, alg.alg_expand, envs; normalize = imaginary_evolution)
+
+        # 2. evolve the (possibly expanded) center tensor forward
         Hac = AC_hamiltonian(i, ψ, H, ψ, envs)
-        ψ.AC[i] = integrate(
+        AC = integrate(
             Hac, ψ.AC[i], t + dt / 2, dt / 2, alg.integrator;
             imaginary_evolution
         )
 
+        # 3. gauge: split AC -> C[i-1], AR[i] and move the center to i-1 (real-time preserves the
+        #    norm; imaginary-time renormalizes)
+        right_gauge!(ψ, i, AC, alg.alg_gauge; normalize = imaginary_evolution)
+
+        # 4. evolve the bond tensor backward
         Hc = C_hamiltonian(i - 1, ψ, H, ψ, envs)
         ψ.C[i - 1] = integrate(
             Hc, ψ.C[i - 1], t + dt, -dt / 2, alg.integrator;
diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl
index a319baede..4b9df6e58 100644
--- a/src/algorithms/toolbox.jl
+++ b/src/algorithms/toolbox.jl
@@ -56,8 +56,7 @@ Concretely, this is the overlap of the current state with the single-site deriva
 function calc_galerkin(pos::Int, below, operator, above, envs)
     AC´ = AC_projection(pos, below, operator, above, envs)
     normalize!(AC´)
-    out = mul!(AC´, below.AL[pos], below.AL[pos]' * AC´, -1, +1)
-    return norm(out)
+    return norm(project_complement!(AC´, below.AL[pos]))
 end
 function calc_galerkin(pos::CartesianIndex{2}, below, operator, above, envs)
     row, col = Tuple(pos)
diff --git a/src/states/orthoview.jl b/src/states/orthoview.jl
index facc91a9c..bb2338340 100644
--- a/src/states/orthoview.jl
+++ b/src/states/orthoview.jl
@@ -60,8 +60,7 @@ function Base.getindex(v::CView{<:FiniteMPS, E}, i::Int)::E where {E}
         end
 
         for j in Iterators.reverse((i + 1):center)
-            v.parent.Cs[j], tmp = lq_compact!(_transpose_tail(v.parent.ACs[j]; copy = true); positive = true)
-            v.parent.ARs[j] = _transpose_front(tmp)
+            v.parent.Cs[j], v.parent.ARs[j] = right_gauge(v.parent.ACs[j])
             if j != i + 1 # last AC not needed
                 v.parent.ACs[j - 1] = _mul_tail(v.parent.ALs[j - 1], v.parent.Cs[j])
             end
@@ -76,7 +75,7 @@ function Base.getindex(v::CView{<:FiniteMPS, E}, i::Int)::E where {E}
         end
 
         for j in center:i
-            v.parent.ALs[j], v.parent.Cs[j + 1] = qr_compact(v.parent.ACs[j]; positive = true)
+            v.parent.ALs[j], v.parent.Cs[j + 1] = left_gauge(v.parent.ACs[j])
             if j != i # last AC not needed
                 v.parent.ACs[j + 1] = _mul_front(v.parent.Cs[j + 1], v.parent.ARs[j + 1])
             end
@@ -89,10 +88,9 @@ end
 function Base.setindex!(v::CView{<:FiniteMPS}, vec, i::Int)
     if ismissing(v.parent.Cs[i + 1])
         if !ismissing(v.parent.ALs[i])
-            v.parent.Cs[i + 1], temp = lq_compact!(_transpose_tail(v.parent.AC[i + 1]; copy = true); positive = true)
-            v.parent.ARs[i + 1] = _transpose_front(temp)
+            v.parent.Cs[i + 1], v.parent.ARs[i + 1] = right_gauge(v.parent.AC[i + 1])
         else
-            v.parent.ALs[i], v.parent.Cs[i + 1] = qr_compact(v.parent.AC[i]; positive = true)
+            v.parent.ALs[i], v.parent.Cs[i + 1] = left_gauge(v.parent.AC[i])
         end
     end
 
@@ -221,3 +219,64 @@ function Base.checkbounds(
         checkbounds(Bool, CView(first(psi.parent.data)), b)
     end
 end
+
+# Gauging routines
+# ----------------
+@doc """
+    left_gauge(AC, [alg]) -> AL, C
+    right_gauge(AC, [alg]) -> C, AR
+
+Factor an updated center MPS tensor `AC` into left- or right-canonical form, `AC ≈ AL * C`
+(left, with `AL` left-isometric) or `AC ≈ C * AR` (right, with `AR` right-isometric), without
+modifying `AC`. `right_gauge` handles the MPS leg permutation internally, so `AR` is returned in
+standard MPS-tensor form.
+
+`alg` selects the factorization and defaults to a (positive) QR/LQ center-move that preserves the virtual bond.
+Passing a [`TruncatedAlgorithm`](@extref MatrixAlgebraKit.TruncatedAlgorithm) instead performs a truncated SVD may shrink the bond.
+"""
+left_gauge
+@doc (@doc left_gauge) right_gauge
+
+left_gauge(AC) = left_gauge(AC, Defaults.alg_orth())
+left_gauge(AC, alg) = left_orth(AC; alg)
+left_gauge(AC, alg::MatrixAlgebraKit.TruncatedAlgorithm) =
+    left_orth(AC; alg = MatrixAlgebraKit.LeftOrthViaSVD(alg))
+
+right_gauge(AC) = right_gauge(AC, Defaults.alg_orth())
+function right_gauge(AC, alg)
+    C, AR = right_orth(_transpose_tail(AC); alg)
+    return C, _transpose_front(AR)
+end
+function right_gauge(AC, alg::MatrixAlgebraKit.TruncatedAlgorithm)
+    C, AR = right_orth(_transpose_tail(AC); alg = MatrixAlgebraKit.RightOrthViaSVD(alg))
+    return C, _transpose_front(AR)
+end
+
+@doc """
+    left_gauge!(ψ, pos, AC, [alg]; normalize = false) -> ψ
+    right_gauge!(ψ, pos, AC, [alg]; normalize = false) -> ψ
+
+Gauge an updated center tensor `AC` at site `pos` and install it into `ψ` in one step: factor
+`AC` with [`left_gauge`](@ref) / [`right_gauge`](@ref) and write the canonical tensors back,
+shifting the gauge center past `pos` (to the right for `left_gauge!`, to the left for
+`right_gauge!`). `alg` is forwarded to [`left_gauge`](@ref) / [`right_gauge`](@ref) and hence may
+be a [`TruncatedAlgorithm`](@extref MatrixAlgebraKit.TruncatedAlgorithm) to truncate the bond.
+
+By default the factors are installed as-is. Pass `normalize = true` to renormalize the bond
+tensor, so `ψ` stays normalized after a local update that changed its norm.
+"""
+left_gauge!
+@doc (@doc left_gauge!) right_gauge!
+
+function left_gauge!(ψ::AbstractFiniteMPS, pos::Int, AC, alg = Defaults.alg_orth(); normalize::Bool = false)
+    AL, C = left_gauge(AC, alg)
+    normalize && normalize!(C)
+    ψ.AC[pos] = (AL, C)
+    return ψ
+end
+function right_gauge!(ψ::AbstractFiniteMPS, pos::Int, AC, alg = Defaults.alg_orth(); normalize::Bool = false)
+    C, AR = right_gauge(AC, alg)
+    normalize && normalize!(C)
+    ψ.AC[pos] = (C, AR)
+    return ψ
+end
diff --git a/src/utility/utility.jl b/src/utility/utility.jl
index 670439523..83f67aadc 100644
--- a/src/utility/utility.jl
+++ b/src/utility/utility.jl
@@ -11,6 +11,26 @@ end
 _mul_front(C, A) = mul_front(C, A) # _transpose_front(C * _transpose_tail(A))
 _mul_tail(A, C) = mul_tail(A, C) # A * C
 
+"""
+    project_complement!(Y, X) -> Y
+
+In-place projection of `Y` onto the orthogonal complement of the range of the left-isometry `X`
+(`X' X = I`): `Y ← (I - X X') Y = Y - X (X' Y)`. `Y` is overwritten and returned.
+
+See also [`project_complement_right!`](@ref).
+"""
+project_complement!(Y, X) = mul!(Y, X, X' * Y, -1, +1)
+
+"""
+    project_complement_right!(Y, X) -> Y
+
+In-place projection of `Y` onto the orthogonal complement of the co-range of the right-isometry
+`X` (`X X' = I`): `Y ← Y (I - X' X) = Y - (Y X') X`. `Y` is overwritten and returned.
+
+See also [`project_complement!`](@ref).
+"""
+project_complement_right!(Y, X) = mul!(Y, Y * X', X, -1, +1)
+
 function _similar_tail(A::AbstractTensorMap)
     cod = _firstspace(A)
     dom = ⊗(dual(_lastspace(A)), dual.(space.(Ref(A), reverse(2:(numind(A) - 1))))...)
diff --git a/test/Project.toml b/test/Project.toml
index 562dd4339..b3c91b9ab 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -12,6 +12,7 @@ MatrixAlgebraKit = "6c742aac-3347-4629-af66-fc926824e5e4"
 ParallelTestRunner = "d3525ed8-44d0-4b2c-a655-542cee43accc"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec"
 TensorKitTensors = "41b62e7d-e9d1-4e23-942c-79a97adf954b"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/test/algorithms/changebonds.jl b/test/algorithms/changebonds.jl
index 0bf95b7cf..a2aea1c08 100644
--- a/test/algorithms/changebonds.jl
+++ b/test/algorithms/changebonds.jl
@@ -9,6 +9,7 @@ using Test, TestExtras
 using MPSKit
 using TensorKit
 using TensorKit: ℙ
+using Random
 
 spacelist = [(ℙ^4, ℙ^3), (Rep[SU₂](1 => 1), Rep[SU₂](0 => 2, 1 => 2, 2 => 1))]
 
@@ -84,11 +85,44 @@ end
     )
     @test dot(state, state_oe) ≈ 1 atol = 1.0e-8
 
+    state_se, _ = changebonds(
+        state, H,
+        SketchedExpand(; trscheme = truncrank(dim(Dspace) * dim(Dspace)), oversampling = 4)
+    )
+    @test dot(state, state_se) ≈ 1 atol = 1.0e-8
+
     state_tr = changebonds(state_oe, SvdCut(; trscheme = truncrank(dim(Dspace))))
 
     @test dim(left_virtualspace(state_tr, 5)) < dim(left_virtualspace(state_oe, 5))
 end
 
+# density-matrix-style MPS: each site carries two physical legs (ket ⊗ bra). The operator-free
+# bond-change algorithms (`RandExpand` expansion, `SvdCut` truncation) must handle the extra
+# physical leg. Operator-based expanders (`OptimalExpand`/`SketchedExpand`) are not covered here
+# because `FiniteMPOHamiltonian`/`FiniteMPO` only support a single physical leg per site.
+@testset "Density-matrix FiniteMPS $(spacetype(pcomp))" for (pcomp, Dspace) in [
+        (ℙ^2 ⊗ (ℙ^2)', ℙ^6),
+        (Rep[SU₂](1 // 2 => 1) ⊗ Rep[SU₂](1 // 2 => 1)', Rep[SU₂](0 => 4, 1 => 3)),
+    ]
+    Random.seed!(2468)
+    L = 8
+    maxbond(ψ) = maximum(i -> dim(left_virtualspace(ψ, i)), 2:length(ψ))
+
+    ψ = FiniteMPS(rand, ComplexF64, fill(pcomp, L), Dspace)
+    @test numind(ψ.AC[L ÷ 2]) == 4    # two physical legs + two virtual legs
+
+    # RandExpand grows the bond while preserving the state (norm-preserving expansion)
+    ψ_re = changebonds(ψ, RandExpand(; trscheme = truncrank(dim(Dspace) * 2)))
+    @test numind(ψ_re.AC[L ÷ 2]) == 4
+    @test abs(dot(ψ, ψ_re)) ≈ 1 atol = 1.0e-8
+    @test maxbond(ψ_re) > maxbond(ψ)
+
+    # SvdCut truncates the enlarged bond back down, leaving a normalized state
+    ψ_tr = changebonds(ψ_re, SvdCut(; trscheme = truncrank(dim(Dspace))))
+    @test maxbond(ψ_tr) < maxbond(ψ_re)
+    @test abs(dot(ψ_tr, ψ_tr)) ≈ 1 atol = 1.0e-8
+end
+
 @testset "MultilineMPS $(spacetype(pspace))" for (pspace, Dspace) in spacelist
     o = rand(ComplexF64, pspace * pspace, pspace * pspace)
     mpo = MultilineMPO(o)
diff --git a/test/algorithms/groundstate.jl b/test/algorithms/groundstate.jl
index cb1faa261..605d45814 100644
--- a/test/algorithms/groundstate.jl
+++ b/test/algorithms/groundstate.jl
@@ -9,6 +9,7 @@ using Test, TestExtras
 using MPSKit
 using TensorKit
 using TensorKit: ℙ
+using Random
 
 verbosity_full = 5
 verbosity_conv = 1
@@ -61,6 +62,59 @@ verbosity_conv = 1
         @test v < 1.0e-2
     end
 
+    @testset "CBEDMRG" begin
+        # start from a small bond so the bond expansion is exercised
+        ψ₀ = FiniteMPS(randn, ComplexF64, L, ℙ^2, ℙ^(D ÷ 2))
+        v₀ = variance(ψ₀, H)
+        expand = OptimalExpand(; trscheme = truncrank(D ÷ 2))
+        trscheme = truncrank(D)
+
+        # test logging
+        ψ, envs, δ = find_groundstate(
+            ψ₀, H, DMRG(; verbosity = verbosity_full, maxiter = 2, alg_expand = expand, trscheme)
+        )
+
+        ψ, envs, δ = find_groundstate(
+            ψ, H, DMRG(; verbosity = verbosity_conv, maxiter = 10, alg_expand = expand, trscheme), envs
+        )
+        v = variance(ψ, H)
+
+        # test using low variance
+        @test sum(δ) ≈ 0 atol = 1.0e-3
+        @test v < v₀
+        @test v < 1.0e-2
+        # the bond should have grown to the truncation target
+        @test dim(left_virtualspace(ψ, L ÷ 2)) == D
+    end
+
+    @testset "CBEDMRG (SketchedExpand)" begin
+        # randomized bond expansion at single-site cost. The sketch is redrawn every sweep, so an
+        # aggressive expansion (a large fraction of the bond) keeps the single-site Galerkin error
+        # noisy; a gentle per-sweep increment lets it converge like the deterministic expanders.
+        Random.seed!(1234)
+        ψ₀ = FiniteMPS(randn, ComplexF64, L, ℙ^2, ℙ^(D ÷ 2))
+        v₀ = variance(ψ₀, H)
+        expand = SketchedExpand(; trscheme = truncrank(2), oversampling = 4)
+        trscheme = truncrank(D)
+
+        # test logging
+        ψ, envs, δ = find_groundstate(
+            ψ₀, H, DMRG(; verbosity = verbosity_full, maxiter = 2, alg_expand = expand, trscheme)
+        )
+
+        ψ, envs, δ = find_groundstate(
+            ψ, H, DMRG(; verbosity = verbosity_conv, maxiter = 15, alg_expand = expand, trscheme), envs
+        )
+        v = variance(ψ, H)
+
+        # test using low variance
+        @test sum(δ) ≈ 0 atol = 1.0e-3
+        @test v < v₀
+        @test v < 1.0e-2
+        # the bond should have grown to the truncation target
+        @test dim(left_virtualspace(ψ, L ÷ 2)) == D
+    end
+
     @testset "GradientGrassmann" begin
         ψ₀ = FiniteMPS(randn, ComplexF64, 10, ℙ^2, ℙ^D)
         v₀ = variance(ψ₀, H)
diff --git a/test/algorithms/timestep.jl b/test/algorithms/timestep.jl
index b6a3eb8dc..a846b2c5c 100644
--- a/test/algorithms/timestep.jl
+++ b/test/algorithms/timestep.jl
@@ -9,6 +9,8 @@ using Test, TestExtras
 using MPSKit
 using TensorKit
 using TensorKit: ℙ
+using LinearAlgebra: norm
+using Random
 
 verbosity_full = 5
 verbosity_conv = 1
@@ -91,6 +93,70 @@ verbosity_conv = 1
     end
 end
 
+@testset "Finite CBE-TDVP" verbose = true begin
+    L = 10
+    H = force_planar(heisenberg_XXX(Float64, Trivial; spin = 1 // 2, L))
+    Dstart, Dcap, dt = 2, 16, 0.05
+
+    # controlled bond expansion lets single-site TDVP grow the bond; the evolution should stay
+    # unitary (norm-preserving) and energy-conserving while tracking the bond-adaptive TDVP2
+    # reference better than fixed-bond single-site TDVP
+    @testset "$(nameof(Exp))" for (Exp, kw) in
+        ((OptimalExpand, (;)), (SketchedExpand, (; oversampling = 4)))
+        Random.seed!(4)
+        ψ₀ = complex(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^Dstart))
+        alg = TDVP(; alg_expand = Exp(; trscheme = truncrank(Dstart), kw...), trscheme = truncrank(Dcap))
+        E₀ = real(expectation_value(ψ₀, H))
+
+        ref, cbe, plain = ψ₀, ψ₀, ψ₀
+        for _ in 1:6
+            ref, = timestep(ref, H, 0.0, dt, TDVP2(; trscheme = truncrank(Dcap)))
+            cbe, = timestep(cbe, H, 0.0, dt, alg)
+            plain, = timestep(plain, H, 0.0, dt, TDVP())   # stuck at Dstart
+        end
+
+        @test norm(cbe) ≈ 1 atol = 1.0e-6
+        @test real(expectation_value(cbe, H)) ≈ E₀ atol = 1.0e-2
+        @test dim(left_virtualspace(cbe, L ÷ 2)) > Dstart
+        @test dim(left_virtualspace(plain, L ÷ 2)) == Dstart
+        @test abs(dot(ref, cbe)) > abs(dot(ref, plain))
+    end
+
+    # the bond truncation must preserve the norm for real-time evolution (the norm reflects the
+    # discarded weight) and only renormalize for imaginary-time evolution
+    @testset "norm handling" begin
+        Random.seed!(6)
+        ψ₀ = complex(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^Dstart))
+        # a deliberately lossy cap so the truncation discards weight every step
+        lossy = TDVP(; alg_expand = OptimalExpand(; trscheme = truncrank(2)), trscheme = truncrank(2))
+
+        ψrt = ψ₀
+        for _ in 1:12
+            ψrt, = timestep(ψrt, H, 0.0, 0.5, lossy)            # real time
+        end
+        @test norm(ψrt) < 1 - 1.0e-3                            # truncation loss is not renormalized away
+
+        ψit = ψ₀
+        for _ in 1:12
+            ψit, = timestep(ψit, H, 0.0, 0.5, lossy; imaginary_evolution = true)  # no external normalize!
+        end
+        @test norm(ψit) ≈ 1 atol = 1.0e-6                       # imaginary-time renormalizes each step
+    end
+
+    @testset "imaginary-time lowers energy" begin
+        Random.seed!(5)
+        ψ₀ = complex(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^Dstart))
+        alg = TDVP(; alg_expand = OptimalExpand(; trscheme = truncrank(Dstart)), trscheme = truncrank(Dcap))
+        E₀ = real(expectation_value(ψ₀, H))
+        ψ = ψ₀
+        for _ in 1:8
+            ψ, = timestep(ψ, H, 0.0, 0.1, alg; imaginary_evolution = true)  # gauge renormalizes
+        end
+        @test real(expectation_value(ψ, H)) < E₀
+        @test dim(left_virtualspace(ψ, L ÷ 2)) > Dstart
+    end
+end
+
 @testset "time_evolve" verbose = true begin
     t_span = 0:0.1:0.1
     algs = [TDVP(), TDVP2(; trscheme = truncrank(10))]