Add docstrings and example for boundary MPS functionality

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])
 ```

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 <peps|peps>
+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 <mps|T|mps>
+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
+# <peps|O|peps>.
+
+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

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

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,39 +33,35 @@ 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."))
 
     Sspaces = adjoint.(circshift(Nspaces, (1, 0, 0)))
     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

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})

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