[WIP] Alternative random expansion scheme

CHANGELOG.md

@@ -0,0 +1,33 @@
+# Changelog
+
+All notable changes to this project will be documented in this file.
+
+The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
+and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
+
+## [Unreleased]
+
+### Added
+
+- `LocalPreferences.toml` file to ensure `TensorOperations` properly precompiles on testing
+  infrastructure
+
+### Fixed
+
+- Dynamic tolerances yielded `NaN` during the initialization stage due to `1 / sqrt(iter)`
+  where `iter = 0`.
+- `InfiniteMPOHamiltonian` environments with low bond dimension and high Krylov dimension now are properly
+  clamped.
+
+### Changed
+
+- The `changebonds(state, ::RandExpand)` algorithm now no longer has to perform a
+  truncated SVD to obtain the desired spaces, and instead sample the space directly
+  and then generates a random isometry. This should be slightly more performant, but
+  otherwise equivalent.
+
+### Deprecated
+
+### Removed
+
+[unreleased]: https://github.com/quantumkithub/pepskit.jl/compare/v0.13.8...HEAD

LocalPreferences.toml

@@ -0,0 +1,2 @@
+[TensorOperations]
+precompile_workload = true

docs/src/man/algorithms.md

@@ -292,7 +292,7 @@ disadvantages:
   and if the bond dimension is grown slow enough, this still obtains a very good expansion
   scheme. Again, The state will remain unchanged and a one-site scheme will now be able to 
   push the optimization further. The subspace used for expansion can be truncated through
-  `trscheme`, which dictates how many singular values will be added.
+  `trscheme`, which dictates how many orthogonal vectors will be added.
 
 * [`VUMPSSvdCut`](@ref): This algorithm is based on the [`VUMPS`](@ref) algorithm, and
   consists of performing a two-site update, and then truncating the state back down. Because

src/MPSKit.jl

@@ -36,7 +36,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

src/algorithms/changebonds/randexpand.jl

@@ -2,8 +2,9 @@
 $(TYPEDEF)
 
 An algorithm that expands the bond dimension by adding random unitary vectors that are
