diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index aad659d0..0adc6bdc 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -21,7 +21,7 @@ jobs: fail-fast: false matrix: version: - - '1.8' # LTS version + - '1.9' - '1' # automatically expands to the latest stable 1.x release of Julia os: - ubuntu-latest diff --git a/Project.toml b/Project.toml index 08904cc6..b5be59ea 100644 --- a/Project.toml +++ b/Project.toml @@ -24,18 +24,18 @@ Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" Accessors = "0.1" ChainRulesCore = "1.0" Compat = "3.46, 4.2" -KrylovKit = "0.6, 0.7, 0.8" +KrylovKit = "0.8" LinearAlgebra = "1" MPSKit = "0.11" OptimKit = "0.3" Printf = "1" Random = "1" Statistics = "1" -TensorKit = "0.12.4" +TensorKit = "0.12.5" TensorOperations = "4" VectorInterface = "0.4" Zygote = "0.6" -julia = "1.8" +julia = "1.9" [extras] ChainRulesTestUtils = "cdddcdb0-9152-4a09-a978-84456f9df70a" diff --git a/README.md b/README.md index a45818d6..d5a41ca5 100644 --- a/README.md +++ b/README.md @@ -36,19 +36,12 @@ For example, in order to obtain the groundstate of the 2D Heisenberg model, we c 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) +H = square_lattice_heisenberg(; Jx=-1, Jy=1, Jz=-1) # sublattice rotation to obtain single-site unit cell # configuring the parameters D = 2 chi = 20 -ctm_alg = CTMRG(; trscheme = truncdim(chi), tol=1e-20, miniter=4, maxiter=100, verbosity=1) +ctm_alg = CTMRG(; tol=1e-10, miniter=4, maxiter=100, verbosity=1, trscheme=truncdim(chi)) opt_alg = PEPSOptimize(; boundary_alg=ctm_alg, optimizer=LBFGS(4; maxiter=100, gradtol=1e-4, verbosity=2), @@ -59,8 +52,8 @@ opt_alg = PEPSOptimize(; # 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) +ctm = leading_boundary(CTMRGEnv(state, ComplexSpace(chi)), state, ctm_alg) +result = fixedpoint(state, H, opt_alg, ctm) @show result.E # -0.6625... ``` diff --git a/docs/src/index.md b/docs/src/index.md index 305f6d84..5ceb20c7 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -21,19 +21,12 @@ For example, in order to obtain the groundstate of the 2D Heisenberg model, we c 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) +H = square_lattice_heisenberg(; Jx=-1, Jy=1, Jz=-1) # sublattice rotation to obtain single-site unit cell # configuring the parameters D = 2 chi = 20 -ctm_alg = CTMRG(; trscheme = truncdim(chi), tol=1e-20, miniter=4, maxiter=100, verbosity=1) +ctm_alg = CTMRG(; tol=1e-10, miniter=4, maxiter=100, verbosity=1, trscheme=truncdim(chi)) opt_alg = PEPSOptimize(; boundary_alg=ctm_alg, optimizer=LBFGS(4; maxiter=100, gradtol=1e-4, verbosity=2), @@ -44,8 +37,8 @@ opt_alg = PEPSOptimize(; # 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) +ctm = leading_boundary(CTMRGEnv(state, ComplexSpace(chi)), state, ctm_alg) +result = fixedpoint(state, H, opt_alg, ctm) @show result.E # -0.6625... ``` diff --git a/examples/boundary_mps.jl b/examples/boundary_mps.jl index abd52e17..eb699b21 100644 --- a/examples/boundary_mps.jl +++ b/examples/boundary_mps.jl @@ -32,9 +32,7 @@ mps, envs, ϵ = leading_boundary(mps, T, VUMPS()) 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)) -) +ctm = leading_boundary(peps, CTMRG(; verbosity=1), CTMRGEnv(peps, ComplexSpace(20))) N´ = abs(norm(peps, ctm)) @show abs(N - N´) / N @@ -55,9 +53,7 @@ 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)) -) +ctm2 = leading_boundary(peps2, CTMRG(; verbosity=1), CTMRGEnv(peps2, ComplexSpace(20))) N2´ = abs(norm(peps2, ctm2)) @show abs(N2 - N2´) / N2 diff --git a/examples/heisenberg.jl b/examples/heisenberg.jl index 0cfd8c06..9dd8dbe3 100644 --- a/examples/heisenberg.jl +++ b/examples/heisenberg.jl @@ -11,9 +11,9 @@ 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, +ctm_alg = CTMRG(; tol=1e-10, miniter=4, maxiter=100, verbosity=1, trscheme=truncdim(χenv)) +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, @@ -27,6 +27,6 @@ alg = PEPSOptimize(; # 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₀) +env₀ = leading_boundary(CTMRGEnv(ψ₀, ℂ^χenv), ψ₀, ctm_alg) +result = fixedpoint(ψ₀, H, opt_alg, env₀) @show result.E diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index fe655d05..93595923 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -11,7 +11,7 @@ using LoggingExtras using MPSKit: loginit!, logiter!, logfinish!, logcancel! include("utility/util.jl") -include("utility/eigsolve.jl") +include("utility/svd.jl") include("utility/rotations.jl") include("utility/hook_pullback.jl") include("utility/autoopt.jl") @@ -61,12 +61,14 @@ module Defaults const fpgrad_tol = 1e-6 end -export CTMRG, CTMRGEnv +export SVDAdjoint, IterSVD, NonTruncSVDAdjoint +export FixedSpaceTruncation, ProjectorAlg, CTMRG, CTMRGEnv export LocalOperator export expectation_value, costfun export leading_boundary export PEPSOptimize, GeomSum, ManualIter, LinSolve export fixedpoint + export InfinitePEPS, InfiniteTransferPEPS export InfinitePEPO, InfiniteTransferPEPO export initializeMPS, initializePEPS diff --git a/src/algorithms/ctmrg.jl b/src/algorithms/ctmrg.jl index bef9c554..a65345d9 100644 --- a/src/algorithms/ctmrg.jl +++ b/src/algorithms/ctmrg.jl @@ -1,26 +1,62 @@ +""" + FixedSpaceTruncation <: TensorKit.TruncationScheme + +CTMRG specific truncation scheme for `tsvd` which keeps the bond space on which the SVD +is performed fixed. Since different environment directions and unit cell entries might +have different spaces, this truncation style is different from `TruncationSpace`. +""" +struct FixedSpaceTruncation <: TensorKit.TruncationScheme end + +# TODO: add option for different projector styles (half-infinite, full-infinite, etc.) +""" + struct ProjectorAlg{S}(; svd_alg = TensorKit.SVD(), trscheme = TensorKit.notrunc(), + fixedspace = false, verbosity = 0) + +Algorithm struct collecting all projector related parameters. The truncation scheme has to be +a `TensorKit.TruncationScheme`, and some SVD algorithms might have further restrictions on what +kind of truncation scheme can be used. If `fixedspace` is true, the truncation scheme is set to +`truncspace(V)` where `V` is the environment bond space, adjusted to the corresponding +environment direction/unit cell entry. +""" +@kwdef struct ProjectorAlg{S<:SVDAdjoint,T} + svd_alg::S = SVDAdjoint() + trscheme::T = FixedSpaceTruncation() + verbosity::Int = 0 +end + # TODO: add abstract Algorithm type? """ - struct CTMRG(; trscheme = TensorKit.notrunc(), tol = Defaults.ctmrg_tol, - maxiter = Defaults.ctmrg_maxiter, miniter = Defaults.ctmrg_miniter, - verbosity = 0, fixedspace = false) + CTMRG(; tol=Defaults.ctmrg_tol, maxiter=Defaults.ctmrg_maxiter, + miniter=Defaults.ctmrg_miniter, verbosity=0, + svd_alg=TensorKit.SVD(), trscheme=FixedSpaceTruncation()) Algorithm struct that represents the CTMRG algorithm for contracting infinite PEPS. -The projector bond dimensions are set via `trscheme` which controls the truncation -properties inside of `TensorKit.tsvd`. Each CTMRG run is converged up to `tol` -where the singular value convergence of the corners as well as the norm is checked. -The maximal and minimal number of CTMRG iterations is set with `maxiter` and `miniter`. -Different levels of output information are printed depending on `verbosity`, where -`0` is silent, `1` only prints warnings, `2` prints start and end of the run, and `3` -displays information each iteration. -Regardless of the truncation scheme, the space can be kept fixed with `fixedspace`. +Each CTMRG run is converged up to `tol` where the singular value convergence of the +corners as well as the norm is checked. The maximal and minimal number of CTMRG iterations +is set with `maxiter` and `miniter`. Different levels of output information are printed +depending on `verbosity`, where `0` suppresses all output, `1` only prints warnings, `2` +gives information at the start and end, and `3` prints information every iteration. +The projectors are computed from `svd_alg` SVDs where the truncation scheme is set via +`trscheme`. """ -@kwdef struct CTMRG - trscheme::TruncationScheme = TensorKit.notrunc() - tol::Float64 = Defaults.ctmrg_tol - maxiter::Int = Defaults.ctmrg_maxiter - miniter::Int = Defaults.ctmrg_miniter - verbosity::Int = 1 - fixedspace::Bool = false +struct CTMRG + tol::Float64 + maxiter::Int + miniter::Int + verbosity::Int + projector_alg::ProjectorAlg +end +function CTMRG(; + tol=Defaults.ctmrg_tol, + maxiter=Defaults.ctmrg_maxiter, + miniter=Defaults.ctmrg_miniter, + verbosity=0, + svd_alg=SVDAdjoint(), + trscheme=FixedSpaceTruncation(), +) + return CTMRG( + tol, maxiter, miniter, verbosity, ProjectorAlg(; svd_alg, trscheme, verbosity) + ) end """ @@ -30,7 +66,7 @@ Contract `state` using CTMRG and return the CTM environment. Per default, a random initial environment is used. """ function MPSKit.leading_boundary(state, alg::CTMRG) - return MPSKit.leading_boundary(CTMRGEnv(state), state, alg) + return MPSKit.leading_boundary(CTMRGEnv(state, oneunit(spacetype(state))), state, alg) end function MPSKit.leading_boundary(envinit, state, alg::CTMRG) CS = map(x -> tsvd(x; alg=TensorKit.SVD())[2], envinit.corners) @@ -48,16 +84,13 @@ function MPSKit.leading_boundary(envinit, state, alg::CTMRG) env, = ctmrg_iter(state, env, alg) # Grow and renormalize in all 4 directions η, CS, TS = calc_convergence(env, CS, TS) N = norm(state, env) - ctmrg_logiter!(log, iter, η, N) (iter > alg.miniter && η <= alg.tol) && break end # Do one final iteration that does not change the spaces - alg_fixed = CTMRG(; - alg.trscheme, alg.tol, alg.maxiter, alg.miniter, alg.verbosity, fixedspace=true - ) + alg_fixed = @set alg.projector_alg.trscheme = FixedSpaceTruncation() env′, = ctmrg_iter(state, env, alg_fixed) envfix = gauge_fix(env, env′) @@ -297,7 +330,7 @@ function ctmrg_iter(state, env::CTMRGEnv{C,T}, alg::CTMRG) where {C,T} ϵ = 0.0 for _ in 1:4 - env, _, _, ϵ₀ = left_move(state, env, alg) + env, _, _, ϵ₀ = left_move(state, env, alg.projector_alg) state = rotate_north(state, EAST) env = rotate_north(env, EAST) ϵ = max(ϵ, ϵ₀) @@ -312,7 +345,7 @@ end Grow, project and renormalize the environment `env` in west direction. Return the updated environment as well as the projectors and truncation error. """ -function left_move(state, env::CTMRGEnv{C,T}, alg::CTMRG) where {C,T} +function left_move(state, env::CTMRGEnv{C,T}, alg::ProjectorAlg) where {C,T} corners::typeof(env.corners) = copy(env.corners) edges::typeof(env.edges) = copy(env.edges) ϵ = 0.0 @@ -320,6 +353,7 @@ function left_move(state, env::CTMRGEnv{C,T}, alg::CTMRG) where {C,T} for col in 1:size(state, 2) cprev = _prev(col, size(state, 2)) + cnext = _next(col, size(state, 2)) # Compute projectors for row in 1:size(state, 1) @@ -341,15 +375,15 @@ function left_move(state, env::CTMRGEnv{C,T}, alg::CTMRG) where {C,T} ) # SVD half-infinite environment - trscheme = if alg.fixedspace == true - truncspace(space(env.edges[WEST, row, col], 1)) + trscheme = if alg.trscheme isa FixedSpaceTruncation + truncspace(space(env.edges[WEST, row, cnext], 1)) else alg.trscheme end @autoopt @tensor QQ[χ_EB D_EBabove D_EBbelow; χ_ET D_ETabove D_ETbelow] := Q_sw[χ_EB D_EBabove D_EBbelow; χ D1 D2] * Q_nw[χ D1 D2; χ_ET D_ETabove D_ETbelow] - U, S, V, ϵ_local = tsvd!(QQ; trunc=trscheme, alg=TensorKit.SVD()) + U, S, V, ϵ_local = PEPSKit.tsvd!(QQ, alg.svd_alg; trunc=trscheme) ϵ = max(ϵ, ϵ_local / norm(S)) # TODO: check if we can just normalize enlarged corners s.t. trunc behaves a bit better diff --git a/src/algorithms/peps_opt.jl b/src/algorithms/peps_opt.jl index f0351b8e..1e757dd1 100644 --- a/src/algorithms/peps_opt.jl +++ b/src/algorithms/peps_opt.jl @@ -62,7 +62,7 @@ in `alg`. The initial environment `env₀` serves as an initial guess for the fi By default, a random initial environment is used. """ function fixedpoint( - ψ₀::InfinitePEPS{T}, H, alg::PEPSOptimize, env₀::CTMRGEnv=CTMRGEnv(ψ₀; Venv=field(T)^20) + ψ₀::InfinitePEPS{T}, H, alg::PEPSOptimize, env₀::CTMRGEnv=CTMRGEnv(ψ₀, field(T)^20) ) where {T} (peps, env), E, ∂E, info = optimize( (ψ₀, env₀), alg.optimizer; retract=my_retract, inner=my_inner diff --git a/src/environments/ctmrgenv.jl b/src/environments/ctmrgenv.jl index fd2bf4c4..cf90e2f0 100644 --- a/src/environments/ctmrgenv.jl +++ b/src/environments/ctmrgenv.jl @@ -2,42 +2,333 @@ struct CTMRGEnv{C,T} Corner transfer-matrix environment containing unit-cell arrays of corner and edge tensors. +The last two indices of the arrays correspond to the row and column indices of the unit +cell, whereas the first index corresponds to the direction of the corner or edge tensor. The +directions are labeled in clockwise direction, starting from the north-west corner and north +edge respectively. + +Given arrays of corners `c` and edges `t`, they connect to the PEPS tensors at site `(r, c)` +in the unit cell as: +``` + c[1,r-1,c-1]---t[1,r-1,c]----c[2,r-1,c+1] + | || | + t[4,r,c-1]=====AA[r,c]=======t[2,r,c+1] + | || | + c[4,r+1,c-1]---t[3,r+1,c]----c[3,r+1,c+1] +``` + +# Fields +- `corners::Array{C,3}`: Array of corner tensors. +- `edges::Array{T,3}`: Array of edge tensors. """ struct CTMRGEnv{C,T} corners::Array{C,3} edges::Array{T,3} end +_spacetype(::Int) = ComplexSpace +_spacetype(::S) where {S<:ElementarySpace} = S + +_to_space(χ::Int) = ℂ^χ +_to_space(χ::ElementarySpace) = χ + +function _corner_tensor( + f, ::Type{T}, left_vspace::S, right_vspace::S=left_vspace +) where {T,S<:Union{Int,ElementarySpace}} + return TensorMap(f, T, _to_space(left_vspace) ← _to_space(right_vspace)) +end + +function _edge_tensor( + f, + ::Type{T}, + left_vspace::S, + top_pspace::S, + bot_pspace::S=top_pspace, + right_vspace::S=left_vspace, +) where {T,S<:Union{Int,ElementarySpace}} + return TensorMap( + f, + T, + _to_space(left_vspace) ⊗ _to_space(top_pspace) ⊗ dual(_to_space(bot_pspace)) ← + _to_space(right_vspace), + ) +end + """ - CTMRGEnv(peps::InfinitePEPS{P}; Venv=oneunit(spacetype(P))) + CTMRGEnv( + [f=randn, ComplexF64], Ds_north, Ds_east::A, chis_north::A, [chis_east::A], [chis_south::A], [chis_west::A] + ) where {A<:AbstractMatrix{<:Union{Int,ElementarySpace}}} + +Construct a CTMRG environment by specifying matrices of north and east virtual spaces of the +corresponding [`InfinitePEPS`](@ref) and the north, east, south and west virtual spaces of +the environment. Each respective matrix entry corresponds to a site in the unit cell. By +default, the virtual environment spaces for all directions are taken to be the same. -Create a random CTMRG environment from a PEPS tensor. The environment bond dimension -defaults to one and can be specified using the `Venv` space. +The environment virtual spaces for each site correspond to the north or east virtual space +of the corresponding edge tensor for each direction. Specifically, for a given site +`(r, c)`, `chis_north[r, c]` corresponds to the east space of the north edge tensor, +`chis_east[r, c]` corresponds to the north space of the east edge tensor, +`chis_south[r, c]` corresponds to the east space of the south edge tensor, and +`chis_west[r, c]` corresponds to the north space of the west edge tensor. """ -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)) +function CTMRGEnv( + Ds_north::A, + Ds_east::A, + chis_north::A, + chis_east::A=chis_north, + chis_south::A=chis_north, + chis_west::A=chis_north, +) where {A<:AbstractMatrix{<:Union{Int,ElementarySpace}}} + return CTMRGEnv( + randn, ComplexF64, Ds_north, Ds_east, chis_north, chis_east, chis_south, chis_west + ) +end +function CTMRGEnv( + f, + T, + Ds_north::A, + Ds_east::A, + chis_north::A, + chis_east::A=chis_north, + chis_south::A=chis_north, + chis_west::A=chis_north, +) where {A<:AbstractMatrix{<:Union{Int,ElementarySpace}}} + Ds_south = adjoint.(circshift(Ds_north, (-1, 0))) + Ds_west = adjoint.(circshift(Ds_east, (0, 1))) + + # do the whole thing + st = _spacetype(first(Ds_north)) + C_type = tensormaptype(st, 1, 1, T) + T_type = tensormaptype(st, 3, 1, T) # 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, 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), - Venv * space(peps[i, j], dir + 1)' * space(peps[i, j], dir + 1), - Venv, + corners = Array{C_type}(undef, 4, size(Ds_north)...) + edges = Array{T_type}(undef, 4, size(Ds_north)...) + + for (r, c) in Iterators.product(axes(Ds_north)...) + edges[NORTH, r, c] = _edge_tensor( + f, + T, + chis_north[r, _prev(c, end)], + Ds_north[_next(r, end), c], + Ds_north[_next(r, end), c], + chis_north[r, c], ) + edges[EAST, r, c] = _edge_tensor( + f, + T, + chis_east[r, c], + Ds_east[r, _prev(c, end)], + Ds_east[r, _prev(c, end)], + chis_east[_next(r, end), c], + ) + edges[SOUTH, r, c] = _edge_tensor( + f, + T, + chis_south[r, c], + Ds_south[_prev(r, end), c], + Ds_south[_prev(r, end), c], + chis_south[r, _prev(c, end)], + ) + edges[WEST, r, c] = _edge_tensor( + f, + T, + chis_west[_next(r, end), c], + Ds_west[r, _next(c, end)], + Ds_west[r, _next(c, end)], + chis_west[r, c], + ) + + corners[NORTHWEST, r, c] = _corner_tensor( + f, T, chis_west[_next(r, end), c], chis_north[r, c] + ) + corners[NORTHEAST, r, c] = _corner_tensor( + f, T, chis_north[r, _prev(c, end)], chis_east[_next(r, end), c] + ) + corners[SOUTHEAST, r, c] = _corner_tensor( + f, T, chis_east[r, c], chis_south[r, _prev(c, end)] + ) + corners[SOUTHWEST, r, c] = _corner_tensor(f, T, chis_south[r, c], chis_west[r, c]) end - @diffset corners[:, :, :] ./= norm.(corners[:, :, :]) - @diffset edges[:, :, :] ./= norm.(edges[:, :, :]) + corners[:, :, :] ./= norm.(corners[:, :, :]) + edges[:, :, :] ./= norm.(edges[:, :, :]) return CTMRGEnv(corners, edges) end +""" + CTMRGEnv( + [f=randn, ComplexF64], D_north::S, D_south::S, chi_north::S, [chi_east::S], [chi_south::S], [chi_west::S]; unitcell::Tuple{Int,Int}=(1, 1), + ) where {S<:Union{Int,ElementarySpace}} + +Construct a CTMRG environment by specifying the north and east virtual spaces of the +corresponding [`InfinitePEPS`](@ref) and the north, east, south and west virtual spaces of +the environment. The PEPS unit cell can be specified by the `unitcell` keyword argument. By +default, the virtual environment spaces for all directions are taken to be the same. + +The environment virtual spaces for each site correspond to virtual space of the +corresponding edge tensor for each direction. +""" +function CTMRGEnv( + D_north::S, + D_south::S, + chi_north::S, + chi_east::S=chi_north, + chi_south::S=chi_north, + chi_west::S=chi_north; + unitcell::Tuple{Int,Int}=(1, 1), +) where {S<:Union{Int,ElementarySpace}} + return CTMRGEnv( + randn, + ComplexF64, + fill(D_north, unitcell), + fill(D_south, unitcell), + fill(chi_north, unitcell), + fill(chi_east, unitcell), + fill(chi_south, unitcell), + fill(chi_west, unitcell), + ) +end +function CTMRGEnv( + f, + T, + D_north::S, + D_south::S, + chi_north::S, + chi_east::S=chi_north, + chi_south::S=chi_north, + chi_west::S=chi_north; + unitcell::Tuple{Int,Int}=(1, 1), +) where {S<:Union{Int,ElementarySpace}} + return CTMRGEnv( + f, + T, + fill(D_north, unitcell), + fill(D_south, unitcell), + fill(chi_north, unitcell), + fill(chi_east, unitcell), + fill(chi_south, unitcell), + fill(chi_west, unitcell), + ) +end + +""" + CTMRGEnv( + [f=randn, T=ComplexF64], peps::InfinitePEPS, chis_north::A, [chis_east::A], [chis_south::A], [chis_west::A] + ) where {A<:AbstractMatrix{<:Union{Int,ElementarySpace}}} + +Construct a CTMRG environment by specifying a corresponding [`InfinitePEPS`](@ref), and the +north, east, south and west virtual spaces of the environment as matrices. Each respective +matrix entry corresponds to a site in the unit cell. By default, the virtual spaces for all +directions are taken to be the same. + +The environment virtual spaces for each site correspond to the north or east virtual space +of the corresponding edge tensor for each direction. Specifically, for a given site +`(r, c)`, `chis_north[r, c]` corresponds to the east space of the north edge tensor, +`chis_east[r, c]` corresponds to the north space of the east edge tensor, +`chis_south[r, c]` corresponds to the east space of the south edge tensor, and +`chis_west[r, c]` corresponds to the north space of the west edge tensor. +""" +function CTMRGEnv( + peps::InfinitePEPS, + chis_north::A, + chis_east::A=chis_north, + chis_south::A=chis_north, + chis_west::A=chis_north, +) where {A<:AbstractMatrix{<:Union{Int,ElementarySpace}}} + Ds_north = map(peps.A) do t + return adjoint(space(t, 2)) + end + Ds_east = map(peps.A) do t + return adjoint(space(t, 3)) + end + return CTMRGEnv( + randn, + ComplexF64, + Ds_north, + Ds_east, + _to_space.(chis_north), + _to_space.(chis_east), + _to_space.(chis_south), + _to_space.(chis_west), + ) +end +function CTMRGEnv( + f, + T, + peps::InfinitePEPS, + chis_north::A, + chis_east::A=chis_north, + chis_south::A=chis_north, + chis_west::A=chis_north, +) where {A<:AbstractMatrix{<:Union{Int,ElementarySpace}}} + Ds_north = map(peps.A) do t + return adjoint(space(t, 2)) + end + Ds_east = map(peps.A) do t + return adjoint(space(t, 3)) + end + return CTMRGEnv( + f, + T, + Ds_north, + Ds_east, + _to_space.(chis_north), + _to_space.(chis_east), + _to_space.(chis_south), + _to_space.(chis_west), + ) +end + +""" + CTMRGEnv( + peps::InfinitePEPS, chi_north::S, [chi_east::S], [chi_south::S], [chi_west::S], + ) where {S<:Union{Int,ElementarySpace}} + +Construct a CTMRG environment by specifying a corresponding [`InfinitePEPS`](@ref), and the +north, east, south and west virtual spaces of the environment. By default, the virtual +spaces for all directions are taken to be the same. + +The environment virtual spaces for each site correspond to virtual space of the +corresponding edge tensor for each direction. +""" +function CTMRGEnv( + peps::InfinitePEPS, + chi_north::S, + chi_east::S=chi_north, + chi_south::S=chi_north, + chi_west::S=chi_north, +) where {S<:Union{Int,ElementarySpace}} + return CTMRGEnv( + peps, + fill(chi_north, size(peps)), + fill(chi_east, size(peps)), + fill(chi_south, size(peps)), + fill(chi_west, size(peps)), + ) +end +function CTMRGEnv( + f, + T, + peps::InfinitePEPS, + chi_north::S, + chi_east::S=chi_north, + chi_south::S=chi_north, + chi_west::S=chi_north, +) where {S<:Union{Int,ElementarySpace}} + return CTMRGEnv( + f, + T, + peps, + fill(chi_north, size(peps)), + fill(chi_east, size(peps)), + fill(chi_south, size(peps)), + fill(chi_west, size(peps)), + ) +end +@non_differentiable CTMRGEnv(peps::InfinitePEPS, args...) + # Custom adjoint for CTMRGEnv constructor, needed for fixed-point differentiation function ChainRulesCore.rrule(::Type{CTMRGEnv}, corners, edges) ctmrgenv_pullback(ē) = NoTangent(), ē.corners, ē.edges diff --git a/src/states/abstractpeps.jl b/src/states/abstractpeps.jl index 467b05c3..93a27142 100644 --- a/src/states/abstractpeps.jl +++ b/src/states/abstractpeps.jl @@ -2,7 +2,8 @@ const PEPSTensor{S} Default type for PEPS tensors with a single physical index, and 4 virtual indices, -conventionally ordered as: T : P ← N ⊗ E ⊗ S ⊗ W. +conventionally ordered as: ``T : P ← N ⊗ E ⊗ S ⊗ W``. Here, ``P``, ``N``, ``E``, ``S`` and +``W`` denote the physics, north, east, south and west spaces, respectively. """ const PEPSTensor{S} = AbstractTensorMap{S,1,4} where {S<:ElementarySpace} diff --git a/src/states/infinitepeps.jl b/src/states/infinitepeps.jl index 4135e8f1..2661c072 100644 --- a/src/states/infinitepeps.jl +++ b/src/states/infinitepeps.jl @@ -21,18 +21,23 @@ end ## Constructors """ - InfinitePEPS(A::AbstractArray{T, 2}) + InfinitePEPS(A::AbstractMatrix{T}) -Allow users to pass in an array of tensors. +Create an `InfinitePEPS` by specifying a matrix containing the PEPS tensors at each site in +the unit cell. """ -function InfinitePEPS(A::AbstractArray{T,2}) where {T<:PEPSTensor} +function InfinitePEPS(A::AbstractMatrix{T}) where {T<:PEPSTensor} return InfinitePEPS(Array(deepcopy(A))) # TODO: find better way to copy end """ - InfinitePEPS(f=randn, T=ComplexF64, Pspaces, Nspaces, Espaces) + InfinitePEPS( + f=randn, T=ComplexF64, Pspaces::A, Nspaces::A, [Espaces::A] + ) where {A<:AbstractMatrix{<:Union{Int,ElementarySpace}}} -Allow users to pass in arrays of spaces. +Create an `InfinitePEPS` by specifying the physical, north virtual and east virtual spaces +of the PEPS tensor at each site in the unit cell as a matrix. Each individual space can be +specified as either an `Int` or an `ElementarySpace`. """ function InfinitePEPS( Pspaces::A, Nspaces::A, Espaces::A @@ -82,8 +87,8 @@ end """ 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`. +Create an InfinitePEPS by specifying its physical, north and east spaces and unit cell. +Spaces can be specified either via `Int` or via `ElementarySpace`. """ function InfinitePEPS( Pspace::S, Nspace::S, Espace::S=Nspace; unitcell::Tuple{Int,Int}=(1, 1) diff --git a/src/utility/eigsolve.jl b/src/utility/eigsolve.jl deleted file mode 100644 index 975adf49..00000000 --- a/src/utility/eigsolve.jl +++ /dev/null @@ -1,251 +0,0 @@ -# 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/svd.jl b/src/utility/svd.jl new file mode 100644 index 00000000..740892ba --- /dev/null +++ b/src/utility/svd.jl @@ -0,0 +1,268 @@ +using TensorKit: + SectorDict, + _tsvd!, + _empty_svdtensors, + _compute_svddata!, + _create_svdtensors, + NoTruncation, + TruncationSpace + +const CRCExt = Base.get_extension(KrylovKit, :KrylovKitChainRulesCoreExt) + +""" + struct SVDAdjoint(; fwd_alg = TensorKit.SVD(), rrule_alg = nothing, + broadening = nothing) + +Wrapper for a SVD algorithm `fwd_alg` with a defined reverse rule `rrule_alg`. +If `isnothing(rrule_alg)`, Zygote differentiates the forward call automatically. +In case of degenerate singular values, one might need a `broadening` scheme which +removes the divergences from the adjoint. +""" +@kwdef struct SVDAdjoint{F,R,B} + fwd_alg::F = TensorKit.SVD() + rrule_alg::R = nothing + broadening::B = nothing +end # Keep truncation algorithm separate to be able to specify CTMRG dependent information + +""" + PEPSKit.tsvd(t::AbstractTensorMap, alg; trunc=notrunc(), p=2) + +Wrapper around `TensorKit.tsvd` which dispatches on the `alg` argument. +This is needed since a custom adjoint for `PEPSKit.tsvd` may be defined, +depending on the algorithm. E.g., for `IterSVD` the adjoint for a truncated +SVD from `KrylovKit.svdsolve` is used. +""" +PEPSKit.tsvd(t::AbstractTensorMap, alg; kwargs...) = PEPSKit.tsvd!(copy(t), alg; kwargs...) +function PEPSKit.tsvd!( + t::AbstractTensorMap, alg::SVDAdjoint; trunc::TruncationScheme=notrunc(), p::Real=2 +) + return TensorKit.tsvd!(t; alg=alg.fwd_alg, trunc, p) +end + +""" + struct IterSVD(; alg = KrylovKit.GKL(), fallback_threshold = Inf) + +Iterative SVD solver based on KrylovKit's GKL algorithm, adapted to (symmmetric) tensors. +The number of targeted singular values is set via the `TruncationSpace` in `ProjectorAlg`. +In particular, this make it possible to specify the targeted singular values block-wise. +In case the symmetry block is too small as compared to the number of singular values, or +the iterative SVD didn't converge, the algorithm falls back to a dense SVD. +""" +@kwdef struct IterSVD + alg::KrylovKit.GKL = KrylovKit.GKL(; tol=1e-14, krylovdim=25) + fallback_threshold::Float64 = Inf +end + +# Compute SVD data block-wise using KrylovKit algorithm +function TensorKit._tsvd!( + t, alg::IterSVD, trunc::Union{NoTruncation,TruncationSpace}, p::Real=2 +) + # early return + if isempty(blocksectors(t)) + truncerr = zero(real(scalartype(t))) + return _empty_svdtensors(t)..., truncerr + end + + Udata, Σdata, Vdata, dims = _compute_svddata!(t, alg, trunc) + U, S, V = _create_svdtensors(t, Udata, Σdata, Vdata, spacetype(t)(dims)) + truncerr = trunc isa NoTruncation ? abs(zero(scalartype(t))) : norm(U * S * V - t, p) + + return U, S, V, truncerr +end +function TensorKit._compute_svddata!( + t::TensorMap, alg::IterSVD, trunc::Union{NoTruncation,TruncationSpace} +) + InnerProductStyle(t) === EuclideanProduct() || throw_invalid_innerproduct(:tsvd!) + I = sectortype(t) + A = storagetype(t) + Udata = SectorDict{I,A}() + Vdata = SectorDict{I,A}() + dims = SectorDict{I,Int}() + local Sdata + for (c, b) in blocks(t) + howmany = trunc isa NoTruncation ? minimum(size(b)) : blockdim(trunc.space, c) + + if howmany / minimum(size(b)) > alg.fallback_threshold # Use dense SVD for small blocks + U, S, V = TensorKit.MatrixAlgebra.svd!(b, TensorKit.SVD()) + Udata[c] = U[:, 1:howmany] + Vdata[c] = V[1:howmany, :] + else + x₀ = randn(eltype(b), size(b, 1)) + S, lvecs, rvecs, info = KrylovKit.svdsolve(b, x₀, howmany, :LR, alg.alg) + if info.converged < howmany # Fall back to dense SVD if not properly converged + @warn "Iterative SVD did not converge for block $c, falling back to dense SVD" + U, S, V = TensorKit.MatrixAlgebra.svd!(b, TensorKit.SVD()) + Udata[c] = U[:, 1:howmany] + Vdata[c] = V[1:howmany, :] + else # Slice in case more values were converged than requested + Udata[c] = stack(view(lvecs, 1:howmany)) + Vdata[c] = stack(conj, view(rvecs, 1:howmany); dims=1) + end + end + + resize!(S, howmany) + if @isdefined Sdata + Sdata[c] = S + else + Sdata = SectorDict(c => S) + end + dims[c] = length(S) + end + return Udata, Sdata, Vdata, dims +end + +function ChainRulesCore.rrule( + ::typeof(PEPSKit.tsvd!), + t::AbstractTensorMap, + alg::SVDAdjoint{IterSVD,R,B}; + trunc::TruncationScheme=notrunc(), + p::Real=2, +) where {R,B} + U, S, V, ϵ = PEPSKit.tsvd(t, alg; trunc, p) + + function tsvd!_itersvd_pullback((ΔU, ΔS, ΔV, Δϵ)) + Δt = similar(t) + for (c, b) in blocks(Δt) + Uc, Sc, Vc = block(U, c), block(S, c), block(V, c) + ΔUc, ΔSc, ΔVc = block(ΔU, c), block(ΔS, c), block(ΔV, c) + Sdc = view(Sc, diagind(Sc)) + ΔSdc = ΔSc isa AbstractZero ? ΔSc : view(ΔSc, diagind(ΔSc)) + + n_vals = length(Sdc) + lvecs = Vector{Vector{scalartype(t)}}(eachcol(Uc)) + rvecs = Vector{Vector{scalartype(t)}}(eachcol(Vc')) + minimal_info = KrylovKit.ConvergenceInfo(length(Sdc), nothing, nothing, -1, -1) # Just supply converged to SVD pullback + + if ΔUc isa AbstractZero && ΔVc isa AbstractZero # Handle ZeroTangent singular vectors + Δlvecs = fill(ZeroTangent(), n_vals) + Δrvecs = fill(ZeroTangent(), n_vals) + else + Δlvecs = Vector{Vector{scalartype(t)}}(eachcol(ΔUc)) + Δrvecs = Vector{Vector{scalartype(t)}}(eachcol(ΔVc')) + end + + xs, ys = CRCExt.compute_svdsolve_pullback_data( + ΔSc isa AbstractZero ? fill(zero(Sc[1]), n_vals) : ΔSdc, + Δlvecs, + Δrvecs, + Sdc, + lvecs, + rvecs, + minimal_info, + block(t, c), + :LR, + alg.fwd_alg.alg, + alg.rrule_alg, + ) + copyto!( + b, + CRCExt.construct∂f_svd(HasReverseMode(), block(t, c), lvecs, rvecs, xs, ys), + ) + end + return NoTangent(), Δt, NoTangent() + end + function tsvd!_itersvd_pullback(::Tuple{ZeroTangent,ZeroTangent,ZeroTangent}) + return NoTangent(), ZeroTangent(), NoTangent() + end + + return (U, S, V, ϵ), tsvd!_itersvd_pullback +end + +""" + struct NonTruncAdjoint(; lorentz_broadening = 0.0) + +Old SVD adjoint that does not account for the truncated part of truncated SVDs. +""" +struct NonTruncSVDAdjoint end + +# Use outdated adjoint in reverse pass (not taking truncated part into account for testing purposes) +function ChainRulesCore.rrule( + ::typeof(PEPSKit.tsvd!), + t::AbstractTensorMap, + alg::SVDAdjoint{F,NonTruncSVDAdjoint,B}; + trunc::TruncationScheme=notrunc(), + p::Real=2, +) where {F,B} + U, S, V, ϵ = PEPSKit.tsvd(t, alg; trunc, p) + + function tsvd!_nontruncsvd_pullback((ΔU, ΔS, ΔV, Δϵ)) + Δt = similar(t) + for (c, b) in blocks(Δt) + Uc, Sc, Vc = block(U, c), block(S, c), block(V, c) + ΔUc, ΔSc, ΔVc = block(ΔU, c), block(ΔS, c), block(ΔV, c) + copyto!( + b, oldsvd_rev(Uc, Sc, Vc, ΔUc, ΔSc, ΔVc; lorentz_broadening=alg.broadening) + ) + end + return NoTangent(), Δt, NoTangent() + end + function tsvd!_nontruncsvd_pullback(::Tuple{ZeroTangent,ZeroTangent,ZeroTangent}) + return NoTangent(), ZeroTangent(), NoTangent() + end + + return (U, S, V, ϵ), tsvd!_nontruncsvd_pullback +end + +function oldsvd_rev( + U::AbstractMatrix, + S::AbstractMatrix, + V::AbstractMatrix, + ΔU, + ΔS, + ΔV; + lorentz_broadening=0, + atol::Real=0, + rtol::Real=atol > 0 ? 0 : eps(scalartype(S))^(3 / 4), +) + tol = atol > 0 ? atol : rtol * S[1, 1] + F = _invert_S²(S, tol, lorentz_broadening) # Includes Lorentzian broadening + S⁻¹ = pinv(S; atol=tol) + + # dS contribution + term = ΔS isa ZeroTangent ? ΔS : Diagonal(diag(ΔS)) + + # dU₁ and dV₁ off-diagonal contribution + J = F .* (U' * ΔU) + term += (J + J') * S + VΔV = (V * ΔV') + K = F .* VΔV + term += S * (K + K') + + # dV₁ diagonal contribution (diagonal of dU₁ is gauged away) + if scalartype(U) <: Complex && !(ΔV isa ZeroTangent) && !(ΔU isa ZeroTangent) + L = Diagonal(diag(VΔV)) + term += 0.5 * S⁻¹ * (L' - L) + end + ΔA = U * term * V + + # Projector contribution for non-square A + UUd = U * U' + VdV = V' * V + Uproj = one(UUd) - UUd + Vproj = one(VdV) - VdV + ΔA += Uproj * ΔU * S⁻¹ * V + U * S⁻¹ * ΔV * Vproj # Wrong truncation contribution + + return ΔA +end + +# Computation of F in SVD adjoint, including Lorentzian broadening +function _invert_S²(S::AbstractMatrix{T}, tol::Real, ε=0) where {T<:Real} + F = similar(S) + @inbounds for i in axes(F, 1), j in axes(F, 2) + F[i, j] = if i == j + zero(T) + else + sᵢ, sⱼ = S[i, i], S[j, j] + Δs = abs(sⱼ - sᵢ) < tol ? tol : sⱼ^2 - sᵢ^2 + ε > 0 && (Δs = _lorentz_broaden(Δs, ε)) + 1 / Δs + end + end + return F +end + +# Lorentzian broadening for SVD adjoint F-singularities +function _lorentz_broaden(x::Real, ε=1e-12) + x′ = 1 / x + return x′ / (x′^2 + ε) +end diff --git a/test/boundarymps/vumps.jl b/test/boundarymps/vumps.jl index bbf69006..a161c073 100644 --- a/test/boundarymps/vumps.jl +++ b/test/boundarymps/vumps.jl @@ -17,9 +17,7 @@ const vumps_alg = VUMPS(; alg_eigsolve=MPSKit.Defaults.alg_eigsolve(; ishermitia mps, envs, ϵ = leading_boundary(mps, T, vumps_alg) N = abs(sum(expectation_value(mps, T))) - ctm = leading_boundary( - CTMRGEnv(psi; Venv=ComplexSpace(20)), psi, CTMRG(; verbosity=1, fixedspace=true) - ) + ctm = leading_boundary(CTMRGEnv(psi, ComplexSpace(20)), psi, CTMRG(; verbosity=1)) N´ = abs(norm(psi, ctm)) @test N ≈ N´ atol = 1e-3 @@ -33,12 +31,10 @@ end mps, envs, ϵ = leading_boundary(mps, T, vumps_alg) N = abs(prod(expectation_value(mps, T))) - ctm = leading_boundary( - CTMRGEnv(psi; Venv=ComplexSpace(20)), psi, CTMRG(; verbosity=1, fixedspace=true) - ) + ctm = leading_boundary(CTMRGEnv(psi, ComplexSpace(20)), psi, CTMRG(; verbosity=1)) N´ = abs(norm(psi, ctm)) - @test N ≈ N´ rtol = 1e-3 + @test N ≈ N´ rtol = 1e-2 end @testset "PEPO runthrough" begin diff --git a/test/ctmrg/gaugefix.jl b/test/ctmrg/gaugefix.jl index 1d26d04a..54bd1058 100644 --- a/test/ctmrg/gaugefix.jl +++ b/test/ctmrg/gaugefix.jl @@ -2,6 +2,7 @@ using Test using Random using PEPSKit using TensorKit +using Accessors using PEPSKit: ctmrg_iter, gauge_fix, calc_elementwise_convergence @@ -43,12 +44,12 @@ end psi = _make_symmetric(psi) Random.seed!(987654321) # Seed RNG to make random environment consistent - ctm = CTMRGEnv(psi; Venv=ctm_space) + ctm = CTMRGEnv(psi, ctm_space) alg = CTMRG(; - trscheme=truncdim(dim(ctm_space)), tol=1e-10, miniter=4, maxiter=400, verbosity + tol=1e-10, miniter=4, maxiter=400, verbosity, trscheme=truncdim(dim(ctm_space)) ) - alg_fixed = CTMRG(; trscheme=truncdim(dim(ctm_space)), verbosity, fixedspace=true) + alg_fixed = @set alg.projector_alg.trscheme = FixedSpaceTruncation() ctm = leading_boundary(ctm, psi, alg) ctm2, = ctmrg_iter(psi, ctm, alg_fixed) @@ -67,12 +68,12 @@ end psi = _make_symmetric(psi) Random.seed!(123456789) # Seed RNG to make random environment consistent - ctm = CTMRGEnv(psi; Venv=ctm_space) + ctm = CTMRGEnv(psi, ctm_space) alg = CTMRG(; - trscheme=truncspace(ctm_space), tol=1e-10, miniter=4, maxiter=400, verbosity + tol=1e-10, miniter=4, maxiter=400, verbosity, trscheme=truncdim(dim(ctm_space)) ) - alg_fixed = CTMRG(; trscheme=truncspace(ctm_space), verbosity, fixedspace=true) + alg_fixed = @set alg.projector_alg.trscheme = FixedSpaceTruncation() ctm = leading_boundary(ctm, psi, alg) ctm2, = ctmrg_iter(psi, ctm, alg_fixed) diff --git a/test/ctmrg/gradients.jl b/test/ctmrg/gradients.jl index 1c81063a..cd4ed24b 100644 --- a/test/ctmrg/gradients.jl +++ b/test/ctmrg/gradients.jl @@ -17,9 +17,7 @@ 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 -) +boundary_alg = CTMRG(; tol=tol, miniter=4, maxiter=100, verbosity=0) gradmodes = [ nothing, GeomSum(; tol), ManualIter(; tol), KrylovKit.GMRES(; tol=tol, maxiter=10) ] @@ -38,7 +36,7 @@ steps = -0.01:0.005:0.01 @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) + env = leading_boundary(CTMRGEnv(psi, Espace), psi, boundary_alg) alphas, fs, dfs1, dfs2 = OptimKit.optimtest( (psi, env), dir; diff --git a/test/ctmrg/gradparts.jl b/test/ctmrg/gradparts.jl index c21a2420..406aa8bc 100644 --- a/test/ctmrg/gradparts.jl +++ b/test/ctmrg/gradparts.jl @@ -30,9 +30,7 @@ Vspaces = [ComplexSpace(χbond), Vect[FermionParity](0 => χbond / 2, 1 => χbon 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 -) +boundary_alg = CTMRG(; tol=tol, miniter=4, maxiter=100, verbosity=0) ## Gauge invariant function of the environment # -------------------------------------------- @@ -51,12 +49,10 @@ end ## Tests # ------ -@testset "Reverse rules for composite parts of the CTMRG fixed point with spacetype $(Vspaces[i])" for i in - eachindex( - Pspaces -) +title = "Reverse rules for composite parts of the CTMRG fixed point with spacetype" +@testset title * "$(Vspaces[i])" for i in eachindex(Pspaces) psi = InfinitePEPS(Pspaces[i], Vspaces[i], Vspaces[i]) - env = CTMRGEnv(psi; Venv=Espaces[i]) + env = CTMRGEnv(psi, Espaces[i]) @testset "$f" for f in functions atol = f == leading_boundary ? sqrt(tol) : tol diff --git a/test/ctmrg/svd_wrapper.jl b/test/ctmrg/svd_wrapper.jl new file mode 100644 index 00000000..3c17dc7d --- /dev/null +++ b/test/ctmrg/svd_wrapper.jl @@ -0,0 +1,93 @@ +using Test +using Random +using LinearAlgebra +using TensorKit +using KrylovKit +using ChainRulesCore, Zygote +using Accessors +using PEPSKit + +# Gauge-invariant loss function +function lossfun(A, alg, R=TensorMap(randn, space(A)), trunc=notrunc()) + U, _, V, = PEPSKit.tsvd(A, alg; trunc) + return real(dot(R, U * V)) # Overlap with random tensor R is gauge-invariant and differentiable, also for m≠n +end + +m, n = 20, 30 +dtype = ComplexF64 +χ = 12 +trunc = truncspace(ℂ^χ) +# lorentz_broadening = 1e-12 +rtol = 1e-9 +r = TensorMap(randn, dtype, ℂ^m, ℂ^n) +R = TensorMap(randn, space(r)) + +full_alg = SVDAdjoint(; fwd_alg=TensorKit.SVD(), rrule_alg=nothing) +old_alg = SVDAdjoint(; + fwd_alg=TensorKit.SVD(), rrule_alg=NonTruncSVDAdjoint(), broadening=0.0 +) +iter_alg = SVDAdjoint(; fwd_alg=IterSVD(), rrule_alg=GMRES(; tol=1e-13)) # Don't make adjoint tolerance too small, g_itersvd will be weird + +@testset "Non-truncacted SVD" begin + l_fullsvd, g_fullsvd = withgradient(A -> lossfun(A, full_alg, R), r) + l_oldsvd, g_oldsvd = withgradient(A -> lossfun(A, old_alg, R), r) + l_itersvd, g_itersvd = withgradient(A -> lossfun(A, iter_alg, R), r) + + @test l_oldsvd ≈ l_itersvd ≈ l_fullsvd + @test g_fullsvd[1] ≈ g_oldsvd[1] rtol = rtol + @test g_fullsvd[1] ≈ g_itersvd[1] rtol = rtol +end + +@testset "Truncated SVD with χ=$χ" begin + l_fullsvd, g_fullsvd = withgradient(A -> lossfun(A, full_alg, R, trunc), r) + l_oldsvd, g_oldsvd = withgradient(A -> lossfun(A, old_alg, R, trunc), r) + l_itersvd, g_itersvd = withgradient(A -> lossfun(A, iter_alg, R, trunc), r) + + @test l_oldsvd ≈ l_itersvd ≈ l_fullsvd + @test !isapprox(g_fullsvd[1], g_oldsvd[1]; rtol) + @test g_fullsvd[1] ≈ g_itersvd[1] rtol = rtol +end + +# TODO: Add when Lorentzian broadening is implemented +# @testset "Truncated SVD with χ=$χ and ε=$lorentz_broadening broadening" begin +# l_fullsvd, g_fullsvd = withgradient( +# A -> lossfun(A, FullSVD(; lorentz_broadening, R; trunc), r +# ) +# l_oldsvd, g_oldsvd = withgradient(A -> lossfun(A, OldSVD(; lorentz_broadening), R; trunc), r) +# l_itersvd, g_itersvd = withgradient( +# A -> lossfun(A, IterSVD(; howmany=χ, lorentz_broadening), R; trunc), r +# ) + +# @test l_oldsvd ≈ l_itersvd ≈ l_fullsvd +# @test norm(g_fullsvd[1] - g_oldsvd[1]) / norm(g_fullsvd[1]) > rtol +# @test norm(g_fullsvd[1] - g_itersvd[1]) / norm(g_fullsvd[1]) < rtol +# end + +symm_m, symm_n = 18, 24 +symm_space = Z2Space(0 => symm_m, 1 => symm_n) +symm_trspace = truncspace(Z2Space(0 => symm_m ÷ 2, 1 => symm_n ÷ 3)) +symm_r = TensorMap(randn, dtype, symm_space, symm_space) +symm_R = TensorMap(randn, dtype, space(symm_r)) + +@testset "IterSVD of symmetric tensors" begin + l_fullsvd, g_fullsvd = withgradient(A -> lossfun(A, full_alg, symm_R), symm_r) + l_itersvd, g_itersvd = withgradient(A -> lossfun(A, iter_alg, symm_R), symm_r) + @test l_itersvd ≈ l_fullsvd + @test g_fullsvd[1] ≈ g_itersvd[1] rtol = rtol + + l_fullsvd_tr, g_fullsvd_tr = withgradient( + A -> lossfun(A, full_alg, symm_R, symm_trspace), symm_r + ) + l_itersvd_tr, g_itersvd_tr = withgradient( + A -> lossfun(A, iter_alg, symm_R, symm_trspace), symm_r + ) + @test l_itersvd_tr ≈ l_fullsvd_tr + @test g_fullsvd_tr[1] ≈ g_itersvd_tr[1] rtol = rtol + + iter_alg_fallback = @set iter_alg.fwd_alg.fallback_threshold = 0.4 # Do dense SVD in one block, sparse SVD in the other + l_itersvd_fb, g_itersvd_fb = withgradient( + A -> lossfun(A, iter_alg_fallback, symm_R, symm_trspace), symm_r + ) + @test l_itersvd_fb ≈ l_fullsvd_tr + @test g_fullsvd_tr[1] ≈ g_itersvd_fb[1] rtol = rtol +end diff --git a/test/ctmrg/unitcell.jl b/test/ctmrg/unitcell.jl index 9d8cad8a..0ffbdc09 100644 --- a/test/ctmrg/unitcell.jl +++ b/test/ctmrg/unitcell.jl @@ -4,48 +4,85 @@ using PEPSKit using PEPSKit: _prev, _next, ctmrg_iter using TensorKit -# initialize (3, 3) PEPS with unique N and E bond spaces +# settings Random.seed!(91283219347) -unitcell = (3, 3) -Pspace = ℂ^2 -Nspace = ℂ^2 -Espace = ℂ^3 stype = ComplexF64 -peps = InfinitePEPS(randn, stype, Pspace, Nspace, Espace; unitcell) - -# initialize (3, 3) environment with unique environment spaces -C_type = tensormaptype(spacetype(peps[1]), 1, 1, storagetype(peps[1])) -T_type = tensormaptype(spacetype(peps[1]), 3, 1, storagetype(peps[1])) -corners = Array{C_type}(undef, 4, unitcell...) -edges = Array{T_type}(undef, 4, unitcell...) -for r in 1:unitcell[1], c in 1:unitcell[2] - corners[1, _prev(r, end), _prev(c, end)] = TensorMap(randn, stype, ℂ^5, ℂ^2) - edges[1, _prev(r, end), c] = TensorMap( - randn, stype, ℂ^2 * space(peps[r, c], 1 + 1)' * space(peps[r, c], 1 + 1), ℂ^2 - ) - corners[2, _prev(r, end), _next(c, end)] = TensorMap(randn, stype, ℂ^2, ℂ^3) - edges[2, r, _next(c, end)] = TensorMap( - randn, stype, ℂ^3 * space(peps[r, c], 2 + 1)' * space(peps[r, c], 2 + 1), ℂ^3 +ctm_alg = CTMRG() + +function test_unitcell( + unitcell, Pspaces, Nspaces, Espaces, chis_north, chis_east, chis_south, chis_west +) + peps = InfinitePEPS(randn, stype, Pspaces, Nspaces, Espaces) + env = CTMRGEnv(randn, stype, peps, chis_north, chis_east, chis_south, chis_west) + + # apply one CTMRG iteration with fixeds + env′, = ctmrg_iter(peps, env, ctm_alg) + + # compute random expecation value to test matching bonds + random_op = LocalOperator( + PEPSKit._to_space.(Pspaces), + [ + (c,) => TensorMap( + randn, + scalartype(peps), + PEPSKit._to_space(Pspaces[c]), + PEPSKit._to_space(Pspaces[c]), + ) for c in CartesianIndices(unitcell) + ]..., ) - corners[3, _next(r, end), _next(c, end)] = TensorMap(randn, stype, ℂ^3, ℂ^4) - edges[3, _next(r, end), c] = TensorMap( - randn, stype, ℂ^4 * space(peps[r, c], 3 + 1)' * space(peps[r, c], 3 + 1), ℂ^4 + @test expectation_value(peps, random_op, env) isa Number + @test expectation_value(peps, random_op, env′) isa Number +end + +function random_dualize!(M::AbstractMatrix{<:ElementarySpace}) + mask = rand([true, false], size(M)) + M[mask] .= adjoint.(M[mask]) + return M +end + +@testset "Integer space specifiers" begin + unitcell = (3, 3) + + Pspaces = rand(2:3, unitcell...) + Nspaces = rand(2:4, unitcell...) + Espaces = rand(2:4, unitcell...) + chis_north = rand(5:10, unitcell...) + chis_east = rand(5:10, unitcell...) + chis_south = rand(5:10, unitcell...) + chis_west = rand(5:10, unitcell...) + + test_unitcell( + unitcell, Pspaces, Nspaces, Espaces, chis_north, chis_east, chis_south, chis_west ) - corners[4, _next(r, end), _prev(c, end)] = TensorMap(randn, stype, ℂ^4, ℂ^5) - edges[4, r, _prev(c, end)] = TensorMap( - randn, stype, ℂ^5 * space(peps[r, c], 4 + 1)' * space(peps[r, c], 4 + 1), ℂ^5 +end + +@testset "Random Cartesian spaces" begin + unitcell = (3, 3) + + Pspaces = random_dualize!(ComplexSpace.(rand(2:3, unitcell...))) + Nspaces = random_dualize!(ComplexSpace.(rand(2:4, unitcell...))) + Espaces = random_dualize!(ComplexSpace.(rand(2:4, unitcell...))) + chis_north = random_dualize!(ComplexSpace.(rand(5:10, unitcell...))) + chis_east = random_dualize!(ComplexSpace.(rand(5:10, unitcell...))) + chis_south = random_dualize!(ComplexSpace.(rand(5:10, unitcell...))) + chis_west = random_dualize!(ComplexSpace.(rand(5:10, unitcell...))) + + test_unitcell( + unitcell, Pspaces, Nspaces, Espaces, chis_north, chis_east, chis_south, chis_west ) end -env = CTMRGEnv(corners, edges) -# apply one CTMRG iteration with fixeds -ctm_alg = CTMRG(; fixedspace=true) -env′, = ctmrg_iter(peps, env, ctm_alg) +@testset "Specific U1 spaces" begin + unitcell = (2, 2) -# compute random expecation value to test matching bonds -random_tensor = TensorMap(randn, scalartype(peps), Pspace, Pspace) -random_op = repeat( - LocalOperator(fill(ℂ^2, 1, 1), (CartesianIndex(1, 1),) => random_tensor), unitcell... -) -@test expectation_value(peps, random_op, env) isa Number -@test expectation_value(peps, random_op, env′) isa Number + PA = U1Space(-1 => 1, 0 => 1) + PB = U1Space(0 => 1, 1 => 1) + Vpeps = U1Space(-1 => 2, 0 => 1, 1 => 2) + Venv = U1Space(-2 => 2, -1 => 3, 0 => 4, 1 => 3, 2 => 2) + + Pspaces = [PA PB; PB PA] + Nspaces = [Vpeps Vpeps'; Vpeps' Vpeps] + chis = [Venv Venv; Venv Venv] + + test_unitcell(unitcell, Pspaces, Nspaces, Nspaces, chis, chis, chis, chis) +end diff --git a/test/heisenberg.jl b/test/heisenberg.jl index b72268ad..a4bf9b4a 100644 --- a/test/heisenberg.jl +++ b/test/heisenberg.jl @@ -8,7 +8,7 @@ using OptimKit # initialize parameters χbond = 2 χenv = 16 -ctm_alg = CTMRG(; trscheme=truncdim(χenv), tol=1e-10, miniter=4, maxiter=100, verbosity=1) +ctm_alg = CTMRG(; tol=1e-10, miniter=4, maxiter=100, verbosity=1, trscheme=truncdim(χenv)) opt_alg = PEPSOptimize(; boundary_alg=ctm_alg, optimizer=LBFGS(4; maxiter=100, gradtol=1e-3, verbosity=2), @@ -21,7 +21,7 @@ opt_alg = PEPSOptimize(; Random.seed!(91283219347) H = square_lattice_heisenberg() psi_init = InfinitePEPS(2, χbond) -env_init = leading_boundary(CTMRGEnv(psi_init; Venv=ComplexSpace(χenv)), psi_init, ctm_alg) +env_init = leading_boundary(CTMRGEnv(psi_init, ComplexSpace(χenv)), psi_init, ctm_alg) # find fixedpoint result = fixedpoint(psi_init, H, opt_alg, env_init) @@ -32,7 +32,7 @@ result = fixedpoint(psi_init, H, opt_alg, env_init) H_2x2 = square_lattice_heisenberg(; unitcell=(2, 2)) psi_init_2x2 = InfinitePEPS(2, χbond; unitcell=(2, 2)) env_init_2x2 = leading_boundary( - CTMRGEnv(psi_init_2x2; Venv=ComplexSpace(χenv)), psi_init_2x2, ctm_alg + CTMRGEnv(psi_init_2x2, ComplexSpace(χenv)), psi_init_2x2, ctm_alg ) result_2x2 = fixedpoint(psi_init_2x2, H_2x2, opt_alg, env_init_2x2) diff --git a/test/pwave.jl b/test/pwave.jl index 0c08c93b..d2f7d3aa 100644 --- a/test/pwave.jl +++ b/test/pwave.jl @@ -10,13 +10,11 @@ unitcell = (2, 2) H = square_lattice_pwave(; unitcell) χbond = 2 χenv = 16 -ctm_alg = CTMRG(; - trscheme=truncdim(χenv), tol=1e-10, miniter=4, maxiter=50, fixedspace=true, verbosity=1 -) +ctm_alg = CTMRG(; tol=1e-10, miniter=4, maxiter=100, 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), + gradient_alg=GMRES(; tol=1e-3, maxiter=2, krylovdim=50), reuse_env=true, verbosity=2, ) @@ -26,7 +24,7 @@ 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) -env_init = leading_boundary(CTMRGEnv(psi_init; Venv=Envspace), psi_init, ctm_alg); +env_init = leading_boundary(CTMRGEnv(psi_init, Envspace), psi_init, ctm_alg); # find fixedpoint result = fixedpoint(psi_init, H, opt_alg, env_init) diff --git a/test/runtests.jl b/test/runtests.jl index 799726b8..9266c8e7 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -17,6 +17,12 @@ end @time @safetestset "Gradients" begin include("ctmrg/gradients.jl") end + @time @safetestset "Unit cell" begin + include("ctmrg/unitcell.jl") + end + @time @safetestset "SVD wrapper" begin + include("ctmrg/svd_wrapper.jl") + end end if GROUP == "ALL" || GROUP == "MPS" @time @safetestset "VUMPS" begin diff --git a/test/tf_ising.jl b/test/tf_ising.jl index 881bdd13..c40b17ba 100644 --- a/test/tf_ising.jl +++ b/test/tf_ising.jl @@ -31,7 +31,7 @@ opt_alg = PEPSOptimize(; # initialize states H = square_lattice_tf_ising(; h) psi_init = InfinitePEPS(2, χbond) -env_init = leading_boundary(CTMRGEnv(psi_init; Venv=ComplexSpace(χenv)), psi_init, ctm_alg) +env_init = leading_boundary(CTMRGEnv(psi_init, ComplexSpace(χenv)), psi_init, ctm_alg) # find fixedpoint result = fixedpoint(psi_init, H, opt_alg, env_init)