Various small fixes and improvements to RandExpand (#335)

* rework random expansion

* update docstring

* add and update changelog

* fix spaces

* fix overexpanding

* add changebonds for MultilineMPS

* ensure spaces are consistent

* add periodicvector

* clamp iteration count to avoid NaN

* rework implementation again

* clamp krylovdim

* add utility fullrank for spaces

* remove stray show

* fix off-by-one error

* discard attempts and go back to original implementation

* update changelog

* update docs

* fix in-place for multiline

* more small fixes and changelog updates

* improve docs

* small fix

* another small fix

Author: lkdvos

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/algorithms/changebonds/randexpand.jl

@@ -2,8 +2,12 @@
 $(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.
+
+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
 
@@ -17,32 +21,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(scalartype(VR), 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 +74,15 @@ function changebonds!(ψ::AbstractFiniteMPS, alg::RandExpand)
 
     return normalize!(ψ)
 end
+
+"""
+    sample_space(V, strategy)
+
+Sample basis states within a given `V::VectorSpace` by creating weights for each state that
+are distributed uniformly, and then truncating according to the given `strategy`.
+"""
+function sample_space(V, strategy)
+    S = TensorKit.SectorDict(c => Random.rand(dim(V, c)) for c in sectors(V))
+    ind = MatrixAlgebraKit.findtruncated(S, strategy)
+    return TensorKit.Factorizations.truncate_space(V, ind)
+end

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