-orthogonal to the existing state. This is achieved by performing a truncated SVD on a random
-two-site MPS tensor, which is made orthogonal to the existing state.
+orthogonal to the existing state. This means that additional directions are added to
+`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.
 
 ## Fields
 
@@ -17,32 +18,35 @@ $(TYPEDFIELDS)
     trscheme::TruncationStrategy
 end
 
-function changebonds(ψ::InfiniteMPS, alg::RandExpand)
+function changebonds!(ψ::InfiniteMPS, alg::RandExpand)
     T = eltype(ψ.AL)
     AL′ = similar(ψ.AL)
     AR′ = similar(ψ.AR, tensormaptype(spacetype(T), 1, numind(T) - 1, storagetype(T)))
-    for i in 1:length(ψ)
-        # determine optimal expansion spaces around bond i
-        AC2 = randomize!(_transpose_front(ψ.AC[i]) * _transpose_tail(ψ.AR[i + 1]))
 
-        # Use the nullspaces and SVD decomposition to determine the optimal expansion space
+    for i in 1:length(ψ)
+        # obtain spaces by sampling the support of both the left and right nullspace
         VL = left_null(ψ.AL[i])
         VR = right_null!(_transpose_tail(ψ.AR[i + 1]; copy = true))
-        intermediate = normalize!(adjoint(VL) * AC2 * adjoint(VR))
-        U, _, Vᴴ = svd_trunc!(intermediate; trunc = alg.trscheme, alg = alg.alg_svd)
+        V = sample_space(infimum(right_virtualspace(VL), space(VR, 1)), alg.trscheme)
 
-        AL′[i] = VL * U
-        AR′[i + 1] = Vᴴ * VR
+        # obtain (orthogonal) directions as isometries in that direction
+        XL = randisometry(scalartype(VL), right_virtualspace(VL) ← V)
+        AL′[i] = VL * XL
+        XR = randisometry(storagetype(VL), space(VR, 1) ← V)
+        AR′[i + 1] = XR * VR
     end
 
-    return _expand(ψ, AL′, AR′)
+    return _expand!(ψ, AL′, AR′)
 end
 
-function changebonds(ψ::MultilineMPS, alg::RandExpand)
-    return Multiline(map(x -> changebonds(x, alg), ψ.data))
+function changebonds!(ψ::MultilineMPS, alg::RandExpand)
+    foreach(Base.Fix2(changebonds!, alg), ψ.data)
+    return ψ
 end
 
-changebonds(ψ::AbstractFiniteMPS, alg::RandExpand) = changebonds!(copy(ψ), alg)
+changebonds(ψ::AbstractMPS, alg::RandExpand) = changebonds!(copy(ψ), alg)
+changebonds(ψ::MultilineMPS, alg::RandExpand) = changebonds!(copy(ψ), alg)
+
 function changebonds!(ψ::AbstractFiniteMPS, alg::RandExpand)
     for i in 1:(length(ψ) - 1)
         AC2 = randomize!(_transpose_front(ψ.AC[i]) * _transpose_tail(ψ.AR[i + 1]))
@@ -67,3 +71,97 @@ function changebonds!(ψ::AbstractFiniteMPS, alg::RandExpand)
 
     return normalize!(ψ)
 end
+
+function sample_space(V, strategy)
+    S = TensorKit.SectorDict(c => sort!(Random.rand(dim(V, c)); rev = true) for c in sectors(V))
+    ind = MatrixAlgebraKit.findtruncated_svd(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!(ψ::InfiniteMPS, 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(ψ::InfiniteMPS, alg::RandPerturbedExpand) = changebonds!(copy(ψ), alg)

src/environments/abstract_envs.jl

@@ -81,7 +81,8 @@ function environment_alg(
         tol = Defaults.tol, maxiter = Defaults.maxiter, krylovdim = Defaults.krylovdim,
         verbosity = Defaults.VERBOSE_NONE
     )
-    return GMRES(; tol, maxiter, krylovdim, verbosity)
+    max_krylovdim = dim(left_virtualspace(above, 1)) * dim(left_virtualspace(below, 1))
+    return GMRES(; tol, maxiter, krylovdim = min(max_krylovdim, krylovdim), verbosity)
 end
 function environment_alg(
         ::Union{InfiniteQP, MultilineQP}, ::Union{InfiniteMPO, MultilineMPO},

src/states/abstractmps.jl

@@ -93,10 +93,11 @@ Determine whether the given tensor is full rank, i.e. whether both the map from
 virtual space and the physical space to the right virtual space, and the map from the right
 virtual space and the physical space to the left virtual space are injective.
 """
-function isfullrank(A::GenericMPSTensor; side = :both)
-    Vₗ = _firstspace(A)
-    Vᵣ = _lastspace(A)
-    P = ⊗(space.(Ref(A), 2:(numind(A) - 1))...)
+isfullrank(A::GenericMPSTensor; kwargs...) = isfullrank(space(A); kwargs...)
+function isfullrank(V::TensorKit.TensorMapSpace; side = :both)
+    Vₗ = V[1]
+    Vᵣ = V[numind(V)]
+    P = ⊗(getindex.(Ref(V), 2:(numind(V) - 1))...)
     return if side === :both
         Vₗ ⊗ P ≿ Vᵣ' && Vₗ' ≾ P ⊗ Vᵣ
     elseif side === :right
@@ -129,6 +130,29 @@ function makefullrank!(A::PeriodicVector{<:GenericMPSTensor}; alg = Defaults.alg
     return A
 end
 
+function makefullrank!(virtualspaces::PeriodicVector{S}, physicalspaces::PeriodicVector{S}) where {S <: ElementarySpace}
+    haschanged = true
+    while haschanged
+        haschanged = false
+        for i in 1:length(virtualspaces)
+            Vmax = fuse(virtualspaces[i - 1], physicalspaces[i - 1])
+            if !(virtualspaces[i] ≾ Vmax)
+                virtualspaces[i] = infimum(virtualspaces[i], Vmax)
+                haschanged = true
+            end
+        end
+        for i in reverse(1:length(virtualspaces))
+            Vmax = fuse(dual(physicalspaces[i - 1]), virtualspaces[i])
+            if !(virtualspaces[i - 1] ≾ Vmax)
+                virtualspaces[i - 1] = infimum(virtualspaces[i - 1], Vmax)
+                haschanged = true
+            end
+        end
+    end
+
+    return virtualspaces
+end
+
 # Tensor accessors
 # ----------------
 @doc """

src/utility/dynamictols.jl

@@ -62,6 +62,7 @@ where the new tolerance is given by
     new_tol = clamp(ϵ * alg.tol_factor / sqrt(iter), alg.tol_min, alg.tol_max)
 """
 function updatetol(alg::DynamicTol, iter::Integer, ϵ::Real)
+    iter = max(iter, one(iter))
     new_tol = clamp(ϵ * alg.tol_factor / sqrt(iter), alg.tol_min, alg.tol_max)
     return _updatetol(alg.alg, new_tol)
 end