From 2af2bb8ac0ca328aa2807a0051fd1aaef6c4b323 Mon Sep 17 00:00:00 2001 From: leburgel Date: Tue, 22 Sep 2026 13:23:51 +0200 Subject: [PATCH 1/2] Add support for local expectation values of tensor product and MPO terms --- .../contractions/local_patch/network_expr.jl | 37 ++++ .../expectation_value/expectation_value.jl | 16 ++ .../expectation_value/patch_contractions.jl | 18 ++ src/operators/localoperator.jl | 183 +++++++++++++++++- test/testsuite/TestSuite.jl | 13 ++ test/testsuite/toolbox/mpo_terms.jl | 168 ++++++++++++++++ test/testsuite/toolbox/tensorproduct_terms.jl | 99 ++++++++++ test/toolbox/mpo_terms.jl | 15 ++ test/toolbox/tensorproduct_terms.jl | 12 ++ 9 files changed, 558 insertions(+), 3 deletions(-) create mode 100644 test/testsuite/toolbox/mpo_terms.jl create mode 100644 test/testsuite/toolbox/tensorproduct_terms.jl create mode 100644 test/toolbox/mpo_terms.jl create mode 100644 test/toolbox/tensorproduct_terms.jl diff --git a/src/algorithms/contractions/local_patch/network_expr.jl b/src/algorithms/contractions/local_patch/network_expr.jl index 6120fa3b3..7e3e39c32 100644 --- a/src/algorithms/contractions/local_patch/network_expr.jl +++ b/src/algorithms/contractions/local_patch/network_expr.jl @@ -260,3 +260,40 @@ function operator_contraction_expr(::Type{<:AbstractTensorMap}, nsites) ), ] end + +# The MPO bond is small - the operator's Schmidt rank, e.g. 3 for Heisenberg XYZ - so it is +# labelled as a physical rather than a virtual dimension, which gives `@autoopt` the right +# order of magnitude when it searches for a contraction order. +mpolabel(args...) = physicallabel(:mpo, args...) + +function operator_contraction_expr(::Type{<:MPOTerm}, nsites) + # one factor per site, linked by a chain of bonds running left to right, rank-2 at the + # ends of the chain and rank-4 in the bulk: + # W₁ : bra₁ ← ket₁ ⊗ b₁ + # Wᵢ : bᵢ₋₁ ⊗ braᵢ ← ketᵢ ⊗ bᵢ + # W_N : b_{N-1} ⊗ bra_N ← ket_N + return map(1:nsites) do i + bra, ket = physicallabel(:O, 2, i), physicallabel(:O, 1, i) + out, ins = if nsites == 1 + (bra,), (ket,) + elseif i == 1 + (bra,), (ket, mpolabel(1)) + elseif i == nsites + (mpolabel(nsites - 1), bra), (ket,) + else + (mpolabel(i - 1), bra), (ket, mpolabel(i)) + end + return tensorexpr(:(operator[$i]), out, ins) + end +end + +# A tensor product term is the special case in which every bond has dimension 1, so its +# factors carry no bonds at all: each simply sits between the ket and bra index of its site. +# This method is more specific than the `MPOTerm` one above, so it takes precedence. +function operator_contraction_expr(::Type{<:TensorProductTerm}, nsites) + return map(1:nsites) do i + return tensorexpr( + :(operator[$i]), (physicallabel(:O, 2, i),), (physicallabel(:O, 1, i),) + ) + end +end diff --git a/src/algorithms/expectation_value/expectation_value.jl b/src/algorithms/expectation_value/expectation_value.jl index d2ea498f6..fe0a3b1a2 100644 --- a/src/algorithms/expectation_value/expectation_value.jl +++ b/src/algorithms/expectation_value/expectation_value.jl @@ -54,6 +54,22 @@ function local_expectation_value(inds, state, operator::AbstractTensorMap, env) return trmul(operator, ρ) end +""" +$(SIGNATURES) + +Compute the contribution of an [`MPOTerm`](@ref) - and hence also of a +[`TensorProductTerm`](@ref) - given as one tensor per site +in `inds`, to the expectation value ⟨bra|O|ket⟩ / ⟨bra|ket⟩. + +Rather than forming the dense operator and tracing it against a reduced density matrix, the +factors are inserted into the patch contraction directly, and the result is divided by the +local norm of the same patch. +""" +function local_expectation_value(inds, bra, operator::MPOTerm, ket, env) + return contract_local_operator(inds, operator, ket, bra, env) / + contract_local_norm(inds, ket, bra, env) +end + # Expectation value of a local partition function tensor # ------------------------------------------------------ diff --git a/src/algorithms/expectation_value/patch_contractions.jl b/src/algorithms/expectation_value/patch_contractions.jl index 324d2ebb5..a1996743b 100644 --- a/src/algorithms/expectation_value/patch_contractions.jl +++ b/src/algorithms/expectation_value/patch_contractions.jl @@ -19,6 +19,24 @@ function contract_local_operator( static_inds = Tuple(Val.(inds)) return _contract_local_operator(static_inds, O, (ket, bra), env) end +function contract_local_operator( + inds::Vector{CartesianIndex{2}}, O::MPOTerm, + ket::InfinitePEPS, bra::InfinitePEPS, env + ) + length(inds) == length(O) || + throw(ArgumentError("Got $(length(inds)) sites but $(length(O)) MPO factors.")) + static_inds = Tuple(Val.(inds)) + return _contract_local_operator(static_inds, O, (ket, bra), env) +end +function contract_local_operator( + inds::Vector{CartesianIndex{2}}, O::TensorProductTerm, + ket::InfinitePEPS, bra::InfinitePEPS, env + ) + length(inds) == length(O) || + throw(ArgumentError("Got $(length(inds)) sites but $(length(O)) product factors.")) + static_inds = Tuple(Val.(inds)) + return _contract_local_operator(static_inds, O, (ket, bra), env) +end function contract_local_operator( inds::Vector{CartesianIndex{2}}, O, state::InfinitePEPO, env ) diff --git a/src/operators/localoperator.jl b/src/operators/localoperator.jl index c35e5ddb9..7f24c47b6 100644 --- a/src/operators/localoperator.jl +++ b/src/operators/localoperator.jl @@ -87,6 +87,152 @@ function add_term!( return operator end +# Tensor product terms +# -------------------- +# A term given as one rank-2 operator per site, i.e. an explicit tensor product which is +# never actually formed. + +""" + TensorProductTerm + +A single term of a [`LocalOperator`](@ref) given as a tensor product of one rank-2 operator +per site acted on, rather than as the rank-`2N` tensor product itself. + +Only operators whose terms are individually tensor products can be represented this way. In +particular several such terms acting on the same set of sites cannot be combined, since a +sum of tensor products is not itself a tensor product. + +The factors are rank-2, carrying no bond indices, which is what distinguishes this from the +more general [`MPOTerm`](@ref): both are vectors of tensors, but a tensor product term is the +narrower type and therefore wins on dispatch wherever it applies. +""" +const TensorProductTerm{T} = AbstractVector{T} where {T <: AbstractTensorMap{<:Any, <:Any, 1, 1}} + +add_term!(operator::LocalOperator, inds::Tuple, term::TensorProductTerm) = + add_term!(operator, collect(inds), term) +add_term!(operator::LocalOperator, inds::Vector, term::TensorProductTerm) = + add_term!(operator, map(CartesianIndex{2}, inds), term) +function add_term!( + operator::LocalOperator, inds::Vector{CartesianIndex{2}}, term::TensorProductTerm; + atol = zero(real(scalartype(first(term)))), + ) + # input checks + length(inds) == length(term) || + throw(ArgumentError("Incompatible number of indices and tensor product factors")) + allunique(inds) || throw(ArgumentError("`inds` should not contain repeated coordinates.")) + for (i, ind) in enumerate(inds) + numin(term[i]) == numout(term[i]) == 1 || + throw(ArgumentError("Tensor product factors should be single-site operators")) + ind_translated = CartesianIndex(mod1.(Tuple(ind), size(operator))) + physicalspace(operator, ind_translated) == domain(term[i])[1] == codomain(term[i])[1] || + throw(SpaceMismatch("Incompatible physical spaces")) + end + prod(norm, term) <= atol && return operator # skip adding negligible terms + + # permute input, which for a product just reorders the factors along with their sites + if !issorted(inds) + I = sortperm(inds) + inds = inds[I] + term = term[I] + end + + # translate coordinates + _shift_into_unitcell!(inds, size(operator)) + + # a sum of tensor products is not a tensor product, so terms cannot be accumulated here + haskey(operator.terms, inds) && throw( + ArgumentError( + "A term acting on $inds is already present. Tensor product terms cannot be \ + summed, since a sum of tensor products is not itself a tensor product." + ) + ) + operator.terms[inds] = collect(term) + + return operator +end + +# MPO terms +# --------- +# A term given as a chain of MPO tensors, one per site, linked by virtual bonds which are +# contracted directly between neighbours rather than being fused into the PEPS bonds. This +# generalizes a tensor product term, which is the special case of bond dimension 1. + +""" + MPOTerm{T} + +A single term of a [`LocalOperator`](@ref) given as a matrix product operator with one +tensor per site acted on, rather than as the dense rank-`2N` tensor. + +The factors are ordered along the chain and follow the usual MPO convention, rank-2 at the +ends and rank-4 in the bulk: + + W₁ : P₁ ← P₁ ⊗ B₁ + Wᵢ : Bᵢ₋₁ ⊗ Pᵢ ← Pᵢ ⊗ Bᵢ + W_N : B_{N-1} ⊗ P_N ← P_N + +so the physical index in the codomain is the bra index and the one in the domain is the ket +index, matching the convention of a dense term, and the bonds run left to right. This is what +[`gate_to_mpo`](@ref) produces, which is the way to obtain an `MPOTerm` from a dense operator. +A single-site term is just a rank-2 operator. + +A [`TensorProductTerm`](@ref) is the special case in which every bond is trivial, so its +factors carry no bond indices and are rank-2 throughout. That is exactly what separates the +two on dispatch: a vector of rank-2 operators is a tensor product term, anything else is an +MPO. +""" +const MPOTerm{T} = AbstractVector{T} where {T <: AbstractTensorMap} + + +add_term!(operator::LocalOperator, inds::Tuple, term::MPOTerm) = + add_term!(operator, collect(inds), term) +add_term!(operator::LocalOperator, inds::Vector, term::MPOTerm) = + add_term!(operator, map(CartesianIndex{2}, inds), term) +function add_term!( + operator::LocalOperator, inds::Vector{CartesianIndex{2}}, term::MPOTerm + ) + # input checks + length(inds) == length(term) || + throw(ArgumentError("Incompatible number of indices and MPO factors")) + allunique(inds) || throw(ArgumentError("`inds` should not contain repeated coordinates.")) + n = length(inds) + for (i, ind) in enumerate(inds) + # a factor carries its physical pair plus a bond towards each neighbour it has, so it + # is rank-2 for a lone site, rank-3 at the ends of a chain and rank-4 in the bulk + nout = (n == 1 || i == 1) ? 1 : 2 + nin = (n == 1 || i == n) ? 1 : 2 + (numout(term[i]) == nout && numin(term[i]) == nin) || throw( + ArgumentError( + "MPO factor $i of $n should have $nout index(es) out and $nin in, got \ + $(numout(term[i])) and $(numin(term[i]))." + ) + ) + # the bra index is the physical one in the codomain: last for a bulk factor, which + # carries the incoming bond first, and only for an end factor + bra = codomain(term[i])[nout] + ind_translated = CartesianIndex(mod1.(Tuple(ind), size(operator))) + physicalspace(operator, ind_translated) == bra || + throw(SpaceMismatch("Incompatible physical spaces")) + end + + # NOTE: `inds` is deliberately *not* sorted here, unlike for dense and tensor product + # terms, since permuting MPOs is not straightforward. + + # translate coordinates + _shift_into_unitcell!(inds, size(operator)) + + # as for tensor products, a sum of MPOs of fixed bond dimension is not one of the same + # bond dimension, so terms are not accumulated here + haskey(operator.terms, inds) && throw( + ArgumentError( + "A term acting on $inds is already present. MPO terms cannot be summed in \ + place; combine the operators before splitting them." + ) + ) + operator.terms[inds] = collect(term) + + return operator +end + """ checklattice(Bool, args...) @@ -153,8 +299,22 @@ Base.eltype(::Type{LocalOperator{O, S}}) where {O, S} = O # Real and imaginary part # ----------------------- +""" +$(SIGNATURES) + +Take the real part of a single term of a [`LocalOperator`](@ref). + +A [`TensorProductTerm`](@ref) is made real factor by factor, so that the real part of a +tensor product is the tensor product of the real parts of its factors. An [`MPOTerm`](@ref) +is treated the same way, factor by factor. +""" +_real_local_term(term) = real(term) +_real_local_term(term::MPOTerm) = map(real, term) + function Base.real(O::LocalOperator) - return LocalOperator(O.lattice, (sites => real(op) for (sites, op) in O.terms)...) + return LocalOperator( + O.lattice, (sites => _real_local_term(op) for (sites, op) in O.terms)... + ) end function Base.imag(O::LocalOperator) return LocalOperator(O.lattice, (sites => imag(op) for (sites, op) in O.terms)...) @@ -162,8 +322,25 @@ end # Linear Algebra # -------------- -Base.:*(α::Number, O::LocalOperator) = - LocalOperator(physicalspace(O), inds => α * operator for (inds, operator) in O.terms) +""" +$(SIGNATURES) + +Scale a single term of a [`LocalOperator`](@ref) by `α`. + +Terms which are not plain tensors need not scale by plain multiplication: an +[`MPOTerm`](@ref) - and hence also a [`TensorProductTerm`](@ref) - is scaled by scaling one +of its factors, since multiplying every factor would scale the term it represents by `α^n`. +""" +_scale_local_term(term, α::Number) = α * term +function _scale_local_term(term::MPOTerm, α::Number) + scaled = collect(term) + scaled[1] = α * scaled[1] + return scaled +end + +Base.:*(α::Number, O::LocalOperator) = LocalOperator( + physicalspace(O), inds => _scale_local_term(operator, α) for (inds, operator) in O.terms +) Base.:*(O::LocalOperator, α::Number) = α * O Base.:/(O::LocalOperator, α::Number) = O * inv(α) diff --git a/test/testsuite/TestSuite.jl b/test/testsuite/TestSuite.jl index 018a68a82..c86d7c290 100644 --- a/test/testsuite/TestSuite.jl +++ b/test/testsuite/TestSuite.jl @@ -239,6 +239,19 @@ module ToolboxDensityMatrices end using .ToolboxDensityMatrices +module ToolboxTensorProductTerms + include("toolbox/tensorproduct_terms.jl") + export toolbox_tensorproduct_ising, toolbox_tensorproduct_bookkeeping +end +using .ToolboxTensorProductTerms + +module ToolboxMPOTerms + include("toolbox/mpo_terms.jl") + export toolbox_mpo_convention, toolbox_mpo_dispatch, toolbox_mpo_terms_dense + export toolbox_mpo_heisenberg, toolbox_mpo_bookkeeping +end +using .ToolboxMPOTerms + # Utility # ------- module UtilityCorrelator diff --git a/test/testsuite/toolbox/mpo_terms.jl b/test/testsuite/toolbox/mpo_terms.jl new file mode 100644 index 000000000..26e821782 --- /dev/null +++ b/test/testsuite/toolbox/mpo_terms.jl @@ -0,0 +1,168 @@ +using Test +using Random +using TensorKit +using PEPSKit, Adapt +using PEPSKit: MPOTerm, TensorProductTerm, gate_to_mpo, add_term! +import TensorKitTensors.SpinOperators as SO + +# `LocalOperator` terms given as a matrix product operator, one tensor per site acted on, +# rather than as the dense rank-2N tensor. `gate_to_mpo` turns a dense operator into one. +# +# A tensor product term is the trivial-bond case, and is the narrower type, so both are +# vectors of tensors and dispatch separates them by whether the factors are rank-2. + +const Dbond, χenv = 2, 8 + +const X = SO.σˣ(ComplexF64, Trivial) +const Z = SO.σᶻ(ComplexF64, Trivial) +const lattice11 = fill(ComplexSpace(2), 1, 1) + +const dense_terms = ( + ("NN horizontal", [(1, 1), (1, 2)], -(X ⊗ X) + Z ⊗ Z), + ("NN vertical", [(1, 1), (2, 1)], -(X ⊗ X) + Z ⊗ Z), + ("NNN diagonal", [(1, 1), (2, 2)], X ⊗ Z), + ("3 site line", [(1, 1), (1, 2), (1, 3)], X ⊗ X ⊗ X + Z ⊗ Z ⊗ Z), + # sites are taken in the order given, not sorted, so an unsorted or reversed order + # has to agree with the dense term written the same way round. `X ⊗ Z` is asymmetric + # under exchange, so this would catch a silently transposed pairing. + ("NNN anti-diagonal", [(1, 2), (2, 1)], X ⊗ Z), + ("reversed sites", [(1, 2), (1, 1)], X ⊗ Z), +) + +""" +Largest bond dimension across the cuts of an MPO term, i.e. the Schmidt rank of the operator +it represents. The outgoing bond of a factor is its last domain index, except for a rank-2 +factor, which carries no bond at all. +""" +function mpobond(term) + length(term) == 1 && return 1 + return maximum(1:(length(term) - 1)) do i + W = term[i] + numind(W) == 2 && return 1 + return dim(space(W, numind(W))) + end +end + +function toolbox_mpo_convention(AT) + return @testset "MPO factors follow the standard convention ($AT)" begin + # rank-2 at the ends of the chain, rank-4 in the bulk + W1, W2 = gate_to_mpo(X ⊗ X + Z ⊗ Z) + @test (numout(W1), numin(W1)) == (1, 2) + @test (numout(W2), numin(W2)) == (2, 1) + + Ws = gate_to_mpo(X ⊗ X ⊗ X + Z ⊗ Z ⊗ Z) + @test length(Ws) == 3 + @test (numout(Ws[1]), numin(Ws[1])) == (1, 2) + @test (numout(Ws[2]), numin(Ws[2])) == (2, 2) + @test (numout(Ws[3]), numin(Ws[3])) == (2, 1) + + # the bond is the Schmidt rank across each cut + @test mpobond(gate_to_mpo(X ⊗ X)) == 1 # a pure product + @test mpobond(gate_to_mpo(X ⊗ X + Z ⊗ Z)) == 2 + @test mpobond(gate_to_mpo(X ⊗ X ⊗ X + Z ⊗ Z ⊗ Z)) == 2 + @test mpobond([X]) == 1 # a lone site has no bond + + # Heisenberg XYZ has Schmidt rank 3 at every spin + for spin in (1 // 2, 1, 3 // 2) + term = SO.S_x_S_x(ComplexF64, Trivial; spin) + + SO.S_y_S_y(ComplexF64, Trivial; spin) + + SO.S_z_S_z(ComplexF64, Trivial; spin) + @test mpobond(gate_to_mpo(term)) == 3 + end + end +end + +function toolbox_mpo_dispatch(AT) + return @testset "Dispatch separates products from MPOs ($AT)" begin + # a vector of rank-2 operators is both, and the narrower type wins where it matters + @test [X, Z] isa TensorProductTerm + @test [X, Z] isa MPOTerm + + # anything carrying a bond is only an MPO term + for O in (X ⊗ Z, X ⊗ X + Z ⊗ Z, X ⊗ X ⊗ X + Z ⊗ Z ⊗ Z) + Ws = gate_to_mpo(O) + @test Ws isa MPOTerm + @test !(Ws isa TensorProductTerm) + end + end +end + +function toolbox_mpo_terms_dense(AT) + return @testset "MPO terms reproduce dense terms ($(name)) ($AT)" for (name, inds, O) in + dense_terms + Random.seed!(2985721) + peps = adapt(AT, InfinitePEPS(ComplexSpace(2), ComplexSpace(Dbond))) + env = CTMRGEnv(randn, ComplexF64, peps, ComplexSpace(χenv)) + + H_dense = LocalOperator(lattice11, inds => O) + H_mpo = LocalOperator(lattice11, inds => gate_to_mpo(O)) + @test expectation_value(peps, H_mpo, env) ≈ expectation_value(peps, H_dense, env) rtol = + 1.0e-9 + end +end + +function toolbox_mpo_heisenberg(AT) + return @testset "MPO terms reproduce the Heisenberg XYZ Hamiltonian ($AT)" begin + Random.seed!(2985721) + H = heisenberg_XYZ(InfiniteSquare()) + H_mpo = LocalOperator( + physicalspace(H), (inds => gate_to_mpo(op) for (inds, op) in H.terms)... + ) + @test all(t -> t isa MPOTerm, values(H_mpo.terms)) + @test all(t -> mpobond(t) == 3, values(H_mpo.terms)) + + peps = adapt(AT, InfinitePEPS(ComplexSpace(2), ComplexSpace(Dbond))) + env, = leading_boundary( + CTMRGEnv(peps, ComplexSpace(χenv)), peps; tol = 1.0e-8, verbosity = 0 + ) + E = expectation_value(peps, H, env) + @test expectation_value(peps, H_mpo, env) ≈ E rtol = 1.0e-9 + + # scaling an MPO term scales the term, not each of its factors + α = 0.37 + @test expectation_value(peps, α * H_mpo, env) ≈ α * E rtol = 1.0e-9 + + # the real part is taken factor by factor + H_real = real(H_mpo) + @test all(t -> t isa MPOTerm, values(H_real.terms)) + @test all(all(W -> norm(imag(W)) < 1.0e-14, t) for t in values(H_real.terms)) + end +end + +function toolbox_mpo_bookkeeping(AT) + return @testset "MPO term bookkeeping ($AT)" begin + # sites are stored in the order given, and are *not* canonicalized by sorting the way + # dense and tensor product terms are: an MPO's bonds are the Schmidt cuts of one + # particular chain ordering, so its factors cannot be permuted after the fact. Factor i + # therefore acts on site inds[i], and the caller owns that pairing. + O_fwd = LocalOperator(lattice11, [(1, 1), (1, 2)] => gate_to_mpo(X ⊗ Z)) + O_rev = LocalOperator(lattice11, [(1, 2), (1, 1)] => gate_to_mpo(X ⊗ Z)) + # the sites are shifted into the unit cell as a block, anchored on the first of them, so + # absolute coordinates depend on which site comes first; what is preserved, and what + # encodes the chain order, is the displacement between consecutive sites + @test diff(only(keys(O_fwd.terms))) == [CartesianIndex(0, 1)] + @test diff(only(keys(O_rev.terms))) == [CartesianIndex(0, -1)] + + # a mismatch between sites and factors is caught + @test_throws ArgumentError LocalOperator( + lattice11, [(1, 1), (1, 2), (1, 3)] => gate_to_mpo(X ⊗ Z) + ) + + # factor ranks must match their position in the chain + @test_throws ArgumentError LocalOperator(lattice11, [(1, 1)] => [X ⊗ X]) + @test_throws ArgumentError LocalOperator(lattice11, [(1, 1)] => gate_to_mpo(X ⊗ Z)) + + # an MPO is not summed in place, since the bond dimension of a sum is not the same + O = LocalOperator(lattice11, [(1, 1), (1, 2)] => gate_to_mpo(X ⊗ Z)) + @test_throws ArgumentError add_term!( + O, [CartesianIndex(1, 1), CartesianIndex(1, 2)], gate_to_mpo(Z ⊗ X) + ) + + # physical spaces are checked + @test_throws SpaceMismatch LocalOperator( + lattice11, [(1, 1), (1, 2)] => gate_to_mpo( + SO.S_x_S_x(ComplexF64, Trivial; spin = 1) + ) + ) + end +end diff --git a/test/testsuite/toolbox/tensorproduct_terms.jl b/test/testsuite/toolbox/tensorproduct_terms.jl new file mode 100644 index 000000000..b651fc349 --- /dev/null +++ b/test/testsuite/toolbox/tensorproduct_terms.jl @@ -0,0 +1,99 @@ +using Test +using Random +using TensorKit +using PEPSKit, Adapt +using PEPSKit: TensorProductTerm, add_term!, _scale_local_term +import TensorKitTensors.SpinOperators as SO + +# `LocalOperator` terms given as an explicit tensor product of single-site operators, one +# per site acted on, rather than as the tensor product formed up front. +# +# The transverse-field Ising model is the natural test case: with trivial symmetry every set +# of sites carries exactly one term, so the whole Hamiltonian is a sum of tensor products, +# a two-factor product per nearest neighbor bond and a one-factor product per site. This +# requires trivial symmetry, since a single-site σᶻ is not symmetric under the Z₂ symmetry +# of the model. + +const Dbond, χenv, g = 2, 8, 3.1 +const unitcells = ((1, 1), (2, 2)) + +""" +Transverse-field Ising Hamiltonian built from tensor product terms, matching the conventions +of [`transverse_field_ising`](@ref). +""" +function transverse_field_ising_tensorproduct( + T::Type{<:Number}, lattice::InfiniteSquare; J = 1.0, g = 1.0 + ) + Z, X = SO.σᶻ(T, Trivial), SO.σˣ(T, Trivial) + spaces = fill(domain(X)[1], (lattice.Nrows, lattice.Ncols)) + return LocalOperator( + spaces, + (neighbor => [-J * Z, copy(Z)] for neighbor in nearest_neighbours(lattice))..., + ([idx] => [(-J * g) * X] for idx in vertices(lattice))..., + ) +end +transverse_field_ising_tensorproduct(lattice::InfiniteSquare; kwargs...) = + transverse_field_ising_tensorproduct(ComplexF64, lattice; kwargs...) + +function toolbox_tensorproduct_ising(AT) + return @testset "Tensor product terms reproduce the two-site Ising Hamiltonian ($uc) ($AT)" for uc in + unitcells + Random.seed!(2985721) + H = transverse_field_ising(InfiniteSquare(uc...); g) + H_prod = transverse_field_ising_tensorproduct(InfiniteSquare(uc...); g) + + # one term per site set: a two-factor product per bond, a one-factor product per site + @test length(H_prod.terms) == length(H.terms) + @test all(t -> t isa TensorProductTerm, values(H_prod.terms)) + @test sort(length.(collect(values(H_prod.terms)))) == + sort(vcat(fill(1, prod(uc)), fill(2, 2 * prod(uc)))) + + peps = adapt(AT, InfinitePEPS(ComplexSpace(2), ComplexSpace(Dbond); unitcell = uc)) + env, = leading_boundary( + CTMRGEnv(peps, ComplexSpace(χenv)), peps; tol = 1.0e-8, verbosity = 0 + ) + + @test expectation_value(peps, H_prod, env) ≈ expectation_value(peps, H, env) rtol = + 1.0e-9 + + # scaling a product term scales the term, not each of its factors + α = 0.37 + @test expectation_value(peps, α * H_prod, env) ≈ + α * expectation_value(peps, H_prod, env) rtol = 1.0e-9 + end +end + +function toolbox_tensorproduct_bookkeeping(AT) + return @testset "Tensor product term bookkeeping ($AT)" begin + lattice = fill(ComplexSpace(2), 1, 1) + X, Z = SO.σˣ(ComplexF64, Trivial), SO.σᶻ(ComplexF64, Trivial) + + # a sum of tensor products is not a tensor product, so terms may not be accumulated + O = LocalOperator(lattice, [(1, 1), (1, 2)] => [X, X]) + @test_throws ArgumentError add_term!(O, [(1, 1), (1, 2)], [Z, Z]) + + # factors are reordered along with their sites + O1 = LocalOperator(lattice, [(1, 1), (1, 2)] => [X, Z]) + O2 = LocalOperator(lattice, [(1, 2), (1, 1)] => [Z, X]) + @test only(keys(O1.terms)) == only(keys(O2.terms)) + @test all(map(≈, only(values(O1.terms)), only(values(O2.terms)))) + + # each factor must be a single-site operator, and match the physical space + @test_throws ArgumentError LocalOperator(lattice, [(1, 1), (1, 2)] => [X]) # arity + @test_throws ArgumentError LocalOperator(lattice, [(1, 1)] => [X ⊗ X]) # not single-site + @test_throws SpaceMismatch LocalOperator( + lattice, [(1, 1)] => [SO.S_x(ComplexF64, Trivial; spin = 1)] + ) + + # scaling touches exactly one factor + scaled = _scale_local_term([X, Z], 3.0) + @test scaled[1] ≈ 3.0 * X + @test scaled[2] ≈ Z + + # the real part of a product is the product of the real parts of its factors + O3 = LocalOperator(lattice, [(1, 1), (1, 2)] => [(2.0 + 3.0im) * Z, copy(Z)]) + re = only(values(real(O3).terms)) + @test re[1] ≈ real((2.0 + 3.0im) * Z) + @test re[2] ≈ real(Z) + end +end diff --git a/test/toolbox/mpo_terms.jl b/test/toolbox/mpo_terms.jl new file mode 100644 index 000000000..5cedc1f2b --- /dev/null +++ b/test/toolbox/mpo_terms.jl @@ -0,0 +1,15 @@ +using Test +using PEPSKit + +@isdefined(TestSuite) || include("../testsuite/TestSuite.jl") +using .TestSuite + +is_buildkite = get(ENV, "BUILDKITE", "false") == "true" + +if !is_buildkite + TestSuite.toolbox_mpo_convention(Vector) + TestSuite.toolbox_mpo_dispatch(Vector) + TestSuite.toolbox_mpo_terms_dense(Vector) + TestSuite.toolbox_mpo_heisenberg(Vector) + TestSuite.toolbox_mpo_bookkeeping(Vector) +end diff --git a/test/toolbox/tensorproduct_terms.jl b/test/toolbox/tensorproduct_terms.jl new file mode 100644 index 000000000..a9b5a8d9f --- /dev/null +++ b/test/toolbox/tensorproduct_terms.jl @@ -0,0 +1,12 @@ +using Test +using PEPSKit + +@isdefined(TestSuite) || include("../testsuite/TestSuite.jl") +using .TestSuite + +is_buildkite = get(ENV, "BUILDKITE", "false") == "true" + +if !is_buildkite + TestSuite.toolbox_tensorproduct_ising(Vector) + TestSuite.toolbox_tensorproduct_bookkeeping(Vector) +end From 9a7b0b578c08112668bb7252d213c1fd73cc7bd9 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Sat, 3 Oct 2026 19:49:55 +0800 Subject: [PATCH 2/2] Fix MPO term validation and arithmetic Validate physical and bond spaces, preserve input containers, and support scalar promotion without losing tensor-product dispatch. Reject unsupported term accumulation and real/imaginary operations, and consolidate regression tests. --- .../contractions/local_patch/network_expr.jl | 2 +- src/operators/localoperator.jl | 123 ++++++---- test/testsuite/TestSuite.jl | 3 +- test/testsuite/toolbox/mpo_terms.jl | 217 +++++++----------- test/testsuite/toolbox/tensorproduct_terms.jl | 100 +++----- test/toolbox/mpo_terms.jl | 4 +- 6 files changed, 193 insertions(+), 256 deletions(-) diff --git a/src/algorithms/contractions/local_patch/network_expr.jl b/src/algorithms/contractions/local_patch/network_expr.jl index 7e3e39c32..5f77d08bc 100644 --- a/src/algorithms/contractions/local_patch/network_expr.jl +++ b/src/algorithms/contractions/local_patch/network_expr.jl @@ -267,7 +267,7 @@ end mpolabel(args...) = physicallabel(:mpo, args...) function operator_contraction_expr(::Type{<:MPOTerm}, nsites) - # one factor per site, linked by a chain of bonds running left to right, rank-2 at the + # one factor per site, linked by a chain of bonds running left to right, rank-3 at the # ends of the chain and rank-4 in the bulk: # W₁ : bra₁ ← ket₁ ⊗ b₁ # Wᵢ : bᵢ₋₁ ⊗ braᵢ ← ketᵢ ⊗ bᵢ diff --git a/src/operators/localoperator.jl b/src/operators/localoperator.jl index 7f24c47b6..e4ac38180 100644 --- a/src/operators/localoperator.jl +++ b/src/operators/localoperator.jl @@ -50,10 +50,8 @@ end # Default to Any for eltype: needs to be abstract anyways so not that much to gain LocalOperator(lattice, terms) = LocalOperator{Any}(lattice, terms) LocalOperator(lattice, terms::Pair...) = LocalOperator(lattice, terms) -# TODO: add terms beyond AbstractTensorMap -# e.g. tensor product of 1-site operators, MPOs -add_term!(operator::LocalOperator, inds::Tuple, term::AbstractTensorMap) = add_term!(operator, collect(inds), term) -add_term!(operator::LocalOperator, inds::Vector, term::AbstractTensorMap) = add_term!(operator, map(CartesianIndex{2}, inds), term) +add_term!(operator::LocalOperator, inds::Tuple, term::AbstractTensorMap; kwargs...) = add_term!(operator, collect(inds), term; kwargs...) +add_term!(operator::LocalOperator, inds::AbstractVector, term::AbstractTensorMap; kwargs...) = add_term!(operator, CartesianIndex{2}[CartesianIndex{2}(ind) for ind in inds], term; kwargs...) function add_term!( operator::LocalOperator, inds::Vector{CartesianIndex{2}}, term::AbstractTensorMap; atol = zero(real(scalartype(term))), @@ -76,9 +74,11 @@ function add_term!( end # translate coordinates - _shift_into_unitcell!(inds, size(operator)) + inds = _shift_into_unitcell!(copy(inds), size(operator)) if haskey(operator.terms, inds) + operator.terms[inds] isa MPOTerm && + throw(ArgumentError("Accumulating terms with the same sites is not implemented for MPO or tensor product terms.")) operator.terms[inds] = VI.add!!(operator.terms[inds], term) else operator.terms[inds] = term @@ -106,17 +106,18 @@ The factors are rank-2, carrying no bond indices, which is what distinguishes th more general [`MPOTerm`](@ref): both are vectors of tensors, but a tensor product term is the narrower type and therefore wins on dispatch wherever it applies. """ -const TensorProductTerm{T} = AbstractVector{T} where {T <: AbstractTensorMap{<:Any, <:Any, 1, 1}} +const TensorProductTerm{T} = AbstractVector{T} where {T <: AbstractTensorMap{<:Number, <:IndexSpace, 1, 1}} -add_term!(operator::LocalOperator, inds::Tuple, term::TensorProductTerm) = - add_term!(operator, collect(inds), term) -add_term!(operator::LocalOperator, inds::Vector, term::TensorProductTerm) = - add_term!(operator, map(CartesianIndex{2}, inds), term) +add_term!(operator::LocalOperator, inds::Tuple, term::TensorProductTerm; kwargs...) = + add_term!(operator, collect(inds), term; kwargs...) +add_term!(operator::LocalOperator, inds::AbstractVector, term::TensorProductTerm; kwargs...) = + add_term!(operator, CartesianIndex{2}[CartesianIndex{2}(ind) for ind in inds], term; kwargs...) function add_term!( operator::LocalOperator, inds::Vector{CartesianIndex{2}}, term::TensorProductTerm; - atol = zero(real(scalartype(first(term)))), + atol = 0, ) # input checks + isempty(term) && throw(ArgumentError("A tensor product term should contain at least one tensor.")) length(inds) == length(term) || throw(ArgumentError("Incompatible number of indices and tensor product factors")) allunique(inds) || throw(ArgumentError("`inds` should not contain repeated coordinates.")) @@ -137,7 +138,7 @@ function add_term!( end # translate coordinates - _shift_into_unitcell!(inds, size(operator)) + inds = _shift_into_unitcell!(copy(inds), size(operator)) # a sum of tensor products is not a tensor product, so terms cannot be accumulated here haskey(operator.terms, inds) && throw( @@ -163,8 +164,7 @@ end A single term of a [`LocalOperator`](@ref) given as a matrix product operator with one tensor per site acted on, rather than as the dense rank-`2N` tensor. -The factors are ordered along the chain and follow the usual MPO convention, rank-2 at the -ends and rank-4 in the bulk: +The factors are ordered along the chain and follow the usual MPO convention, rank-3 at the ends and rank-4 in the bulk: W₁ : P₁ ← P₁ ⊗ B₁ Wᵢ : Bᵢ₋₁ ⊗ Pᵢ ← Pᵢ ⊗ Bᵢ @@ -183,14 +183,11 @@ MPO. const MPOTerm{T} = AbstractVector{T} where {T <: AbstractTensorMap} -add_term!(operator::LocalOperator, inds::Tuple, term::MPOTerm) = - add_term!(operator, collect(inds), term) -add_term!(operator::LocalOperator, inds::Vector, term::MPOTerm) = - add_term!(operator, map(CartesianIndex{2}, inds), term) -function add_term!( - operator::LocalOperator, inds::Vector{CartesianIndex{2}}, term::MPOTerm - ) - # input checks +""" +Validate an ordered MPO's sites, tensor partitions, physical spaces, and adjacent bond spaces. +""" +function _validate_mpo_term(inds, term::MPOTerm, lattice::AbstractMatrix{<:ElementarySpace}) + isempty(term) && throw(ArgumentError("An MPO term should contain at least one tensor.")) length(inds) == length(term) || throw(ArgumentError("Incompatible number of indices and MPO factors")) allunique(inds) || throw(ArgumentError("`inds` should not contain repeated coordinates.")) @@ -209,16 +206,32 @@ function add_term!( # the bra index is the physical one in the codomain: last for a bulk factor, which # carries the incoming bond first, and only for an end factor bra = codomain(term[i])[nout] - ind_translated = CartesianIndex(mod1.(Tuple(ind), size(operator))) - physicalspace(operator, ind_translated) == bra || + ind_translated = CartesianIndex(mod1.(Tuple(ind), size(lattice))) + lattice[ind_translated] == bra == domain(term[i])[1] || throw(SpaceMismatch("Incompatible physical spaces")) end + for i in 1:(n - 1) + domain(term[i])[numin(term[i])] == codomain(term[i + 1])[1] || + throw(SpaceMismatch("Incompatible MPO bond spaces between tensors $i and $(i + 1).")) + end + return nothing +end + +add_term!(operator::LocalOperator, inds::Tuple, term::MPOTerm) = + add_term!(operator, collect(inds), term) +add_term!(operator::LocalOperator, inds::AbstractVector, term::MPOTerm) = + add_term!(operator, CartesianIndex{2}[CartesianIndex{2}(ind) for ind in inds], term) +function add_term!( + operator::LocalOperator, inds::Vector{CartesianIndex{2}}, term::MPOTerm + ) + _validate_mpo_term(inds, term, physicalspace(operator)) + _local_term_iszero(term) && return operator # NOTE: `inds` is deliberately *not* sorted here, unlike for dense and tensor product # terms, since permuting MPOs is not straightforward. # translate coordinates - _shift_into_unitcell!(inds, size(operator)) + inds = _shift_into_unitcell!(copy(inds), size(operator)) # as for tensor products, a sum of MPOs of fixed bond dimension is not one of the same # bond dimension, so terms are not accumulated here @@ -233,6 +246,12 @@ function add_term!( return operator end +""" +Detect an identically zero dense tensor or an MPO containing a zero factor. +""" +_local_term_iszero(term::AbstractTensorMap) = iszero(norm(term)) +_local_term_iszero(term::MPOTerm) = any(_local_term_iszero, term) + """ checklattice(Bool, args...) @@ -299,43 +318,33 @@ Base.eltype(::Type{LocalOperator{O, S}}) where {O, S} = O # Real and imaginary part # ----------------------- -""" -$(SIGNATURES) - -Take the real part of a single term of a [`LocalOperator`](@ref). - -A [`TensorProductTerm`](@ref) is made real factor by factor, so that the real part of a -tensor product is the tensor product of the real parts of its factors. An [`MPOTerm`](@ref) -is treated the same way, factor by factor. -""" -_real_local_term(term) = real(term) -_real_local_term(term::MPOTerm) = map(real, term) - function Base.real(O::LocalOperator) - return LocalOperator( - O.lattice, (sites => _real_local_term(op) for (sites, op) in O.terms)... - ) + any(term -> term isa MPOTerm, values(O.terms)) && + throw(ArgumentError("Taking the real part is not implemented for MPO or tensor product terms.")) + return LocalOperator(O.lattice, (sites => real(op) for (sites, op) in O.terms)...) end function Base.imag(O::LocalOperator) + any(term -> term isa MPOTerm, values(O.terms)) && + throw(ArgumentError("Taking the imaginary part is not implemented for MPO or tensor product terms.")) return LocalOperator(O.lattice, (sites => imag(op) for (sites, op) in O.terms)...) end # Linear Algebra # -------------- """ -$(SIGNATURES) - Scale a single term of a [`LocalOperator`](@ref) by `α`. -Terms which are not plain tensors need not scale by plain multiplication: an -[`MPOTerm`](@ref) - and hence also a [`TensorProductTerm`](@ref) - is scaled by scaling one -of its factors, since multiplying every factor would scale the term it represents by `α^n`. +An [`MPOTerm`](@ref), including a [`TensorProductTerm`](@ref), is scaled by scaling its first factor only. +The factor container may widen its element type when the scalar type changes, while retaining tensor-product dispatch. """ _scale_local_term(term, α::Number) = α * term function _scale_local_term(term::MPOTerm, α::Number) - scaled = collect(term) - scaled[1] = α * scaled[1] - return scaled + return AbstractTensorMap[i == 1 ? α * tensor : tensor for (i, tensor) in enumerate(term)] +end +function _scale_local_term(term::TensorProductTerm, α::Number) + return AbstractTensorMap{<:Number, <:IndexSpace, 1, 1}[ + i == 1 ? α * tensor : tensor for (i, tensor) in enumerate(term) + ] end Base.:*(α::Number, O::LocalOperator) = LocalOperator( @@ -348,7 +357,16 @@ Base.:\(α::Number, O::LocalOperator) = inv(α) * O function Base.:+(O1::LocalOperator, O2::LocalOperator) checklattice(O1, O2) - return LocalOperator(physicalspace(O1), mergewith(VI.add, O1.terms, O2.terms)) + return LocalOperator(physicalspace(O1), mergewith(_add_local_terms, O1.terms, O2.terms)) +end + +""" +Accumulate dense terms while rejecting addition involving an MPO or tensor product at the same sites. +""" +function _add_local_terms(term1, term2) + (term1 isa MPOTerm || term2 isa MPOTerm) && + throw(ArgumentError("Accumulating terms with the same sites is not implemented for MPO or tensor product terms.")) + return VI.add(term1, term2) end Base.:-(O::LocalOperator) = -1 * O @@ -359,9 +377,14 @@ Base.:-(O1::LocalOperator, O2::LocalOperator) = O1 + (-O2) # Since we allow abstract types in T, value and type domain might not match function VI.scalartype(operator::LocalOperator) - return promote_type((scalartype(term[2]) for term in operator.terms)...) + return promote_type((_local_term_scalartype(term) for term in values(operator.terms))...) end +""" +Return the promoted scalar type of a dense tensor or the factors of an MPO. +""" +_local_term_scalartype(term::AbstractTensorMap) = scalartype(term) +_local_term_scalartype(term::MPOTerm) = promote_type((scalartype(tensor) for tensor in term)...) # Equivalence # ----------- diff --git a/test/testsuite/TestSuite.jl b/test/testsuite/TestSuite.jl index c86d7c290..ccd1e9fa0 100644 --- a/test/testsuite/TestSuite.jl +++ b/test/testsuite/TestSuite.jl @@ -247,8 +247,7 @@ using .ToolboxTensorProductTerms module ToolboxMPOTerms include("toolbox/mpo_terms.jl") - export toolbox_mpo_convention, toolbox_mpo_dispatch, toolbox_mpo_terms_dense - export toolbox_mpo_heisenberg, toolbox_mpo_bookkeeping + export toolbox_mpo_terms_dense, toolbox_mpo_bookkeeping, toolbox_mpo_validation end using .ToolboxMPOTerms diff --git a/test/testsuite/toolbox/mpo_terms.jl b/test/testsuite/toolbox/mpo_terms.jl index 26e821782..58697db42 100644 --- a/test/testsuite/toolbox/mpo_terms.jl +++ b/test/testsuite/toolbox/mpo_terms.jl @@ -2,19 +2,13 @@ using Test using Random using TensorKit using PEPSKit, Adapt -using PEPSKit: MPOTerm, TensorProductTerm, gate_to_mpo, add_term! +using PEPSKit: gate_to_mpo, add_term! import TensorKitTensors.SpinOperators as SO -# `LocalOperator` terms given as a matrix product operator, one tensor per site acted on, -# rather than as the dense rank-2N tensor. `gate_to_mpo` turns a dense operator into one. -# -# A tensor product term is the trivial-bond case, and is the narrower type, so both are -# vectors of tensors and dispatch separates them by whether the factors are rank-2. - const Dbond, χenv = 2, 8 -const X = SO.σˣ(ComplexF64, Trivial) -const Z = SO.σᶻ(ComplexF64, Trivial) +const X = SO.σˣ(Float64, Trivial) +const Z = SO.σᶻ(Float64, Trivial) const lattice11 = fill(ComplexSpace(2), 1, 1) const dense_terms = ( @@ -22,147 +16,98 @@ const dense_terms = ( ("NN vertical", [(1, 1), (2, 1)], -(X ⊗ X) + Z ⊗ Z), ("NNN diagonal", [(1, 1), (2, 2)], X ⊗ Z), ("3 site line", [(1, 1), (1, 2), (1, 3)], X ⊗ X ⊗ X + Z ⊗ Z ⊗ Z), - # sites are taken in the order given, not sorted, so an unsorted or reversed order - # has to agree with the dense term written the same way round. `X ⊗ Z` is asymmetric - # under exchange, so this would catch a silently transposed pairing. + # Asymmetric operators expose accidental reordering of sites or factors. ("NNN anti-diagonal", [(1, 2), (2, 1)], X ⊗ Z), ("reversed sites", [(1, 2), (1, 1)], X ⊗ Z), ) -""" -Largest bond dimension across the cuts of an MPO term, i.e. the Schmidt rank of the operator -it represents. The outgoing bond of a factor is its last domain index, except for a rank-2 -factor, which carries no bond at all. -""" -function mpobond(term) - length(term) == 1 && return 1 - return maximum(1:(length(term) - 1)) do i - W = term[i] - numind(W) == 2 && return 1 - return dim(space(W, numind(W))) - end -end - -function toolbox_mpo_convention(AT) - return @testset "MPO factors follow the standard convention ($AT)" begin - # rank-2 at the ends of the chain, rank-4 in the bulk - W1, W2 = gate_to_mpo(X ⊗ X + Z ⊗ Z) - @test (numout(W1), numin(W1)) == (1, 2) - @test (numout(W2), numin(W2)) == (2, 1) - - Ws = gate_to_mpo(X ⊗ X ⊗ X + Z ⊗ Z ⊗ Z) - @test length(Ws) == 3 - @test (numout(Ws[1]), numin(Ws[1])) == (1, 2) - @test (numout(Ws[2]), numin(Ws[2])) == (2, 2) - @test (numout(Ws[3]), numin(Ws[3])) == (2, 1) - - # the bond is the Schmidt rank across each cut - @test mpobond(gate_to_mpo(X ⊗ X)) == 1 # a pure product - @test mpobond(gate_to_mpo(X ⊗ X + Z ⊗ Z)) == 2 - @test mpobond(gate_to_mpo(X ⊗ X ⊗ X + Z ⊗ Z ⊗ Z)) == 2 - @test mpobond([X]) == 1 # a lone site has no bond - - # Heisenberg XYZ has Schmidt rank 3 at every spin - for spin in (1 // 2, 1, 3 // 2) - term = SO.S_x_S_x(ComplexF64, Trivial; spin) + - SO.S_y_S_y(ComplexF64, Trivial; spin) + - SO.S_z_S_z(ComplexF64, Trivial; spin) - @test mpobond(gate_to_mpo(term)) == 3 - end - end -end - -function toolbox_mpo_dispatch(AT) - return @testset "Dispatch separates products from MPOs ($AT)" begin - # a vector of rank-2 operators is both, and the narrower type wins where it matters - @test [X, Z] isa TensorProductTerm - @test [X, Z] isa MPOTerm - - # anything carrying a bond is only an MPO term - for O in (X ⊗ Z, X ⊗ X + Z ⊗ Z, X ⊗ X ⊗ X + Z ⊗ Z ⊗ Z) - Ws = gate_to_mpo(O) - @test Ws isa MPOTerm - @test !(Ws isa TensorProductTerm) - end - end -end - function toolbox_mpo_terms_dense(AT) - return @testset "MPO terms reproduce dense terms ($(name)) ($AT)" for (name, inds, O) in - dense_terms - Random.seed!(2985721) + Random.seed!(2985721) + return @testset "MPO contractions ($AT)" begin peps = adapt(AT, InfinitePEPS(ComplexSpace(2), ComplexSpace(Dbond))) env = CTMRGEnv(randn, ComplexF64, peps, ComplexSpace(χenv)) + @testset "$name" for (name, inds, O) in dense_terms + O = adapt(AT, O) + H_dense = LocalOperator(lattice11, inds => O) + H_mpo = LocalOperator(lattice11, inds => gate_to_mpo(O)) + @test expectation_value(peps, H_mpo, env) ≈ expectation_value(peps, H_dense, env) rtol = 1.0e-9 + end - H_dense = LocalOperator(lattice11, inds => O) + # Complex scaling must promote real factors without scaling the coefficient twice. + _, inds, O = first(dense_terms) + O = adapt(AT, O) H_mpo = LocalOperator(lattice11, inds => gate_to_mpo(O)) - @test expectation_value(peps, H_mpo, env) ≈ expectation_value(peps, H_dense, env) rtol = - 1.0e-9 + α = 2 + 3im + @test expectation_value(peps, α * H_mpo, env) ≈ + α * expectation_value(peps, LocalOperator(lattice11, inds => O), env) rtol = 1.0e-9 end end -function toolbox_mpo_heisenberg(AT) - return @testset "MPO terms reproduce the Heisenberg XYZ Hamiltonian ($AT)" begin - Random.seed!(2985721) - H = heisenberg_XYZ(InfiniteSquare()) - H_mpo = LocalOperator( - physicalspace(H), (inds => gate_to_mpo(op) for (inds, op) in H.terms)... - ) - @test all(t -> t isa MPOTerm, values(H_mpo.terms)) - @test all(t -> mpobond(t) == 3, values(H_mpo.terms)) - - peps = adapt(AT, InfinitePEPS(ComplexSpace(2), ComplexSpace(Dbond))) - env, = leading_boundary( - CTMRGEnv(peps, ComplexSpace(χenv)), peps; tol = 1.0e-8, verbosity = 0 - ) - E = expectation_value(peps, H, env) - @test expectation_value(peps, H_mpo, env) ≈ E rtol = 1.0e-9 - - # scaling an MPO term scales the term, not each of its factors - α = 0.37 - @test expectation_value(peps, α * H_mpo, env) ≈ α * E rtol = 1.0e-9 - - # the real part is taken factor by factor - H_real = real(H_mpo) - @test all(t -> t isa MPOTerm, values(H_real.terms)) - @test all(all(W -> norm(imag(W)) < 1.0e-14, t) for t in values(H_real.terms)) +function toolbox_mpo_bookkeeping(AT) + Random.seed!(2985721) + return @testset "MPO bookkeeping ($AT)" begin + d = ℂ^2 + lattice = fill(d, 2, 2) + dense = adapt(AT, randn(Float64, d ⊗ d ← d ⊗ d)) + mpo = gate_to_mpo(dense; trunc = notrunc()) + sites = CartesianIndex.([(3, 4), (3, 3)]) + original_sites, original_mpo = copy(sites), copy(mpo) + operator = LocalOperator(lattice, sites => mpo) + @test sites == original_sites + reverse!(sites) + reverse!(mpo) + @test only(operator.terms) == (CartesianIndex.([(1, 2), (1, 1)]) => original_mpo) + + snapshot = deepcopy(operator) + scaled = (2 + 3im) * operator + @test operator == snapshot + onsite = adapt(AT, randn(ComplexF64, d ← d)) + @test scalartype(operator + LocalOperator(lattice, ((2, 1),) => onsite)) == ComplexF64 + @test scalartype(scaled) == ComplexF64 + @test_throws ArgumentError real(scaled) + @test_throws ArgumentError imag(scaled) + + # Both insertion and addition must reject accumulation involving an MPO. + inds = CartesianIndex.([(1, 1), (1, 2)]) + dense_operator = LocalOperator(lattice, inds => dense) + mpo_operator = LocalOperator(lattice, inds => original_mpo) + for (left, right) in ((mpo_operator, mpo_operator), (mpo_operator, dense_operator), (dense_operator, mpo_operator)) + @test_throws ArgumentError add_term!(deepcopy(left), inds, only(values(right.terms))) + @test_throws ArgumentError left + right + end + @test isempty(LocalOperator(lattice, inds => [zero(first(original_mpo)), last(original_mpo)]).terms) + + # Unequal physical spaces expose factor reordering during coordinate transformations. + lattice = [ℂ^2 ℂ^3 ℂ^4; ℂ^5 ℂ^6 ℂ^7] + dense = adapt(AT, randn(Float64, ℂ^6 ⊗ ℂ^2 ← ℂ^6 ⊗ ℂ^2)) + operator = LocalOperator(lattice, ((2, 2), (1, 1)) => gate_to_mpo(dense; trunc = notrunc())) + @test rotr90(rotl90(operator)) == operator + sites, term = only(operator.terms) + @test repeat(operator, 2, 1).terms == Dict(sites => term, (sites .+ CartesianIndex(2, 0)) => term) end end -function toolbox_mpo_bookkeeping(AT) - return @testset "MPO term bookkeeping ($AT)" begin - # sites are stored in the order given, and are *not* canonicalized by sorting the way - # dense and tensor product terms are: an MPO's bonds are the Schmidt cuts of one - # particular chain ordering, so its factors cannot be permuted after the fact. Factor i - # therefore acts on site inds[i], and the caller owns that pairing. - O_fwd = LocalOperator(lattice11, [(1, 1), (1, 2)] => gate_to_mpo(X ⊗ Z)) - O_rev = LocalOperator(lattice11, [(1, 2), (1, 1)] => gate_to_mpo(X ⊗ Z)) - # the sites are shifted into the unit cell as a block, anchored on the first of them, so - # absolute coordinates depend on which site comes first; what is preserved, and what - # encodes the chain order, is the displacement between consecutive sites - @test diff(only(keys(O_fwd.terms))) == [CartesianIndex(0, 1)] - @test diff(only(keys(O_rev.terms))) == [CartesianIndex(0, -1)] - - # a mismatch between sites and factors is caught - @test_throws ArgumentError LocalOperator( - lattice11, [(1, 1), (1, 2), (1, 3)] => gate_to_mpo(X ⊗ Z) - ) - - # factor ranks must match their position in the chain - @test_throws ArgumentError LocalOperator(lattice11, [(1, 1)] => [X ⊗ X]) - @test_throws ArgumentError LocalOperator(lattice11, [(1, 1)] => gate_to_mpo(X ⊗ Z)) - - # an MPO is not summed in place, since the bond dimension of a sum is not the same - O = LocalOperator(lattice11, [(1, 1), (1, 2)] => gate_to_mpo(X ⊗ Z)) - @test_throws ArgumentError add_term!( - O, [CartesianIndex(1, 1), CartesianIndex(1, 2)], gate_to_mpo(Z ⊗ X) - ) - - # physical spaces are checked - @test_throws SpaceMismatch LocalOperator( - lattice11, [(1, 1), (1, 2)] => gate_to_mpo( - SO.S_x_S_x(ComplexF64, Trivial; spin = 1) - ) - ) +function toolbox_mpo_validation(AT) + Random.seed!(2985721) + return @testset "MPO validation ($AT)" begin + d, b1, b2 = ℂ^2, ℂ^3, ℂ^4 + lattice = fill(d, 2, 2) + first_tensor = adapt(AT, randn(Float64, d ← d ⊗ b1)) + last_tensor = adapt(AT, randn(Float64, b1 ⊗ d ← d)) + mpo = [first_tensor, last_tensor] + sites = ((1, 1), (1, 2)) + + @test_throws ArgumentError LocalOperator(lattice, CartesianIndex{2}[] => AbstractTensorMap[]) + @test_throws ArgumentError LocalOperator(lattice, ((1, 1),) => mpo) + @test_throws ArgumentError LocalOperator(lattice, ((1, 1), (1, 1)) => mpo) + @test_throws ArgumentError LocalOperator(lattice, sites => reverse(mpo)) + + # Bond matching and the two physical legs are independent constraints. + wrong_bond = adapt(AT, randn(Float64, b2 ⊗ d ← d)) + wrong_input = adapt(AT, randn(Float64, b1 ⊗ d ← ℂ^3)) + wrong_output = adapt(AT, randn(Float64, b1 ⊗ ℂ^3 ← d)) + @test_throws SpaceMismatch LocalOperator(lattice, sites => [first_tensor, wrong_bond]) + @test_throws SpaceMismatch LocalOperator(lattice, sites => [first_tensor, wrong_input]) + @test_throws SpaceMismatch LocalOperator(lattice, sites => [first_tensor, wrong_output]) end end diff --git a/test/testsuite/toolbox/tensorproduct_terms.jl b/test/testsuite/toolbox/tensorproduct_terms.jl index b651fc349..58d7e3616 100644 --- a/test/testsuite/toolbox/tensorproduct_terms.jl +++ b/test/testsuite/toolbox/tensorproduct_terms.jl @@ -2,98 +2,70 @@ using Test using Random using TensorKit using PEPSKit, Adapt -using PEPSKit: TensorProductTerm, add_term!, _scale_local_term +using PEPSKit: add_term! import TensorKitTensors.SpinOperators as SO -# `LocalOperator` terms given as an explicit tensor product of single-site operators, one -# per site acted on, rather than as the tensor product formed up front. -# -# The transverse-field Ising model is the natural test case: with trivial symmetry every set -# of sites carries exactly one term, so the whole Hamiltonian is a sum of tensor products, -# a two-factor product per nearest neighbor bond and a one-factor product per site. This -# requires trivial symmetry, since a single-site σᶻ is not symmetric under the Z₂ symmetry -# of the model. - const Dbond, χenv, g = 2, 8, 3.1 const unitcells = ((1, 1), (2, 2)) """ -Transverse-field Ising Hamiltonian built from tensor product terms, matching the conventions -of [`transverse_field_ising`](@ref). +Build the transverse-field Ising Hamiltonian from real single-site factors. """ -function transverse_field_ising_tensorproduct( - T::Type{<:Number}, lattice::InfiniteSquare; J = 1.0, g = 1.0 - ) - Z, X = SO.σᶻ(T, Trivial), SO.σˣ(T, Trivial) +function transverse_field_ising_tensorproduct(AT, lattice::InfiniteSquare; g = 1.0) + Z, X = adapt.(Ref(AT), (SO.σᶻ(Float64, Trivial), SO.σˣ(Float64, Trivial))) spaces = fill(domain(X)[1], (lattice.Nrows, lattice.Ncols)) return LocalOperator( spaces, - (neighbor => [-J * Z, copy(Z)] for neighbor in nearest_neighbours(lattice))..., - ([idx] => [(-J * g) * X] for idx in vertices(lattice))..., + (neighbor => [-Z, Z] for neighbor in nearest_neighbours(lattice))..., + ([idx] => [-g * X] for idx in vertices(lattice))..., ) end -transverse_field_ising_tensorproduct(lattice::InfiniteSquare; kwargs...) = - transverse_field_ising_tensorproduct(ComplexF64, lattice; kwargs...) function toolbox_tensorproduct_ising(AT) - return @testset "Tensor product terms reproduce the two-site Ising Hamiltonian ($uc) ($AT)" for uc in - unitcells + return @testset "Tensor product Ising Hamiltonian ($uc) ($AT)" for uc in unitcells Random.seed!(2985721) H = transverse_field_ising(InfiniteSquare(uc...); g) - H_prod = transverse_field_ising_tensorproduct(InfiniteSquare(uc...); g) - - # one term per site set: a two-factor product per bond, a one-factor product per site - @test length(H_prod.terms) == length(H.terms) - @test all(t -> t isa TensorProductTerm, values(H_prod.terms)) - @test sort(length.(collect(values(H_prod.terms)))) == - sort(vcat(fill(1, prod(uc)), fill(2, 2 * prod(uc)))) - + H_prod = transverse_field_ising_tensorproduct(AT, InfiniteSquare(uc...); g) peps = adapt(AT, InfinitePEPS(ComplexSpace(2), ComplexSpace(Dbond); unitcell = uc)) env, = leading_boundary( CTMRGEnv(peps, ComplexSpace(χenv)), peps; tol = 1.0e-8, verbosity = 0 ) - - @test expectation_value(peps, H_prod, env) ≈ expectation_value(peps, H, env) rtol = - 1.0e-9 - - # scaling a product term scales the term, not each of its factors - α = 0.37 - @test expectation_value(peps, α * H_prod, env) ≈ - α * expectation_value(peps, H_prod, env) rtol = 1.0e-9 + E = expectation_value(peps, H, env) + @test expectation_value(peps, H_prod, env) ≈ E rtol = 1.0e-9 + # Real factors must retain tensor-product dispatch after complex scaling. + α = 2 + 3im + @test expectation_value(peps, α * H_prod, env) ≈ α * E rtol = 1.0e-9 end end function toolbox_tensorproduct_bookkeeping(AT) - return @testset "Tensor product term bookkeeping ($AT)" begin + return @testset "Tensor product bookkeeping ($AT)" begin lattice = fill(ComplexSpace(2), 1, 1) - X, Z = SO.σˣ(ComplexF64, Trivial), SO.σᶻ(ComplexF64, Trivial) - - # a sum of tensor products is not a tensor product, so terms may not be accumulated - O = LocalOperator(lattice, [(1, 1), (1, 2)] => [X, X]) - @test_throws ArgumentError add_term!(O, [(1, 1), (1, 2)], [Z, Z]) + X, Z = adapt.(Ref(AT), (SO.σˣ(Float64, Trivial), SO.σᶻ(Float64, Trivial))) + sites = CartesianIndex.([(3, 4), (3, 5)]) + factors = [X, Z] + product = LocalOperator(lattice, sites => factors) + @test sites == CartesianIndex.([(3, 4), (3, 5)]) + reverse!(sites) + reverse!(factors) + @test only(product.terms) == (CartesianIndex.([(1, 1), (1, 2)]) => [X, Z]) + @test product == LocalOperator(lattice, [(1, 2), (1, 1)] => [Z, X]) - # factors are reordered along with their sites - O1 = LocalOperator(lattice, [(1, 1), (1, 2)] => [X, Z]) - O2 = LocalOperator(lattice, [(1, 2), (1, 1)] => [Z, X]) - @test only(keys(O1.terms)) == only(keys(O2.terms)) - @test all(map(≈, only(values(O1.terms)), only(values(O2.terms)))) + snapshot = deepcopy(product) + scaled = (2 + 3im) * product + @test product == snapshot + @test scalartype(scaled) == ComplexF64 + @test_throws ArgumentError add_term!(product, [(1, 1), (1, 2)], [Z, X]) + @test_throws ArgumentError product + product - # each factor must be a single-site operator, and match the physical space - @test_throws ArgumentError LocalOperator(lattice, [(1, 1), (1, 2)] => [X]) # arity - @test_throws ArgumentError LocalOperator(lattice, [(1, 1)] => [X ⊗ X]) # not single-site + # Factorwise real/imaginary parts are incorrect even for tensor products. + imaginary = LocalOperator(lattice, [(1, 1), (1, 2)] => [im * Z, im * Z]) + @test_throws ArgumentError real(imaginary) + @test_throws ArgumentError imag(imaginary) + @test_throws ArgumentError LocalOperator(lattice, CartesianIndex{2}[] => typeof(X)[]) + @test_throws ArgumentError LocalOperator(lattice, [(1, 1), (1, 2)] => [X]) @test_throws SpaceMismatch LocalOperator( - lattice, [(1, 1)] => [SO.S_x(ComplexF64, Trivial; spin = 1)] + lattice, [(1, 1)] => [SO.S_x(Float64, Trivial; spin = 1)] ) - - # scaling touches exactly one factor - scaled = _scale_local_term([X, Z], 3.0) - @test scaled[1] ≈ 3.0 * X - @test scaled[2] ≈ Z - - # the real part of a product is the product of the real parts of its factors - O3 = LocalOperator(lattice, [(1, 1), (1, 2)] => [(2.0 + 3.0im) * Z, copy(Z)]) - re = only(values(real(O3).terms)) - @test re[1] ≈ real((2.0 + 3.0im) * Z) - @test re[2] ≈ real(Z) end end diff --git a/test/toolbox/mpo_terms.jl b/test/toolbox/mpo_terms.jl index 5cedc1f2b..e29ac44ea 100644 --- a/test/toolbox/mpo_terms.jl +++ b/test/toolbox/mpo_terms.jl @@ -7,9 +7,7 @@ using .TestSuite is_buildkite = get(ENV, "BUILDKITE", "false") == "true" if !is_buildkite - TestSuite.toolbox_mpo_convention(Vector) - TestSuite.toolbox_mpo_dispatch(Vector) TestSuite.toolbox_mpo_terms_dense(Vector) - TestSuite.toolbox_mpo_heisenberg(Vector) TestSuite.toolbox_mpo_bookkeeping(Vector) + TestSuite.toolbox_mpo_validation(Vector) end