From f9d4507fe965877fe9c7b65ff319ca332639ca52 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 10 Sep 2025 14:10:59 -0400 Subject: [PATCH 1/5] Support swap gate for MPOs --- src/operators/mpo.jl | 32 ++++++++++++++++++++++++++++++++ 1 file changed, 32 insertions(+) diff --git a/src/operators/mpo.jl b/src/operators/mpo.jl index c152c578a..63f24bdf5 100644 --- a/src/operators/mpo.jl +++ b/src/operators/mpo.jl @@ -390,3 +390,35 @@ function Base.isapprox( # don't take square roots to avoid precision loss return norm₁₂² ≤ max(atol^2, rtol^2 * max(norm₁², norm₂²)) end + +""" + swap(mpo::FiniteMPO, i::Integer; inv::Bool=false, alg=SDD(), trscheme) + swap!(mpo::FiniteMPO, i::Integer; inv::Bool=false, alg=SDD(), trscheme) + +Compose the mpo with a swap gate applied to indices `i` and `i + 1`, effectively creating an +operator that acts on the Hilbert spaces with those factors swapped. +The keyword arguments `alg` and `trscheme` can be used to control how the resulting tensor +is truncated again. +""" +swap(mpo::FiniteMPO{<:MPOTensor}, i::Integer; kwargs...) = swap!(copy(mpo), i; kwargs...) +function swap!( + mpo::FiniteMPO{<:MPOTensor}, i::Integer; + inv::Bool = false, + alg = SDD(), trscheme = truncbelow(eps(real(scalartype(mpo)))^(4 / 5)) + ) + O₁, O₂ = mpo[i], mpo[i + 1] + + if inv + @plansor O₂₁[-1 -2 -3; -4 -5 -6] := + τ'[-3 -6; 4 5] * O₁[-2 4; 2 1] * O₂[1 5; 3 -5] * τ[2 3; -1 -4] + else + @plansor O₂₁[-1 -2 -3; -4 -5 -6] := + τ[-3 -6; 4 5] * O₁[-2 4; 2 1] * O₂[1 5; 3 -5] * τ'[2 3; -1 -4] + end + + U, S, Vᴴ, = tsvd!(O₂₁; alg, trunc = trscheme) + sqrtS = sqrt(S) + @plansor mpo[i][-1 -2; -3 -4] := U[-3 -1 -2; 1] * sqrtS[1; -4] + @plansor mpo[i + 1][-1 -2; -3 -4] := sqrtS[-1; 1] * Vᴴ[1; -3 -4 -2] + return mpo +end From e9b2b1d8bba2233cf301a5e568852649affd2d7b Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 10 Sep 2025 17:21:29 -0400 Subject: [PATCH 2/5] support for non-sorted indices --- src/operators/mpohamiltonian.jl | 18 +++++++++++++++++- 1 file changed, 17 insertions(+), 1 deletion(-) diff --git a/src/operators/mpohamiltonian.jl b/src/operators/mpohamiltonian.jl index 476983658..53f26df73 100644 --- a/src/operators/mpohamiltonian.jl +++ b/src/operators/mpohamiltonian.jl @@ -313,10 +313,16 @@ function instantiate_operator(lattice::AbstractArray{<:VectorSpace}, (inds′, O mpo = O isa FiniteMPO ? O : FiniteMPO(O) # convert to linear index type - operators = parent(mpo) indices = map(inds) do I return Base._to_linear_index(lattice, Tuple(I)...) # this should mean all inds are valid... end + + # sort indices and deduplicate + indices, mpo = canonicalize_indices(indices, mpo) + operators = parent(mpo) + + @assert allunique(indices) && issorted(indices) "From here on we require unique and ascending indices\n$indices" + T = eltype(mpo) local_mpo = Union{T, scalartype(T)}[] sites = Int[] @@ -339,6 +345,16 @@ function instantiate_operator(lattice::AbstractArray{<:VectorSpace}, (inds′, O return sites => local_mpo end +function canonicalize_indices(indices, mpo) + issorted(indices) && return indices, mpo + mpo = copy(mpo) + I = TensorKit.TupleTools.sortperm(indices) + for s in TensorKit.permutation2swaps(I) + swap!(mpo, s) + end + return TensorKit.TupleTools.getindices(indices, I), mpo +end + # yields the promoted tensortype of all tensors function _find_tensortype(nonzero_operators::AbstractArray) return mapreduce(promote_type, nonzero_operators) do x From d6f13103d676fbd7326632382c47166576d19c0b Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sun, 21 Sep 2025 20:17:52 +0200 Subject: [PATCH 3/5] handle repeated indices --- src/operators/mpo.jl | 8 ++++++++ src/operators/mpohamiltonian.jl | 34 +++++++++++++++++++++------------ 2 files changed, 30 insertions(+), 12 deletions(-) diff --git a/src/operators/mpo.jl b/src/operators/mpo.jl index 63f24bdf5..9277318d0 100644 --- a/src/operators/mpo.jl +++ b/src/operators/mpo.jl @@ -422,3 +422,11 @@ function swap!( @plansor mpo[i + 1][-1 -2; -3 -4] := sqrtS[-1; 1] * Vᴴ[1; -3 -4 -2] return mpo end + +function multiply_neighbours!(mpo::FiniteMPO{<:MPOTensor}, i::Integer) + 1 <= i < length(mpo) || throw(BoundsError(mpo, i)) + O₁ = mpo[i] + O₂ = popat!(parent(mpo), i + 1) + @plansor mpo[i][-1 -2; -3 -4] := O₁[-1 -2; 1 2] * τ[1 2; 3 4] * O₂[3 4; -3 -4] + return mpo +end diff --git a/src/operators/mpohamiltonian.jl b/src/operators/mpohamiltonian.jl index 53f26df73..3e04ed5fe 100644 --- a/src/operators/mpohamiltonian.jl +++ b/src/operators/mpohamiltonian.jl @@ -309,16 +309,17 @@ function instantiate_operator(state::AbstractMPS, O::Pair) return instantiate_operator(physicalspace(state), O) end function instantiate_operator(lattice::AbstractArray{<:VectorSpace}, (inds′, O)::Pair) - inds = inds′ isa Int ? tuple(inds′) : inds′ - mpo = O isa FiniteMPO ? O : FiniteMPO(O) + inds = inds′ isa Int ? [inds′] : inds′ + mpo = O isa FiniteMPO ? copy(O) : FiniteMPO(O) # convert to linear index type - indices = map(inds) do I - return Base._to_linear_index(lattice, Tuple(I)...) # this should mean all inds are valid... + indices = Vector{Int}(undef, length(inds)) + for i in eachindex(indices) + indices[i] = Base._to_linear_index(lattice, Tuple(inds[i])...) # this should mean all inds are valid... end # sort indices and deduplicate - indices, mpo = canonicalize_indices(indices, mpo) + indices, mpo = canonicalize_indices!(indices, mpo) operators = parent(mpo) @assert allunique(indices) && issorted(indices) "From here on we require unique and ascending indices\n$indices" @@ -345,14 +346,23 @@ function instantiate_operator(lattice::AbstractArray{<:VectorSpace}, (inds′, O return sites => local_mpo end -function canonicalize_indices(indices, mpo) - issorted(indices) && return indices, mpo - mpo = copy(mpo) - I = TensorKit.TupleTools.sortperm(indices) - for s in TensorKit.permutation2swaps(I) - swap!(mpo, s) +function canonicalize_indices!(indices, mpo) + # swap non-sorted entries + for i in 2:length(indices) + for j in reverse(i:length(indices)) + if indices[j] < indices[j - 1] + swap!(mpo, j - 1) + indices[j - 1], indices[j] = indices[j], indices[j - 1] + end + end + end + for i in length(indices):-1:2 + if indices[i] == indices[i - 1] + multiply_neighbours!(mpo, i - 1) + popat!(indices, i) + end end - return TensorKit.TupleTools.getindices(indices, I), mpo + return indices, mpo end # yields the promoted tensortype of all tensors From 8d8154138c07216325bc7b8e7e81c8cdf5d359d0 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Mon, 22 Sep 2025 08:33:50 +0200 Subject: [PATCH 4/5] update docstrings --- src/operators/mpo.jl | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/src/operators/mpo.jl b/src/operators/mpo.jl index 9277318d0..3fb28fa29 100644 --- a/src/operators/mpo.jl +++ b/src/operators/mpo.jl @@ -391,7 +391,7 @@ function Base.isapprox( return norm₁₂² ≤ max(atol^2, rtol^2 * max(norm₁², norm₂²)) end -""" +@doc """ swap(mpo::FiniteMPO, i::Integer; inv::Bool=false, alg=SDD(), trscheme) swap!(mpo::FiniteMPO, i::Integer; inv::Bool=false, alg=SDD(), trscheme) @@ -399,8 +399,9 @@ Compose the mpo with a swap gate applied to indices `i` and `i + 1`, effectively operator that acts on the Hilbert spaces with those factors swapped. The keyword arguments `alg` and `trscheme` can be used to control how the resulting tensor is truncated again. -""" -swap(mpo::FiniteMPO{<:MPOTensor}, i::Integer; kwargs...) = swap!(copy(mpo), i; kwargs...) +""" swap, swap! + +swap(mpo::FiniteMPO, i::Integer; kwargs...) = swap!(copy(mpo), i; kwargs...) function swap!( mpo::FiniteMPO{<:MPOTensor}, i::Integer; inv::Bool = false, @@ -423,6 +424,15 @@ function swap!( return mpo end +@doc """ + multiply_neighbours(mpo::FiniteMPO, i::Integer) + multiply_neighbours!(mpo::FiniteMPO, i::Integer) + +Construct the mpo of length `length(mpo) - 1` which is formed by multiplying the operators +on site `i` and `i + 1`. +""" multiply_neighbours, multiply_neighbours! + +multiply_neighbours(mpo::FiniteMPO, i::Integer) = multiply_neighbours!(copy(mpo), i) function multiply_neighbours!(mpo::FiniteMPO{<:MPOTensor}, i::Integer) 1 <= i < length(mpo) || throw(BoundsError(mpo, i)) O₁ = mpo[i] From 6bd55828f699c7bdb23b0d3660c511f60e0b05c1 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 25 Sep 2025 17:28:51 -0400 Subject: [PATCH 5/5] Add tests repeated indices --- test/operators.jl | 37 +++++++++++++++++++++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/test/operators.jl b/test/operators.jl index 61e26c67d..2b8bb7f7c 100644 --- a/test/operators.jl +++ b/test/operators.jl @@ -218,6 +218,43 @@ module TestOperators end end + @testset "Finite MPOHamiltonian repeated indices" begin + X = randn(ComplexF64, ℂ^2, ℂ^2) + X += X' + Y = randn(ComplexF64, ℂ^2, ℂ^2) + Y += Y' + L = 4 + chain = fill(space(X, 1), 4) + + H1 = FiniteMPOHamiltonian(chain, (1,) => (X * X * Y * Y)) + H2 = FiniteMPOHamiltonian(chain, (1, 1, 1, 1) => (X ⊗ X ⊗ Y ⊗ Y)) + @test convert(TensorMap, H1) ≈ convert(TensorMap, H2) + + H1 = FiniteMPOHamiltonian(chain, (1, 2) => ((X * Y) ⊗ (X * Y))) + H2 = FiniteMPOHamiltonian(chain, (1, 2, 1, 2) => (X ⊗ X ⊗ Y ⊗ Y)) + @test convert(TensorMap, H1) ≈ convert(TensorMap, H2) + + H1 = FiniteMPOHamiltonian(chain, (1, 2) => ((X * X * Y) ⊗ Y)) + H2 = FiniteMPOHamiltonian(chain, (1, 1, 1, 2) => (X ⊗ X ⊗ Y ⊗ Y)) + @test convert(TensorMap, H1) ≈ convert(TensorMap, H2) + + H1 = FiniteMPOHamiltonian(chain, (1, 2) => ((X * Y * Y) ⊗ X)) + H2 = FiniteMPOHamiltonian(chain, (1, 2, 1, 1) => (X ⊗ X ⊗ Y ⊗ Y)) + @test convert(TensorMap, H1) ≈ convert(TensorMap, H2) + + H1 = FiniteMPOHamiltonian(chain, (1, 2, 3) => FiniteMPO((X * X) ⊗ Y ⊗ Y)) + H2 = FiniteMPOHamiltonian(chain, (1, 1, 2, 3) => FiniteMPO(X ⊗ X ⊗ Y ⊗ Y)) + @test convert(TensorMap, H1) ≈ convert(TensorMap, H2) + + H1 = FiniteMPOHamiltonian(chain, (1, 2, 3) => FiniteMPO((Y * Y) ⊗ X ⊗ X)) + H2 = FiniteMPOHamiltonian(chain, (2, 3, 1, 1) => FiniteMPO(X ⊗ X ⊗ Y ⊗ Y)) + @test convert(TensorMap, H1) ≈ convert(TensorMap, H2) + + H1 = FiniteMPOHamiltonian(chain, (1, 2, 3) => FiniteMPO(X ⊗ (X * Y) ⊗ Y)) + H2 = FiniteMPOHamiltonian(chain, (1, 2, 2, 3) => FiniteMPO(X ⊗ X ⊗ Y ⊗ Y)) + @test convert(TensorMap, H1) ≈ convert(TensorMap, H2) + end + @testset "InfiniteMPOHamiltonian $(sectortype(pspace))" for (pspace, Dspace) in zip(pspaces, vspaces) # generate a 1-2-3 body interaction operators = ntuple(3) do i