From 683acb1d2a3aa5181d7c4c73bccc6d57285a12bb Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Tue, 21 Oct 2025 11:01:32 -0400 Subject: [PATCH 01/21] rework random expansion --- src/algorithms/changebonds/randexpand.jl | 37 ++++++++++++------------ 1 file changed, 18 insertions(+), 19 deletions(-) diff --git a/src/algorithms/changebonds/randexpand.jl b/src/algorithms/changebonds/randexpand.jl index 82cb684cf..79442e9ee 100644 --- a/src/algorithms/changebonds/randexpand.jl +++ b/src/algorithms/changebonds/randexpand.jl @@ -17,32 +17,31 @@ $(TYPEDFIELDS) trscheme::TruncationStrategy end -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 - 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) +function changebonds!(ψ::InfiniteMPS, alg::RandExpand) + AL′ = map(ψ.AL) do A + # find random orthogonal vectors + A_perp = randn!(similar(A)) + add!(A_perp, A * (A' * A_perp), -1) + A′, _, _ = svd_trunc!(A_perp; alg = alg.alg_svd, trunc = alg.trscheme) + return A′ + end - AL′[i] = VL * U - AR′[i + 1] = Vᴴ * VR + AR′ = map(ψ.AR) do A + At = _transpose_tail(A) + A_perp = randn!(similar(At)) + add!(A_perp, (A_perp * At') * At, -1) + _, _, A′ = svd_trunc!(A_perp; alg = alg.alg_svd, trunc = alg.trscheme) + return A′ 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) + return Multiline(map(x -> changebonds!(x, alg), ψ.data)) end -changebonds(ψ::AbstractFiniteMPS, alg::RandExpand) = changebonds!(copy(ψ), alg) +changebonds(ψ::AbstractMPS, 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])) From dabfae7be513d6c20e2c520de019fcaa90651bf5 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Tue, 21 Oct 2025 11:02:40 -0400 Subject: [PATCH 02/21] update docstring --- docs/src/man/algorithms.md | 4 +++- src/algorithms/changebonds/randexpand.jl | 3 +-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/docs/src/man/algorithms.md b/docs/src/man/algorithms.md index 378bb1af0..384295d10 100644 --- a/docs/src/man/algorithms.md +++ b/docs/src/man/algorithms.md @@ -292,7 +292,9 @@ 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 singular values will be added. Additionally, this + method does not project onto the local two-site basis, which might overcome symmetry-based + obstructions that can be encountered in the other methods. * [`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/algorithms/changebonds/randexpand.jl b/src/algorithms/changebonds/randexpand.jl index 79442e9ee..084c555c2 100644 --- a/src/algorithms/changebonds/randexpand.jl +++ b/src/algorithms/changebonds/randexpand.jl @@ -2,8 +2,7 @@ $(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. ## Fields From 96e968f71c4785af4cd39843f4825407df269171 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Tue, 21 Oct 2025 11:11:37 -0400 Subject: [PATCH 03/21] add and update changelog --- CHANGELOG.md | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) create mode 100644 CHANGELOG.md diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 000000000..810c8d927 --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,25 @@ +# 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 + +### Fixed + +### Changed + +- The `changebonds(state, ::RandExpand)` algorithm now no longer projects onto the nullspace + on both sides. This ensures that the expanded symmetry sectors can be selected beyond what + is allowed by two-site updates, which can be relevant for certain systems that have + symmetry-related obstructions. + +### Deprecated + +### Removed + +[unreleased]: https://github.com/quantumkithub/pepskit.jl/compare/v0.13.8...HEAD From 5326ff806668177061a5f89e0fb74a60d0367ab9 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Tue, 21 Oct 2025 13:49:37 -0400 Subject: [PATCH 04/21] fix spaces --- src/algorithms/changebonds/randexpand.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/algorithms/changebonds/randexpand.jl b/src/algorithms/changebonds/randexpand.jl index 084c555c2..0a40b27ac 100644 --- a/src/algorithms/changebonds/randexpand.jl +++ b/src/algorithms/changebonds/randexpand.jl @@ -19,7 +19,7 @@ end function changebonds!(ψ::InfiniteMPS, alg::RandExpand) AL′ = map(ψ.AL) do A # find random orthogonal vectors - A_perp = randn!(similar(A)) + A_perp = randn!(similar(A, codomain(A) ← fuse(codomain(A)))) add!(A_perp, A * (A' * A_perp), -1) A′, _, _ = svd_trunc!(A_perp; alg = alg.alg_svd, trunc = alg.trscheme) return A′ @@ -27,7 +27,7 @@ function changebonds!(ψ::InfiniteMPS, alg::RandExpand) AR′ = map(ψ.AR) do A At = _transpose_tail(A) - A_perp = randn!(similar(At)) + A_perp = randn!(similar(At, fuse(domain(At)) ← domain(At))) add!(A_perp, (A_perp * At') * At, -1) _, _, A′ = svd_trunc!(A_perp; alg = alg.alg_svd, trunc = alg.trscheme) return A′ From 5b1e29c8fef5b9e0ff64afe2bc6914a522b874ab Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Tue, 21 Oct 2025 13:52:43 -0400 Subject: [PATCH 05/21] fix overexpanding --- src/algorithms/changebonds/randexpand.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/algorithms/changebonds/randexpand.jl b/src/algorithms/changebonds/randexpand.jl index 0a40b27ac..ed854fdd3 100644 --- a/src/algorithms/changebonds/randexpand.jl +++ b/src/algorithms/changebonds/randexpand.jl @@ -19,16 +19,16 @@ end function changebonds!(ψ::InfiniteMPS, alg::RandExpand) AL′ = map(ψ.AL) do A # find random orthogonal vectors - A_perp = randn!(similar(A, codomain(A) ← fuse(codomain(A)))) - add!(A_perp, A * (A' * A_perp), -1) + A_perp = randn!(similar(A, codomain(A) ← fuse(codomain(A)) ⊖ right_virtualspace(A))) + normalize!(add!(A_perp, A * (A' * A_perp), -1)) A′, _, _ = svd_trunc!(A_perp; alg = alg.alg_svd, trunc = alg.trscheme) return A′ end AR′ = map(ψ.AR) do A At = _transpose_tail(A) - A_perp = randn!(similar(At, fuse(domain(At)) ← domain(At))) - add!(A_perp, (A_perp * At') * At, -1) + A_perp = randn!(similar(At, fuse(domain(At)) ⊖ left_virtualspace(A) ← domain(At))) + normalize!(add!(A_perp, (A_perp * At') * At, -1)) _, _, A′ = svd_trunc!(A_perp; alg = alg.alg_svd, trunc = alg.trscheme) return A′ end From 2ad7ba57176dae70265eac9cd8e00bbd500175c7 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 22 Oct 2025 09:13:23 -0400 Subject: [PATCH 06/21] add changebonds for MultilineMPS --- src/algorithms/changebonds/randexpand.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/algorithms/changebonds/randexpand.jl b/src/algorithms/changebonds/randexpand.jl index ed854fdd3..2c0e6e0e4 100644 --- a/src/algorithms/changebonds/randexpand.jl +++ b/src/algorithms/changebonds/randexpand.jl @@ -41,6 +41,8 @@ function changebonds!(ψ::MultilineMPS, alg::RandExpand) end 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])) From 9d8faf3cb3bf86c8b83ade2a3eabdc55cf1dbc5f Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 22 Oct 2025 09:19:14 -0400 Subject: [PATCH 07/21] ensure spaces are consistent --- src/algorithms/changebonds/randexpand.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/algorithms/changebonds/randexpand.jl b/src/algorithms/changebonds/randexpand.jl index 2c0e6e0e4..2b5122ddf 100644 --- a/src/algorithms/changebonds/randexpand.jl +++ b/src/algorithms/changebonds/randexpand.jl @@ -25,11 +25,12 @@ function changebonds!(ψ::InfiniteMPS, alg::RandExpand) return A′ end - AR′ = map(ψ.AR) do A + AR′ = map(enumerate(ψ.AR)) do (i, A) At = _transpose_tail(A) A_perp = randn!(similar(At, fuse(domain(At)) ⊖ left_virtualspace(A) ← domain(At))) normalize!(add!(A_perp, (A_perp * At') * At, -1)) - _, _, A′ = svd_trunc!(A_perp; alg = alg.alg_svd, trunc = alg.trscheme) + trunc = truncspace(right_virtualspace(AL′[i - 1])) + _, _, A′ = svd_trunc!(A_perp; alg = alg.alg_svd, trunc) return A′ end From c85e30efb4817c65e61257405750c506aad3db80 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 22 Oct 2025 10:25:14 -0400 Subject: [PATCH 08/21] add periodicvector --- src/algorithms/changebonds/randexpand.jl | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/src/algorithms/changebonds/randexpand.jl b/src/algorithms/changebonds/randexpand.jl index 2b5122ddf..9761e07a9 100644 --- a/src/algorithms/changebonds/randexpand.jl +++ b/src/algorithms/changebonds/randexpand.jl @@ -25,14 +25,16 @@ function changebonds!(ψ::InfiniteMPS, alg::RandExpand) return A′ end - AR′ = map(enumerate(ψ.AR)) do (i, A) - At = _transpose_tail(A) - A_perp = randn!(similar(At, fuse(domain(At)) ⊖ left_virtualspace(A) ← domain(At))) - normalize!(add!(A_perp, (A_perp * At') * At, -1)) - trunc = truncspace(right_virtualspace(AL′[i - 1])) - _, _, A′ = svd_trunc!(A_perp; alg = alg.alg_svd, trunc) - return A′ - end + AR′ = PeriodicVector( + map(enumerate(ψ.AR)) do (i, A) + At = _transpose_tail(A) + A_perp = randn!(similar(At, fuse(domain(At)) ⊖ left_virtualspace(A) ← domain(At))) + normalize!(add!(A_perp, (A_perp * At') * At, -1)) + trunc = truncspace(right_virtualspace(AL′[i - 1])) + _, _, A′ = svd_trunc!(A_perp; alg = alg.alg_svd, trunc) + return A′ + end + ) return _expand!(ψ, AL′, AR′) end From 5c5185fbadf30e1326ae9173cda2b666244d316f Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 22 Oct 2025 13:10:19 -0400 Subject: [PATCH 09/21] clamp iteration count to avoid NaN --- src/utility/dynamictols.jl | 1 + 1 file changed, 1 insertion(+) 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 From cce516ddcb19243b69ccec6d59b120fcba84fdff Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 22 Oct 2025 15:06:45 -0400 Subject: [PATCH 10/21] rework implementation again --- src/algorithms/changebonds/randexpand.jl | 52 ++++++++++++++++-------- 1 file changed, 36 insertions(+), 16 deletions(-) diff --git a/src/algorithms/changebonds/randexpand.jl b/src/algorithms/changebonds/randexpand.jl index 9761e07a9..edf3b11d0 100644 --- a/src/algorithms/changebonds/randexpand.jl +++ b/src/algorithms/changebonds/randexpand.jl @@ -14,28 +14,42 @@ $(TYPEDFIELDS) "algorithm used for [truncation](@extref MatrixAlgebraKit.TruncationStrategy] the expanded space" trscheme::TruncationStrategy + + "expansion range that is considered for selecting the orthogonal subspace" + range::Int = 1 end function changebonds!(ψ::InfiniteMPS, alg::RandExpand) - AL′ = map(ψ.AL) do A - # find random orthogonal vectors - A_perp = randn!(similar(A, codomain(A) ← fuse(codomain(A)) ⊖ right_virtualspace(A))) - normalize!(add!(A_perp, A * (A' * A_perp), -1)) - A′, _, _ = svd_trunc!(A_perp; alg = alg.alg_svd, trunc = alg.trscheme) - return A′ - end - - AR′ = PeriodicVector( - map(enumerate(ψ.AR)) do (i, A) - At = _transpose_tail(A) - A_perp = randn!(similar(At, fuse(domain(At)) ⊖ left_virtualspace(A) ← domain(At))) - normalize!(add!(A_perp, (A_perp * At') * At, -1)) - trunc = truncspace(right_virtualspace(AL′[i - 1])) - _, _, A′ = svd_trunc!(A_perp; alg = alg.alg_svd, trunc) - return A′ + # obtain the spaces for the expanded directions by sampling + virtualspaces = PeriodicVector( + map(1:length(ψ)) do i + Vmax = mapreduce(fuse, reverse(0:alg.range)) do j + j == alg.range ? left_virtualspace(ψ, i - j) : physicalspace(ψ, i - j - 1) + end ⊖ left_virtualspace(ψ, i) + return sample_space(Vmax, alg.trscheme) end ) + # ensure the resulting tensors are full rank (space-wise) + makefullrank!(virtualspaces, physicalspace(ψ)) + + # add vectors orthogonal to the current state + AL′ = similar(ψ.AL) + T = eltype(AL′) + AR′ = similar(ψ.AR, tensormaptype(spacetype(T), 1, numind(T) - 1, storagetype(T))) + + for i in 1:length(ψ) + VL = left_null(ψ.AL[i]) + XL = similar(VL, right_virtualspace(VL) ← virtualspaces[i + 1]) + foreach(((c, b),) -> TensorKit.one!(b), blocks(XL)) + AL′[i] = VL * XL + + VR = right_null!(_transpose_tail(ψ.AR[i])) + XR = similar(VR, virtualspaces[i] ← space(VR, 1)) + foreach(((c, b),) -> TensorKit.one!(b), blocks(XR)) + AR′[i] = XR * VR + end + return _expand!(ψ, AL′, AR′) end @@ -70,3 +84,9 @@ function changebonds!(ψ::AbstractFiniteMPS, alg::RandExpand) return normalize!(ψ) end + +function sample_space(V, strategy) + S = TensorKit.SectorDict(c => Random.randexp(dim(V, c)) for c in sectors(V)) + ind = MatrixAlgebraKit.findtruncated(S, strategy) + return TensorKit.Factorizations.truncate_space(V, ind) +end From ee8277bc8ab010cd1bbedefa96f064acc8991e3e Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 22 Oct 2025 15:06:55 -0400 Subject: [PATCH 11/21] clamp krylovdim --- src/environments/abstract_envs.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/environments/abstract_envs.jl b/src/environments/abstract_envs.jl index cbdcafee1..54c20dee4 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 = @show(min(max_krylovdim, krylovdim)), verbosity) end function environment_alg( ::Union{InfiniteQP, MultilineQP}, ::Union{InfiniteMPO, MultilineMPO}, From d88293610cab9afb0f41336546e8ff69fe3f9fea Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 22 Oct 2025 15:07:07 -0400 Subject: [PATCH 12/21] add utility fullrank for spaces --- src/states/abstractmps.jl | 32 ++++++++++++++++++++++++++++---- 1 file changed, 28 insertions(+), 4 deletions(-) diff --git a/src/states/abstractmps.jl b/src/states/abstractmps.jl index dd756b2f4..30adb58a3 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(physicalspaces[i - 1], virtualspaces[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]), 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 """ From a620ed4100802382faa9b0f7ca4a41bd02a55ad0 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 23 Oct 2025 14:19:41 -0400 Subject: [PATCH 13/21] remove stray show --- src/environments/abstract_envs.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/environments/abstract_envs.jl b/src/environments/abstract_envs.jl index 54c20dee4..1ec89a5ee 100644 --- a/src/environments/abstract_envs.jl +++ b/src/environments/abstract_envs.jl @@ -82,7 +82,7 @@ function environment_alg( verbosity = Defaults.VERBOSE_NONE ) max_krylovdim = dim(left_virtualspace(above, 1)) * dim(left_virtualspace(below, 1)) - return GMRES(; tol, maxiter, krylovdim = @show(min(max_krylovdim, krylovdim)), verbosity) + return GMRES(; tol, maxiter, krylovdim = min(max_krylovdim, krylovdim), verbosity) end function environment_alg( ::Union{InfiniteQP, MultilineQP}, ::Union{InfiniteMPO, MultilineMPO}, From 75b0c7a92466656752440255a3323a21fa769f5c Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 23 Oct 2025 14:19:48 -0400 Subject: [PATCH 14/21] fix off-by-one error --- src/states/abstractmps.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/states/abstractmps.jl b/src/states/abstractmps.jl index 30adb58a3..6dd8a1341 100644 --- a/src/states/abstractmps.jl +++ b/src/states/abstractmps.jl @@ -135,14 +135,14 @@ function makefullrank!(virtualspaces::PeriodicVector{S}, physicalspaces::Periodi while haschanged haschanged = false for i in 1:length(virtualspaces) - Vmax = fuse(physicalspaces[i - 1], virtualspaces[i - 1]) + 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]), virtualspaces[i]) + Vmax = fuse(dual(physicalspaces[i - 1]), virtualspaces[i]) if !(virtualspaces[i - 1] ≾ Vmax) virtualspaces[i - 1] = infimum(virtualspaces[i - 1], Vmax) haschanged = true From 5928d5434c25d1434a551f958df51b9d5df35ad1 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sun, 26 Oct 2025 17:45:33 -0400 Subject: [PATCH 15/21] discard attempts and go back to original implementation --- src/algorithms/changebonds/randexpand.jl | 34 +++++++----------------- 1 file changed, 9 insertions(+), 25 deletions(-) diff --git a/src/algorithms/changebonds/randexpand.jl b/src/algorithms/changebonds/randexpand.jl index edf3b11d0..fae5fc4e8 100644 --- a/src/algorithms/changebonds/randexpand.jl +++ b/src/algorithms/changebonds/randexpand.jl @@ -14,40 +14,24 @@ $(TYPEDFIELDS) "algorithm used for [truncation](@extref MatrixAlgebraKit.TruncationStrategy] the expanded space" trscheme::TruncationStrategy - - "expansion range that is considered for selecting the orthogonal subspace" - range::Int = 1 end function changebonds!(ψ::InfiniteMPS, alg::RandExpand) - # obtain the spaces for the expanded directions by sampling - virtualspaces = PeriodicVector( - map(1:length(ψ)) do i - Vmax = mapreduce(fuse, reverse(0:alg.range)) do j - j == alg.range ? left_virtualspace(ψ, i - j) : physicalspace(ψ, i - j - 1) - end ⊖ left_virtualspace(ψ, i) - return sample_space(Vmax, alg.trscheme) - end - ) - - # ensure the resulting tensors are full rank (space-wise) - makefullrank!(virtualspaces, physicalspace(ψ)) - - # add vectors orthogonal to the current state + T = eltype(ψ.AL) AL′ = similar(ψ.AL) - T = eltype(AL′) AR′ = similar(ψ.AR, tensormaptype(spacetype(T), 1, numind(T) - 1, storagetype(T))) for i in 1:length(ψ) + # obtain spaces by sampling the support of both the left and right nullspace VL = left_null(ψ.AL[i]) - XL = similar(VL, right_virtualspace(VL) ← virtualspaces[i + 1]) - foreach(((c, b),) -> TensorKit.one!(b), blocks(XL)) - AL′[i] = VL * XL + VR = right_null!(_transpose_tail(ψ.AR[i + 1]; copy = true)) + V = sample_space(infimum(right_virtualspace(VL), space(VR, 1)), alg.trscheme) - VR = right_null!(_transpose_tail(ψ.AR[i])) - XR = similar(VR, virtualspaces[i] ← space(VR, 1)) - foreach(((c, b),) -> TensorKit.one!(b), blocks(XR)) - AR′[i] = XR * VR + # obtain (orthogonal) directions as isometries in that direction + XL = randisometry(storagetype(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′) From 47b9044de42ea374f87c73a7285620af8800c045 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sun, 26 Oct 2025 17:50:26 -0400 Subject: [PATCH 16/21] update changelog --- CHANGELOG.md | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 810c8d927..4923ebceb 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,10 +13,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Changed -- The `changebonds(state, ::RandExpand)` algorithm now no longer projects onto the nullspace - on both sides. This ensures that the expanded symmetry sectors can be selected beyond what - is allowed by two-site updates, which can be relevant for certain systems that have - symmetry-related obstructions. +- 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 From 61eb366cf685fceb18b8ff1ddf1f17f1524fa65f Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sun, 26 Oct 2025 17:53:39 -0400 Subject: [PATCH 17/21] update docs --- docs/src/man/algorithms.md | 4 +--- src/algorithms/changebonds/randexpand.jl | 4 +++- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/src/man/algorithms.md b/docs/src/man/algorithms.md index 384295d10..261bf133d 100644 --- a/docs/src/man/algorithms.md +++ b/docs/src/man/algorithms.md @@ -292,9 +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. Additionally, this - method does not project onto the local two-site basis, which might overcome symmetry-based - obstructions that can be encountered in the other methods. + `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/algorithms/changebonds/randexpand.jl b/src/algorithms/changebonds/randexpand.jl index fae5fc4e8..813397eff 100644 --- a/src/algorithms/changebonds/randexpand.jl +++ b/src/algorithms/changebonds/randexpand.jl @@ -2,7 +2,9 @@ $(TYPEDEF) An algorithm that expands the bond dimension by adding random unitary vectors that are -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 From 13232ba1b64294304a6a41d3a87d4d7e3ff8d8e8 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sun, 26 Oct 2025 17:54:59 -0400 Subject: [PATCH 18/21] fix in-place for multiline --- src/algorithms/changebonds/randexpand.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/algorithms/changebonds/randexpand.jl b/src/algorithms/changebonds/randexpand.jl index 813397eff..4364b0c15 100644 --- a/src/algorithms/changebonds/randexpand.jl +++ b/src/algorithms/changebonds/randexpand.jl @@ -40,7 +40,8 @@ function changebonds!(ψ::InfiniteMPS, alg::RandExpand) end function changebonds!(ψ::MultilineMPS, alg::RandExpand) - return Multiline(map(x -> changebonds!(x, alg), ψ.data)) + foreach(Base.Fix2(changebonds!, alg), ψ.data) + return ψ end changebonds(ψ::AbstractMPS, alg::RandExpand) = changebonds!(copy(ψ), alg) From 11257fbf5e1fe83896002bcaeb976becff711b6e Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sun, 26 Oct 2025 18:02:14 -0400 Subject: [PATCH 19/21] more small fixes and changelog updates --- CHANGELOG.md | 8 ++++++++ LocalPreferences.toml | 2 ++ 2 files changed, 10 insertions(+) create mode 100644 LocalPreferences.toml diff --git a/CHANGELOG.md b/CHANGELOG.md index 4923ebceb..26cc89842 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,8 +9,16 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### 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 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 From 24624a1a80b8e47bb357d23a486caaa58d11bd44 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sun, 26 Oct 2025 18:41:19 -0400 Subject: [PATCH 20/21] add new expansion algorithm --- src/MPSKit.jl | 2 +- src/algorithms/changebonds/randexpand.jl | 57 ++++++++++++++++++++++++ 2 files changed, 58 insertions(+), 1 deletion(-) 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 4364b0c15..71edc0ed9 100644 --- a/src/algorithms/changebonds/randexpand.jl +++ b/src/algorithms/changebonds/randexpand.jl @@ -77,3 +77,60 @@ 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!(ψ::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(storagetype(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: + ψ.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 = :LR, alg.alg_gauge.maxiter, alg.alg_gauge.tol) + mul!.(ψ.AC, ψ.AL, ψ.C) + + return ψ +end From edba5a4ea057b99e3098b1ef75c4146226e14cea Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sun, 26 Oct 2025 20:20:06 -0400 Subject: [PATCH 21/21] fixes and improvements --- src/algorithms/changebonds/randexpand.jl | 43 ++++++++++++++++++++---- 1 file changed, 37 insertions(+), 6 deletions(-) diff --git a/src/algorithms/changebonds/randexpand.jl b/src/algorithms/changebonds/randexpand.jl index 71edc0ed9..94c39eb7b 100644 --- a/src/algorithms/changebonds/randexpand.jl +++ b/src/algorithms/changebonds/randexpand.jl @@ -30,7 +30,7 @@ function changebonds!(ψ::InfiniteMPS, alg::RandExpand) V = sample_space(infimum(right_virtualspace(VL), space(VR, 1)), alg.trscheme) # obtain (orthogonal) directions as isometries in that direction - XL = randisometry(storagetype(VL), right_virtualspace(VL) ← V) + 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 @@ -73,8 +73,8 @@ function changebonds!(ψ::AbstractFiniteMPS, alg::RandExpand) end function sample_space(V, strategy) - S = TensorKit.SectorDict(c => Random.randexp(dim(V, c)) for c in sectors(V)) - ind = MatrixAlgebraKit.findtruncated(S, 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 @@ -110,7 +110,7 @@ function changebonds!(ψ::InfiniteMPS, alg::RandPerturbedExpand) # add (orthogonal) directions as isometries in that direction VL = left_null(ψ.AL[i]) V = sample_space(right_virtualspace(VL), alg.trscheme) - XL = randisometry(storagetype(VL), right_virtualspace(VL) ← V) + 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 @@ -122,15 +122,46 @@ function changebonds!(ψ::InfiniteMPS, alg::RandPerturbedExpand) end # properly regauge the state: + makefullrank!(ψ.AL) ψ.AR .= similar.(ψ.AL) - ψ.AC .= 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 = :LR, alg.alg_gauge.maxiter, alg.alg_gauge.tol) + 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)