diff --git a/.github/dependabot.yml b/.github/dependabot.yml new file mode 100644 index 00000000..ff6499d6 --- /dev/null +++ b/.github/dependabot.yml @@ -0,0 +1,7 @@ +# https://docs.github.com/github/administering-a-repository/configuration-options-for-dependency-updates +version: 2 +updates: + - package-ecosystem: "github-actions" + directory: "/" # Location of package manifests + schedule: + interval: "weekly" \ No newline at end of file diff --git a/.github/workflows/ci.yml b/.github/workflows/CI.yml similarity index 85% rename from .github/workflows/ci.yml rename to .github/workflows/CI.yml index 17148087..4d0a1f65 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/CI.yml @@ -31,17 +31,20 @@ jobs: - x64 steps: - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v1 + - uses: julia-actions/setup-julia@v2 with: version: ${{ matrix.version }} arch: ${{ matrix.arch }} - - uses: julia-actions/cache@v1 + - uses: julia-actions/cache@v2 - uses: julia-actions/julia-buildpkg@latest - uses: julia-actions/julia-runtest@latest env: JULIA_NUM_THREADS: 4 - uses: julia-actions/julia-processcoverage@v1 - - uses: codecov/codecov-action@v3 + - uses: codecov/codecov-action@v4 + if: always() + env: + CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }} with: file: lcov.info test-nightly: @@ -61,11 +64,11 @@ jobs: - x64 steps: - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v1 + - uses: julia-actions/setup-julia@v2 with: version: ${{ matrix.version }} arch: ${{ matrix.arch }} - - uses: julia-actions/cache@v1 + - uses: julia-actions/cache@v2 - uses: julia-actions/julia-buildpkg@latest - uses: julia-actions/julia-runtest@latest env: diff --git a/.github/workflows/CompatHelper.yml b/.github/workflows/CompatHelper.yml new file mode 100644 index 00000000..06033911 --- /dev/null +++ b/.github/workflows/CompatHelper.yml @@ -0,0 +1,43 @@ +name: CompatHelper +on: + schedule: + - cron: 0 0 * * * + workflow_dispatch: +permissions: + contents: write + pull-requests: write +jobs: + CompatHelper: + runs-on: ubuntu-latest + steps: + - name: Check if Julia is already available in the PATH + id: julia_in_path + run: which julia + continue-on-error: true + - name: Install Julia, but only if it is not already available in the PATH + uses: julia-actions/setup-julia@v2 + with: + version: '1' + arch: ${{ runner.arch }} + if: steps.julia_in_path.outcome != 'success' + - name: "Add the General registry via Git" + run: | + import Pkg + ENV["JULIA_PKG_SERVER"] = "" + Pkg.Registry.add("General") + shell: julia --color=yes {0} + - name: "Install CompatHelper" + run: | + import Pkg + name = "CompatHelper" + uuid = "aa819f21-2bde-4658-8897-bab36330d9b7" + version = "3" + Pkg.add(; name, uuid, version) + shell: julia --color=yes {0} + - name: "Run CompatHelper" + run: | + import CompatHelper + CompatHelper.main() + shell: julia --color=yes {0} + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} diff --git a/.github/workflows/docs.yml b/.github/workflows/Documentation.yml similarity index 90% rename from .github/workflows/docs.yml rename to .github/workflows/Documentation.yml index f3269903..d19c2f5e 100644 --- a/.github/workflows/docs.yml +++ b/.github/workflows/Documentation.yml @@ -1,5 +1,9 @@ name: Documentation +permissions: + contents: write + pages: write + on: push: branches: @@ -8,6 +12,7 @@ on: - 'release-' tags: '*' pull_request: + workflow_dispatch: jobs: build: @@ -31,5 +36,4 @@ jobs: - name: Build and deploy env: GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # For authentication with GitHub Actions token - DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # For authentication with SSH deploy key run: julia --project=docs/ docs/make.jl \ No newline at end of file diff --git a/.github/workflows/formatter.yml b/.github/workflows/FormatCheck.yml similarity index 60% rename from .github/workflows/formatter.yml rename to .github/workflows/FormatCheck.yml index a57dbf68..1f0d0515 100644 --- a/.github/workflows/formatter.yml +++ b/.github/workflows/FormatCheck.yml @@ -1,4 +1,4 @@ -name: format-check +name: FormatCheck on: push: @@ -35,7 +35,6 @@ jobs: julia -e 'using Pkg; Pkg.add(PackageSpec(name="JuliaFormatter"))' julia -e 'using JuliaFormatter; format(".", verbose=true)' - name: Format check - id: format run: | julia -e ' out = Cmd(`git diff --name-only`) |> read |> String @@ -45,23 +44,4 @@ jobs: @error "Some files have not been formatted !!!" write(stdout, out) exit(1) - end' - - - name: Create pull request - if: ${{ failure() && steps.format.conclusion == 'failure' }} - id: cpr - uses: peter-evans/create-pull-request@v3 - with: - token: ${{ secrets.GITHUB_TOKEN }} - commit-message: Format .jl files - title: 'Automatic JuliaFormatter.jl run' - base: ${{ github.head_ref }} - branch: auto-juliaformatter-pr - delete-branch: true - labels: formatting, automated pr, no changelog - - - name: Check outputs - if: ${{ success() && steps.cpr.conclusion == 'success' }} - run: | - echo "Pull Request Number - ${{ steps.cpr.outputs.pull-request-number }}" - echo "Pull Request URL - ${{ steps.cpr.outputs.pull-request-url }}" + end' \ No newline at end of file diff --git a/.github/workflows/TagBot.yml b/.github/workflows/TagBot.yml new file mode 100644 index 00000000..c466a764 --- /dev/null +++ b/.github/workflows/TagBot.yml @@ -0,0 +1,34 @@ +name: TagBot + +on: + issue_comment: + types: + - created + workflow_dispatch: + inputs: + lookback: + default: 3 + +permissions: + actions: read + checks: read + contents: write + deployments: read + issues: read + discussions: read + packages: read + pages: read + pull-requests: read + repository-projects: read + security-events: read + statuses: read + +jobs: + TagBot: + if: github.event_name == 'workflow_dispatch' || github.actor == 'JuliaTagBot' + runs-on: ubuntu-latest + steps: + - uses: JuliaRegistries/TagBot@v1 + with: + token: ${{ secrets.GITHUB_TOKEN }} + ssh: ${{ secrets.DOCUMENTER_KEY }} \ No newline at end of file diff --git a/.gitignore b/.gitignore index 3b29f370..c7634e80 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,4 @@ docs/build Manifest.toml examples/.ipynb_checkpoints +.vscode/ \ No newline at end of file diff --git a/LICENSE b/LICENSE new file mode 100644 index 00000000..c0dcc7ad --- /dev/null +++ b/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2024 QuantumGroup@UGent + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/Project.toml b/Project.toml index 2826518e..f8b9f357 100644 --- a/Project.toml +++ b/Project.toml @@ -1,25 +1,43 @@ name = "PEPSKit" uuid = "52969e89-939e-4361-9b68-9bc7cde4bdeb" -authors = ["maarten "] +authors = ["Maarten Vandamme", "Paul Brehmer", "Lander Burgelman", "Rui-Zhen Huang", "Daan Maertens", "Lukas Devos - - PEPSKit.jl logo + + PEPSKit.jl logo -# PepsKit.jl +# PEPSKit.jl [![docs][docs-dev-img]][docs-dev-url] ![CI][ci-url] [![codecov][codecov-img]][codecov-url] [docs-dev-img]: https://img.shields.io/badge/docs-dev-blue.svg -[docs-dev-url]: https://quantumghent.github.io/PEPSKit.jl/dev/ +[docs-dev-url]: https://QuantumKitHub.github.io/PEPSKit.jl/dev/ -[codecov-img]: https://codecov.io/gh/quantumghent/PEPSKit.jl/branch/master/graph/badge.svg -[codecov-url]: https://codecov.io/gh/quantumghent/PEPSKit.jl +[codecov-img]: https://codecov.io/gh/QuantumKitHub/PEPSKit.jl/graph/badge.svg?token=1OBDY03SUP +[codecov-url]: https://codecov.io/gh/QuantumKitHub/PEPSKit.jl -[ci-url]: https://github.com/quantumghent/PEPSKit.jl/workflows/CI/badge.svg +[ci-url]: https://github.com/QuantumKitHub/PEPSKit.jl/workflows/CI/badge.svg **Tools for working with projected entangled-pair states** -It contracts, it optimizes, it may be broken at any point. +It contracts, it optimizes, it may break. + +## Installation + +The package can be installed through the Julia general registry, via the package manager: + +```julia-repl +pkg> add PEPSKit +``` + +## Quickstart + +After following the installation process, it should now be possible to load the packages and start simulating. +For example, in order to obtain the groundstate of the 2D Heisenberg model, we can use the following code: + +```julia +using TensorKit, PEPSKit, KrylovKit, OptimKit + +# constructing the Hamiltonian: +Jx, Jy, Jz = (-1, 1, -1) # sublattice rotation to obtain single-site unit cell +physical_space = ComplexSpace(2) +T = ComplexF64 +σx = TensorMap(T[0 1; 1 0], physical_space, physical_space) +σy = TensorMap(T[0 im; -im 0], physical_space, physical_space) +σz = TensorMap(T[1 0; 0 -1], physical_space, physical_space) +H = (Jx * σx ⊗ σx) + (Jy * σy ⊗ σy) + (Jz * σz ⊗ σz) +Heisenberg_hamiltonian = NLocalOperator{NearestNeighbor}(H / 4) + +# configuring the parameters +D = 2 +chi = 20 +ctm_alg = CTMRG(; trscheme = truncdim(chi), tol=1e-20, miniter=4, maxiter=100, verbosity=1) +opt_alg = PEPSOptimize(; + boundary_alg=ctm_alg, + optimizer=LBFGS(4; maxiter=100, gradtol=1e-4, verbosity=2), + gradient_alg=GMRES(; tol=1e-6, maxiter=100), + reuse_env=true, + verbosity=2, +) + +# ground state search +state = InfinitePEPS(2, D) +ctm = leading_boundary(CTMRGEnv(state; Venv=ComplexSpace(chi)), state, ctm_alg) +result = fixedpoint(state, Heisenberg_hamiltonian, opt_alg, ctm) + +@show result.E # -0.6625... +``` + diff --git a/docs/Project.toml b/docs/Project.toml index 3a52a5db..818c9dd4 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,5 +1,6 @@ [deps] Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +PEPSKit = "52969e89-939e-4361-9b68-9bc7cde4bdeb" [compat] Documenter = "0.27" diff --git a/docs/make.jl b/docs/make.jl index 143b0c4b..1814ee54 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -25,4 +25,4 @@ makedocs(; ], ) -deploydocs(; repo="https://github.com/quantumghent/PEPSKit.jl.git") +deploydocs(; repo="https://github.com/QuantumKitHub/PEPSKit.jl.git") diff --git a/docs/src/examples/index.md b/docs/src/examples/index.md index fcc39fe6..f33fdb82 100644 --- a/docs/src/examples/index.md +++ b/docs/src/examples/index.md @@ -1 +1 @@ -Coming soon. \ No newline at end of file +For now, refer to the [examples folder](https://github.com/QuantumKitHub/PEPSKit.jl/tree/master/examples) on GitHub. \ No newline at end of file diff --git a/docs/src/index.md b/docs/src/index.md index 9da83b7a..305f6d84 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -2,4 +2,50 @@ **Tools for working with projected entangled-pair states** -It contracts, it optimizes, it may be broken at any point. +It contracts, it optimizes, it may break. + +## Installation + +The package can be installed through the Julia general registry, via the package manager: + +```julia-repl +pkg> add PEPSKit +``` + +## Quickstart + +After following the installation process, it should now be possible to load the packages and start simulating. +For example, in order to obtain the groundstate of the 2D Heisenberg model, we can use the following code: + +```julia +using TensorKit, PEPSKit, KrylovKit, OptimKit + +# constructing the Hamiltonian: +Jx, Jy, Jz = (-1, 1, -1) # sublattice rotation to obtain single-site unit cell +physical_space = ComplexSpace(2) +T = ComplexF64 +σx = TensorMap(T[0 1; 1 0], physical_space, physical_space) +σy = TensorMap(T[0 im; -im 0], physical_space, physical_space) +σz = TensorMap(T[1 0; 0 -1], physical_space, physical_space) +H = (Jx * σx ⊗ σx) + (Jy * σy ⊗ σy) + (Jz * σz ⊗ σz) +Heisenberg_hamiltonian = NLocalOperator{NearestNeighbor}(H / 4) + +# configuring the parameters +D = 2 +chi = 20 +ctm_alg = CTMRG(; trscheme = truncdim(chi), tol=1e-20, miniter=4, maxiter=100, verbosity=1) +opt_alg = PEPSOptimize(; + boundary_alg=ctm_alg, + optimizer=LBFGS(4; maxiter=100, gradtol=1e-4, verbosity=2), + gradient_alg=GMRES(; tol=1e-6, maxiter=100), + reuse_env=true, + verbosity=2, +) + +# ground state search +state = InfinitePEPS(2, D) +ctm = leading_boundary(CTMRGEnv(state; Venv=ComplexSpace(chi)), state, ctm_alg) +result = fixedpoint(state, Heisenberg_hamiltonian, opt_alg, ctm) + +@show result.E # -0.6625... +``` diff --git a/docs/src/lib/lib.md b/docs/src/lib/lib.md index fcc39fe6..b77bf0a4 100644 --- a/docs/src/lib/lib.md +++ b/docs/src/lib/lib.md @@ -1 +1,6 @@ -Coming soon. \ No newline at end of file +# Library + +```@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 00000000..abd52e17 --- /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/heisenberg.jl b/examples/heisenberg.jl index 7e8581d7..0cfd8c06 100644 --- a/examples/heisenberg.jl +++ b/examples/heisenberg.jl @@ -1,140 +1,32 @@ -using Revise, PEPSKit, TensorKit, Zygote, MPSKit -using MPSKitModels, LinearAlgebra, OptimKit -using PEPSKit: - NORTH, SOUTH, WEST, EAST, NORTHWEST, NORTHEAST, SOUTHEAST, SOUTHWEST, @diffset -using JLD2, ChainRulesCore - -function two_site_rho(r::Int, c::Int, ψ::InfinitePEPS, env::PEPSKit.CTMRGEnv) - cp = mod1(c + 1, size(ψ, 2)) - @tensor ρ[-11, -20; -12, -18] := - env.corners[NORTHWEST, r, c][1, 3] * - env.edges[WEST, r, c][2, 7, 9, 1] * - env.corners[SOUTHWEST, r, c][4, 2] * - env.edges[NORTH, r, c][3, 5, 8, 13] * - env.edges[SOUTH, r, c][14, 6, 10, 4] * - ψ[r, c][-12, 5, 15, 6, 7] * - conj(ψ[r, c][-11, 8, 19, 10, 9]) * - env.edges[NORTH, r, cp][13, 16, 22, 23] * - env.edges[SOUTH, r, cp][28, 17, 21, 14] * - ψ[r, cp][-18, 16, 25, 17, 15] * - conj(ψ[r, cp][-20, 22, 26, 21, 19]) * - env.corners[NORTHEAST, r, cp][23, 24] * - env.edges[EAST, r, cp][24, 25, 26, 27] * - env.corners[SOUTHEAST, r, cp][27, 28] - return ρ -end - -function iCtmGsEh( - ψ::InfinitePEPS, env::PEPSKit.CTMRGEnv, H::AbstractTensorMap{S,2,2} -) where {S} - #Es = Matrix{eltype(H)}(undef,size(ψ,1),size(ψ,2)) - E = 0.0 - for r in 1:size(ψ, 1), c in 1:size(ψ, 2) - ρ = two_site_rho(r, c, ψ, env) - @tensor nn = ρ[1 2; 1 2] - @tensor Eh = H[1 2; 3 4] * ρ[1 2; 3 4] - Eh = Eh / nn - E = E + Eh - #@diffset Es[r,c] = Eh; - end - return real(E) -end - -function H_expectation_value( - ψ::InfinitePEPS, env::PEPSKit.CTMRGEnv, H::AbstractTensorMap{S,2,2} -) where {S} - Eh = iCtmGsEh(ψ, env, H) - ψ1 = rotl90(ψ) - env1 = PEPSKit.rotate_north(env, EAST) - Ev = iCtmGsEh(ψ1, env1, H) - E = real(Eh + Ev) - return E -end - -function SqLatHeisenberg() - Sx, Sy, Sz, _ = spinmatrices(1//2) - - Dphys = ComplexSpace(2) - σx = TensorMap(Sx, Dphys, Dphys) - σy = TensorMap(Sy, Dphys, Dphys) - σz = TensorMap(Sz, Dphys, Dphys) - - @tensor H[-1 -3; -2 -4] := - -σx[-1, -2] * σx[-3, -4] + σy[-1, -2] * σy[-3, -4] + -σz[-1, -2] * σz[-3, -4] - - return H -end - -H = SqLatHeisenberg() - -function cfun(x) - (ψ, env) = x - - function fun(peps) - env = leading_boundary(peps, alg_ctm, env) - x = H_expectation_value(peps, env, H) - return x - end - env = leading_boundary(ψ, alg_ctm, env) - E = H_expectation_value(ψ, env, H) - ∂E = fun'(ψ) - - @assert !isnan(norm(∂E)) - return E, ∂E -end - -# my_retract is not an in place function which should not change x -function my_retract(x, dx, α::Number) - (ϕ, env0) = x - ψ = deepcopy(ϕ) - env = deepcopy(env0) - ψ.A .+= dx.A .* α - #env = leading_boundary(ψ, alg_ctm,env) - return (ψ, env), dx -end - -my_inner(x, dx1, dx2) = real(dot(dx1, dx2)) - -function my_add!(Y, X, a) - Y.A .+= X.A .* a - return Y -end - -function my_scale!(η, β) - rmul!(η.A, β) - return η -end - -function init_psi(d::Int, D::Int, Lx::Int, Ly::Int) - Pspaces = fill(ℂ^d, Lx, Ly) - Nspaces = fill(ℂ^D, Lx, Ly) - Espaces = fill(ℂ^D, Lx, Ly) - - Sspaces = adjoint.(circshift(Nspaces, (1, 0))) - Wspaces = adjoint.(circshift(Espaces, (0, -1))) - - A = map(Pspaces, Nspaces, Espaces, Sspaces, Wspaces) do P, N, E, S, W - return TensorMap(rand, ComplexF64, P ← N * E * S * W) - end - - return InfinitePEPS(A) -end - -alg_ctm = CTMRG(; verbose=1, tol=1e-4, trscheme=truncdim(10), miniter=4, maxiter=200) - -function main(; d=2, D=2, Lx=1, Ly=1) - ψ = init_psi(d, D, Lx, Ly) - env = leading_boundary(ψ, alg_ctm) - optimize( - cfun, - (ψ, env), - ConjugateGradient(; verbosity=2); - inner=my_inner, - retract=my_retract, - (scale!)=my_scale!, - (add!)=my_add!, - ) - return ψ -end - -main() +using LinearAlgebra +using TensorKit, OptimKit +using PEPSKit, KrylovKit + +# Square lattice Heisenberg Hamiltonian +# We use the parameters (J₁, J₂, J₃) = (-1, 1, -1) by default to capture +# the ground state in a single-site unit cell. This can be seen from +# sublattice rotating H from parameters (1, 1, 1) to (-1, 1, -1). +H = square_lattice_heisenberg(; Jx=-1, Jy=1, Jz=-1) + +# Parameters +χbond = 2 +χenv = 20 +ctmalg = CTMRG(; trscheme=truncdim(χenv), tol=1e-10, miniter=4, maxiter=100, verbosity=1) +alg = PEPSOptimize(; + boundary_alg=ctmalg, + optimizer=LBFGS(4; maxiter=100, gradtol=1e-4, verbosity=2), + gradient_alg=GMRES(; tol=1e-6, maxiter=100), + reuse_env=true, + verbosity=2, +) + +# Ground state search +# We initialize a random PEPS with bond dimension χbond and from that converge +# a CTMRG environment with dimension χenv on the environment bonds before +# starting the optimization. The ground-state energy should approximately approach +# E/N = −0.6694421, which is a QMC estimate from https://arxiv.org/abs/1101.3281. +# Of course there is a noticable bias for small χbond and χenv. +ψ₀ = InfinitePEPS(2, χbond) +env₀ = leading_boundary(CTMRGEnv(ψ₀; Venv=ℂ^χenv), ψ₀, ctmalg) +result = fixedpoint(ψ₀, H, alg, env₀) +@show result.E diff --git a/examples/ising_pepo.jl b/examples/ising_pepo.jl new file mode 100644 index 00000000..4936e602 --- /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/examples/test_pepo.jl b/examples/test_pepo.jl deleted file mode 100644 index 46db0e95..00000000 --- a/examples/test_pepo.jl +++ /dev/null @@ -1,129 +0,0 @@ -# this should run once the compatocalypse has been resolved - -using Revise -using TensorOperations -using TensorKit -using MPSKit -using PEPSKit -using OptimKit - -## Set some options - -optim_method = ConjugateGradient -optim_tol = 1e-5 -optim_maxiter = 100 -verbosity = 2 # >= 5 for debugging -boundary_maxiter = 50 -tol_min = 1e-12 -tol_max = 1e-4 -tol_factor = 1e-3 -symm = Full() -hermitian = true -boundary_method = VUMPS(; - tol_galerkin=tol_max, - maxiter=boundary_maxiter, - dynamical_tols=true, - eigs_tolfactor=1e-3, - envs_tolfactor=1e-3, - gauge_tolfactor=1e-6, - tol_max=1e-4, - verbose=verbosity >= 5, -) -# boundary_method = GradientGrassmann - -## Set model parameters - -# 3D Ising temperature -beta = 1 / 3 # should give f = -1.0219139... @ D = 3 - -# bond dimensions -D = 3 -χ = 30 - -# pick size of unit cell -depth = 1 -width = 1 -height = 1 - -## PEPO defintion: 3D classical Ising model - -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] - -M = zeros(2, 2, 2, 2, 2, 2) -M[1, 1, 1, 1, 1, 1] = 1 -M[2, 2, 2, 2, 2, 2] = -1 -@tensor m[-1 -2; -3 -4 -5 -6] := - M[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)') -m = TensorMap(m, ℂ^2 ⊗ (ℂ^2)' ← ℂ^2 ⊗ ℂ^2 ⊗ (ℂ^2)' ⊗ (ℂ^2)') - -O = InfinitePEPO(repeat([o], depth, width, height)) - -## Optimization algorithm - -pepo_alg = PEPOOptimize(; - optim_method, - optim_tol, - optim_maxiter, - verbosity, - boundary_method, - boundary_maxiter, - tol_min, - tol_max, - tol_factor, - symm, - hermitian, -) - -## Initialize PEPS and MPS fixed points - -peps = symmetrize(initializePEPS(O, ℂ^D), symm) -envs = pepo_opt_environments( - peps, O, pepo_alg.boundary_method; vspaces=[ℂ^χ], hermitian=pepo_alg.hermitian -) -normalize!(peps, envs.peps_boundary) - -## Perform optimization - -x, f, normgrad = leading_boundary(peps, O, pepo_alg, envs) - -## Unpack result and check magnetization - -peps = x.state -above = x.envs.pepo_boundary.boundaries[1] -below = x.envs.pepo_boundary.boundaries[2] -lw, rw = x.envs.pepo_boundary.envs[3] - -pos = (1, 1) # pick a site - -if height == 1 # no crazy contractions for now... - row, col = pos - fliprow = size(peps, 1) - row + 1 # below starts counting from below - mag = @tensor lw[row, col][13 8 10 18; 1] * - above.AC[row, col][1 9 11 15; 12] * - peps[row, col][5; 9 3 4 8] * - m[14 5; 11 6 7 10] * - rw[row, col][12 3 6 16; 2] * - conj(below.AC[fliprow, col][13 4 7 17; 2]) * - conj(peps[row, col][14; 15 16 17 18]) - - λ = @tensor lw[row, col][13 8 10 18; 1] * - above.AC[row, col][1 9 11 15; 12] * - peps[row, col][5; 9 3 4 8] * - o[14 5; 11 6 7 10] * - rw[row, col][12 3 6 16; 2] * - conj(below.AC[fliprow, col][13 4 7 17; 2]) * - conj(peps[row, col][14; 15 16 17 18]) - - println("m = $(mag / λ)") # should be ~0.945 @ beta = 1/3 and D = 3 -end - -nothing diff --git a/examples/test_vumps.jl b/examples/test_vumps.jl deleted file mode 100644 index b3492c2c..00000000 --- a/examples/test_vumps.jl +++ /dev/null @@ -1,8 +0,0 @@ -using Revise, PEPSKit, TensorKit, Zygote, MPSKit - -p = InfinitePEPS(fill(ℂ^2, 1, 1), fill(ℂ^2, 1, 1)); - -trans = PEPSKit.InfiniteTransferPEPS(p, 1, 1); -mps = PEPSKit.initializeMPS(trans, [ℂ^5]); - -leading_boundary(mps, trans, VUMPS()) diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index f195b882..f0cbc7fb 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -1,18 +1,18 @@ module PEPSKit +using LinearAlgebra, Statistics, Base.Threads, Base.Iterators, Printf +using Base: @kwdef +using Compat using Accessors using VectorInterface -using TensorKit, - KrylovKit, MPSKit, OptimKit, Base.Threads, Base.Iterators, Parameters, Printf -using ChainRulesCore - -using LinearAlgebra - -export CTMRG, CTMRG2 -export leading_boundary +using TensorKit, KrylovKit, MPSKit, OptimKit +using ChainRulesCore, Zygote include("utility/util.jl") include("utility/svd.jl") +include("utility/eigsolve.jl") +include("utility/rotations.jl") +include("utility/hook_pullback.jl") include("states/abstractpeps.jl") include("states/infinitepeps.jl") @@ -26,28 +26,52 @@ include("mpskit_glue/transferpeps_environments.jl") include("mpskit_glue/transferpepo_environments.jl") include("environments/ctmrgenv.jl") -include("environments/boundarympsenv.jl") +include("operators/localoperator.jl") +include("operators/models.jl") include("algorithms/ctmrg.jl") -include("algorithms/expval.jl") +include("algorithms/peps_opt.jl") include("utility/symmetrization.jl") -include("algorithms/pepo_opt.jl") -include("utility/rotations.jl") +""" + module Defaults + const ctmrg_maxiter = 100 + const ctmrg_miniter = 4 + const ctmrg_tol = 1e-12 + const fpgrad_maxiter = 100 + const fpgrad_tol = 1e-6 + end + +Module containing default values that represent typical algorithm parameters. -#default settings +- `ctmrg_maxiter = 100`: Maximal number of CTMRG iterations per run +- `ctmrg_miniter = 4`: Minimal number of CTMRG carried out +- `ctmrg_tol = 1e-12`: Tolerance checking singular value and norm convergence +- `fpgrad_maxiter = 100`: Maximal number of iterations for computing the CTMRG fixed-point gradient +- `fpgrad_tol = 1e-6`: Convergence tolerance for the fixed-point gradient iteration +""" module Defaults - const maxiter = 100 - const tol = 1e-12 + const ctmrg_maxiter = 100 + const ctmrg_miniter = 4 + const ctmrg_tol = 1e-12 + const fpgrad_maxiter = 100 + const fpgrad_tol = 1e-6 end -export FullSVD, IterSVD, OldSVD, svdwrap +export FullSVD, IterSVD, OldSVD, svdwrap, itersvd +export CTMRG, CTMRGEnv +export NLocalOperator, AnisotropicNNOperator, OnSite, NearestNeighbor +export expectation_value, costfun +export leading_boundary +export PEPSOptimize, GeomSum, ManualIter, LinSolve +export fixedpoint + export InfinitePEPS, InfiniteTransferPEPS export InfinitePEPO, InfiniteTransferPEPO export initializeMPS, initializePEPS -export PEPOOptimize, pepo_opt_environments export symmetrize, None, Depth, Full -export itersvd +export showtypeofgrad +export square_lattice_heisenberg, square_lattice_pwave end # module diff --git a/src/algorithms/ctmrg.jl b/src/algorithms/ctmrg.jl index 65312e0c..a69d6da9 100644 --- a/src/algorithms/ctmrg.jl +++ b/src/algorithms/ctmrg.jl @@ -41,29 +41,30 @@ function MPSKit.leading_boundary(envinit, state, alg::CTMRG) for i in 1:(alg.maxiter) env, ϵ = ctmrg_iter(state, env, alg) # Grow and renormalize in all 4 directions - # Compute convergence criteria and take max (TODO: How should we handle logging all of this?) - Δϵ = abs((ϵold - ϵ) / ϵold) - normnew = norm(state, env) - Δnorm = abs(normold - normnew) / abs(normold) - CSnew = map(c -> tsvd(c; alg=TensorKit.SVD())[2], env.corners) - ΔCS = maximum(zip(CSold, CSnew)) do (c_old, c_new) - # only compute the difference on the smallest part of the spaces - smallest = infimum(MPSKit._firstspace(c_old), MPSKit._firstspace(c_new)) - e_old = isometry(MPSKit._firstspace(c_old), smallest) - e_new = isometry(MPSKit._firstspace(c_new), smallest) - return norm(e_new' * c_new * e_new - e_old' * c_old * e_old) - end - TSnew = map(t -> tsvd(t; alg=TensorKit.SVD())[2], env.edges) + conv_condition, normold, CSold, TSold, ϵ = ignore_derivatives() do + # Compute convergence criteria and take max (TODO: How should we handle logging all of this?) + Δϵ = abs((ϵold - ϵ) / ϵold) + normnew = norm(state, env) + Δnorm = abs(normold - normnew) / abs(normold) + CSnew = map(c -> tsvd(c; alg=TensorKit.SVD())[2], env.corners) + ΔCS = maximum(zip(CSold, CSnew)) do (c_old, c_new) + # only compute the difference on the smallest part of the spaces + smallest = infimum(MPSKit._firstspace(c_old), MPSKit._firstspace(c_new)) + e_old = isometry(MPSKit._firstspace(c_old), smallest) + e_new = isometry(MPSKit._firstspace(c_new), smallest) + return norm(e_new' * c_new * e_new - e_old' * c_old * e_old) + end + TSnew = map(t -> tsvd(t; alg=TensorKit.SVD())[2], env.edges) - ΔTS = maximum(zip(TSold, TSnew)) do (t_old, t_new) - MPSKit._firstspace(t_old) == MPSKit._firstspace(t_new) || - return scalartype(t_old)(Inf) - # TODO: implement when spaces aren't the same - return norm(t_new - t_old) - end + ΔTS = maximum(zip(TSold, TSnew)) do (t_old, t_new) + MPSKit._firstspace(t_old) == MPSKit._firstspace(t_new) || + return scalartype(t_old)(Inf) + # TODO: implement when spaces aren't the same + return norm(t_new - t_old) + end + + conv_condition = max(Δnorm, ΔCS, ΔTS) < alg.tol && i > alg.miniter - conv_condition = max(Δnorm, ΔCS, ΔTS) < alg.tol && i > alg.miniter - ignore_derivatives() do # Print verbose info if alg.verbosity > 1 || (alg.verbosity == 1 && (i == 1 || conv_condition)) @printf( "CTMRG iter: %3d norm: %.2e Δnorm: %.2e ΔCS: %.2e ΔTS: %.2e ϵ: %.2e Δϵ: %.2e\n", @@ -81,14 +82,9 @@ function MPSKit.leading_boundary(envinit, state, alg::CTMRG) @warn( "CTMRG reached maximal number of iterations at (Δnorm=$Δnorm, ΔCS=$ΔCS, ΔTS=$ΔTS)" ) + return conv_condition, normnew, CSnew, TSnew, ϵ end conv_condition && break # Converge if maximal Δ falls below tolerance - - # Update convergence criteria - normold = normnew - CSold = CSnew - TSold = TSnew - ϵold = ϵ end # Do one final iteration that does not change the spaces @@ -97,7 +93,7 @@ function MPSKit.leading_boundary(envinit, state, alg::CTMRG) ) env′, = ctmrg_iter(state, env, alg_fixed) envfix = gauge_fix(env, env′) - check_elementwise_convergence(env, envfix; atol=alg.tol^(3 / 4)) || + check_elementwise_convergence(env, envfix; atol=alg.tol^(1 / 2)) || @warn "CTMRG did not converge elementwise." return envfix end @@ -146,13 +142,14 @@ function gauge_fix(envprev::CTMRGEnv{C,T}, envfinal::CTMRGEnv{C,T}) where {C,T} scalartype(T), MPSKit._lastspace(Tsfinal[end])' ← MPSKit._lastspace(M[end])', ) - ρprev = eigsolve(TransferMatrix(Tsprev, M), ρinit, 1, :LM)[2][1] - ρfinal = eigsolve(TransferMatrix(Tsfinal, M), ρinit, 1, :LM)[2][1] + + ρ_prev = transfermatrix_fixedpoint(Tsprev, M, ρinit) + ρ_final = transfermatrix_fixedpoint(Tsfinal, M, ρinit) # Decompose and multiply - Up, _, Vp = tsvd(ρprev) - Uf, _, Vf = tsvd(ρfinal) + Up, _, Vp = tsvd!(ρ_prev) Qprev = Up * Vp + Uf, _, Vf = tsvd!(ρ_final) Qfinal = Uf * Vf σ = Qprev * Qfinal' @@ -162,17 +159,25 @@ function gauge_fix(envprev::CTMRGEnv{C,T}, envfinal::CTMRGEnv{C,T}) where {C,T} cornersfix, edgesfix = fix_relative_phases(envfinal, signs) # Fix global phase - cornersgfix = map(zip(envprev.corners, cornersfix)) do (Cprev, Cfix) - φ = dot(Cprev, Cfix) - φ' * Cfix + cornersgfix = map(envprev.corners, cornersfix) do Cprev, Cfix + return dot(Cfix, Cprev) * Cfix end - edgesgfix = map(zip(envprev.edges, edgesfix)) do (Tprev, Tfix) - φ = dot(Tprev, Tfix) - φ' * Tfix + edgesgfix = map(envprev.edges, edgesfix) do Tprev, Tfix + return dot(Tfix, Tprev) * Tfix end - envfix = CTMRGEnv(cornersgfix, edgesgfix) + return CTMRGEnv(cornersgfix, edgesgfix) +end - return envfix +# this is a bit of a hack to get the fixed point of the mixed transfer matrix +# because MPSKit is not compatible with AD +function transfermatrix_fixedpoint(tops, bottoms, ρinit) + _, vecs, info = eigsolve(ρinit, 1, :LM, Arnoldi()) do ρ + return foldr(zip(tops, bottoms); init=ρ) do (top, bottom), ρ + return @tensor ρ′[-1; -2] := top[-1 4 3; 1] * conj(bottom[-2 4 3; 2]) * ρ[1; 2] + end + end + info.converged > 0 || @warn "eigsolve did not converge" + return first(vecs) end # Explicit fixing of relative phases (doing this compactly in a loop is annoying) @@ -330,7 +335,8 @@ function left_move(state, env::CTMRGEnv{C,T}, alg::CTMRG) where {C,T} else alg.trscheme end - U, S, V, ϵ_local = svdwrap(Q_sw * Q_nw, alg.svdalg; trunc=trscheme) + @tensor QQ[-1 -2 -3; -4 -5 -6] := Q_sw[-1 -2 -3; 1 2 3] * Q_nw[1 2 3; -4 -5 -6] + U, S, V, ϵ_local = svdwrap(QQ, alg.svdalg; trunc=trscheme) ϵ = max(ϵ, ϵ_local / norm(S)) # TODO: check if we can just normalize enlarged corners s.t. trunc behaves a bit better @@ -410,8 +416,8 @@ function build_projectors( U::AbstractTensorMap{E,3,1}, S, V::AbstractTensorMap{E,1,3}, Q_sw, Q_nw ) where {E<:ElementarySpace} isqS = sdiag_inv_sqrt(S) - @tensor Pl[-1 -2 -3; -4] := Q_nw[-1 -2 -3; 1 2 3] * conj(V[4; 1 2 3]) * isqS[4; -4] - @tensor Pr[-1; -2 -3 -4] := isqS[-1; 1] * conj(U[2 3 4; 1]) * Q_sw[2 3 4; -2 -3 -4] + Pl = Q_nw * V' * isqS + Pr = isqS * U' * Q_sw return Pl, Pr end @@ -449,12 +455,10 @@ function LinearAlgebra.norm(peps::InfinitePEPS, env::CTMRGEnv) peps[r, c][17; 6 10 14 2] * conj(peps[r, c][17; 7 11 15 3]) - total *= tr( - env.corners[NORTHWEST, r, c] * - env.corners[NORTHEAST, r, mod1(c - 1, end)] * - env.corners[SOUTHEAST, mod1(r - 1, end), mod1(c - 1, end)] * - env.corners[SOUTHWEST, mod1(r - 1, end), c], - ) + total *= @tensor env.corners[NORTHWEST, r, c][1; 2] * + env.corners[NORTHEAST, r, mod1(c - 1, end)][2; 3] * + env.corners[SOUTHEAST, mod1(r - 1, end), mod1(c - 1, end)][3; 4] * + env.corners[SOUTHWEST, mod1(r - 1, end), c][4; 1] total /= @tensor env.edges[WEST, r, c][1 10 11; 4] * env.corners[NORTHWEST, r, c][4; 5] * diff --git a/src/algorithms/expval.jl b/src/algorithms/expval.jl deleted file mode 100644 index 99a7b1a5..00000000 --- a/src/algorithms/expval.jl +++ /dev/null @@ -1,34 +0,0 @@ -function MPSKit.expectation_value(envs::CTMRGEnv, ham::AbstractTensorMap{S,1,1}) where {S} - @assert envs.peps_above == envs.peps_below - peps = envs.peps_above - result = Matrix{eltype(ham)}(undef, size(peps, 1), size(peps, 2)) - - for r in 1:size(peps, 1), c in 1:size(peps, 2) - e = @tensor envs.edges[WEST, r, c][1 2 3; 4] * - envs.corners[NORTHWEST, r, c][4; 5] * - envs.edges[NORTH, r, c][5 6 7; 8] * - envs.corners[NORTHEAST, r, c][8; 9] * - envs.edges[EAST, r, c][9 10 11; 12] * - envs.corners[SOUTHEAST, r, c][12; 13] * - envs.edges[SOUTH, r, c][13 14 15; 16] * - envs.corners[SOUTHWEST, r, c][16; 1] * - peps[r, c][17; 6 10 14 2] * - conj(peps[r, c][18; 7 11 15 3]) * - ham[18; 17] - - n = @tensor envs.edges[WEST, r, c][1 2 3; 4] * - envs.corners[NORTHWEST, r, c][4; 5] * - envs.edges[NORTH, r, c][5 6 7; 8] * - envs.corners[NORTHEAST, r, c][8; 9] * - envs.edges[EAST, r, c][9 10 11; 12] * - envs.corners[SOUTHEAST, r, c][12; 13] * - envs.edges[SOUTH, r, c][13 14 15; 16] * - envs.corners[SOUTHWEST, r, c][16; 1] * - peps[r, c][17; 6 10 14 2] * - conj(peps[r, c][17; 7 11 15 3]) - - @diffset result[r, c] = e / n - end - - return result -end diff --git a/src/algorithms/pepo_opt.jl b/src/algorithms/pepo_opt.jl deleted file mode 100644 index 7182ccd5..00000000 --- a/src/algorithms/pepo_opt.jl +++ /dev/null @@ -1,263 +0,0 @@ -# First go at PEPO fixed point optimization using PEPSKit - -## Characterize PEPO optimization manifold - -# point on optimization manifold contains PEPS, its 'boundary' -# the norm of the gradient is for convenience for now, but should be removed -mutable struct PEPOOptPoint{T,E,F} - state::T - envs::E - normgrad::F -end -PEPOOptPoint(state, envs) = PEPOOptPoint(state, envs, Inf) - -# the gradient is just a 2D array of PEPS tensors, is this sufficient? -function _retract(x::PEPOOptPoint, dx, α::Number) - # move the state, keep the rest? - new_state = copy(x.state) - new_state.A .+= dx .* α - return PEPOOptPoint(new_state, x.envs, x.normgrad), dx -end - -_inner(x, dx1, dx2) = real(dot(dx1, dx2)) - -function _add!(Y, X, a) - return add!(Y, X, a) -end - -function _scale!(η, β) - return scale!(η, β) -end - -## Utility; TODO: remove all of this - -function _nthroot(x::Real, n::Integer) - return if isodd(n) || x ≥ 0 - copysign(abs(x)^(1//n), x) - else - throw( - DomainError( - "Exponentiation yielding a complex result requires a complex argument. Replace nthroot(x, n) with complex(x)^(1//n).", - ), - ) - end -end -function _hacked_root(x, n) - m = abs(x) - p = angle(x) - mY = _nthroot(m, n) - pYs = p / n .+ 2 * pi / n .* (0:(n - 1)) - I = argmin(minimum([abs.(pYs .- p); abs.(2 * pi .- pYs .+ p)]; dims=1)) # find phase with minimal angular distance - return mY * exp(pYs[I] * 1im) -end - -isverbose(alg::VUMPS) = alg.verbose -isverbose(alg::GradientGrassmann) = alg.method.verbosity >= 0 - -## Characterize environments for PEPO optimization - -mutable struct PEPOOptEnv{T,O,F} <: Cache - peps_boundary::T - pepo_boundary::O - alg::F -end - -# TODO: think about name; just picked something to avoid conflict with MPSKit.environments -function pepo_opt_environments( - peps::InfinitePEPS, - pepo::InfinitePEPO, - boundary_alg; - vspaces=[oneunit(space(peps, 1, 1))], - hermitian=false, - kwargs..., -) - # boundary_alg handles everything - peps_boundary = environments(peps, boundary_alg; vspaces, hermitian) - pepo_boundary = environments(peps, pepo, boundary_alg; vspaces, hermitian) - - return PEPOOptEnv(peps_boundary, pepo_boundary, boundary_alg) -end - -# I did overload recalculate!, this seemed to make sense -function MPSKit.recalculate!( - envs::PEPOOptEnv, - peps::InfinitePEPS, - pepo::InfinitePEPO; - tol=algtol(envs.alg), - hermitian=false, - kwargs..., -) - recalculate!(envs.peps_boundary, peps; tol, hermitian, kwargs...) - recalculate!(envs.pepo_boundary, peps, pepo; tol, hermitian, kwargs...) - - return envs -end - -## PEPO fixed point optimization algorithm - -# first attempt at a PEPO-fixed-point-optimization algorithm, bit of a mess... -struct PEPOOptimize{A} - optim_method::OptimKit.OptimizationAlgorithm - optim_finalize!::Function - boundary_method::A - tol_min::Float64 # some basic tolerance scaling - tol_max::Float64 - tol_factor::Float64 - symm::SymmetrizationStyle # symmetrization - hermitian::Bool - - function PEPOOptimize(; - optim_method=ConjugateGradient, - (optim_finalize!)=OptimKit._finalize!, - optim_tol=Defaults.tol, - optim_maxiter=Defaults.maxiter, - verbosity=2, - boundary_method=VUMPS, - boundary_maxiter=Defaults.maxiter, - boundary_finalize=MPSKit.Defaults._finalize, - tol_min=1e-12, - tol_max=1e-5, - tol_factor=1e-3, - symm=None(), - hermitian=false, - ) - if isa(optim_method, OptimKit.OptimizationAlgorithm) - # We were given an optimisation method, just use it. - m = optim_method - elseif optim_method <: OptimKit.OptimizationAlgorithm - # We were given an optimisation method type, construct an instance of it. - m = optim_method(; - gradtol=optim_tol, maxiter=optim_maxiter, verbosity=verbosity - ) - else - msg = "optim_method should be either an instance or a subtype of `OptimKit.OptimizationAlgorithm`." - throw(ArgumentError(msg)) - end - if isa(boundary_method, Union{<:VUMPS,CTMRG,GradientGrassmann}) - bm = boundary_method - elseif boundary_method <: Union{Type{CTMRG},MPSKit.Algorithm} - # total syntax clusterfuck, need to clean this up - if boundary_method <: VUMPS - bm = boundary_method(; - tol_galerkin=tol_max, - verbose=verbosity >= 5, - maxiter=boundary_maxiter, - finalize=boundary_finalize, - ) - elseif boundary_method <: GradientGrassmann - bm = boundary_method(; - tol=tol_max, verbosity=verbosity - 4, maxiter=boundary_maxiter - ) - elseif boundary_method <: CTMRG - bm = method(; tol=tol_max, verbose=verbosity >= 5, maxiter=boundary_maxiter) - else - msg = "Unknown boundary contraction method." - throw(ArgumentError(msg)) - end - else - msg = "boundary_method should be a valid boundary contraction algorithm." - throw(ArgumentError(msg)) - end - return new{typeof(bm)}( - m, optim_finalize!, bm, tol_min, tol_max, tol_factor, symm, hermitian - ) - end -end - -# default PEPO optimization cost function for given PEPO and optimization algorithm -function pepo_opt_costfun(op::InfinitePEPO, alg::PEPOOptimize) - D, W, H = size(op) - nrm = D * W * H - function pepo_opt_fg(x::PEPOOptPoint) - # unpack state - peps = symmetrize(x.state, alg.symm) - envs = x.envs - ng = x.normgrad - - # recompute environment with scaled tolerance - boundary_tol = min(max(alg.tol_min, ng * alg.tol_factor), alg.tol_max) - recalculate!(envs, peps, op; tol=boundary_tol, hermitian=alg.hermitian) - - # compute cost function - lambdas_peps = expectation_value(peps, envs.peps_boundary) - lambdas_pepo = expectation_value(peps, op, envs.pepo_boundary) - f = -log(real(prod(lambdas_pepo) / prod(lambdas_peps))) / nrm - - # compute gradient - ∂p_peps = ∂∂peps(peps, envs.peps_boundary) - ∂p_pepo = ∂∂peps(peps, op, envs.pepo_boundary) - grad = -(2 / nrm) .* ∂p_pepo ./ lambdas_pepo .+ (2 / nrm) .* ∂p_peps ./ lambdas_peps - grad = symmetrize(grad, alg.symm) - # TODO: test if gradient is actually correct - - # TODO: decide whether we want to: - # Actually update the state in place after symmetrization? - # Update the environments and norm of gradient after each cost function - # evaluation, or move this to the optim_finalize to only update after line search - # has terminated? - x.state = peps # is this sensible? - x.envs = envs - x.normgrad = sum(norm.(grad)) - - # some temporary debugging verbosity, to be removed - if isverbose(alg.boundary_method) - lambda_peps = prod(lambdas_peps) - lambda_pepo = prod(lambdas_pepo) - - lambda_pepo_root = _hacked_root(lambda_pepo, D * W * H) - rel_im_pepo = abs(imag(lambda_pepo_root) / abs(lambda_pepo_root)) - @printf( - "\n\tCurrent lambda_pepo: %f + %fim;\tRelative im. part: %e;\n", - real(lambda_pepo_root), - imag(lambda_pepo_root), - rel_im_pepo, - ) - - lambda_peps_root = _hacked_root(lambda_peps, D * W) - rel_im_peps = abs(imag(lambda_peps_root) / abs(lambda_peps_root)) - @printf( - "\tCurrent lambda_peps: %f + %fim;\tRelative im. part: %e;\n\n", - real(lambda_peps_root), - imag(lambda_peps_root), - rel_im_peps, - ) - - @printf("\tCurrent f: %f\n\n", f) - - lambdas_peps_bis = map(∂p_peps, peps.A) do ∂p, p - @tensor ∂p[1; 2 3 4 5] * conj(p[1; 2 3 4 5]) - end - lambdas_pepo_bis = map(∂p_pepo, peps.A) do ∂p, p - @tensor ∂p[1; 2 3 4 5] * conj(p[1; 2 3 4 5]) - end - @show diff_peps = maximum(abs.(lambdas_peps .- lambdas_peps_bis)) - @show diff_pepo = maximum(abs.(lambdas_pepo .- lambdas_pepo_bis)) - end - - return f, grad - end - return pepo_opt_fg -end - -## The actual leading boundary routine - -function MPSKit.leading_boundary( - state::InfinitePEPS, - op::InfinitePEPO, - alg::PEPOOptimize, - envs=pepo_opt_environments(state, op, alg.boundary_method; hermitian=alg.hermitian), - fg=pepo_opt_costfun(op, alg), -) - res = optimize( - fg, - PEPOOptPoint(state, envs), - alg.optim_method; - retract=_retract, - inner=_inner, - (scale!)=_scale!, - (add!)=_add!, - (finalize!)=alg.optim_finalize!, - ) - (x, fx, gx, numfg, normgradhistory) = res - return x, fx, normgradhistory[end] -end diff --git a/src/algorithms/peps_opt.jl b/src/algorithms/peps_opt.jl new file mode 100644 index 00000000..cdee687f --- /dev/null +++ b/src/algorithms/peps_opt.jl @@ -0,0 +1,208 @@ +abstract type GradMode end + +""" + struct NaiveAD <: GradMode + +Gradient mode for CTMRG using AD. +""" +struct NaiveAD <: GradMode end + +""" + struct GeomSum <: GradMode + +Gradient mode for CTMRG using explicit evaluation of the geometric sum. +""" +@kwdef struct GeomSum <: GradMode + maxiter::Int = Defaults.fpgrad_maxiter + tol::Real = Defaults.fpgrad_tol + verbosity::Int = 0 +end + +""" + struct ManualIter <: GradMode + +Gradient mode for CTMRG using manual iteration to solve the linear problem. +""" +@kwdef struct ManualIter <: GradMode + maxiter::Int = Defaults.fpgrad_maxiter + tol::Real = Defaults.fpgrad_tol + verbosity::Int = 0 +end + +""" + PEPSOptimize{G}(; boundary_alg = CTMRG(), optimizer::OptimKit.OptimizationAlgorithm = LBFGS() + reuse_env::Bool = true, gradient_alg::G, verbosity::Int = 0) + +Algorithm struct that represent PEPS ground-state optimization using AD. +Set the algorithm to contract the infinite PEPS in `boundary_alg`; +currently only `CTMRG` is supported. The `optimizer` computes the gradient directions +based on the CTMRG gradient and updates the PEPS parameters. In this optimization, +the CTMRG runs can be started on the converged environments of the previous optimizer +step by setting `reuse_env` to true. Otherwise a random environment is used at each +step. The CTMRG gradient itself is computed using the `gradient_alg` algorithm. +Different levels of output verbosity can be activated using `verbosity` (0, 1 or 2). +""" +@kwdef struct PEPSOptimize{G} + boundary_alg::CTMRG = CTMRG() # Algorithm to find boundary environment + optimizer::OptimKit.OptimizationAlgorithm = LBFGS( + 4; maxiter=100, gradtol=1e-4, verbosity=2 + ) + reuse_env::Bool = true # Reuse environment of previous optimization as initial guess for next + gradient_alg::G = GeomSum() # Algorithm to solve gradient linear problem + verbosity::Int = 0 +end + +""" + fixedpoint(ψ₀::InfinitePEPS{T}, H, alg::PEPSOptimize, [env₀::CTMRGEnv]) where {T} + +Optimize `ψ₀` with respect to the Hamiltonian `H` according to the parameters supplied +in `alg`. The initial environment `env₀` serves as an initial guess for the first CTMRG run. +By default, a random initial environment is used. +""" +function fixedpoint( + ψ₀::InfinitePEPS{T}, H, alg::PEPSOptimize, env₀::CTMRGEnv=CTMRGEnv(ψ₀; Venv=field(T)^20) +) where {T} + (peps, env), E, ∂E, info = optimize( + (ψ₀, env₀), alg.optimizer; retract=my_retract, inner=my_inner + ) do (peps, envs) + E, g = withgradient(peps) do ψ + envs´ = hook_pullback( + leading_boundary, + envs, + ψ, + alg.boundary_alg; + alg_rrule=alg.gradient_alg, + ) + ignore_derivatives() do + alg.reuse_env && update!(envs, envs´) + end + return costfun(ψ, envs´, H) + end + # withgradient returns tuple of gradients `g` + return E, only(g) + end + return (; peps, env, E, ∂E, info) +end + +# Update PEPS unit cell in non-mutating way +# Note: Both x and η are InfinitePEPS during optimization +function my_retract(x, η, α) + peps = deepcopy(x[1]) + peps.A .+= η.A .* α + env = deepcopy(x[2]) + return (peps, env), η +end + +# Take real valued part of dot product +my_inner(_, η₁, η₂) = real(dot(η₁, η₂)) + +#= +Evaluating the gradient of the cost function for CTMRG: +- The gradient of the cost function for CTMRG can be computed using automatic differentiation (AD) or explicit evaluation of the geometric sum. +- With AD, the gradient is computed by differentiating the cost function with respect to the PEPS tensors, including computing the environment tensors. +- With explicit evaluation of the geometric sum, the gradient is computed by differentiating the cost function with the environment kept fixed, and then manually adding the gradient contributions from the environments. +=# + +function _rrule( + gradmode::Union{GradMode,KrylovKit.LinearSolver}, + ::RuleConfig, + ::typeof(MPSKit.leading_boundary), + envinit, + state, + alg::CTMRG, +) + envs = leading_boundary(envinit, state, alg) + #TODO: fixed space for unit cells + + function leading_boundary_pullback(Δenvs′) + Δenvs = unthunk(Δenvs′) + + # find partial gradients of gauge_fixed single CTMRG iteration + # TODO: make this rrule_via_ad so it's zygote-agnostic + _, env_vjp = pullback(state, envs) do A, x + return gauge_fix(x, ctmrg_iter(A, x, alg)[1]) + end + + # evaluate the geometric sum + ∂f∂A(x)::typeof(state) = env_vjp(x)[1] + ∂f∂x(x)::typeof(envs) = env_vjp(x)[2] + ∂F∂envs = fpgrad(Δenvs, ∂f∂x, ∂f∂A, Δenvs, gradmode) + + return NoTangent(), ZeroTangent(), ∂F∂envs, NoTangent() + end + + return envs, leading_boundary_pullback +end + +@doc """ + fpgrad(∂F∂x, ∂f∂x, ∂f∂A, y0, alg) + +Compute the gradient of the cost function for CTMRG by solving the following equation: + +dx = ∑ₙ (∂f∂x)ⁿ ∂f∂A dA = (1 - ∂f∂x)⁻¹ ∂f∂A dA + +where `∂F∂x` is the gradient of the cost function with respect to the PEPS tensors, `∂f∂x` +is the partial gradient of the CTMRG iteration with respect to the environment tensors, +`∂f∂A` is the partial gradient of the CTMRG iteration with respect to the PEPS tensors, and +`y0` is the initial guess for the fixed-point iteration. The function returns the gradient +`dx` of the fixed-point iteration. +""" +fpgrad + +# TODO: can we construct an implementation that does not need to evaluate the vjp +# twice if both ∂f∂A and ∂f∂x are needed? +function fpgrad(∂F∂x, ∂f∂x, ∂f∂A, _, alg::GeomSum) + g = ∂F∂x + dx = ∂f∂A(g) # n = 0 term: ∂F∂x ∂f∂A + ϵ = 2 * alg.tol + for i in 1:(alg.maxiter) + g = ∂f∂x(g) + Σₙ = ∂f∂A(g) + dx += Σₙ + ϵnew = norm(Σₙ) # TODO: normalize this error? + Δϵ = ϵ - ϵnew + alg.verbosity > 1 && + @printf("Gradient iter: %3d ‖Σₙ‖: %.2e Δ‖Σₙ‖: %.2e\n", i, ϵnew, Δϵ) + ϵ = ϵnew + + ϵ < alg.tol && break + if alg.verbosity > 0 && i == alg.maxiter + @warn "gradient fixed-point iteration reached maximal number of iterations at ‖Σₙ‖ = $ϵ" + end + end + return dx +end + +function fpgrad(∂F∂x, ∂f∂x, ∂f∂A, y₀, alg::ManualIter) + y = deepcopy(y₀) # Do not mutate y₀ + dx = ∂f∂A(y) + ϵ = 1.0 + for i in 1:(alg.maxiter) + y′ = ∂F∂x + ∂f∂x(y) + + dxnew = ∂f∂A(y′) + ϵnew = norm(dxnew - dx) + Δϵ = ϵ - ϵnew + alg.verbosity > 1 && @printf( + "Gradient iter: %3d ‖Cᵢ₊₁-Cᵢ‖/N: %.2e Δ‖Cᵢ₊₁-Cᵢ‖/N: %.2e\n", i, ϵnew, Δϵ + ) + y = y′ + dx = dxnew + ϵ = ϵnew + + ϵ < alg.tol && break + if alg.verbosity > 0 && i == alg.maxiter + @warn "gradient fixed-point iteration reached maximal number of iterations at ‖Cᵢ₊₁-Cᵢ‖ = $ϵ" + end + end + return dx +end + +function fpgrad(∂F∂x, ∂f∂x, ∂f∂A, y₀, alg::KrylovKit.LinearSolver) + y, info = linsolve(∂f∂x, ∂F∂x, y₀, alg, 1, -1) + if alg.verbosity > 0 && info.converged != 1 + @warn("gradient fixed-point iteration reached maximal number of iterations:", info) + end + + return ∂f∂A(y) +end diff --git a/src/environments/boundarympsenv.jl b/src/environments/boundarympsenv.jl deleted file mode 100644 index 245ce7bc..00000000 --- a/src/environments/boundarympsenv.jl +++ /dev/null @@ -1,242 +0,0 @@ -# Some form of boundary MPS environments for infinite PEPS and PEPO routines - -## Utility - -algtol(alg::VUMPS) = alg.tol_galerkin -algtol(alg::GradientGrassmann) = alg.method.gradtol -update_tol(alg::VUMPS, tol) = @set alg.tol_galerkin = tol -function update_tol(alg::GradientGrassmann, tol) # annoying disparity between typedef and actual constructor... - m = alg.method - m = @set m.gradtol = tol - return GradientGrassmann(; method=m, (finalize!)=alg.finalize!) -end - -## Boundary MPS environment manager - -mutable struct BoundaryMPSEnv{A,E,F} <: Cache - boundaries::A - envs::E - alg::F -end - -## PEPS boundary MPS - -function MPSKit.environments( - peps::InfinitePEPS, - alg::A=VUMPS(); - vspaces=[oneunit(space(peps, 1, 1))], - hermitian=false, - kwargs..., -) where {A<:Union{VUMPS,GradientGrassmann}} - tr = TransferPEPSMultiline(peps, 1) - return environments(tr, alg; vspaces, hermitian, kwargs...) -end - -function MPSKit.recalculate!( - envs::BoundaryMPSEnv, - peps::InfinitePEPS; - tol=algtol(envs.alg), - hermitian=false, - kwargs..., -) - tr = TransferPEPSMultiline(peps, 1) - return recalculate!(envs, tr; tol, hermitian, kwargs...) -end - -# this probably needs a different name? -# computes norm-per-site of PEPS-PEPS overlap using above, below and mixed environments -function MPSKit.expectation_value(peps::InfinitePEPS, ca::BoundaryMPSEnv) - retval = PeriodicArray{scalartype(peps),2}(undef, size(peps)...) - above = ca.boundaries[1] - below = ca.boundaries[2] - (lw, rw) = ca.envs[3] - for (row, col) in product(1:size(peps, 1), 1:size(peps, 2)) - fliprow = size(peps, 1) - row + 1 # below starts counting from below - retval[row, col] = @tensor lw[row, col][1 2 4; 7] * - conj(below.AC[fliprow, col][1 3 6; 13]) * - peps[row, col][5; 8 11 3 2] * - conj(peps[row, col][5; 9 12 6 4]) * - above.AC[row, col][7 8 9; 10] * - rw[row, col][10 11 12; 13] - end - return retval -end - -function MPSKit.normalize!(peps::InfinitePEPS, envs::BoundaryMPSEnv) - norm_per_site = expectation_value(peps, envs) - for (row, col) in product(1:size(peps, 1), 1:size(peps, 2)) - scale!(peps[row, col], 1 / sqrt(norm_per_site[row, col])) - end -end - -# naming convention? this does not return an effective operator, so naming is probably misleading -function ∂∂peps(peps::InfinitePEPS{T}, ca::BoundaryMPSEnv) where {T<:PEPSTensor} - retval = PeriodicArray{T,2}(undef, size(peps)...) - above = ca.boundaries[1] - below = ca.boundaries[2] - (lw, rw) = ca.envs[3] - for (row, col) in product(1:size(peps, 1), 1:size(peps, 2)) - fliprow = size(peps, 1) - row + 1 # below starts counting from below - retval[row, col] = ∂peps( - above.AC[row, col], - below.AC[fliprow, col], - peps[row, col], - lw[row, col], - rw[row, col], - ) - end - return retval -end - -## PEPO boundary MPS - -function MPSKit.environments( - peps::InfinitePEPS, - pepo::InfinitePEPO, - alg::A=VUMPS(); - vspaces=[oneunit(space(peps, 1, 1))], - hermitian=false, - kwargs..., -) where {A<:Union{VUMPS,GradientGrassmann}} - tr = TransferPEPOMultiline(peps, pepo, 1) - return environments(tr, alg; vspaces, hermitian, kwargs...) -end - -function MPSKit.recalculate!( - envs::BoundaryMPSEnv, - peps::InfinitePEPS, - pepo::InfinitePEPO; - tol=algtol(envs.alg), - hermitian=false, - kwargs..., -) - tr = TransferPEPOMultiline(peps, pepo, 1) - return recalculate!(envs, tr; tol, hermitian, kwargs...) -end - -# this probably needs a different name? -# computes norm-per-site of PEPS-PEPO-PEPS sandwich using above, below and mixed environments -function MPSKit.expectation_value( - peps::InfinitePEPS, pepo::InfinitePEPO, ca::BoundaryMPSEnv -) - retval = PeriodicArray{scalartype(peps),2}(undef, size(peps)...) - opp = TransferPEPOMultiline(peps, pepo, 1) # for convenience... - N = height(opp[1]) + 4 - above = ca.boundaries[1] - below = ca.boundaries[2] - (lw, rw) = ca.envs[3] - for (row, col) in product(1:size(peps, 1), 1:size(peps, 2)) - fliprow = size(peps, 1) - row + 1 # below starts counting from below - O_rc = opp[row, col] - GL´ = transfer_left(lw[row, col], O_rc, above.AC[row, col], below.AC[fliprow, col]) - retval[row, col] = TensorKit.TensorOperations.tensorscalar( - ncon([GL´, rw[row, col]], [[N, (2:(N - 1))..., 1], [(1:N)...]]) - ) - end - return retval -end - -# naming convention? this does not return an effective operator, so naming is probably misleading -function ∂∂peps( - peps::InfinitePEPS{T}, pepo::InfinitePEPO, ca::BoundaryMPSEnv -) where {T<:PEPSTensor} - retval = PeriodicArray{T,2}(undef, size(peps)...) - opp = TransferPEPOMultiline(peps, pepo, 1) # just for convenience - above = ca.boundaries[1] - below = ca.boundaries[2] - (lw, rw) = ca.envs[3] - for (row, col) in product(1:size(peps, 1), 1:size(peps, 2)) - fliprow = size(peps, 1) - row + 1 # below starts counting from below - O_rc = opp[row, col] - retval[row, col] = ∂peps( - above.AC[row, col], - below.AC[fliprow, col], - (O_rc[1], O_rc[3]), - lw[row, col], - rw[row, col], - ) - end - return retval -end - -## The actual routines - -# because below starts counting from below -_flippedy(M::MPSMultiline) = MPSKit.Multiline(circshift(reverse(M.data), 1)) - -function MPSKit.environments( - tr::Union{TransferPEPSMultiline,TransferPEPOMultiline}, - alg::A=VUMPS(); - vspaces=[oneunit(space(tr, 1, 1))], # defaults to trivial - hermitian=false, -) where {A<:Union{VUMPS,GradientGrassmann}} - # above boundary - above = initializeMPS(tr, vspaces) - envs_above = environments(above, tr) - - # maybe below boundary - if hermitian # literally the same - below = above - envs_below = envs_above - envs_mixed = (envs_above.lw, envs_above.rw) - else - error("not implemented yet") - tr_dag = dagger(tr) - below = initializeMPS(tr_dag, vspaces) - envs_below = envirlnments(below, tr_dag) - - # mixed environments - envs_mixed = MPSKit.mixed_fixpoints(above, tr, _flippedy(below)) - end - - # collect - boundaries = [above, below] - envs = [envs_above, envs_below, envs_mixed] - - return recalculate!(BoundaryMPSEnv(boundaries, envs, alg), tr; hermitian) -end - -function MPSKit.recalculate!( - envs::BoundaryMPSEnv, - tr::Union{TransferPEPSMultiline,TransferPEPOMultiline}; - tol=algtol(envs.alg), - hermitian=false, -) - alg = algtol(envs.alg) == tol ? envs.alg : update_tol(envs.alg, tol) - - # above boundary - envs.envs[1].opp = tr # because pre-supplied environments only ever use their own operator... - envs.boundaries[1], envs.envs[1], err = leading_boundary( - envs.boundaries[1], tr, alg, envs.envs[1] - ) - - # maybe below boundary - if hermitian # literally the same - envs.boundaries[2], envs.envs[2] = envs.boundaries[1], envs.envs[1] - envs.envs[3] = (envs.envs[1].lw, envs.envs[1].rw) - else - error("not implemented yet") - tr_dag = dagger(tr) # TODO: actually define this... - envs.envs[2].opp = tr_dag # because pre-supplied environments only ever use their own operator... - envs.boundaries[2], envs.envs[2], err = leading_boundary( - envs.boundaries[2], tr_dag, alg, envs.envs[2] - ) - # mixed environments - envs.envs[3] = MPSKit.mixed_fixpoints( - envs.boundaries[1], tr, _flippedy(envs.boundaries[2]), envs.envs[3] - ) - end - - envs.alg = alg - - return envs -end - -## Channel MPS environments - -# TODO: actually try this again at some point? -mutable struct ChannelMPSEnv{A,C,F} <: Cache - boundaries::A - corners::C - alg::F -end diff --git a/src/environments/ctmrgenv.jl b/src/environments/ctmrgenv.jl index a00627e4..fd2bf4c4 100644 --- a/src/environments/ctmrgenv.jl +++ b/src/environments/ctmrgenv.jl @@ -1,52 +1,194 @@ -struct CTMRGEnv{P,C,T} - peps_above::InfinitePEPS{P} - peps_below::InfinitePEPS{P} +""" + struct CTMRGEnv{C,T} + +Corner transfer-matrix environment containing unit-cell arrays of corner and edge tensors. +""" +struct CTMRGEnv{C,T} corners::Array{C,3} edges::Array{T,3} end -# initialize ctmrg environments with some random tensors -CTMRGEnv(peps::InfinitePEPS) = CTMRGEnv(peps, peps); - -function CTMRGEnv(peps_above::InfinitePEPS{P}, peps_below::InfinitePEPS{P}) where {P} - ou = oneunit(space(peps_above, 1, 1)) # the bogus space +""" + CTMRGEnv(peps::InfinitePEPS{P}; Venv=oneunit(spacetype(P))) +Create a random CTMRG environment from a PEPS tensor. The environment bond dimension +defaults to one and can be specified using the `Venv` space. +""" +function CTMRGEnv(peps::InfinitePEPS{P}; Venv=oneunit(spacetype(P))) where {P} C_type = tensormaptype(spacetype(P), 1, 1, storagetype(P)) - T_type = tensormaptype(spacetype(P), 3, 1, storagetype(P)) # debatable how we should do the legs? + T_type = tensormaptype(spacetype(P), 3, 1, storagetype(P)) - #first index is de - corners = Array{C_type}(undef, 4, size(peps_above)...) - edges = Array{T_type}(undef, 4, size(peps_above)...) + # First index is direction + corners = Array{C_type}(undef, 4, size(peps)...) + edges = Array{T_type}(undef, 4, size(peps)...) - for dir in 1:4, i in 1:size(peps_above, 1), j in 1:size(peps_above, 2) - @diffset corners[dir, i, j] = TensorMap(randn, scalartype(P), ou, ou) + for dir in 1:4, i in 1:size(peps, 1), j in 1:size(peps, 2) + @diffset corners[dir, i, j] = TensorMap(randn, scalartype(P), Venv, Venv) @diffset edges[dir, i, j] = TensorMap( randn, scalartype(P), - ou * space(peps_above[i, j], dir + 1)' * space(peps_below[i, j], dir + 1), - ou, + Venv * space(peps[i, j], dir + 1)' * space(peps[i, j], dir + 1), + Venv, ) end @diffset corners[:, :, :] ./= norm.(corners[:, :, :]) @diffset edges[:, :, :] ./= norm.(edges[:, :, :]) - return CTMRGEnv(peps_above, peps_below, corners, edges) + return CTMRGEnv(corners, edges) +end + +# Custom adjoint for CTMRGEnv constructor, needed for fixed-point differentiation +function ChainRulesCore.rrule(::Type{CTMRGEnv}, corners, edges) + ctmrgenv_pullback(ē) = NoTangent(), ē.corners, ē.edges + return CTMRGEnv(corners, edges), ctmrgenv_pullback end -function Base.rotl90(envs::CTMRGEnv{P,C,T}) where {P,C,T} - n_peps_above = rotl90(envs.peps_above) - n_peps_below = rotl90(envs.peps_below) - n_corners = Array{C,3}(undef, size(envs.corners)...) - n_edges = Array{T,3}(undef, size(envs.edges)...) +# Custom adjoint for CTMRGEnv getproperty, to avoid creating named tuples in backward pass +function ChainRulesCore.rrule(::typeof(getproperty), e::CTMRGEnv, name::Symbol) + result = getproperty(e, name) + if name === :corners + function corner_pullback(Δcorners) + return NoTangent(), CTMRGEnv(Δcorners, zerovector.(e.edges)), NoTangent() + end + return result, corner_pullback + elseif name === :edges + function edge_pullback(Δedges) + return NoTangent(), CTMRGEnv(zerovector.(e.corners), Δedges), NoTangent() + end + return result, edge_pullback + else + # this should never happen because already errored in forwards pass + throw(ArgumentError("No rrule for getproperty of $name")) + end +end + +# Rotate corners & edges counter-clockwise +function Base.rotl90(env::CTMRGEnv{C,T}) where {C,T} + # Initialize rotated corners & edges with rotated sizes + corners′ = Zygote.Buffer( + Array{C,3}(undef, 4, size(env.corners, 3), size(env.corners, 2)) + ) + edges′ = Zygote.Buffer(Array{T,3}(undef, 4, size(env.edges, 3), size(env.edges, 2))) for dir in 1:4 - dirm = mod1(dir - 1, 4) - @diffset n_corners[dirm, :, :] .= rotl90(envs.corners[dir, :, :]) - @diffset n_edges[dirm, :, :] .= rotl90(envs.edges[dir, :, :]) + corners′[_prev(dir, 4), :, :] = rotl90(env.corners[dir, :, :]) + edges′[_prev(dir, 4), :, :] = rotl90(env.edges[dir, :, :]) end - return CTMRGEnv(n_peps_above, n_peps_below, n_corners, n_edges) + return CTMRGEnv(copy(corners′), copy(edges′)) +end + +Base.eltype(env::CTMRGEnv) = eltype(env.corners[1]) + +# In-place update of environment +function update!(env::CTMRGEnv{C,T}, env´::CTMRGEnv{C,T}) where {C,T} + env.corners .= env´.corners + env.edges .= env´.edges + return env +end + +# Functions used for FP differentiation and by KrylovKit.linsolve +function Base.:+(e₁::CTMRGEnv, e₂::CTMRGEnv) + return CTMRGEnv(e₁.corners + e₂.corners, e₁.edges + e₂.edges) +end + +function Base.:-(e₁::CTMRGEnv, e₂::CTMRGEnv) + return CTMRGEnv(e₁.corners - e₂.corners, e₁.edges - e₂.edges) +end + +Base.:*(α::Number, e::CTMRGEnv) = CTMRGEnv(α * e.corners, α * e.edges) +Base.similar(e::CTMRGEnv) = CTMRGEnv(similar(e.corners), similar(e.edges)) + +function LinearAlgebra.mul!(edst::CTMRGEnv, esrc::CTMRGEnv, α::Number) + edst.corners .= α * esrc.corners + edst.edges .= α * esrc.edges + return edst +end + +function LinearAlgebra.rmul!(e::CTMRGEnv, α::Number) + rmul!(e.corners, α) + rmul!(e.edges, α) + return e +end + +function LinearAlgebra.axpy!(α::Number, e₁::CTMRGEnv, e₂::CTMRGEnv) + e₂.corners .+= α * e₁.corners + e₂.edges .+= α * e₁.edges + return e₂ end -Base.eltype(envs::CTMRGEnv) = eltype(envs.corners[1]) +function LinearAlgebra.axpby!(α::Number, e₁::CTMRGEnv, β::Number, e₂::CTMRGEnv) + e₂.corners .= α * e₁.corners + β * e₂.corners + e₂.edges .= α * e₁.edges + β * e₂.edges + return e₂ +end + +function LinearAlgebra.dot(e₁::CTMRGEnv, e₂::CTMRGEnv) + return dot(e₁.corners, e₂.corners) + dot(e₁.edges, e₂.edges) +end + +# VectorInterface +# --------------- + +# Note: the following methods consider the environment tensors as separate components of one +# big vector. In other words, the associated vector space is not the natural one associated +# to the original (physical) system, and addition, scaling, etc. are performed element-wise. + +import VectorInterface as VI + +function VI.scalartype(::Type{CTMRGEnv{C,T}}) where {C,T} + S₁ = scalartype(C) + S₂ = scalartype(T) + return promote_type(S₁, S₂) +end + +function VI.zerovector(env::CTMRGEnv, ::Type{S}) where {S<:Number} + _zerovector = Base.Fix2(zerovector, S) + return CTMRGEnv(map(_zerovector, env.corners), map(_zerovector, env.edges)) +end +function VI.zerovector!(env::CTMRGEnv) + foreach(zerovector!, env.corners) + foreach(zerovector!, env.edges) + return env +end +VI.zerovector!!(env::CTMRGEnv) = zerovector!(env) + +function VI.scale(env::CTMRGEnv, α::Number) + _scale = Base.Fix2(scale, α) + return CTMRGEnv(map(_scale, env.corners), map(_scale, env.edges)) +end +function VI.scale!(env::CTMRGEnv, α::Number) + _scale! = Base.Fix2(scale!, α) + foreach(_scale!, env.corners) + foreach(_scale!, env.edges) + return env +end +function VI.scale!(env₁::CTMRGEnv, env₂::CTMRGEnv, α::Number) + _scale!(x, y) = scale!(x, y, α) + foreach(_scale!, env₁.corners, env₂.corners) + foreach(_scale!, env₁.edges, env₂.edges) + return env₁ +end +VI.scale!!(env::CTMRGEnv, α::Number) = scale!(env, α) +VI.scale!!(env₁::CTMRGEnv, env₂::CTMRGEnv, α::Number) = scale!(env₁, env₂, α) + +function VI.add(env₁::CTMRGEnv, env₂::CTMRGEnv, α::Number, β::Number) + _add(x, y) = add(x, y, α, β) + return CTMRGEnv( + map(_add, env₁.corners, env₂.corners), map(_add, env₁.edges, env₂.edges) + ) +end +function VI.add!(env₁::CTMRGEnv, env₂::CTMRGEnv, α::Number, β::Number) + _add!(x, y) = add!(x, y, α, β) + foreach(_add!, env₁.corners, env₂.corners) + foreach(_add!, env₁.edges, env₂.edges) + return env₁ +end +VI.add!!(env₁::CTMRGEnv, env₂::CTMRGEnv, α::Number, β::Number) = add!(env₁, env₂, α, β) + +# Exploiting the fact that VectorInterface works for tuples: +function VI.inner(env₁::CTMRGEnv, env₂::CTMRGEnv) + return inner((env₁.corners, env₁.edges), (env₂.corners, env₂.edges)) +end +VI.norm(env::CTMRGEnv) = norm((env.corners, env.edges)) diff --git a/src/operators/infinitepepo.jl b/src/operators/infinitepepo.jl index c8ffc9f0..d8e24e00 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) + InfinitePEPO(A; unitcell=(1, 1, 1)) -Allow users to pass in single tensor. +Create an InfinitePEPO by specifying a tensor and unit cell. """ -function InfinitePEPO(A::T) where {T<:PEPOTensor} - As = Array{T,3}(undef, (1, 1, 1)) - As[1, 1, 1] = A - return InfinitePEPO(As) +function InfinitePEPO(A::T; unitcell::Tuple{Int,Int,Int}=(1, 1, 1)) where {T<:PEPOTensor} + return InfinitePEPO(fill(A, unitcell)) 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(f=randn, T=ComplexF64, Pspace, Nspace, [Espace]; unitcell=(1,1,1)) +Create an InfinitePEPO by specifying its spaces and unit cell. """ - InfinitePEPO(d, D) - -Allow users to pass in integers. -""" -function InfinitePEPO(d::Integer, D::Integer) - T = TensorMap(rand, ComplexF64, ℂ^d ⊗ (ℂ^d)' ← ℂ^D ⊗ ℂ^D ⊗ (ℂ^D)' ⊗ (ℂ^D)') - return InfinitePEPO(T) -end - -""" -InfinitePEPO(d, D, L) -InfinitePEPO(d, D, (Lx, Ly, Lz))) - -Allow users to pass in integers and specify 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 @@ -142,8 +121,6 @@ Base.getindex(T::InfinitePEPO, args...) = Base.getindex(T.A, args...) Base.axes(T::InfinitePEPO, args...) = axes(T.A, args...) TensorKit.space(T::InfinitePEPO, i, j) = space(T[i, j, end], 1) -Base.rotl90(T::InfinitePEPO) = InfinitePEPO(rotl90(rotl90.(T.A))); - function initializePEPS( T::InfinitePEPO{<:PEPOTensor{S}}, vspace::S ) where {S<:ElementarySpace} diff --git a/src/operators/localoperator.jl b/src/operators/localoperator.jl new file mode 100644 index 00000000..386869f8 --- /dev/null +++ b/src/operators/localoperator.jl @@ -0,0 +1,219 @@ +# TODO: change this implementation to a type-stable one + +abstract type AbstractInteraction end + +""" + struct OnSite <: AbstractInteraction + +Trivial interaction representing terms that act on one isolated site. +""" +struct OnSite <: AbstractInteraction end + +""" + struct NearestNeighbor <: AbstractInteraction + +Interaction representing nearest neighbor terms that act on two adjacent sites. +""" +struct NearestNeighbor <: AbstractInteraction end + +""" + struct NLocalOperator{I<:AbstractInteraction} + +Operator in form of a `AbstractTensorMap` which is parametrized by an interaction type. +Mostly, this is used to define Hamiltonian terms and observables. +""" +struct NLocalOperator{I<:AbstractInteraction} + op::AbstractTensorMap +end + +""" + struct AnisotropicNNOperator{I<:AbstractInteraction} + +Operator which includes an on-site term and two nearest-neighbor terms, vertical and horizontal. +""" +struct AnisotropicNNOperator + h0::NLocalOperator{OnSite} + hx::NLocalOperator{NearestNeighbor} + hy::NLocalOperator{NearestNeighbor} +end +function AnisotropicNNOperator( + h0::AbstractTensorMap{S,1,1}, + hx::AbstractTensorMap{S,2,2}, + hy::AbstractTensorMap{S,2,2}=hx, +) where {S} + return AnisotropicNNOperator( + NLocalOperator{OnSite}(h0), + NLocalOperator{NearestNeighbor}(hx), + NLocalOperator{NearestNeighbor}(hy), + ) +end +# TODO: include the on-site term in the two-site terms, to reduce number of contractions. + +@doc """ + operator_env(peps::InfinitePEPS, env::CTMRGEnv, ::AbstractInteraction) + +Contract a PEPS and a CTMRG environment to form an operator environment. +The open bonds correspond to the indices of an operator with the specified +`AbstractInteraction` type. +""" +operator_env + +function operator_env(peps::InfinitePEPS, env::CTMRGEnv, ::OnSite) + return map(Iterators.product(axes(env.corners, 2), axes(env.corners, 3))) do (r, c) + @tensor opt = true ρ[-1; -2] := + env.corners[NORTHWEST, r, c][1; 2] * + env.edges[NORTH, r, c][2 3 4; 5] * + env.corners[NORTHEAST, r, c][5; 6] * + env.edges[EAST, r, c][6 7 8; 9] * + env.corners[SOUTHEAST, r, c][9; 10] * + env.edges[SOUTH, r, c][10 11 12; 13] * + env.corners[SOUTHWEST, r, c][13; 14] * + env.edges[WEST, r, c][14 15 16; 1] * + peps[r, c][-1; 3 7 11 15] * + conj(peps[r, c][-2; 4 8 12 16]) + end +end + +function operator_env(peps::InfinitePEPS, env::CTMRGEnv, ::NearestNeighbor) + return map(Iterators.product(axes(env.corners, 2), axes(env.corners, 3))) do (r, c) + cnext = _next(c, size(peps, 2)) + @tensor opt = true ρ[-12 -18; -11 -20] := + env.corners[NORTHWEST, r, c][1; 3] * + env.edges[NORTH, r, c][3 5 8; 13] * + env.edges[NORTH, r, cnext][13 16 22; 23] * + env.corners[NORTHEAST, r, cnext][23; 24] * + env.edges[EAST, r, cnext][24 25 26; 27] * + env.corners[SOUTHEAST, r, cnext][27; 28] * + env.edges[SOUTH, r, cnext][28 17 21; 14] * + env.edges[SOUTH, r, c][14 6 10; 4] * + env.corners[SOUTHWEST, r, c][4; 2] * + env.edges[WEST, r, c][2 7 9; 1] * + peps[r, c][-12; 5 15 6 7] * + conj(peps[r, c][-11; 8 19 10 9]) * + peps[r, cnext][-18; 16 25 17 15] * + conj(peps[r, cnext][-20; 22 26 21 19]) + end +end + +@doc """ + MPSKit.expectation_value(peps::InfinitePEPS, env, O::NLocalOperator) + +Evaluate the expectation value of any `NLocalOperator` on each unit-cell entry +of `peps` and `env`. +""" MPSKit.expectation_value(::InfinitePEPS, ::Any, ::NLocalOperator) + +# Optimal contraction order is obtained by manually trying out some space sizes and using costcheck = warn +# in principle, we would like to use opt = true, but this does not give optimal results without also supplying costs +# However, due to a bug in tensoroperations this is currently not possible with integer labels. +# Thus, this is a workaround until the bug is fixed. (otherwise we'd need to rewrite all the labels to symbols...) + +# 1-site operator expectation values on unit cell +function MPSKit.expectation_value( + peps::InfinitePEPS, env::CTMRGEnv, O::NLocalOperator{OnSite} +) + return map(Iterators.product(axes(env.corners, 2), axes(env.corners, 3))) do (r, c) + o = @tensor order = (6, 2, 5, 10, 14, 13, 11, 15, 7, 9, 1, 3, 4, 8, 12, 16, 18, 17) begin + env.corners[NORTHWEST, r, c][1; 2] * + env.edges[NORTH, r, c][2 3 4; 5] * + env.corners[NORTHEAST, r, c][5; 6] * + env.edges[EAST, r, c][6 7 8; 9] * + env.corners[SOUTHEAST, r, c][9; 10] * + env.edges[SOUTH, r, c][10 11 12; 13] * + env.corners[SOUTHWEST, r, c][13; 14] * + env.edges[WEST, r, c][14 15 16; 1] * + peps[r, c][17; 3 7 11 15] * + conj(peps[r, c][18; 4 8 12 16]) * + O.op[18; 17] + end + n = @tensor order = (9, 13, 10, 5, 1, 2, 4, 16, 6, 8, 14, 12, 17, 3, 7, 11, 15) begin + env.corners[NORTHWEST, r, c][1; 2] * + env.edges[NORTH, r, c][2 3 4; 5] * + env.corners[NORTHEAST, r, c][5; 6] * + env.edges[EAST, r, c][6 7 8; 9] * + env.corners[SOUTHEAST, r, c][9; 10] * + env.edges[SOUTH, r, c][10 11 12; 13] * + env.corners[SOUTHWEST, r, c][13; 14] * + env.edges[WEST, r, c][14 15 16; 1] * + peps[r, c][17; 3 7 11 15] * + conj(peps[r, c][17; 4 8 12 16]) + end + o / n + end +end + +#! format: off +function MPSKit.expectation_value( + peps::InfinitePEPS, env, O::NLocalOperator{NearestNeighbor} +) + return map(Iterators.product(axes(env.corners, 2), axes(env.corners, 3))) do (r, c) + cnext = _next(c, size(peps, 2)) + o = @tensor order = ( + 28, 24, 23, 16, 25, 22, 26, 27, 17, 21, 4, 1, 3, 5, 7, 8, 9, 2, 6, 10, 14, 19, + 15, 13, 31, 32, 29, 30, + ) begin # physical spaces + env.corners[NORTHWEST, r, c][1; 3] * + env.edges[NORTH, r, c][3 5 8; 13] * + env.edges[NORTH, r, cnext][13 16 22; 23] * + env.corners[NORTHEAST, r, cnext][23; 24] * + env.edges[EAST, r, cnext][24 25 26; 27] * + env.corners[SOUTHEAST, r, cnext][27; 28] * + env.edges[SOUTH, r, cnext][28 17 21; 14] * + env.edges[SOUTH, r, c][14 6 10; 4] * + env.corners[SOUTHWEST, r, c][4; 2] * + env.edges[WEST, r, c][2 7 9; 1] * + peps[r, c][29; 5 15 6 7] * + conj(peps[r, c][31; 8 19 10 9]) * + peps[r, cnext][30; 16 25 17 15] * + conj(peps[r, cnext][32; 22 26 21 19]) * + O.op[31 32; 29 30] + end + + n = @tensor order = ( + 2, 3, 1, 5, 7, 28, 24, 23, 16, 25, 30, 22, 26, 27, 17, 21, 14, 15, 6, 4, 13, 29, + 8, 19, 10, 9, + ) begin + env.corners[NORTHWEST, r, c][1; 3] * + env.edges[NORTH, r, c][3 5 8; 13] * + env.edges[NORTH, r, cnext][13 16 22; 23] * + env.corners[NORTHEAST, r, cnext][23; 24] * + env.edges[EAST, r, cnext][24 25 26; 27] * + env.corners[SOUTHEAST, r, cnext][27; 28] * + env.edges[SOUTH, r, cnext][28 17 21; 14] * + env.edges[SOUTH, r, c][14 6 10; 4] * + env.corners[SOUTHWEST, r, c][4; 2] * + env.edges[WEST, r, c][2 7 9; 1] * + peps[r, c][29; 5 15 6 7] * + conj(peps[r, c][29; 8 19 10 9]) * + peps[r, cnext][30; 16 25 17 15] * + conj(peps[r, cnext][30; 22 26 21 19]) + end + o / n + end +end +#! format: on + +""" + costfun(peps::InfinitePEPS, env, op::NLocalOperator{NearestNeighbor}) + +Compute the expectation value of a nearest-neighbor operator. +This is used to evaluate and differentiate the energy in ground-state PEPS optimizations. +""" +function costfun(peps::InfinitePEPS, env, op::NLocalOperator{NearestNeighbor}) + oh = sum(expectation_value(peps, env, op)) + ov = sum(expectation_value(rotl90(peps), rotl90(env), op)) + return real(oh + ov) +end + +""" + costfun(peps::InfinitePEPS, env, op::AnisotropicNNOperator) + +Compute the expectation value of an on-site and an anisotropic nearest-neighbor operator. +This is used to evaluate and differentiate the energy in ground-state PEPS optimizations. +""" +function costfun(peps::InfinitePEPS, env, op::AnisotropicNNOperator) + oos = sum(expectation_value(peps, env, op.h0)) + oh = sum(expectation_value(peps, env, op.hx)) + ov = sum(expectation_value(rotr90(peps), rotate_north(env, WEST), op.hy)) + #ov = sum(expectation_value(rotl90(peps), rotl90(env), op.hy)) + return real(oos + oh + ov) +end diff --git a/src/operators/models.jl b/src/operators/models.jl new file mode 100644 index 00000000..0c36d7ad --- /dev/null +++ b/src/operators/models.jl @@ -0,0 +1,45 @@ +## Model Hamiltonians +# ------------------- + +""" + square_lattice_heisenberg(; Jx=-1, Jy=1, Jz=-1) + +Square lattice Heisenberg model. +By default, this implements a single site unit cell via a sublattice rotation. +""" +function square_lattice_heisenberg( + ::Type{T}=ComplexF64; Jx=-1, Jy=1, Jz=-1 +) where {T<:Number} + physical_space = ComplexSpace(2) + σx = TensorMap(T[0 1; 1 0], physical_space, physical_space) + σy = TensorMap(T[0 im; -im 0], physical_space, physical_space) + σz = TensorMap(T[1 0; 0 -1], physical_space, physical_space) + H = (Jx * σx ⊗ σx) + (Jy * σy ⊗ σy) + (Jz * σz ⊗ σz) + return NLocalOperator{NearestNeighbor}(H / 4) +end + +""" + square_lattice_pwave(; t=1, μ=2, Δ=1) + +Square lattice p-wave superconductor model. +""" +function square_lattice_pwave( + ::Type{T}=ComplexF64; t::Number=1, μ::Number=2, Δ::Number=1 +) where {T<:Number} + physical_space = Vect[FermionParity](0 => 1, 1 => 1) + # on-site + h0 = TensorMap(zeros, T, physical_space ← physical_space) + block(h0, FermionParity(1)) .= -μ + + # two-site (x-direction) + hx = TensorMap(zeros, T, physical_space^2 ← physical_space^2) + block(hx, FermionParity(0)) .= [0 -Δ; -Δ 0] + block(hx, FermionParity(1)) .= [0 -t; -t 0] + + # two-site (y-direction) + hy = TensorMap(zeros, T, physical_space^2 ← physical_space^2) + block(hy, FermionParity(0)) .= [0 Δ*im; -Δ*im 0] + block(hy, FermionParity(1)) .= [0 -t; -t 0] + + return AnisotropicNNOperator(h0, hx, hy) +end diff --git a/src/operators/transferpepo.jl b/src/operators/transferpepo.jl index a05143c7..eae23fd4 100644 --- a/src/operators/transferpepo.jl +++ b/src/operators/transferpepo.jl @@ -1,11 +1,26 @@ +""" + 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} - bot::PeriodicArray{T,1} # this is assumed to be conjed by default, but do we want this? + bot::PeriodicArray{T,1} 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) @@ -15,7 +30,7 @@ end Base.size(transfer::InfiniteTransferPEPO) = size(transfer.top) Base.size(transfer::InfiniteTransferPEPO, args...) = size(transfer.top, args...) Base.length(transfer::InfiniteTransferPEPO) = size(transfer, 1) -height(transfer::InfiniteTransferPEPO) = size(transfer.mid, 2) # will I ever need this? +height(transfer::InfiniteTransferPEPO) = size(transfer.mid, 2) Base.getindex(O::InfiniteTransferPEPO, i) = (O.top[i], O.bot[i], Tuple(O.mid[i, :])) # TODO: not too sure about this Base.iterate(O::InfiniteTransferPEPO, i=1) = i > length(O) ? nothing : (O[i], i + 1) @@ -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 @@ -77,7 +105,7 @@ function MPSKit.transfer_left( A[1 9 6 3; -5] end -# it actually works, somehow +# general case function MPSKit.transfer_left( GL::GenericMPSTensor{S,N}, O::Tuple{T,T,Tuple{Vararg{P,H}}}, @@ -134,6 +162,7 @@ function MPSKit.transfer_right( A[-1 11 12 13; 10] end +# general case function MPSKit.transfer_right( GR::GenericMPSTensor{S,N}, O::Tuple{T,T,Tuple{Vararg{P,H}}}, @@ -194,7 +223,7 @@ function MPSKit.expectation_value( for (i, j) in product(1:size(st, 1), 1:size(st, 2)) O_ij = opp[i, j] N = height(opp[1]) + 4 - # just reuse left environment contraction? + # just reuse left environment contraction GL´ = transfer_left(leftenv(ca, i, j, st), O_ij, st.AC[i, j], st.AC[i + 1, j]) retval[i, j] = TensorKit.TensorOperations.tensorscalar( ncon([GL´, rightenv(ca, i, j, st)], [[N, (2:(N - 1))..., 1], [(1:N)...]]) diff --git a/src/operators/transferpeps.jl b/src/operators/transferpeps.jl index 05859bbd..f1597647 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,24 +88,15 @@ 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) end function initializeMPS(O::MPSKit.Multiline, virtualspaces::AbstractArray{S,1}) where {S} - return initializeMPS(O, repeat(virtualspaces', length(O), 1)) + return initializeMPS(O, repeat(virtualspaces, length(O), 1)) +end +function initializeMPS(O::MPSKit.Multiline, V::ElementarySpace) + return initializeMPS(O, repeat([V], length(O), length(O[1]))) end function initializeMPS(O::MPSKit.Multiline, χ::Int) return initializeMPS(O, repeat([ℂ^χ], length(O), length(O[1]))) @@ -90,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) @@ -119,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/src/states/abstractpeps.jl b/src/states/abstractpeps.jl index 490f4695..c324e711 100644 --- a/src/states/abstractpeps.jl +++ b/src/states/abstractpeps.jl @@ -1,9 +1,3 @@ -# abstractpeps.jl - -########################### -## Abstract tensor types ## -########################### - """ const PEPSTensor{S} @@ -12,6 +6,16 @@ conventionally ordered as: T : P ← N ⊗ E ⊗ S ⊗ W. """ const PEPSTensor{S} = AbstractTensorMap{S,1,4} where {S<:ElementarySpace} +""" + PEPSTensor(f, ::Type{T}, Pspace::S, Nspace::S, + [Espace::S], [Sspace::S], [Wspace::S]) where {T,S<:ElementarySpace} + PEPSTensor(f, ::Type{T}, Pspace::Int, Nspace::Int, + [Espace::Int], [Sspace::Int], [Wspace::Int]) where {T} + +Construct a PEPS tensor based on the physical, north, east, west and south spaces. +Alternatively, only the space dimensions can be provided and ℂ is assumed as the field. +The tensor elements are generated based on `f` and the element type is specified in `T`. +""" function PEPSTensor( f, ::Type{T}, @@ -39,27 +43,20 @@ end const PEPOTensor{S} Default type for PEPO tensors with a single incoming and outgoing physical index, and 4 - virtual indices, conventionally ordered as: O : P ⊗ P' ← N ⊗ E ⊗ S ⊗ W. +virtual indices, conventionally ordered as: O : P ⊗ P' ← N ⊗ E ⊗ S ⊗ W. """ const PEPOTensor{S} = AbstractTensorMap{S,2,4} where {S<:ElementarySpace} -########################## -## Abstract state types ## -########################## - """ abstract type AbstractPEPS end -Abstract supertype for a 2D projected entangled pairs state. +Abstract supertype for a 2D projected entangled-pair state. """ abstract type AbstractPEPS end """ abstract type AbstractPEPO end -Abstract supertype for a 2D projected entangled pairs operator. +Abstract supertype for a 2D projected entangled-pair operator. """ abstract type AbstractPEPO end - -Base.rotl90(t::PEPSTensor) = permute(t, ((1,), (3, 4, 5, 2))) -Base.rotl90(t::PEPOTensor) = permute(t, ((1, 2), (4, 5, 6, 3))); diff --git a/src/states/infinitepeps.jl b/src/states/infinitepeps.jl index dc717656..48285bc7 100644 --- a/src/states/infinitepeps.jl +++ b/src/states/infinitepeps.jl @@ -1,15 +1,11 @@ -# not everything is a PeriodicArray anymore -_next(i, total) = mod1(i + 1, total) -_prev(i, total) = mod1(i - 1, total) - """ struct InfinitePEPS{T<:PEPSTensor} Represents an infinite projected entangled-pair state on a 2D square lattice. """ struct InfinitePEPS{T<:PEPSTensor} <: AbstractPEPS - A::Array{T,2} # TODO: switch back to PeriodicArray? - + A::Matrix{T} + InfinitePEPS{T}(A::Matrix{T}) where {T<:PEPSTensor} = new{T}(A) function InfinitePEPS(A::Array{T,2}) where {T<:PEPSTensor} for (d, w) in Tuple.(CartesianIndices(A)) space(A[d, w], 2) == space(A[_prev(d, end), w], 4)' || throw( @@ -17,6 +13,7 @@ struct InfinitePEPS{T<:PEPSTensor} <: AbstractPEPS ) space(A[d, w], 3) == space(A[d, _next(w, end)], 5)' || throw(SpaceMismatch("East virtual space at site $((d, w)) does not match.")) + dim(space(A[d, w])) > 0 || @warn "no fusion channels at site ($d, $w)" end return new{T}(A) end @@ -33,15 +30,18 @@ function InfinitePEPS(A::AbstractArray{T,2}) where {T<:PEPSTensor} end """ - InfinitePEPS(Pspaces, Nspaces, Espaces) + InfinitePEPS(f=randn, T=ComplexF64, Pspaces, Nspaces, Espaces) Allow users to pass in arrays of spaces. """ function InfinitePEPS( - Pspaces::AbstractArray{S,2}, - Nspaces::AbstractArray{S,2}, - Espaces::AbstractArray{S,2}=Nspaces, -) where {S<:ElementarySpace} + Pspaces::A, Nspaces::A, Espaces::A +) where {A<:AbstractMatrix{<:Union{Int,ElementarySpace}}} + return InfinitePEPS(randn, ComplexF64, Pspaces, Nspaces, Espaces) +end +function InfinitePEPS( + f, T, Pspaces::M, Nspaces::M, Espaces::M=Nspaces +) where {M<:AbstractMatrix{<:Union{Int,ElementarySpace}}} size(Pspaces) == size(Nspaces) == size(Espaces) || throw(ArgumentError("Input spaces should have equal sizes.")) @@ -49,7 +49,7 @@ function InfinitePEPS( Wspaces = adjoint.(circshift(Espaces, (0, -1))) A = map(Pspaces, Nspaces, Espaces, Sspaces, Wspaces) do P, N, E, S, W - return PEPSTensor(randn, ComplexF64, P, N, E, S, W) + return PEPSTensor(f, T, P, N, E, S, W) end return InfinitePEPS(A) @@ -65,7 +65,7 @@ function InfinitePEPS(A::T; unitcell::Tuple{Int,Int}=(1, 1)) where {T<:PEPSTenso end """ - InfinitePEPS(Pspace, Nspace, [Espace]; unitcell=(1,1)) + InfinitePEPS(f=randn, T=ComplexF64, Pspace, Nspace, [Espace]; unitcell=(1,1)) Create an InfinitePEPS by specifying its spaces and unit cell. Spaces can be specified either via `Int` or via `ElementarySpace`. @@ -74,7 +74,18 @@ function InfinitePEPS( Pspace::S, Nspace::S, Espace::S=Nspace; unitcell::Tuple{Int,Int}=(1, 1) ) where {S<:Union{ElementarySpace,Int}} return InfinitePEPS( - fill(Pspace, unitcell), fill(Nspace, unitcell), fill(Espace, unitcell) + randn, + ComplexF64, + fill(Pspace, unitcell), + fill(Nspace, unitcell), + fill(Espace, unitcell), + ) +end +function InfinitePEPS( + f, T, Pspace::S, Nspace::S, Espace::S=Nspace; unitcell::Tuple{Int,Int}=(1, 1) +) where {S<:Union{ElementarySpace,Int}} + return InfinitePEPS( + f, T, fill(Pspace, unitcell), fill(Nspace, unitcell), fill(Espace, unitcell) ) end @@ -87,7 +98,7 @@ VectorInterface.scalartype(T::InfinitePEPS) = scalartype(T.A) ## Copy Base.copy(T::InfinitePEPS) = InfinitePEPS(copy(T.A)) -Base.similar(T::InfinitePEPS) = InfinitePEPS(similar(T.A)) +# Base.similar(T::InfinitePEPS) = InfinitePEPS(similar(T.A)) # TODO: This is incompatible with inner constructor Base.repeat(T::InfinitePEPS, counts...) = InfinitePEPS(repeat(T.A, counts...)) Base.getindex(T::InfinitePEPS, args...) = Base.getindex(T.A, args...) @@ -95,4 +106,62 @@ Base.setindex!(T::InfinitePEPS, args...) = (Base.setindex!(T.A, args...); T) Base.axes(T::InfinitePEPS, args...) = axes(T.A, args...) TensorKit.space(t::InfinitePEPS, i, j) = space(t[i, j], 1) -Base.rotl90(t::InfinitePEPS) = InfinitePEPS(rotl90(rotl90.(t.A))); +## Math +Base.:+(ψ₁::InfinitePEPS, ψ₂::InfinitePEPS) = InfinitePEPS(ψ₁.A + ψ₂.A) +Base.:-(ψ₁::InfinitePEPS, ψ₂::InfinitePEPS) = InfinitePEPS(ψ₁.A - ψ₂.A) +Base.:*(α::Number, ψ::InfinitePEPS) = InfinitePEPS(α * ψ.A) +LinearAlgebra.dot(ψ₁::InfinitePEPS, ψ₂::InfinitePEPS) = dot(ψ₁.A, ψ₂.A) +LinearAlgebra.norm(ψ::InfinitePEPS) = norm(ψ.A) + +# Used in _scale during OptimKit.optimize +function LinearAlgebra.rmul!(ψ::InfinitePEPS, α::Number) + rmul!(ψ.A, α) + return ψ +end + +# Used in _add during OptimKit.optimize +function LinearAlgebra.axpy!(α::Number, ψ₁::InfinitePEPS, ψ₂::InfinitePEPS) + ψ₂.A .+= α * ψ₁.A + return ψ₂ +end + +# VectorInterface +VectorInterface.zerovector(x::InfinitePEPS) = InfinitePEPS(zerovector(x.A)) + +# Chainrules +function ChainRulesCore.rrule( + ::typeof(Base.getindex), state::InfinitePEPS, row::Int, col::Int +) + pepstensor = state[row, col] + + function getindex_pullback(Δpepstensor) + Δstate = zerovector(state) + Δstate[row, col] = Δpepstensor + return NoTangent(), Δstate, NoTangent(), NoTangent() + end + return pepstensor, getindex_pullback +end + +function ChainRulesCore.rrule(::Type{<:InfinitePEPS}, A::Matrix{T}) where {T<:PEPSTensor} + peps = InfinitePEPS(A) + function InfinitePEPS_pullback(Δpeps) + return NoTangent(), Δpeps.A + end + return peps, InfinitePEPS_pullback +end + +function ChainRulesCore.rrule(::typeof(rotl90), peps::InfinitePEPS) + peps′ = rotl90(peps) + function rotl90_pullback(Δpeps) + return NoTangent(), rotr90(Δpeps) + end + return peps′, rotl90_pullback +end + +function ChainRulesCore.rrule(::typeof(rotr90), peps::InfinitePEPS) + peps′ = rotr90(peps) + function rotr90_pullback(Δpeps) + return NoTangent(), rotl90(Δpeps) + end + return peps′, rotr90_pullback +end diff --git a/src/utility/eigsolve.jl b/src/utility/eigsolve.jl new file mode 100644 index 00000000..975adf49 --- /dev/null +++ b/src/utility/eigsolve.jl @@ -0,0 +1,251 @@ +# Copied from Jutho/KrylovKit.jl/pull/56, with minor tweaks + +function ChainRulesCore.rrule( + ::typeof(eigsolve), A::AbstractMatrix, x₀, howmany, which, alg +) + (vals, vecs, info) = eigsolve(A, x₀, howmany, which, alg) + project_A = ProjectTo(A) + T = eltype(vecs[1]) # will be real for real symmetric problems and complex otherwise + + function eigsolve_pullback(ΔX) + _Δvals = unthunk(ΔX[1]) + _Δvecs = unthunk(ΔX[2]) + + ∂self = NoTangent() + ∂x₀ = ZeroTangent() + ∂howmany = NoTangent() + ∂which = NoTangent() + ∂alg = NoTangent() + if _Δvals isa AbstractZero && _Δvecs isa AbstractZero + ∂A = ZeroTangent() + return ∂self, ∂A, ∂x₀, ∂howmany, ∂which, ∂alg + end + + if _Δvals isa AbstractZero + Δvals = fill(NoTangent(), length(Δvecs)) + else + Δvals = _Δvals + end + if _Δvecs isa AbstractZero + Δvecs = fill(NoTangent(), length(Δvals)) + else + Δvecs = _Δvecs + end + + @assert length(Δvals) == length(Δvecs) + @assert length(Δvals) <= length(vals) + + # Determine algorithm to solve linear problem + # TODO: Is there a better choice? Should we make this user configurable? + linalg = GMRES(; + tol=alg.tol, krylovdim=alg.krylovdim, maxiter=alg.maxiter, orth=alg.orth + ) + + ws = similar(vecs, length(Δvecs)) + for i in 1:length(Δvecs) + Δλ = Δvals[i] + Δv = Δvecs[i] + λ = vals[i] + v = vecs[i] + + # First threat special cases + if isa(Δv, AbstractZero) && isa(Δλ, AbstractZero) # no contribution + ws[i] = Δv # some kind of zero + continue + end + if isa(Δv, AbstractZero) && isa(alg, Lanczos) # simple contribution + ws[i] = Δλ * v + continue + end + + # General case : + if isa(Δv, AbstractZero) + b = RecursiveVec(zero(T) * v, T[Δλ]) + else + @assert isa(Δv, typeof(v)) + b = RecursiveVec(Δv, T[Δλ]) + end + + if i > 1 && + eltype(A) <: Real && + vals[i] == conj(vals[i - 1]) && + Δvals[i] == conj(Δvals[i - 1]) && + vecs[i] == conj(vecs[i - 1]) && + Δvecs[i] == conj(Δvecs[i - 1]) + ws[i] = conj(ws[i - 1]) + continue + end + + w, reverse_info = let λ = λ, v = v, Aᴴ = A' + linsolve(b, zero(T) * b, linalg) do x + x1, x2 = x + γ = 1 + # γ can be chosen freely and does not affect the solution theoretically + # The current choice guarantees that the extended matrix is Hermitian if A is + # TODO: is this the best choice in all cases? + y1 = axpy!(-γ * x2[], v, axpy!(-conj(λ), x1, A' * x1)) + y2 = T[-dot(v, x1)] + return RecursiveVec(y1, y2) + end + end + if info.converged >= i && reverse_info.converged == 0 + @warn "The cotangent linear problem did not converge, whereas the primal eigenvalue problem did." + end + ws[i] = w[1] + end + + if A isa StridedMatrix + ∂A = InplaceableThunk( + Ā -> _buildĀ!(Ā, ws, vecs), @thunk(_buildĀ!(zero(A), ws, vecs)) + ) + else + ∂A = @thunk(project_A(_buildĀ!(zero(A), ws, vecs))) + end + return ∂self, ∂A, ∂x₀, ∂howmany, ∂which, ∂alg + end + return (vals, vecs, info), eigsolve_pullback +end + +function _buildĀ!(Ā, ws, vs) + for i in 1:length(ws) + w = ws[i] + v = vs[i] + if !(w isa AbstractZero) + if eltype(Ā) <: Real && eltype(w) <: Complex + mul!(Ā, _realview(w), _realview(v)', -1, 1) + mul!(Ā, _imagview(w), _imagview(v)', -1, 1) + else + mul!(Ā, w, v', -1, 1) + end + end + end + return Ā +end +function _realview(v::AbstractVector{Complex{T}}) where {T} + return view(reinterpret(T, v), 2 * (1:length(v)) .- 1) +end +function _imagview(v::AbstractVector{Complex{T}}) where {T} + return view(reinterpret(T, v), 2 * (1:length(v))) +end + +function ChainRulesCore.rrule( + config::RuleConfig{>:HasReverseMode}, + ::typeof(eigsolve), + A::AbstractMatrix, + x₀, + howmany, + which, + alg, +) + return ChainRulesCore.rrule(eigsolve, A, x₀, howmany, which, alg) +end + +function ChainRulesCore.rrule( + config::RuleConfig{>:HasReverseMode}, ::typeof(eigsolve), f, x₀, howmany, which, alg +) + (vals, vecs, info) = eigsolve(f, x₀, howmany, which, alg) + resize!(vecs, howmany) + resize!(vals, howmany) + T = typeof(dot(vecs[1], vecs[1])) + f_pullbacks = map(x -> rrule_via_ad(config, f, x)[2], vecs) + function eigsolve_pullback(ΔX) + _Δvals = unthunk(ΔX[1]) + _Δvecs = unthunk(ΔX[2]) + + ∂self = NoTangent() + ∂x₀ = ZeroTangent() + ∂howmany = NoTangent() + ∂which = NoTangent() + ∂alg = NoTangent() + if _Δvals isa AbstractZero && _Δvecs isa AbstractZero + ∂A = ZeroTangent() + return (∂self, ∂A, ∂x₀, ∂howmany, ∂which, ∂alg) + end + + if _Δvals isa AbstractZero + Δvals = fill(NoTangent(), howmany) + else + Δvals = _Δvals + end + if _Δvecs isa AbstractZero + Δvecs = fill(NoTangent(), howmany) + else + Δvecs = _Δvecs + end + + # filter ZeroTangents, added compared to Jutho/KrylovKit.jl/pull/56 + Δvecs = filter(x -> !(x isa AbstractZero), Δvecs) + @assert length(Δvals) == length(Δvecs) + + # Determine algorithm to solve linear problem + # TODO: Is there a better choice? Should we make this user configurable? + linalg = GMRES(; + tol=alg.tol, + krylovdim=alg.krylovdim + 10, + maxiter=alg.maxiter * 10, + orth=alg.orth, + ) + # linalg = BiCGStab(; + # tol = alg.tol, + # maxiter = alg.maxiter*alg.krylovdim, + # ) + + ws = similar(Δvecs) + for i in 1:length(Δvecs) + Δλ = Δvals[i] + Δv = Δvecs[i] + λ = vals[i] + v = vecs[i] + + # First threat special cases + if isa(Δv, AbstractZero) && isa(Δλ, AbstractZero) # no contribution + ws[i] = 0 * v # some kind of zero + continue + end + if isa(Δv, AbstractZero) && isa(alg, Lanczos) # simple contribution + ws[i] = Δλ * v + continue + end + + # General case : + if isa(Δv, AbstractZero) + b = RecursiveVec(zero(T) * v, T[-Δλ]) + else + @assert isa(Δv, typeof(v)) + b = RecursiveVec(-Δv, T[-Δλ]) + end + + # TODO: is there any analogy to this for general vector-like user types + # if i > 1 && eltype(A) <: Real && + # vals[i] == conj(vals[i-1]) && Δvals[i] == conj(Δvals[i-1]) && + # vecs[i] == conj(vecs[i-1]) && Δvecs[i] == conj(Δvecs[i-1]) + # + # ws[i] = conj(ws[i-1]) + # continue + # end + + w, reverse_info = let λ = λ, v = v, fᴴ = x -> f_pullbacks[i](x)[2] + linsolve(b, zero(T) * b, linalg) do x + x1, x2 = x + γ = 1 + # γ can be chosen freely and does not affect the solution theoretically + # The current choice guarantees that the extended matrix is Hermitian if A is + # TODO: is this the best choice in all cases? + y1 = axpy!(-γ * x2[], v, axpy!(-conj(λ), x1, fᴴ(x1))) + y2 = T[-dot(v, x1)] + return RecursiveVec(y1, y2) + end + end + if info.converged >= i && reverse_info.converged == 0 + @warn "The cotangent linear problem ($i) did not converge, whereas the primal eigenvalue problem did." reverse_info b + end + ws[i] = w[1] + end + ∂f = f_pullbacks[1](ws[1])[1] + for i in 2:length(ws) + ∂f = VectorInterface.add!!(∂f, f_pullbacks[i](ws[i])[1]) + end + return ∂self, ∂f, ∂x₀, ∂howmany, ∂which, ∂alg + end + return (vals, vecs, info), eigsolve_pullback +end diff --git a/src/utility/hook_pullback.jl b/src/utility/hook_pullback.jl new file mode 100644 index 00000000..7b5e2307 --- /dev/null +++ b/src/utility/hook_pullback.jl @@ -0,0 +1,48 @@ +#= +In order to use the machinery of AD but still have the option of controlling the algorithm +for computing gradients, it is possible to hook into the pullback of a function and define +a custom rrule. This is achieved via a wrapper function that takes an optional keyword +argument `alg_rrule` which is consequently used to customize the pullback. +In order to be able to specialize on this, we wrap the rrule in a function `_rrule` which +receives this as its first positional argument. +=# + +""" + hook_pullback(f, args...; alg_rrule=nothing, kwargs...) + +Wrapper function to customize the pullback of a function `f`. This function is equivalent to +`f(args...; kwargs...)`, but the pullback can be customized by implementing the following +function: + + _rrule(alg_rrule, config, f, args...; kwargs...) -> NoTangent(), ∂f, ∂args... + +This function can specialize on its first argument in order to customize the pullback. If no +specialization is needed, the default `alg_rrule=nothing` results in the default AD +pullback. + +See also [`_rrule`](@ref). +""" +function hook_pullback(@nospecialize(f), args...; alg_rrule=nothing, kwargs...) + return f(args...; kwargs...) +end + +function ChainRulesCore.rrule( + config::RuleConfig, ::typeof(hook_pullback), f, args...; alg_rrule=nothing, kwargs... +) + # Need to add ∂hook_pullback + y, f_pullback = _rrule(alg_rrule, config, f, args...; kwargs...) + return y, Δ -> (NoTangent(), f_pullback(Δ)...) +end + +""" + _rrule(alg_rrule, config, f, args...; kwargs...) -> ∂f, ∂args... + +Customize the pullback of a function `f`. This function can specialize on its first argument +in order to have multiple implementations for a pullback. If no specialization is needed, +the default `alg_rrule=nothing` results in the default AD pullback. + +!!! warning + No tangent is expected for the `alg_rrule` argument +""" +_rrule(::Nothing, config::RuleConfig, f, args...; kwargs...) = + rrule_via_ad(config, f, args...; kwargs...) diff --git a/src/utility/rotations.jl b/src/utility/rotations.jl index 58d96c3d..634e6199 100644 --- a/src/utility/rotations.jl +++ b/src/utility/rotations.jl @@ -1,13 +1,16 @@ -const NORTH = 1; -const EAST = 2; -const SOUTH = 3; -const WEST = 4; +const NORTH = 1 +const EAST = 2 +const SOUTH = 3 +const WEST = 4 -const NORTHWEST = 1; -const NORTHEAST = 2; -const SOUTHEAST = 3; -const SOUTHWEST = 4; +const NORTHWEST = 1 +const NORTHEAST = 2 +const SOUTHEAST = 3 +const SOUTHWEST = 4 -rotate_north(t, dir) = mod1(dir, 4) == NORTH ? t : rotate_north(rotl90(t), dir - 1) +""" + rotate_north(t, dir) -#Base.rotl90(t::Tuple) = rotl90.(t); +Rotate north direction of `t` to `dir` by successive applications of `rotl90`. +""" +rotate_north(t, dir) = mod1(dir, 4) == NORTH ? t : rotate_north(rotl90(t), dir - 1) diff --git a/src/utility/symmetrization.jl b/src/utility/symmetrization.jl index 88cc6be1..55a97820 100644 --- a/src/utility/symmetrization.jl +++ b/src/utility/symmetrization.jl @@ -1,5 +1,4 @@ -# an initial attempt at some basic symmetrization routines for PEPS -# it looks horrible, but it works, maybe +# some basic symmetrization routines for PEPS abstract type SymmetrizationStyle end @@ -18,7 +17,6 @@ function herm_depth(x::PEPOTensor) end function herm_width(x::PEPSTensor) - x = Permute(Conj(x), [1, 4, 3, 2, 5:(x.legs)]) return permute(x', ((5,), (1, 4, 3, 2))) end function herm_width(x::PEPOTensor) @@ -26,7 +24,6 @@ function herm_width(x::PEPOTensor) end function herm_height(x::PEPOTensor) - x = Permute(Conj(x), [1:4, 6, 5]) return permute(x', ((6, 5), (1, 2, 3, 4))) end @@ -63,10 +60,22 @@ function herm_height_inv(x::Union{PEPSTensor,PEPOTensor}) end # rotation invariance - +Base.rotl90(t::PEPSTensor) = permute(t, ((1,), (3, 4, 5, 2))) Base.rotr90(t::PEPSTensor) = permute(t, ((1,), (5, 2, 3, 4))) Base.rot180(t::PEPSTensor) = permute(t, ((1,), (4, 5, 2, 3))) +Base.rotl90(t::PEPOTensor) = permute(t, ((1, 2), (4, 5, 6, 3))) +Base.rotr90(t::PEPOTensor) = permute(t, ((1, 2), (6, 3, 4, 5))) +Base.rot180(t::PEPOTensor) = permute(t, ((1, 2), (5, 6, 3, 4))) + +Base.rotl90(t::InfinitePEPS) = InfinitePEPS(rotl90(rotl90.(t.A))) +Base.rotr90(t::InfinitePEPS) = InfinitePEPS(rotr90(rotr90.(t.A))) +Base.rot180(t::InfinitePEPS) = InfinitePEPS(rot180(rot180.(t.A))) + +Base.rotl90(T::InfinitePEPO) = InfinitePEPO(stack(rotl90, eachslice(T.A; dims=3))) +Base.rotr90(T::InfinitePEPO) = InfinitePEPO(stack(rotr90, eachslice(T.A; dims=3))) +Base.rot180(T::InfinitePEPO) = InfinitePEPO(stack(rot180, eachslice(T.A; dims=3))) + function rot_inv(x) return 0.25 * ( x + diff --git a/src/utility/util.jl b/src/utility/util.jl index 5ebf8d4e..0797a412 100644 --- a/src/utility/util.jl +++ b/src/utility/util.jl @@ -1,36 +1,71 @@ +# Get next and previous directional CTM enviroment index, respecting periodicity +_next(i, total) = mod1(i + 1, total) +_prev(i, total) = mod1(i - 1, total) + +# Element-wise multiplication of TensorMaps respecting block structure +function _elementwise_mult(a::AbstractTensorMap, b::AbstractTensorMap) + dst = similar(a) + for (k, block) in blocks(dst) + copyto!(block, blocks(a)[k] .* blocks(b)[k]) + end + return dst +end + +# Compute √S⁻¹ for diagonal TensorMaps function sdiag_inv_sqrt(S::AbstractTensorMap) - toret = similar(S) + invsq = similar(S) if sectortype(S) == Trivial - copyto!(toret.data, LinearAlgebra.diagm(LinearAlgebra.diag(S.data) .^ (-1 / 2))) + copyto!(invsq.data, LinearAlgebra.diagm(LinearAlgebra.diag(S.data) .^ (-1 / 2))) else for (k, b) in blocks(S) copyto!( - blocks(toret)[k], LinearAlgebra.diagm(LinearAlgebra.diag(b) .^ (-1 / 2)) + blocks(invsq)[k], LinearAlgebra.diagm(LinearAlgebra.diag(b) .^ (-1 / 2)) ) end end - return toret + return invsq end + function ChainRulesCore.rrule(::typeof(sdiag_inv_sqrt), S::AbstractTensorMap) - toret = sdiag_inv_sqrt(S) - return toret, - c̄ -> (ChainRulesCore.NoTangent(), -1 / 2 * _elementwise_mult(c̄, toret'^3)) + invsq = sdiag_inv_sqrt(S) + function sdiag_inv_sqrt_pullback(c̄) + return (ChainRulesCore.NoTangent(), -1 / 2 * _elementwise_mult(c̄, invsq'^3)) + end + return invsq, sdiag_inv_sqrt_pullback end -function _elementwise_mult(a::AbstractTensorMap, b::AbstractTensorMap) - dst = similar(a) - for (k, block) in blocks(dst) - copyto!(block, blocks(a)[k] .* blocks(b)[k]) + +# Check whether diagonals contain degenerate values up to absolute or relative tolerance +function is_degenerate_spectrum( + S; atol::Real=0, rtol::Real=atol > 0 ? 0 : sqrt(eps(scalartype(S))) +) + for (_, b) in blocks(S) + s = real(diag(b)) + for i in 1:(length(s) - 1) + isapprox(s[i], s[i + 1]; atol, rtol) && return true + end end - return dst + return false end -#rotl90 appeared to lose PeriodicArray'ness, which in turn caused zygote problems -#Base.rotl90(a::Array) = Array(rotl90(a)); -#Base.rotr90(a::Array) = Array(rotr90(a)); +""" + projector_type(T::DataType, size) + +Create two arrays of specified `size` that contain undefined tensors representing +left and right acting projectors, respectively. The projector types are inferred +from the TensorMap type `T` which avoids having to recompute transpose tensors. +""" +function projector_type(T::DataType, size) + Pleft = Array{T,length(size)}(undef, size) + Prtype = tensormaptype(spacetype(T), numin(T), numout(T), storagetype(T)) + Pright = Array{Prtype,length(size)}(undef, size) + return Pleft, Pright +end + +# There are no rrules for rotl90 and rotr90 in ChainRules.jl function ChainRulesCore.rrule(::typeof(rotl90), a::AbstractMatrix) - function pb(x) + function rotl90_pullback(x) if !iszero(x) x = if x isa Tangent ChainRulesCore.construct(typeof(a), ChainRulesCore.backing(x)) @@ -40,13 +75,28 @@ function ChainRulesCore.rrule(::typeof(rotl90), a::AbstractMatrix) x = rotr90(x) end - return (ZeroTangent(), x) + return NoTangent(), x end - return rotl90(a), pb + return rotl90(a), rotl90_pullback end -structure(t) = codomain(t) ← domain(t); +function ChainRulesCore.rrule(::typeof(rotr90), a::AbstractMatrix) + function rotr90_pullback(x) + if !iszero(x) + x = if x isa Tangent + ChainRulesCore.construct(typeof(a), ChainRulesCore.backing(x)) + else + x + end + x = rotl90(x) + end + + return NoTangent(), x + end + return rotr90(a), rotr90_pullback +end +# Differentiable setindex! alternative function _setindex(a::AbstractArray, v, args...) b::typeof(a) = copy(a) b[args...] = v @@ -56,7 +106,7 @@ end function ChainRulesCore.rrule(::typeof(_setindex), a::AbstractArray, tv, args...) t = _setindex(a, tv, args...) - function toret(v) + function _setindex_pullback(v) if iszero(v) backwards_tv = ZeroTangent() backwards_a = ZeroTangent() @@ -66,6 +116,7 @@ function ChainRulesCore.rrule(::typeof(_setindex), a::AbstractArray, tv, args... else v end + # TODO: Fix this for ZeroTangents v = typeof(v) != typeof(a) ? convert(typeof(a), v) : v #v = convert(typeof(a),v); backwards_tv = v[args...] @@ -80,9 +131,18 @@ function ChainRulesCore.rrule(::typeof(_setindex), a::AbstractArray, tv, args... NoTangent(), backwards_a, backwards_tv, fill(ZeroTangent(), length(args))... ) end - return t, toret + return t, _setindex_pullback end +""" + @diffset assign + +Helper macro which allows in-place operations in the forward-pass of Zygote, but +resorts to non-mutating operations in the backwards-pass. The expression `assign` +should assign an object to an pre-existing `AbstractArray` and the use of updating +operators is also possible. This is especially needed when in-place assigning +tensors to unit-cell arrays of environments. +""" macro diffset(ex) return esc(parse_ex(ex)) end @@ -121,3 +181,19 @@ end is_indexing(ex) = false is_indexing(ex::Expr) = ex.head == :ref + +""" + @showtypeofgrad(x) + +Macro utility to show to type of the gradient that is about to accumulate for `x`. + +See also [`Zygote.@showgrad`](@ref). +""" +macro showtypeofgrad(x) + return :( + Zygote.hook($(esc(x))) do x̄ + println($"∂($x) = ", repr(typeof(x̄))) + x̄ + end + ) +end diff --git a/test/Project.toml b/test/Project.toml deleted file mode 100644 index b9240c28..00000000 --- a/test/Project.toml +++ /dev/null @@ -1,5 +0,0 @@ -[deps] -MPSKit = "bb1c41ca-d63c-52ed-829e-0820dda26502" -OptimKit = "77e91f04-9b3b-57a6-a776-40b61faaebe0" -TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/test/boundarymps/vumps.jl b/test/boundarymps/vumps.jl new file mode 100644 index 00000000..260f40b3 --- /dev/null +++ b/test/boundarymps/vumps.jl @@ -0,0 +1,72 @@ +using Test +using Random +using PEPSKit +using TensorKit +using MPSKit +using LinearAlgebra + +Random.seed!(29384293742893) + +@testset "(1, 1) PEPS" begin + psi = InfinitePEPS(ComplexSpace(2), ComplexSpace(2)) + + T = PEPSKit.InfiniteTransferPEPS(psi, 1, 1) + mps = PEPSKit.initializeMPS(T, [ComplexSpace(20)]) + + mps, envs, ϵ = leading_boundary(mps, T, VUMPS()) + N = abs(sum(expectation_value(mps, T))) + + ctm = leading_boundary( + CTMRGEnv(psi; Venv=ComplexSpace(20)), psi, CTMRG(; verbosity=1, fixedspace=true) + ) + 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( + CTMRGEnv(psi; Venv=ComplexSpace(20)), psi, CTMRG(; verbosity=1, fixedspace=true) + ) + 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 diff --git a/test/ctmrg/gaugefix.jl b/test/ctmrg/gaugefix.jl new file mode 100644 index 00000000..45bee5d5 --- /dev/null +++ b/test/ctmrg/gaugefix.jl @@ -0,0 +1,82 @@ +using Test +using Random +using PEPSKit +using TensorKit + +using PEPSKit: ctmrg_iter, gauge_fix, check_elementwise_convergence + +scalartypes = [Float64, ComplexF64] +unitcells = [(1, 1), (2, 2), (3, 2)] +χ = 6 + +function _make_symmetric(psi) + if ==(size(psi)...) + return PEPSKit.symmetrize(psi, PEPSKit.Full()) + else + return PEPSKit.symmetrize(PEPSKit.symmetrize(psi, PEPSKit.Depth()), PEPSKit.Width()) + end +end + +# If I can't make the rng seed behave, I'll just randomly define a peps somehow +function semi_random_peps!(psi::InfinitePEPS) + i = 0 + A′ = map(psi.A) do a + for (_, b) in blocks(a) + l = length(b) + b .= reshape(collect((1:l) .+ i), size(b)) + i += l + end + return a + end + return InfinitePEPS(A′) +end + +@testset "Trivial symmetry ($T) - ($unitcell)" for (T, unitcell) in + Iterators.product(scalartypes, unitcells) + physical_space = ComplexSpace(2) + peps_space = ComplexSpace(2) + ctm_space = ComplexSpace(χ) + + psi = InfinitePEPS(undef, T, physical_space, peps_space; unitcell) + semi_random_peps!(psi) + psi = _make_symmetric(psi) + + Random.seed!(987654321) # Seed RNG to make random environment consistent + ctm = CTMRGEnv(psi; Venv=ctm_space) + + verbosity = 1 + alg = CTMRG(; + trscheme=truncdim(dim(ctm_space)), tol=1e-10, miniter=4, maxiter=400, verbosity + ) + alg_fixed = CTMRG(; trscheme=truncdim(dim(ctm_space)), verbosity, fixedspace=true) + + ctm = leading_boundary(ctm, psi, alg) + ctm2, = ctmrg_iter(psi, ctm, alg_fixed) + ctm_fixed = gauge_fix(ctm, ctm2) + @test PEPSKit.check_elementwise_convergence(ctm, ctm_fixed; atol=1e-4) +end + +@testset "Z2 symmetry ($T) - ($unitcell)" for (T, unitcell) in + Iterators.product(scalartypes, unitcells) + physical_space = Z2Space(0 => 1, 1 => 1) + peps_space = Z2Space(0 => 1, 1 => 1) + ctm_space = Z2Space(0 => χ ÷ 2, 1 => χ ÷ 2) + + psi = InfinitePEPS(undef, T, physical_space, peps_space; unitcell) + semi_random_peps!(psi) + psi = _make_symmetric(psi) + + Random.seed!(123456789) # Seed RNG to make random environment consistent + ctm = CTMRGEnv(psi; Venv=ctm_space) + + verbosity = 1 + alg = CTMRG(; + trscheme=truncspace(ctm_space), tol=1e-10, miniter=4, maxiter=400, verbosity + ) + alg_fixed = CTMRG(; trscheme=truncspace(ctm_space), verbosity, fixedspace=true) + + ctm = leading_boundary(ctm, psi, alg) + ctm2, = ctmrg_iter(psi, ctm, alg_fixed) + ctm_fixed = gauge_fix(ctm, ctm2) + @test PEPSKit.check_elementwise_convergence(ctm, ctm_fixed; atol=1e-4) +end diff --git a/test/ctmrg/gradients.jl b/test/ctmrg/gradients.jl new file mode 100644 index 00000000..2c9b86f9 --- /dev/null +++ b/test/ctmrg/gradients.jl @@ -0,0 +1,57 @@ +using Test +using Random +using PEPSKit +using TensorKit +using Zygote +using OptimKit +using KrylovKit + +## Test models, gradmodes and CTMRG algorithm +# ------------------------------------------- +χbond = 2 +χenv = 4 +Pspaces = [ComplexSpace(2), Vect[FermionParity](0 => 1, 1 => 1)] +Vspaces = [ComplexSpace(χbond), Vect[FermionParity](0 => χbond / 2, 1 => χbond / 2)] +Espaces = [ComplexSpace(χenv), Vect[FermionParity](0 => χenv / 2, 1 => χenv / 2)] +models = [square_lattice_heisenberg(), square_lattice_pwave()] +names = ["Heisenberg", "p-wave superconductor"] +Random.seed!(42039482030) +tol = 1e-8 +boundary_alg = CTMRG(; + trscheme=truncdim(χenv), tol=tol, miniter=4, maxiter=100, fixedspace=true, verbosity=0 +) +gradmodes = [ + nothing, GeomSum(; tol), ManualIter(; tol), KrylovKit.GMRES(; tol=tol, maxiter=10) +] +steps = -0.01:0.005:0.01 + +## Tests +# ------ +@testset "AD CTMRG energy gradients for $(names[i]) model" for i in eachindex(models) + Pspace = Pspaces[i] + Vspace = Pspaces[i] + Espace = Espaces[i] + psi_init = InfinitePEPS(Pspace, Vspace, Vspace) + @testset "$alg_rrule" for alg_rrule in gradmodes + dir = InfinitePEPS(Pspace, Vspace, Vspace) + psi = InfinitePEPS(Pspace, Vspace, Vspace) + env = leading_boundary(CTMRGEnv(psi; Venv=Espace), psi, boundary_alg) + alphas, fs, dfs1, dfs2 = OptimKit.optimtest( + (psi, env), + dir; + alpha=steps, + retract=PEPSKit.my_retract, + inner=PEPSKit.my_inner, + ) do (peps, envs) + E, g = Zygote.withgradient(peps) do psi + envs2 = PEPSKit.hook_pullback( + leading_boundary, envs, psi, boundary_alg, ; alg_rrule + ) + return costfun(psi, envs2, models[i]) + end + + return E, only(g) + end + @test dfs1 ≈ dfs2 atol = 1e-2 + end +end diff --git a/test/ctmrg/gradparts.jl b/test/ctmrg/gradparts.jl new file mode 100644 index 00000000..c21a2420 --- /dev/null +++ b/test/ctmrg/gradparts.jl @@ -0,0 +1,84 @@ +using Test +using Random +using PEPSKit +using TensorKit +using PEPSKit: + NORTH, + SOUTH, + WEST, + EAST, + NORTHWEST, + NORTHEAST, + SOUTHEAST, + SOUTHWEST, + rotate_north, + left_move, + ctmrg_iter + +using Zygote +using ChainRulesCore +using ChainRulesTestUtils + +include(joinpath("..", "utility.jl")) + +## Test spaces, tested functions and CTMRG algorithm +# -------------------------------------------------- +χbond = 2 +χenv = 4 +Pspaces = [ComplexSpace(2), Vect[FermionParity](0 => 1, 1 => 1)] +Vspaces = [ComplexSpace(χbond), Vect[FermionParity](0 => χbond / 2, 1 => χbond / 2)] +Espaces = [ComplexSpace(χenv), Vect[FermionParity](0 => χenv / 2, 1 => χenv / 2)] +functions = [left_move, ctmrg_iter, leading_boundary] +tol = 1e-8 +boundary_alg = CTMRG(; + trscheme=truncdim(χenv), tol=tol, miniter=4, maxiter=100, fixedspace=true, verbosity=0 +) + +## Gauge invariant function of the environment +# -------------------------------------------- +function rho(env) + @tensor ρ[-1 -2 -3 -4 -5 -6 -7 -8] := + env.edges[WEST, 1, 1][1 -1 -2; 4] * + env.corners[NORTHWEST, 1, 1][4; 5] * + env.edges[NORTH, 1, 1][5 -3 -4; 8] * + env.corners[NORTHEAST, 1, 1][8; 9] * + env.edges[EAST, 1, 1][9 -5 -6; 12] * + env.corners[SOUTHEAST, 1, 1][12; 13] * + env.edges[SOUTH, 1, 1][13 -7 -8; 16] * + env.corners[SOUTHWEST, 1, 1][16; 1] + return ρ +end + +## Tests +# ------ +@testset "Reverse rules for composite parts of the CTMRG fixed point with spacetype $(Vspaces[i])" for i in + eachindex( + Pspaces +) + psi = InfinitePEPS(Pspaces[i], Vspaces[i], Vspaces[i]) + env = CTMRGEnv(psi; Venv=Espaces[i]) + + @testset "$f" for f in functions + atol = f == leading_boundary ? sqrt(tol) : tol + g = if f == leading_boundary + function (state, env) + return rho(f(env, state, boundary_alg)) + end + else + function (state, env) + return rho(f(state, env, boundary_alg)[1]) + end + end + + # use rrule testing functionality but compute rrule via Zygote + test_rrule( + Zygote.ZygoteRuleConfig(), + g, + psi, + env; + check_inferred=false, + atol, + rrule_f=rrule_via_ad, + ) + end +end diff --git a/test/heisenberg.jl b/test/heisenberg.jl new file mode 100644 index 00000000..35e12ef7 --- /dev/null +++ b/test/heisenberg.jl @@ -0,0 +1,29 @@ +using Test +using Random +using PEPSKit +using TensorKit +using KrylovKit +using OptimKit + +# Initialize parameters +H = square_lattice_heisenberg() +χbond = 2 +χenv = 16 +ctm_alg = CTMRG(; trscheme=truncdim(χenv), tol=1e-10, miniter=4, maxiter=100, verbosity=1) +opt_alg = PEPSOptimize(; + boundary_alg=ctm_alg, + optimizer=LBFGS(4; maxiter=100, gradtol=1e-3, verbosity=2), + gradient_alg=GMRES(; tol=1e-6, maxiter=100), + reuse_env=true, + verbosity=2, +) + +# initialize states +Random.seed!(91283219347) +psi_init = InfinitePEPS(2, χbond) +env_init = leading_boundary(CTMRGEnv(psi_init; Venv=ComplexSpace(χenv)), psi_init, ctm_alg) + +# find fixedpoint +result = fixedpoint(psi_init, H, opt_alg, env_init) + +@test result.E ≈ -0.6694421 atol = 1e-2 diff --git a/test/pwave.jl b/test/pwave.jl new file mode 100644 index 00000000..4ebb9c85 --- /dev/null +++ b/test/pwave.jl @@ -0,0 +1,35 @@ +using Test +using Random +using PEPSKit +using TensorKit +using KrylovKit +using OptimKit + +# Initialize parameters +H = square_lattice_pwave() +χbond = 2 +χenv = 16 +ctm_alg = CTMRG(; + trscheme=truncdim(χenv), tol=1e-10, miniter=4, maxiter=50, fixedspace=true, verbosity=1 +) +opt_alg = PEPSOptimize(; + boundary_alg=ctm_alg, + optimizer=LBFGS(4; maxiter=10, gradtol=1e-3, verbosity=2), + gradient_alg=GMRES(; tol=1e-3, maxiter=2, krylovdim=50, verbosity=2), + reuse_env=true, + verbosity=2, +) + +# initialize states +Random.seed!(96678827397) +Pspace = Vect[FermionParity](0 => 1, 1 => 1) +Vspace = Vect[FermionParity](0 => χbond ÷ 2, 1 => χbond ÷ 2) +Envspace = Vect[FermionParity](0 => χenv ÷ 2, 1 => χenv ÷ 2) +psi_init = InfinitePEPS(Pspace, Vspace, Vspace; unitcell=(2, 2)) +env_init = leading_boundary(CTMRGEnv(psi_init; Venv=Envspace), psi_init, ctm_alg); + +# find fixedpoint +result = fixedpoint(psi_init, H, opt_alg, env_init) + +# comparison with Gaussian PEPS minimum at D=2 on 1000x1000 square lattice with aPBC +@test result.E / prod(size(psi_init)) ≈ -2.6241 atol = 5e-2 diff --git a/test/runtests.jl b/test/runtests.jl index 3255ddb4..af2023e1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,19 +1,27 @@ -using Test, PEPSKit, MPSKit, TensorKit +using SafeTestsets -@testset "Boundary MPS" begin - peps = InfinitePEPS(2, 3) - tpeps = InfiniteTransferPEPS(peps, 1, 1) +const GROUP = get(ENV, "GROUP", "All") - mps = initializeMPS(tpeps, 4) - - mps, _, _ = leading_boundary(mps, tpeps, VUMPS()) -end - -@testset "CTMRG" begin - peps = InfinitePEPS(2, 3) - tpeps = InfiniteTransferPEPS(peps, 1, 1) - - mps = initializeMPS(tpeps, 4) - - env = leading_boundary(peps, CTMRG(; tol=1e-10)) +@time begin + if GROUP == "All" || GROUP == "CTMRG" + @time @safetestset "Gauge Fixing" begin + include("ctmrg/gaugefix.jl") + end + @time @safetestset "Gradients" begin + include("ctmrg/gradients.jl") + end + end + if GROUP == "All" || GROUP == "MPS" + @time @safetestset "VUMPS" begin + include("boundarymps/vumps.jl") + end + end + if GROUP == "All" || GROUP == "EXAMPLES" + @time @safetestset "Heisenberg model" begin + include("heisenberg.jl") + end + @time @safetestset "P-wave superconductor" begin + include("pwave.jl") + end + end end diff --git a/test/utility.jl b/test/utility.jl new file mode 100644 index 00000000..383113da --- /dev/null +++ b/test/utility.jl @@ -0,0 +1,78 @@ +using TensorKit +using ChainRulesTestUtils + +using TensorKit: sqrtdim, isqrtdim +using VectorInterface: scale! +using FiniteDifferences + +## Test utility +# ------------- +function ChainRulesTestUtils.rand_tangent(rng::AbstractRNG, x::AbstractTensorMap) + return TensorMap(randn, scalartype(x), space(x)) +end +function ChainRulesTestUtils.rand_tangent(rng::AbstractRNG, x::CTMRGEnv) + Ctans = x.corners + Etans = x.edges + for i in eachindex(x.corners) + Ctans[i] = rand_tangent(rng, x.corners[i]) + end + for i in eachindex(x.edges) + Etans[i] = rand_tangent(rng, x.edges[i]) + end + return CTMRGEnv(Ctans, Etans) +end +function ChainRulesTestUtils.test_approx( + actual::AbstractTensorMap, expected::AbstractTensorMap, msg=""; kwargs... +) + for (c, b) in blocks(actual) + ChainRulesTestUtils.@test_msg msg isapprox(b, block(expected, c); kwargs...) + end +end +function ChainRulesTestUtils.test_approx( + actual::InfinitePEPS, expected::InfinitePEPS, msg=""; kwargs... +) + for i in eachindex(size(actual, 1)) + for j in eachindex(size(actual, 2)) + ChainRulesTestUtils.@test_msg msg isapprox( + actual[i, j], expected[i, j]; kwargs... + ) + end + end +end +function ChainRulesTestUtils.test_approx( + actual::CTMRGEnv, expected::CTMRGEnv, msg=""; kwargs... +) + for i in eachindex(actual.corners) + ChainRulesTestUtils.@test_msg msg isapprox( + actual.corners[i], expected.corners[i]; kwargs... + ) + end + for i in eachindex(actual.edges) + ChainRulesTestUtils.@test_msg msg isapprox( + actual.edges[i], expected.edges[i]; kwargs... + ) + end +end + +# TODO: remove these functions once TensorKit is updated +function FiniteDifferences.to_vec(t::T) where {T<:TensorKit.TrivialTensorMap} + vec, from_vec = to_vec(t.data) + return vec, x -> T(from_vec(x), codomain(t), domain(t)) +end +function FiniteDifferences.to_vec(t::AbstractTensorMap) + # convert to vector of vectors to make use of existing functionality + vec_of_vecs = [b * sqrtdim(c) for (c, b) in blocks(t)] + vec, back = FiniteDifferences.to_vec(vec_of_vecs) + + function from_vec(x) + t′ = similar(t) + xvec_of_vecs = back(x) + for (i, (c, b)) in enumerate(blocks(t′)) + scale!(b, xvec_of_vecs[i], isqrtdim(c)) + end + return t′ + end + + return vec, from_vec +end +FiniteDifferences.to_vec(t::TensorKit.AdjointTensorMap) = to_vec(copy(t))