Controlled Bond Expansion: randomized SVD approach

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")

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

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

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)