perf: DMRG(2) with smarter adaptive tolerances

src/MPSKit.jl

@@ -52,7 +52,7 @@ export r_LL, l_LL, r_RR, l_RR, r_RL, r_LR, l_RL, l_LR # TODO: rename
 
 # unexported
 using Compat: @compat
-@compat public DynamicTols
+@compat public DynamicTols, DynamicTol, AdaptiveKrylov
 @compat public VERBOSE_NONE, VERBOSE_WARN, VERBOSE_CONV, VERBOSE_ITER, VERBOSE_ALL
 @compat public IterLog, loginit!, logiter!, logfinish!, logcancel!
 

src/algorithms/approximate/fvomps.jl

@@ -8,7 +8,7 @@ function approximate!(ψ::AbstractFiniteMPS, Oϕ, alg::DMRG2, envs = environment
             ϵ = 0.0
             for pos in [1:(length(ψ) - 1); (length(ψ) - 2):-1:1]
                 AC2′ = AC2_projection(pos, ψ, Oϕ, envs)
-                al, c, ar = svd_trunc!(AC2′; trunc = alg.trscheme, alg = alg.alg_svd)
+                al, c, ar, = svd_trunc!(AC2′, alg.alg_gauge)
 
                 AC2 = ψ.AC[pos] * _transpose_tail(ψ.AR[pos + 1])
                 ϵ = max(ϵ, norm(al * c * ar - AC2) / norm(AC2))

src/algorithms/approximate/idmrg.jl

@@ -52,7 +52,7 @@ function approximate!(
     end
 
     # TODO: immediately compute in-place
-    alg_gauge = updatetol(alg.alg_gauge, iter, ϵ)
+    alg_gauge = adapt_solver(alg.alg_gauge; iter, g_global = ϵ)
     ψ′ = MultilineMPS(map(x -> x, ψ.AR); alg_gauge.tol, alg_gauge.maxiter)
     copy!(ψ, ψ′) # ensure output destination is unchanged
 
@@ -179,7 +179,7 @@ function approximate!(
     end
 
     # TODO: immediately compute in-place
-    alg_gauge = updatetol(alg.alg_gauge, iter, ϵ)
+    alg_gauge = adapt_solver(alg.alg_gauge; iter, g_global = ϵ)
     ψ′ = MultilineMPS(map(identity, ψ.AR); alg_gauge.tol, alg_gauge.maxiter)
     copy!(ψ, ψ′) # ensure output destination is unchanged
 

src/algorithms/approximate/vomps.jl

@@ -17,7 +17,7 @@ function approximate(
     log = IterLog("VOMPS")
     iter = 0
     ϵ = calc_galerkin(mps, toapprox..., envs)
-    alg_environments = updatetol(alg.alg_environments, iter, ϵ)
+    alg_environments = adapt_solver(alg.alg_environments; iter, g_global = ϵ)
     recalculate!(envs, mps, toapprox..., alg_environments)
 
     state = VOMPSState(mps, toapprox, envs, iter, ϵ)
@@ -65,7 +65,7 @@ end
 function localupdate_step!(
         it::IterativeSolver{<:VOMPS}, state::VOMPSState{<:Any, <:Tuple}, ::SerialScheduler
     )
-    alg_gauge = updatetol(it.alg_gauge, state.iter, state.ϵ)
+    alg_gauge = adapt_solver(it.alg_gauge; iter = state.iter, g_global = state.ϵ)
     alg_orth = alg_gauge.alg_orth
 
     ACs = similar(state.mps.AC)
@@ -89,7 +89,7 @@ end
 function localupdate_step!(
         it::IterativeSolver{<:VOMPS}, state::VOMPSState{<:Any, <:Tuple}, scheduler
     )
-    alg_gauge = updatetol(it.alg_gauge, state.iter, state.ϵ)
+    alg_gauge = adapt_solver(it.alg_gauge; iter = state.iter, g_global = state.ϵ)
     alg_orth = alg_gauge.alg_orth
 
     ACs = similar(state.mps.AC)
@@ -119,6 +119,6 @@ function localupdate_step!(
 end
 
 function envs_step!(it::IterativeSolver{<:VOMPS}, state::VOMPSState{<:Any, <:Tuple}, mps)
-    alg_environments = updatetol(it.alg_environments, state.iter, state.ϵ)
+    alg_environments = adapt_solver(it.alg_environments; iter = state.iter, g_global = state.ϵ)
     return recalculate!(state.envs, mps, state.operator..., alg_environments)
 end

src/algorithms/changebonds/optimalexpand.jl

@@ -114,7 +114,7 @@ function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::Op
     # 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))
+    nal, nc, _ = left_gauge(absorb!(zerovector!(similar(left, nal_space)), left))
     nar = _transpose_front(catcodomain(_transpose_tail(right), ar_re))
 
     normalize && normalize!(nc)

src/algorithms/changebonds/randexpand.jl

@@ -98,7 +98,7 @@ function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::Ra
     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))
+    nal, nc, _ = left_gauge(absorb!(zerovector!(similar(left, nal_space)), left))
     nar = _transpose_front(catcodomain(_transpose_tail(right), ar_re))
 
     normalize && normalize!(nc)

src/algorithms/changebonds/sketchedexpand.jl

@@ -53,7 +53,7 @@ end
 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
+    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)
@@ -81,7 +81,7 @@ function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::Sk
     # 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))
+    nal, nc, _ = left_gauge(absorb!(zerovector!(similar(left, nal_space)), left))
     nar = _transpose_front(catcodomain(_transpose_tail(right), ar_re))
 
     normalize && normalize!(nc)

src/algorithms/excitation/dmrgexcitation.jl

@@ -18,12 +18,21 @@ $(TYPEDFIELDS)
     weight::Float64 = 10.0
 end
 
+# Initialize excited state by perturbing the current eigenvector to avoid local minima
+function _perturbed_state(ψ::FiniteMPS; atol = 1.0e-2)
+    return FiniteMPS(
+        map(1:length(ψ)) do i
+            A = i == length(ψ) ? ψ.AC[i] : ψ.AL[i]
+            noise = randomize!(similar(A))
+            return add!(noise, A, 1, atol / norm(noise))
+        end
+    )
+end
+
 function excitations(
         H::FiniteMPOHamiltonian, alg::FiniteExcited,
         states::Tuple{T, Vararg{T}};
-        init = FiniteMPS(
-            [ copy(first(states).AC[i]) for i in 1:length(first(states)) ]
-        ), num = 1
+        init = _perturbed_state(first(states)), num = 1
     ) where {T <: FiniteMPS}
     num == 0 && return (scalartype(T)[], T[])
 

src/algorithms/fixedpoint.jl

@@ -1,32 +1,23 @@
 # wrapper around KrylovKit.jl's eigsolve function
 
+function fixedpoint(A, x₀, which::Symbol; kwargs...)
+    alg = KrylovKit.eigselector(A, scalartype(x₀); kwargs...)
+    return fixedpoint(A, x₀, which, alg)
+end
+
 """
-    fixedpoint(A, x₀, which::Symbol; kwargs...) -> val, vec
-    fixedpoint(A, x₀, which::Symbol, alg) -> val, vec
+    fixedpoint(A, x₀, which::Symbol; kwargs...) -> val, vec, info
+    fixedpoint(A, x₀, which::Symbol, alg) -> val, vec, info
 
-Compute the fixed point of a linear operator `A` using the specified eigensolver `alg`. The
-fixedpoint is assumed to be unique.
+Compute the fixed point of a given linear operator `A` with initial guess `x₀`.
+The dominant eigenvector is assumed to be unique.
 """
 function fixedpoint(A, x₀, which::Symbol, alg::Lanczos)
     vals, vecs, info = eigsolve(A, x₀, 1, which, alg)
-
-    info.converged == 0 &&
-        @warnv 1 "fixed point not converged after $(info.numiter) iterations: normres = $(info.normres[1])"
-
-    return vals[1], vecs[1]
+    return vals[1], vecs[1], info
 end
-
 function fixedpoint(A, x₀, which::Symbol, alg::Arnoldi)
     TT, vecs, vals, info = schursolve(A, x₀, 1, which, alg)
-
-    info.converged == 0 &&
-        @warnv 1 "fixed point not converged after $(info.numiter) iterations: normres = $(info.normres[1])"
     size(TT, 2) > 1 && TT[2, 1] != 0 && @warnv 1 "non-unique fixed point detected"
-
-    return vals[1], vecs[1]
-end
-
-function fixedpoint(A, x₀, which::Symbol; kwargs...)
-    alg = KrylovKit.eigselector(A, scalartype(x₀); kwargs...)
-    return fixedpoint(A, x₀, which, alg)
+    return vals[1], vecs[1], info
 end