styles, find_groundstate restructure, local subspace expansion, randperturbexpand

Author: AFeuerpfeil

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, RandPerturbedExpand
 export propagator
 export DynamicalDMRG, NaiveInvert, Jeckelmann
 export exact_diagonalization, fidelity_susceptibility
@@ -155,6 +155,7 @@ include("algorithms/changebonds/optimalexpand.jl")
 include("algorithms/changebonds/vumpssvd.jl")
 include("algorithms/changebonds/svdcut.jl")
 include("algorithms/changebonds/randexpand.jl")
+include("algorithms/changebonds/localexpand.jl")
 
 include("algorithms/timestep/tdvp.jl")
 include("algorithms/timestep/taylorcluster.jl")

src/algorithms/approximate/idmrg.jl

@@ -120,12 +120,9 @@ function approximate!(
 
             # update error
             ϵ = sum(zip(C_current, ψ.C[:, 0])) do (c1, c2)
-                smallest = infimum(_firstspace(c1), _firstspace(c2))
-                e1 = isometry(_firstspace(c1), smallest)
-                e2 = isometry(_firstspace(c2), smallest)
-                return norm(e2' * c2 * e2 - e1' * c1 * e1)
+                return bond_error(c1, c2)
             end
-
+            
             if ϵ < alg.tol
                 @infov 2 logfinish!(log, iter, ϵ)
                 break

src/algorithms/changebonds/changebonds.jl

@@ -20,7 +20,7 @@ function changebonds(
 end
 
 _expand(ψ, AL′, AR′) = _expand!(copy(ψ), AL′, AR′)
-function _expand!(ψ::InfiniteMPS, AL′::PeriodicVector, AR′::PeriodicVector)
+function _expand!(ψ, AL′::AbstractVector, AR′::AbstractVector)
     for i in 1:length(ψ)
         # update AL: add vectors, make room for new vectors:
         # AL -> [AL expansion; 0 0]
@@ -46,7 +46,7 @@ function _expand!(ψ::InfiniteMPS, AL′::PeriodicVector, AR′::PeriodicVector)
     end
     return normalize!(ψ)
 end
-function _expand!(ψ::MultilineMPS, AL′::PeriodicMatrix, AR′::PeriodicMatrix)
+function _expand!(ψ::MultilineMPS, AL′::AbstractMatrix, AR′::AbstractMatrix)
     for i in 1:size(ψ, 1)
         _expand!(ψ[i], AL′[i, :], AR′[i, :])
     end

src/algorithms/changebonds/localexpand.jl

@@ -0,0 +1,76 @@
+struct NoExpand <: Algorithm end
+
+function changebonds_left(AL, C, alg::RandPerturbedExpand)
+    AL = _expand_leftisometry(AL, alg)
+    C′ = _embed_left_space(AL, C, alg)
+    return AL, C′
+end
+function changebonds_left(AL, C1, C2, alg::RandPerturbedExpand)
+    AL = _expand_leftisometry(AL, alg)
+    C1′ = _embed_left_space(AL, C1, alg)
+    C2′ = _embed_left_space(AL, C2, alg)
+    return AL, C1′, C2′
+end
+function changebonds_right(C, AR, alg::RandPerturbedExpand)
+    AR = _expand_rightisometry(AR, alg)
+    C′ = _embed_right_space(C, AR, alg)
+    return C′, AR
+end
+function changebonds_right(C1, C2, AR, alg::RandPerturbedExpand)
+    AR = _expand_rightisometry(AR, alg)
+    C1′ = _embed_right_space(C1, AR, alg)
+    C2′ = _embed_right_space(C2, AR, alg)
+    return C1′, C2′, AR
+end
+
+function changebonds(AL, C, AR, alg::RandPerturbedExpand)
+    AL = _expand_leftisometry(AL, alg)
+    C = _embed_left_space(AL, C, alg)
+    AR = _expand_rightisometry(AR, alg)
+    C = _embed_right_space(C, AR, alg)
+    return AL, C, AR
+end
+
+function _expand_leftisometry(A::MPSTensor, alg)
+    VL = left_null(A)
+    V = sample_space(right_virtualspace(VL), alg.trscheme)
+    XL = randisometry(scalartype(VL), right_virtualspace(VL) ← V)
+    A = catdomain(A, VL * XL)
+    return A
+end
+function _expand_rightisometry(A::MPSTensor, alg)
+    return _transpose_front(_expand_rightisometry(_transpose_tail(A; copy = true), alg))
+end
+function _expand_rightisometry(AR_tail::AbstractTensorMap, alg)
+    VR = right_null(AR_tail)
+    V = sample_space(space(VR, 1), alg.trscheme)
+    XR = randisometry(scalartype(VR), space(VR, 1) ← V)
+    AR_tail = catcodomain(AR_tail, XR' * VR)
+    return AR_tail
+end
+
+function _embed_left_space(A::MPSTensor, C::MPSBondTensor, alg)
+    C′ = similar(C, right_virtualspace(A) ← right_virtualspace(C))
+    scale!(randn!(C′), alg.noisefactor)
+    C′ = TensorKit.absorb!(C′, C)
+    return C′
+end
+function _embed_left_space(A::MPSTensor, Anext::MPSTensor, alg)
+    Anext′ = similar(Anext, right_virtualspace(A) ⊗ physicalspace(Anext) ← right_virtualspace(Anext))
+    scale!(randn!(Anext′), alg.noisefactor)
+    Anext′ = TensorKit.absorb!(Anext′, Anext)
+    return Anext′
+end
+
+function _embed_right_space(C::MPSBondTensor, A::AbstractTensorMap, alg)
+    C′ = similar(C, left_virtualspace(C) ← space(A, 1))
+    scale!(randn!(C′), alg.noisefactor)
+    C′ = TensorKit.absorb!(C′, C)
+    return C′
+end
+function _embed_right_space(Anext::MPSTensor, A::AbstractTensorMap, alg)
+    Anext′ = similar(Anext, left_virtualspace(Anext) ⊗ physicalspace(Anext) ← space(A, 1))
+    scale!(randn!(Anext′), alg.noisefactor)
+    Anext′ = TensorKit.absorb!(Anext′, Anext)
+    return Anext′
+end

src/algorithms/changebonds/randexpand.jl

@@ -1,3 +1,4 @@
+
 """
 $(TYPEDEF)
 
@@ -6,9 +7,6 @@ orthogonal to the existing state. This means that additional directions are adde
 `AL` and `AR` that are contained in the nullspace of both. Note that this is happens in
 parallel, and therefore the expansion will never go beyond the local two-site subspace.
 
-The truncation strategy dictates the number of expanded states, by generating uniformly
-distributed weights for each state in the two-site space and truncating that.
-
 ## Fields
 
 $(TYPEDFIELDS)
@@ -35,8 +33,8 @@ function changebonds!(ψ::InfiniteMPS, alg::RandExpand)
         # obtain (orthogonal) directions as isometries in that direction
         XL = randisometry(scalartype(VL), right_virtualspace(VL) ← V)
         AL′[i] = VL * XL
-        XR = randisometry(scalartype(VR), space(VR, 1) ← V)
-        AR′[i + 1] = XR' * VR
+        XR = randisometry(storagetype(VL), space(VR, 1) ← V)
+        AR′[i + 1] = XR * VR
     end
 
     return _expand!(ψ, AL′, AR′)
@@ -87,3 +85,92 @@ function sample_space(V, strategy)
     ind = MatrixAlgebraKit.findtruncated(S, strategy)
     return TensorKit.Factorizations.truncate_space(V, ind)
 end
+
+
+
+"""
+$(TYPEDEF)
+
+An algorithm that expands the bond dimension by adding random unitary vectors that are
+orthogonal to the existing state, in a sweeping fashion. Additionally, some random noise
+is added to the state in order for it to remain gauge-able.
+
+## Fields
+
+$(TYPEDFIELDS)
+"""
+@kwdef struct RandPerturbedExpand{S} <: Algorithm
+    "algorithm used for the singular value decomposition"
+    alg_svd::S = Defaults.alg_svd()
+
+    "algorithm used for [truncation](@extref MatrixAlgebraKit.TruncationStrategy] the expanded space"
+    trscheme::TruncationStrategy
+
+    "amount of noise that is added to the current state"
+    noisefactor::Float64 = eps()^(3 / 4)
+
+    "algorithm used for gauging the state"
+    alg_gauge = Defaults.alg_gauge(; dynamic_tols = false)
+end
+
+function changebonds!(ψ, alg::RandPerturbedExpand)
+    for i in 1:length(ψ)
+        # obtain space by sampling the support of left nullspace
+        # add (orthogonal) directions as isometries in that direction
+        VL = left_null(ψ.AL[i])
+        V = sample_space(right_virtualspace(VL), alg.trscheme)
+        XL = randisometry(scalartype(VL), right_virtualspace(VL) ← V)
+        ψ.AL[i] = catdomain(ψ.AL[i], VL * XL)
+
+        # make sure the next site fits, by "absorbing" into a larger tensor
+        # with some random noise to ensure state is still gauge-able
+        AL = ψ.AL[i + 1]
+        AL′ = similar(AL, right_virtualspace(ψ.AL[i]) ⊗ physicalspace(AL) ← right_virtualspace(AL))
+        scale!(randn!(AL′), alg.noisefactor)
+        ψ.AL[i + 1] = TensorKit.absorb!(AL′, AL)
+    end
+
+    # properly regauge the state:
+    makefullrank!(ψ.AL)
+    ψ.AR .= similar.(ψ.AL)
+    # ψ.AC .= similar.(ψ.AL)
+
+    # initial guess for gauge is embedded original C
+    C₀ = similar(ψ.C[0], right_virtualspace(ψ.AL[end]) ← left_virtualspace(ψ.AL[1]))
+    absorb!(id!(C₀), ψ.C[0])
+
+    gaugefix!(ψ, ψ.AL, C₀; order = :R, alg.alg_gauge.maxiter, alg.alg_gauge.tol)
+
+    for i in reverse(1:length(ψ))
+        # obtain space by sampling the support of left nullspace
+        # add (orthogonal) directions as isometries in that direction
+        AR_tail = _transpose_tail(ψ.AR[i]; copy = true)
+        VR = right_null(AR_tail)
+        V = sample_space(space(VR, 1), alg.trscheme)
+        XR = randisometry(scalartype(VR), space(VR, 1) ← V)
+        ψ.AR[i] = _transpose_front(catcodomain(AR_tail, XR' * VR))
+
+        # make sure the next site fits, by "absorbing" into a larger tensor
+        # with some random noise to ensure state is still gauge-able
+        AR = ψ.AR[i - 1]
+        AR′ = similar(AR, left_virtualspace(AR) ⊗ physicalspace(AR) ← left_virtualspace(ψ.AR[i]))
+        scale!(randn!(AR′), alg.noisefactor)
+        ψ.AR[i - 1] = TensorKit.absorb!(AR′, AR)
+    end
+
+    # properly regauge the state:
+    makefullrank!(ψ.AR)
+    ψ.AL .= similar.(ψ.AR)
+    ψ.AC .= similar.(ψ.AR)
+
+    # initial guess for gauge is embedded original C
+    C₀ = similar(ψ.C[0], right_virtualspace(ψ.AR[end]) ← left_virtualspace(ψ.AR[1]))
+    absorb!(id!(C₀), ψ.C[0])
+
+    gaugefix!(ψ, ψ.AR, C₀; order = :LR, alg.alg_gauge.maxiter, alg.alg_gauge.tol)
+    mul!.(ψ.AC, ψ.AL, ψ.C)
+
+    return ψ
+end
+
+changebonds(ψ, alg::RandPerturbedExpand) = changebonds!(copy(ψ), alg)

src/algorithms/groundstate/dmrg.jl

@@ -14,25 +14,30 @@ struct DMRG{A, F} <: Algorithm
     "maximal amount of iterations"
     maxiter::Int
 
+    miniter::Int
+
     "setting for how much information is displayed"
     verbosity::Int
 
     "algorithm used for the eigenvalue solvers"
     alg_eigsolve::A
 
+    expscheme::Algorithm
+
     "callback function applied after each iteration, of signature `finalize(iter, ψ, H, envs) -> ψ, envs`"
     finalize::F
 end
 function DMRG(;
         tol = Defaults.tol, maxiter = Defaults.maxiter, alg_eigsolve = (;),
-        verbosity = Defaults.verbosity, finalize = Defaults._finalize
+        verbosity = Defaults.verbosity, finalize = Defaults._finalize,
+        miniter = 0, expscheme = NoExpand()
     )
     alg_eigsolve′ = alg_eigsolve isa NamedTuple ? Defaults.alg_eigsolve(; alg_eigsolve...) :
         alg_eigsolve
-    return DMRG(tol, maxiter, verbosity, alg_eigsolve′, finalize)
+    return DMRG(tol, maxiter, miniter, verbosity, alg_eigsolve′, expscheme, finalize)
 end
 
-function find_groundstate!(ψ::AbstractFiniteMPS, H, alg::DMRG, envs = environments(ψ, H))
+function find_groundstate!(::FiniteChainStyle, ψ, H, alg::DMRG, envs = environments(ψ, H))
     ϵs = map(pos -> calc_galerkin(pos, ψ, H, ψ, envs), 1:length(ψ))
     ϵ = maximum(ϵs)
     log = IterLog("DMRG")
@@ -43,17 +48,29 @@ function find_groundstate!(ψ::AbstractFiniteMPS, H, alg::DMRG, envs = environme
             alg_eigsolve = updatetol(alg.alg_eigsolve, iter, ϵ)
 
             zerovector!(ϵs)
+            dir = 1
             for pos in [1:(length(ψ) - 1); length(ψ):-1:2]
+                pos == length(ψ) && (dir = -1)
                 h = AC_hamiltonian(pos, ψ, H, ψ, envs)
                 _, vec = fixedpoint(h, ψ.AC[pos], :SR, alg_eigsolve)
                 ϵs[pos] = max(ϵs[pos], calc_galerkin(pos, ψ, H, ψ, envs))
-                ψ.AC[pos] = vec
+                if alg.expscheme isa NoExpand 
+                    ψ.AC[pos] = vec 
+                elseif dir == 1
+                    AL, C = left_orth!(vec; positive = true)
+                    AL, C, ψ.AC[pos + 1] = changebonds_left(AL, C, ψ.AC[pos + 1], alg.expscheme)
+                    ψ.AC[pos] = (AL, C)
+                elseif dir == -1 
+                    C, temp = right_orth!(_transpose_tail(ψ.AC[pos]; copy = true); positive = true)
+                    C, ψ.AC[pos - 1], temp = changebonds_right(C, ψ.AC[pos - 1], temp, alg.expscheme)
+                    ψ.AC[pos] = (C, _transpose_front(temp))
+                end
             end
             ϵ = maximum(ϵs)
 
             ψ, envs = alg.finalize(iter, ψ, H, envs)::Tuple{typeof(ψ), typeof(envs)}
 
-            if ϵ <= alg.tol
+            if ϵ <= alg.tol && iter > alg.miniter
                 @infov 2 logfinish!(log, iter, ϵ, expectation_value(ψ, H, envs))
                 break
             end
@@ -83,6 +100,8 @@ struct DMRG2{A, S, F} <: Algorithm
     "maximal amount of iterations"
     maxiter::Int
 
+    miniter::Int
+
     "setting for how much information is displayed"
     verbosity::Int
 
@@ -95,21 +114,23 @@ struct DMRG2{A, S, F} <: Algorithm
     "algorithm used for [truncation](@extref MatrixAlgebraKit.TruncationStrategy) of the two-site update"
     trscheme::TruncationStrategy
 
+    expscheme::Algorithm
+
     "callback function applied after each iteration, of signature `finalize(iter, ψ, H, envs) -> ψ, envs`"
     finalize::F
 end
 # TODO: find better default truncation
 function DMRG2(;
         tol = Defaults.tol, maxiter = Defaults.maxiter, verbosity = Defaults.verbosity,
-        alg_eigsolve = (;), alg_svd = Defaults.alg_svd(), trscheme,
-        finalize = Defaults._finalize
+        miniter = 0, alg_eigsolve = (;), alg_svd = Defaults.alg_svd(), trscheme,
+        expscheme = NoExpand(), finalize = Defaults._finalize
     )
     alg_eigsolve′ = alg_eigsolve isa NamedTuple ? Defaults.alg_eigsolve(; alg_eigsolve...) :
         alg_eigsolve
-    return DMRG2(tol, maxiter, verbosity, alg_eigsolve′, alg_svd, trscheme, finalize)
+    return DMRG2(tol, maxiter, miniter, verbosity, alg_eigsolve′, alg_svd, trscheme, expscheme, finalize)
 end
 
-function find_groundstate!(ψ::AbstractFiniteMPS, H, alg::DMRG2, envs = environments(ψ, H))
+function find_groundstate!(::FiniteChainStyle, ψ, H, alg::DMRG2, envs = environments(ψ, H))
     ϵs = map(pos -> calc_galerkin(pos, ψ, H, ψ, envs), 1:length(ψ))
     ϵ = maximum(ϵs)
     log = IterLog("DMRG2")
@@ -126,6 +147,7 @@ function find_groundstate!(ψ::AbstractFiniteMPS, H, alg::DMRG2, envs = environm
                 _, newA2center = fixedpoint(Hac2, ac2, :SR, alg_eigsolve)
 
                 al, c, ar = svd_trunc!(newA2center; trunc = alg.trscheme, alg = alg.alg_svd)
+                al, c = changebonds_left(al, c, alg.expscheme)
                 normalize!(c)
                 v = @plansor ac2[1 2; 3 4] * conj(al[1 2; 5]) * conj(c[5; 6]) * conj(ar[6; 3 4])
                 ϵs[pos] = max(ϵs[pos], abs(1 - abs(v)))
@@ -141,6 +163,7 @@ function find_groundstate!(ψ::AbstractFiniteMPS, H, alg::DMRG2, envs = environm
                 _, newA2center = fixedpoint(Hac2, ac2, :SR, alg_eigsolve)
 
                 al, c, ar = svd_trunc!(newA2center; trunc = alg.trscheme, alg = alg.alg_svd)
+                c, ar = changebonds_right(c, ar, alg.expscheme)
                 normalize!(c)
                 v = @plansor ac2[1 2; 3 4] * conj(al[1 2; 5]) * conj(c[5; 6]) * conj(ar[6; 3 4])
                 ϵs[pos] = max(ϵs[pos], abs(1 - abs(v)))
@@ -152,7 +175,7 @@ function find_groundstate!(ψ::AbstractFiniteMPS, H, alg::DMRG2, envs = environm
             ϵ = maximum(ϵs)
             ψ, envs = alg.finalize(iter, ψ, H, envs)::Tuple{typeof(ψ), typeof(envs)}
 
-            if ϵ <= alg.tol
+            if ϵ <= alg.tol && iter > alg.miniter
                 @infov 2 logfinish!(log, iter, ϵ, expectation_value(ψ, H, envs))
                 break
             end
@@ -165,7 +188,3 @@ function find_groundstate!(ψ::AbstractFiniteMPS, H, alg::DMRG2, envs = environm
     end
     return ψ, envs, ϵ
 end
-
-function find_groundstate(ψ, H, alg::Union{DMRG, DMRG2}, envs...; kwargs...)
-    return find_groundstate!(copy(ψ), H, alg, envs...; kwargs...)
-end