diff --git a/docs/src/lib/lib.md b/docs/src/lib/lib.md index 592894e..b77bf0a 100644 --- a/docs/src/lib/lib.md +++ b/docs/src/lib/lib.md @@ -2,4 +2,5 @@ ```@autodocs Modules = [PEPSKit, PEPSKit.Defaults] +Filter = t -> !(t in [PEPSKit.PEPS_∂∂AC, PEPSKit.PEPS_∂∂C, PEPSKit.PEPO_∂∂AC, PEPSKit.PEPO_∂∂C]) ``` diff --git a/examples/boundary_mps.jl b/examples/boundary_mps.jl new file mode 100644 index 0000000..abd52e1 --- /dev/null +++ b/examples/boundary_mps.jl @@ -0,0 +1,86 @@ +using Random +using PEPSKit +using TensorKit +using MPSKit +using LinearAlgebra + +include("ising_pepo.jl") + +# This example demonstrates some boundary-MPS methods for working with 2D projected +# entangled-pair states and operators. + +## Computing a PEPS norm + +# We start by initializing a random initial infinite PEPS +Random.seed!(29384293742893) +peps = InfinitePEPS(ComplexSpace(2), ComplexSpace(2)) + +# To compute its norm, we start by constructing the transfer operator corresponding to +# the partition function representing the overlap +T = InfiniteTransferPEPS(peps, 1, 1) + +# We then find its leading boundary MPS fixed point, where the corresponding eigenvalue +# encodes the norm of the state + +# Fist we build an initial guess for the boundary MPS, choosing a bond dimension of 20 +mps = PEPSKit.initializeMPS(T, [ComplexSpace(20)]) + +# We then find the leading boundary MPS fixed point using the VUMPS algorithm +mps, envs, ϵ = leading_boundary(mps, T, VUMPS()) + +# The norm of the state per unit cell is then given by the expectation value +N = abs(prod(expectation_value(mps, T))) + +# This can be compared to the result obtained using the CTMRG algorithm +ctm = leading_boundary( + peps, CTMRG(; verbosity=1, fixedspace=true), CTMRGEnv(peps; Venv=ComplexSpace(20)) +) +N´ = abs(norm(peps, ctm)) + +@show abs(N - N´) / N + +## Working with unit cells + +# For PEPS with non-trivial unit cells, the principle is exactly the same. +# The only difference is that now the transfer operator of the PEPS norm partition function +# has multiple lines, each of which can be represented by an `InfiniteTransferPEPS` object. +# Such a multi-line transfer operator is represented by a `TransferPEPSMultiline` object. +# In this case, the boundary MPS is an `MPSMultiline` object, which should be initialized +# by specifying a virtual space for each site in the partition function unit cell. + +peps2 = InfinitePEPS(ComplexSpace(2), ComplexSpace(2); unitcell=(2, 2)) +T2 = PEPSKit.TransferPEPSMultiline(peps2, 1) + +mps2 = PEPSKit.initializeMPS(T2, fill(ComplexSpace(20), 2, 2)) +mps2, envs2, ϵ = leading_boundary(mps2, T2, VUMPS()) +N2 = abs(prod(expectation_value(mps2, T2))) + +ctm2 = leading_boundary( + peps2, CTMRG(; verbosity=1, fixedspace=true), CTMRGEnv(peps2; Venv=ComplexSpace(20)) +) +N2´ = abs(norm(peps2, ctm2)) + +@show abs(N2 - N2´) / N2 + +# Note that for larger unit cells and non-Hermitian PEPS the VUMPS algorithm may become +# unstable, in which case the CTMRG algorithm is recommended. + +## Contracting PEPO overlaps + +# Using exactly the same machinery, we can contract partition functions which encode the +# expectation value of a PEPO for a given PEPS state. +# As an example, we can consider the overlap of the PEPO correponding to the partition +# function of 3D classical ising model with our random PEPS from before and evaluate +# . + +pepo = ising_pepo(1) +T3 = InfiniteTransferPEPO(peps, pepo, 1, 1) + +mps3 = PEPSKit.initializeMPS(T3, [ComplexSpace(20)]) +mps3, envs3, ϵ = leading_boundary(mps3, T3, VUMPS()) +@show N3 = abs(prod(expectation_value(mps3, T3))) + +# These objects and routines can be used to optimize PEPS fixed points of 3D partition +# functions, see for example https://arxiv.org/abs/1805.10598 + +nothing diff --git a/examples/ising_pepo.jl b/examples/ising_pepo.jl new file mode 100644 index 0000000..4936e60 --- /dev/null +++ b/examples/ising_pepo.jl @@ -0,0 +1,23 @@ +using TensorOperations +using TensorKit + +""" + ising_pepo(beta; unitcell=(1, 1, 1)) + +Return the PEPO tensor for partition function of the 3D classical Ising model at inverse +temperature `beta`. +""" +function ising_pepo(beta; unitcell=(1, 1, 1)) + t = ComplexF64[exp(beta) exp(-beta); exp(-beta) exp(beta)] + q = sqrt(t) + + O = zeros(2, 2, 2, 2, 2, 2) + O[1, 1, 1, 1, 1, 1] = 1 + O[2, 2, 2, 2, 2, 2] = 1 + @tensor o[-1 -2; -3 -4 -5 -6] := + O[1 2; 3 4 5 6] * q[-1; 1] * q[-2; 2] * q[-3; 3] * q[-4; 4] * q[-5; 5] * q[-6; 6] + + O = TensorMap(o, ℂ^2 ⊗ (ℂ^2)' ← ℂ^2 ⊗ ℂ^2 ⊗ (ℂ^2)' ⊗ (ℂ^2)') + + return InfinitePEPO(O; unitcell) +end diff --git a/src/operators/infinitepepo.jl b/src/operators/infinitepepo.jl index afd90bb..5633cff 100644 --- a/src/operators/infinitepepo.jl +++ b/src/operators/infinitepepo.jl @@ -1,7 +1,7 @@ """ struct InfinitePEPO{T<:PEPOTensor} -Represents an infinte projected entangled-pair operator (PEPO) on a 3D cubic lattice. +Represents an infinite projected entangled-pair operator (PEPO) on a 3D cubic lattice. """ struct InfinitePEPO{T<:PEPOTensor} <: AbstractPEPO A::Array{T,3} @@ -24,7 +24,7 @@ end ## Constructors """ - InfinitePEPO(A::AbstractArray{T, 2}) + InfinitePEPO(A::AbstractArray{T, 3}) Allow users to pass in an array of tensors. """ @@ -33,15 +33,18 @@ function InfinitePEPO(A::AbstractArray{T,3}) where {T<:PEPOTensor} end """ - InfinitePEPO(Pspaces, Nspaces, Espaces) + InfinitePEPO(f=randn, T=ComplexF64, Pspaces, Nspaces, Espaces) Allow users to pass in arrays of spaces. """ function InfinitePEPO( - Pspaces::AbstractArray{S,3}, - Nspaces::AbstractArray{S,3}, - Espaces::AbstractArray{S,3}=Nspaces, -) where {S<:ElementarySpace} + Pspaces::A, Nspaces::A, Espaces::A=Nspaces +) where {A<:AbstractArray{<:ElementarySpace,3}} + return InfinitePEPO(randn, ComplexF64, Pspaces, Nspaces, Espaces) +end +function InfinitePEPO( + f, T, Pspaces::A, Nspaces::A, Espaces::A=Nspaces +) where {A<:AbstractArray{<:ElementarySpace,3}} size(Pspaces) == size(Nspaces) == size(Espaces) || throw(ArgumentError("Input spaces should have equal sizes.")) @@ -49,23 +52,16 @@ function InfinitePEPO( Wspaces = adjoint.(circshift(Espaces, (0, -1, 0))) Ppspaces = adjoint.(circshift(Pspaces, (0, 0, -1))) - A = map(Pspaces, Ppspaces, Nspaces, Espaces, Sspaces, Wspaces) do P, Pp, N, E, S, W - return TensorMap(rand, ComplexF64, P * Pp ← N * E * S * W) + P = map(Pspaces, Ppspaces, Nspaces, Espaces, Sspaces, Wspaces) do P, Pp, N, E, S, W + return TensorMap(f, T, P * Pp ← N * E * S * W) end - return InfinitePEPO(A) + return InfinitePEPO(P) end -""" - InfinitePEPO(Pspaces, Nspaces, Espaces) - -Allow users to pass in arrays of spaces, single layer special case. -""" function InfinitePEPO( - Pspaces::AbstractArray{S,2}, - Nspaces::AbstractArray{S,2}, - Espaces::AbstractArray{S,2}=Nspaces, -) where {S<:ElementarySpace} + Pspaces::A, Nspaces::A, Espaces::A=Nspaces +) where {A<:AbstractArray{<:ElementarySpace,2}} size(Pspaces) == size(Nspaces) == size(Espaces) || throw(ArgumentError("Input spaces should have equal sizes.")) @@ -77,53 +73,36 @@ function InfinitePEPO( end """ - InfinitePEPO(A) - -Allow users to pass in single tensor. -""" -function InfinitePEPO(A::T) where {T<:PEPOTensor} - As = Array{T,3}(undef, (1, 1, 1)) - As[1, 1, 1] = A - return InfinitePEPO(As) -end - -""" - InfinitePEPO(Pspace, Nspace, Espace) - -Allow users to pass in single space. -""" -function InfinitePEPO(Pspace::S, Nspace::S, Espace::S=Nspace) where {S<:ElementarySpace} - Pspaces = Array{S,3}(undef, (1, 1, 1)) - Pspaces[1, 1] = Pspace - Nspaces = Array{S,3}(undef, (1, 1, 1)) - Nspaces[1, 1] = Nspace - Espaces = Array{S,3}(undef, (1, 1, 1)) - Espaces[1, 1] = Espace - return InfinitePEPO(Pspaces, Nspaces, Espaces) -end - -""" - InfinitePEPO(d, D) + InfinitePEPO(A; unitcell=(1, 1, 1)) -Allow users to pass in integers. +Create an InfinitePEPO by specifying a tensor and unit cell. """ -function InfinitePEPO(d::Integer, D::Integer) - T = TensorMap(rand, ComplexF64, ℂ^d ⊗ (ℂ^d)' ← ℂ^D ⊗ ℂ^D ⊗ (ℂ^D)' ⊗ (ℂ^D)') - return InfinitePEPO(T) +function InfinitePEPO(A::T; unitcell::Tuple{Int,Int,Int}=(1, 1, 1)) where {T<:PEPOTensor} + return InfinitePEPO(fill(A, unitcell)) end """ - InfinitePEPO(d, D, L) - InfinitePEPO(d, D, (Lx, Ly, Lz))) + InfinitePEPO(f=randn, T=ComplexF64, Pspace, Nspace, [Espace]; unitcell=(1,1,1)) -Allow users to pass in integers and specify unit cell. +Create an InfinitePEPO by specifying its spaces and unit cell. """ -function InfinitePEPO(d::Integer, D::Integer, L::Integer) - return InfinitePEPO(d, D, (L, L, L)) +function InfinitePEPO( + Pspace::S, Nspace::S, Espace::S=Nspace; unitcell::Tuple{Int,Int,Int}=(1, 1, 1) +) where {S<:ElementarySpace} + return InfinitePEPO( + randn, + ComplexF64, + fill(Pspace, unitcell), + fill(Nspace, unitcell), + fill(Espace, unitcell), + ) end -function InfinitePEPO(d::Integer, D::Integer, Ls::NTuple{3,Integer}) - T = [TensorMap(rand, ComplexF64, ℂ^d ⊗ (ℂ^d)' ← ℂ^D ⊗ ℂ^D ⊗ (ℂ^D)' ⊗ (ℂ^D)')] - return InfinitePEPO(Array(repeat(T, Ls...))) +function InfinitePEPO( + f, T, Pspace::S, Nspace::S, Espace::S=Nspace; unitcell::Tuple{Int,Int,Int}=(1, 1, 1) +) where {S<:ElementarySpace} + return InfinitePEPO( + f, T, fill(Pspace, unitcell), fill(Nspace, unitcell), fill(Espace, unitcell) + ) end ## Shape and size diff --git a/src/operators/localoperator.jl b/src/operators/localoperator.jl index 53834a4..5f8e57b 100644 --- a/src/operators/localoperator.jl +++ b/src/operators/localoperator.jl @@ -75,8 +75,7 @@ end Evaluate the expectation value of any `NLocalOperator` on each unit-cell entry of `peps` and `env`. -""" -MPSKit.expectation_value +""" MPSKit.expectation_value(::InfinitePEPS, ::Any, ::NLocalOperator) # 1-site operator expectation values on unit cell function MPSKit.expectation_value(peps::InfinitePEPS, env, O::NLocalOperator{OnSite}) diff --git a/src/operators/transferpepo.jl b/src/operators/transferpepo.jl index a331479..eae23fd 100644 --- a/src/operators/transferpepo.jl +++ b/src/operators/transferpepo.jl @@ -1,3 +1,10 @@ +""" + InfiniteTransferPEPO{T,O} + +Represents an infinite transfer operator corresponding to a single row of a partition +function which corresponds to the expectation value of an `InfinitePEPO` between 'ket' and +'bra' `InfinitePEPS` states. +""" struct InfiniteTransferPEPO{T,O} top::PeriodicArray{T,1} mid::PeriodicArray{O,2} @@ -6,6 +13,14 @@ end InfiniteTransferPEPO(top, mid) = InfiniteTransferPEPO(top, mid, top) +""" + InfiniteTransferPEPO(T::InfinitePEPS, O::InfinitePEPO, dir, row) + +Constructs a transfer operator corresponding to a single row of a partition function +representing the expectation value of `O` for the state `T`. The partition function is first +rotated such that the direction `dir` faces north, after which its `row`th row from the +north is selected. +""" function InfiniteTransferPEPO(T::InfinitePEPS, O::InfinitePEPO, dir, row) T = rotate_north(T, dir) O = rotate_north(O, dir) @@ -51,12 +66,25 @@ function initializeMPS(O::InfiniteTransferPEPO, χ::Int) ]) end +""" + const TransferPEPOMultiline = MPSKit.Multiline{<:InfiniteTransferPEPO} + +Type that represents a multi-line transfer operator, where each line each corresponds to a +row of a partition function encoding the overlap of an `InfinitePEPO` between 'ket' and +'bra' `InfinitePEPS` states. +""" const TransferPEPOMultiline = MPSKit.Multiline{<:InfiniteTransferPEPO} Base.convert(::Type{TransferPEPOMultiline}, O::InfiniteTransferPEPO) = MPSKit.Multiline([O]) Base.getindex(t::TransferPEPOMultiline, i::Colon, j::Int) = Base.getindex.(t.data[i], j) Base.getindex(t::TransferPEPOMultiline, i::Int, j) = Base.getindex(t.data[i], j) -# multiline patch +""" + TransferPEPOMultiline(T::InfinitePEPS, O::InfinitePEPO, dir) + +Construct a multi-row transfer operator corresponding to the partition function representing +the expectation value of `O` for the state `T`. The partition function is first rotated such +that the direction `dir` faces north. +""" function TransferPEPOMultiline(T::InfinitePEPS, O::InfinitePEPO, dir) return MPSKit.Multiline(map(cr -> InfiniteTransferPEPO(T, O, dir, cr), 1:size(T, 1))) end diff --git a/src/operators/transferpeps.jl b/src/operators/transferpeps.jl index e7375ca..f159764 100644 --- a/src/operators/transferpeps.jl +++ b/src/operators/transferpeps.jl @@ -1,4 +1,9 @@ +""" + InfiniteTransferPEPS{T} +Represents an infinite transfer operator corresponding to a single row of a partition +function which corresponds to the overlap between 'ket' and 'bra' `InfinitePEPS` states. +""" struct InfiniteTransferPEPS{T} top::PeriodicArray{T,1} bot::PeriodicArray{T,1} @@ -6,6 +11,13 @@ end InfiniteTransferPEPS(top) = InfiniteTransferPEPS(top, top) +""" + InfiniteTransferPEPS(T::InfinitePEPS, dir, row) + +Constructs a transfer operator corresponding to a single row of a partition function +representing the norm of the state `T`. The partition function is first rotated such that +the direction `dir` faces north, after which its `row`th row from the north is selected. +""" function InfiniteTransferPEPS(T::InfinitePEPS, dir, row) T = rotate_north(T, dir) return InfiniteTransferPEPS(PeriodicArray(T[row, :])) @@ -18,6 +30,44 @@ Base.getindex(O::InfiniteTransferPEPS, i) = (O.top[i], O.bot[i]) Base.iterate(O::InfiniteTransferPEPS, i=1) = i > length(O) ? nothing : (O[i], i + 1) +import MPSKit.GenericMPSTensor + +""" + const TransferPEPSMultiline = MPSKit.Multiline{<:InfiniteTransferPEPS} + +Type that represents a multi-line transfer operator, where each line each corresponds to a +row of a partition function encoding the overlap between 'ket' and 'bra' `InfinitePEPS` +states. +""" +const TransferPEPSMultiline = MPSKit.Multiline{<:InfiniteTransferPEPS} +Base.convert(::Type{TransferPEPSMultiline}, O::InfiniteTransferPEPS) = MPSKit.Multiline([O]) +Base.getindex(t::TransferPEPSMultiline, i::Colon, j::Int) = Base.getindex.(t.data[i], j) +Base.getindex(t::TransferPEPSMultiline, i::Int, j) = Base.getindex(t.data[i], j) + +""" + TransferPEPSMultiline(T::InfinitePEPS, dir) + +Construct a multi-row transfer operator corresponding to the partition function representing +the norm of the state `T`. The partition function is first rotated such +that the direction `dir` faces north. +""" +function TransferPEPSMultiline(T::InfinitePEPS, dir) + return MPSKit.Multiline(map(cr -> InfiniteTransferPEPS(T, dir, cr), 1:size(T, 1))) +end + +""" + initializeMPS( + O::Union{InfiniteTransferPEPS,InfiniteTransferPEPO}, + virtualspaces::AbstractArray{<:ElementarySpace,1} + ) + initializeMPS( + O::Union{TransferPEPSMultiline,TransferPEPOMultiline}, + virtualspaces::AbstractArray{<:ElementarySpace,2} + ) + +Inialize a boundary MPS for the transfer operator `O` by specifying an array of virtual +spaces consistent with the unit cell. +""" function initializeMPS(O::InfiniteTransferPEPS, virtualspaces::AbstractArray{S,1}) where {S} return InfiniteMPS([ TensorMap( @@ -28,7 +78,6 @@ function initializeMPS(O::InfiniteTransferPEPS, virtualspaces::AbstractArray{S,1 ) for i in 1:length(O) ]) end - function initializeMPS(O::InfiniteTransferPEPS, χ::Int) return InfiniteMPS([ TensorMap( @@ -39,18 +88,6 @@ function initializeMPS(O::InfiniteTransferPEPS, χ::Int) ) for i in 1:length(O) ]) end - -import MPSKit.GenericMPSTensor - -const TransferPEPSMultiline = MPSKit.Multiline{<:InfiniteTransferPEPS} -Base.convert(::Type{TransferPEPSMultiline}, O::InfiniteTransferPEPS) = MPSKit.Multiline([O]) -Base.getindex(t::TransferPEPSMultiline, i::Colon, j::Int) = Base.getindex.(t.data[i], j) -Base.getindex(t::TransferPEPSMultiline, i::Int, j) = Base.getindex(t.data[i], j) - -# multiline patch -function TransferPEPSMultiline(T::InfinitePEPS, dir) - return MPSKit.Multiline(map(cr -> InfiniteTransferPEPS(T, dir, cr), 1:size(T, 1))) -end function initializeMPS(O::MPSKit.Multiline, virtualspaces::AbstractArray{S,2}) where {S} mpss = map(cr -> initializeMPS(O[cr], virtualspaces[cr, :]), 1:size(O, 1)) return MPSKit.Multiline(mpss) @@ -93,6 +130,14 @@ function MPSKit.transfer_right( A[-1 9 8 7] end +@doc """ + MPSKit.expectation_value(st::InfiniteMPS, op::Union{InfiniteTransferPEPS,InfiniteTransferPEPO}) + MPSKit.expectation_value(st::MPSMultiline, op::Union{TransferPEPSMultiline,TransferPEPOMultiline}) + +Compute expectation value of the transfer operator `op` for the state `st` for each site in +the unit cell. +""" MPSKit.expectation_value(st, op) + function MPSKit.expectation_value(st::InfiniteMPS, transfer::InfiniteTransferPEPS) return expectation_value( convert(MPSMultiline, st), convert(TransferPEPSMultiline, transfer) @@ -122,5 +167,14 @@ function MPSKit.expectation_value( return retval end -# PEPS derivatives -# ---------------- +@doc """ + MPSKit.leading_boundary( + st::InfiniteMPS, op::Union{InfiniteTransferPEPS,InfiniteTransferPEPO}, alg, [envs] + ) + MPSKit.leading_boundary( + st::MPSMulitline, op::Union{TransferPEPSMultiline,TransferPEPOMultiline}, alg, [envs] + ) + +Approximate the leading boundary MPS eigenvector for the transfer operator `op` using `st` +as initial guess. +""" MPSKit.leading_boundary(st, op, alg) diff --git a/test/boundarymps/vumps.jl b/test/boundarymps/vumps.jl index 9943cac..aa0b99c 100644 --- a/test/boundarymps/vumps.jl +++ b/test/boundarymps/vumps.jl @@ -6,17 +6,67 @@ using MPSKit using LinearAlgebra Random.seed!(29384293742893) -psi = InfinitePEPS(ComplexSpace(2), ComplexSpace(2)) -T = PEPSKit.InfiniteTransferPEPS(psi, 1, 1) -mps = PEPSKit.initializeMPS(T, [ComplexSpace(20)]) +@testset "(1, 1) PEPS" begin + psi = InfinitePEPS(ComplexSpace(2), ComplexSpace(2)) -mps, envs, ϵ = leading_boundary(mps, T, VUMPS()) -N = abs(sum(expectation_value(mps, T))) + T = PEPSKit.InfiniteTransferPEPS(psi, 1, 1) + mps = PEPSKit.initializeMPS(T, [ComplexSpace(20)]) -ctm = leading_boundary( - psi, CTMRG(; verbosity=1, fixedspace=true), CTMRGEnv(psi; Venv=ComplexSpace(20)) -) -N2 = abs(norm(psi, ctm)) + mps, envs, ϵ = leading_boundary(mps, T, VUMPS()) + N = abs(sum(expectation_value(mps, T))) -@test N ≈ N2 atol = 1e-3 + ctm = leading_boundary( + psi, CTMRG(; verbosity=1, fixedspace=true), CTMRGEnv(psi; Venv=ComplexSpace(20)) + ) + N´ = abs(norm(psi, ctm)) + + @test N ≈ N´ atol = 1e-3 +end + +@testset "(2, 2) PEPS" begin + psi = InfinitePEPS(ComplexSpace(2), ComplexSpace(2); unitcell=(2, 2)) + T = PEPSKit.TransferPEPSMultiline(psi, 1) + + mps = PEPSKit.initializeMPS(T, fill(ComplexSpace(20), 2, 2)) + mps, envs, ϵ = leading_boundary(mps, T, VUMPS()) + N = abs(prod(expectation_value(mps, T))) + + ctm = leading_boundary( + psi, CTMRG(; verbosity=1, fixedspace=true), CTMRGEnv(psi; Venv=ComplexSpace(20)) + ) + N´ = abs(norm(psi, ctm)) + + @test N ≈ N´ rtol = 1e-3 +end + +@testset "PEPO runthrough" begin + function ising_pepo(beta; unitcell=(1, 1, 1)) + t = ComplexF64[exp(beta) exp(-beta); exp(-beta) exp(beta)] + q = sqrt(t) + + O = zeros(2, 2, 2, 2, 2, 2) + O[1, 1, 1, 1, 1, 1] = 1 + O[2, 2, 2, 2, 2, 2] = 1 + @tensor o[-1 -2; -3 -4 -5 -6] := + O[1 2; 3 4 5 6] * + q[-1; 1] * + q[-2; 2] * + q[-3; 3] * + q[-4; 4] * + q[-5; 5] * + q[-6; 6] + + O = TensorMap(o, ℂ^2 ⊗ (ℂ^2)' ← ℂ^2 ⊗ ℂ^2 ⊗ (ℂ^2)' ⊗ (ℂ^2)') + + return InfinitePEPO(O; unitcell) + end + + psi = InfinitePEPS(ComplexSpace(2), ComplexSpace(2)) + O = ising_pepo(1) + T = InfiniteTransferPEPO(psi, O, 1, 1) + + mps = PEPSKit.initializeMPS(T, [ComplexSpace(10)]) + mps, envs, ϵ = leading_boundary(mps, T, VUMPS()) + f = abs(prod(expectation_value(mps, T))) +end