diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 000000000..26cc89842 --- /dev/null +++ b/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 diff --git a/LocalPreferences.toml b/LocalPreferences.toml new file mode 100644 index 000000000..cdb275e21 --- /dev/null +++ b/LocalPreferences.toml @@ -0,0 +1,2 @@ +[TensorOperations] +precompile_workload = true diff --git a/docs/src/man/algorithms.md b/docs/src/man/algorithms.md index 378bb1af0..261bf133d 100644 --- a/docs/src/man/algorithms.md +++ b/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 diff --git a/src/MPSKit.jl b/src/MPSKit.jl index a6c586144..342094cbd 100644 --- a/src/MPSKit.jl +++ b/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 diff --git a/src/algorithms/changebonds/randexpand.jl b/src/algorithms/changebonds/randexpand.jl index 82cb684cf..94c39eb7b 100644 --- a/src/algorithms/changebonds/randexpand.jl +++ b/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) diff --git a/src/environments/abstract_envs.jl b/src/environments/abstract_envs.jl index cbdcafee1..1ec89a5ee 100644 --- a/src/environments/abstract_envs.jl +++ b/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}, diff --git a/src/states/abstractmps.jl b/src/states/abstractmps.jl index dd756b2f4..6dd8a1341 100644 --- a/src/states/abstractmps.jl +++ b/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 """ diff --git a/src/utility/dynamictols.jl b/src/utility/dynamictols.jl index dc2a15781..75bea4d2f 100644 --- a/src/utility/dynamictols.jl +++ b/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