Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Enhancements for AD-VUMPS support #82

Open
wants to merge 25 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
25 commits
Select commit Hold shift + click to select a range
e44eb72
InfiniteTransferPEPS(PeriodicArray is so weird!
XingyuZhang2018 Oct 10, 2024
fd4b42d
::Type for constructer!
XingyuZhang2018 Oct 11, 2024
a193ea8
backup
XingyuZhang2018 Oct 14, 2024
1c490d8
why does vumps fix gauge not work?
XingyuZhang2018 Oct 14, 2024
bef7ce0
backup
XingyuZhang2018 Oct 18, 2024
85b1983
backup
XingyuZhang2018 Oct 20, 2024
4d32282
start own VUMPS without MPSKit
XingyuZhang2018 Oct 21, 2024
b4be238
to do: leftorth and rightorth
XingyuZhang2018 Oct 21, 2024
ad1ad05
finish canonical form, start FLFR
XingyuZhang2018 Oct 22, 2024
8dafeb2
finish leftenv
XingyuZhang2018 Oct 22, 2024
7c274df
finished vumps forward
XingyuZhang2018 Oct 23, 2024
ab980ae
add test
XingyuZhang2018 Oct 23, 2024
2f25070
add updown
XingyuZhang2018 Oct 24, 2024
0eb4515
to do: fix updown index
XingyuZhang2018 Oct 24, 2024
2732121
finish vumps forward
XingyuZhang2018 Oct 25, 2024
7152ba9
AD for leftenv and rightenv
XingyuZhang2018 Oct 25, 2024
d7e146c
finish ad one unit cell vumps
XingyuZhang2018 Oct 26, 2024
3c2620f
two side vumps ad is not stable
XingyuZhang2018 Oct 27, 2024
d708320
another way to calculate Z; but still incorrect two side AD
XingyuZhang2018 Oct 28, 2024
ff76337
C4 iPEPS opt with ad-vumps
XingyuZhang2018 Oct 30, 2024
1dd02a1
delete some things
XingyuZhang2018 Oct 30, 2024
da2adff
two side vumps
XingyuZhang2018 Oct 30, 2024
3e542c3
correct large unit cell index
XingyuZhang2018 Nov 2, 2024
03c0e82
backup
XingyuZhang2018 Nov 7, 2024
14e97ed
alg_rrule=KrylovKit.GMRES() for degeneracy
XingyuZhang2018 Nov 8, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -7,12 +7,14 @@ version = "0.2.1"
Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If I'm not mistaken, this is not being used right?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I didn't realize I added it. Maybe it can be removed if other packages don't rely on it.

KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
LoggingExtras = "e6f89c97-d47a-5376-807f-9c37f3926c36"
MPSKit = "bb1c41ca-d63c-52ed-829e-0820dda26502"
MPSKitModels = "ca635005-6f8c-4cd1-b51d-8491250ef2ab"
OptimKit = "77e91f04-9b3b-57a6-a776-40b61faaebe0"
Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not so familiar with this package, as far as I know this was mostly used for the @kwdef, which is now a part of Base julia.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have employed the @unpack macro to manage constructors, a practice I find convenient and concise.

Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
Expand Down
55 changes: 9 additions & 46 deletions src/PEPSKit.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@ module PEPSKit

using LinearAlgebra, Statistics, Base.Threads, Base.Iterators, Printf
using Base: @kwdef
using Parameters
using Compat
using Accessors
using VectorInterface
Expand All @@ -11,89 +12,51 @@ using LoggingExtras
using MPSKit: loginit!, logiter!, logfinish!, logcancel!
using MPSKitModels


include("utility/util.jl")
include("utility/svd.jl")
include("utility/rotations.jl")
include("utility/diffset.jl")
include("utility/hook_pullback.jl")
include("utility/autoopt.jl")
include("utility/loginfo.jl")

include("states/abstractpeps.jl")
include("states/infinitepeps.jl")

include("operators/transferpeps.jl")
include("operators/infinitepepo.jl")
include("operators/transferpepo.jl")
include("operators/derivatives.jl")
include("operators/localoperator.jl")
include("operators/lattices/squarelattice.jl")
include("operators/models.jl")

include("environments/ctmrg_environments.jl")
include("environments/transferpeps_environments.jl")
include("environments/transferpepo_environments.jl")
include("environments/vumps_environments.jl")

include("algorithms/contractions/localoperator.jl")
include("algorithms/contractions/ctmrg_contractions.jl")
include("algorithms/contractions/vumps_contractions.jl")

include("algorithms/ctmrg/ctmrg.jl")
include("algorithms/ctmrg/gaugefix.jl")

include("algorithms/vumps/vumps.jl")

include("algorithms/toolbox.jl")

include("algorithms/peps_opt.jl")

include("utility/symmetrization.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.

- `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
using TensorKit, KrylovKit, OptimKit
using PEPSKit: LinSolver, FixedSpaceTruncation, SVDAdjoint
const ctmrg_maxiter = 100
const ctmrg_miniter = 4
const ctmrg_tol = 1e-8
const fpgrad_maxiter = 30
const fpgrad_tol = 1e-6
const ctmrgscheme = :simultaneous
const reuse_env = true
const trscheme = FixedSpaceTruncation()
const iterscheme = :fixed
const fwd_alg = TensorKit.SVD()
const rrule_alg = GMRES(; tol=1e1ctmrg_tol)
const svd_alg = SVDAdjoint(; fwd_alg, rrule_alg)
const optimizer = LBFGS(32; maxiter=100, gradtol=1e-4, verbosity=2)
const gradient_linsolver = KrylovKit.BiCGStab(;
maxiter=Defaults.fpgrad_maxiter, tol=Defaults.fpgrad_tol
)
const gradient_alg = LinSolver(; solver=gradient_linsolver, iterscheme)
end
include("utility/defaults.jl")

export SVDAdjoint, IterSVD, NonTruncSVDAdjoint
export FixedSpaceTruncation, ProjectorAlg, CTMRG, CTMRGEnv, correlation_length
export FixedSpaceTruncation, ProjectorAlg, CTMRG, CTMRGEnv, VUMPS, VUMPSRuntime, VUMPSEnv, correlation_length
export LocalOperator
export expectation_value, costfun, product_peps
export leading_boundary
export PEPSOptimize, GeomSum, ManualIter, LinSolver
export fixedpoint

export InfinitePEPS, InfiniteTransferPEPS
export InfinitePEPO, InfiniteTransferPEPO
export initializeMPS, initializePEPS
export ReflectDepth, ReflectWidth, Rotate, RotateReflect
export symmetrize!, symmetrize_retract_and_finalize!
Expand Down
234 changes: 234 additions & 0 deletions src/algorithms/contractions/vumps_contractions.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,234 @@

"""
ρ = ρmap(ρ::Matrix{<:AbstractTensorMap}, A::Matrix{<:AbstractTensorMap})
````
┌─ Aᵢⱼ─ ┌─
ρᵢⱼ │ = ρⱼ₊₁
└─ Aᵢⱼ─ └─
````
"""
function ρmap(ρ::Matrix{<:AbstractTensorMap}, A::Matrix{<:AbstractTensorMap})
Ni, Nj = size(ρ)
ρ = deepcopy(ρ)
@inbounds for j in 1:Nj, i in 1:Ni
jr = mod1(j + 1, Nj)
@tensor ρ[i,jr][-1; -2] = ρ[i,j][4; 1] * A[i,j][1 2 3; -2] * conj(A[i,j][4 2 3; -1])
end
return ρ
end

"""
C = LRtoC(L::Matrix{<:AbstractTensorMap}, R::Matrix{<:AbstractTensorMap})

```
── Cᵢⱼ ── = ── Lᵢⱼ ── Rᵢⱼ₊₁ ──
```
"""
function LRtoC(L::Matrix{<:AbstractTensorMap}, R::Matrix{<:AbstractTensorMap})
Rijr = circshift(R, (0, -1))
return L .* Rijr
end

"""
AC = ALCtoAC(AL::Matrix{<:AbstractTensorMap}, C::Matrix{<:AbstractTensorMap})

```
── ACᵢⱼ ── = ── ALᵢⱼ ── Cᵢⱼ ──
| |
```
"""
function ALCtoAC(AL::Matrix{<:AbstractTensorMap}, C::Matrix{<:AbstractTensorMap})
return AL .* C
end

"""
FLm = FLmap(FLi::Vector{<:AbstractTensorMap},
ALui::Vector{<:AbstractTensorMap},
ALdir::Vector{<:AbstractTensorMap},
Ati::Vector{<:AbstractTensorMap},
Abi::Vector{<:AbstractTensorMap})

```
┌── ┌── ALuᵢⱼ ──
│ │ │
FLᵢⱼ₊₁ = FLᵢⱼ ─ Oᵢⱼ ──
│ │ │
└── └── ALdᵢᵣⱼ ─
```
"""
function FLmap(FLi::Vector{<:AbstractTensorMap},
ALui::Vector{<:AbstractTensorMap},
ALdir::Vector{<:AbstractTensorMap},
Ati::Vector{<:AbstractTensorMap},
Abi::Vector{<:AbstractTensorMap})
FLm = [@tensoropt FL[-1 -2 -3; -4] := FL[6 5 4; 1] * ALu[1 2 3; -4] * At[9; 2 -2 8 5] *
Ab[3 -3 7 4; 9] * ALd[-1; 6 8 7] for (FL, ALu, ALd, At, Ab) in zip(FLi, ALui, ALdir, Ati, Abi)]

return circshift(FLm, 1)
end

"""
```
┌── ALuᵢⱼ ── ┌──
Lᵢⱼ | = Lᵢⱼ₊₁
└── ALdᵢᵣⱼ ── └──
```
"""
function Lmap(Li::Vector{<:AbstractTensorMap},
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this not the same as ρmap defined above?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Initially, I created ρmap, but later I realized that the index of Ad might vary in certain scenarios. As a result, I will eliminate ρmap and substitute it with Lmap.

ALui::Vector{<:AbstractTensorMap},
ALdir::Vector{<:AbstractTensorMap})
Lm = [@tensoropt L[-6; -4] := ALu[1 2 3; -4] * L[5; 1] * ALd[-6; 5 2 3] for (L, ALu, ALd) in zip(Li, ALui, ALdir)]

return circshift(Lm, 1)
end

"""
FRm = FRmap(FRi::Vector{<:AbstractTensorMap},
ARui::Vector{<:AbstractTensorMap},
ARdir::Vector{<:AbstractTensorMap},
Ati::Vector{<:AbstractTensorMap},
Abi::Vector{<:AbstractTensorMap})

```
── ARuᵢⱼ ──┐ ──┐
│ │ │
── Oᵢⱼ ──FRᵢⱼ = ──FRᵢⱼ₋₁
│ │ │
── ARdᵢᵣⱼ ──┘ ──┘
```
"""
function FRmap(FRi::Vector{<:AbstractTensorMap},
ARui::Vector{<:AbstractTensorMap},
ARdir::Vector{<:AbstractTensorMap},
Ati::Vector{<:AbstractTensorMap},
Abi::Vector{<:AbstractTensorMap})
FRm = [@tensoropt FR[-1 -2 -3; -4] := ARu[-1 1 2; 3] * FR[3 4 5; 8] * At[9; 1 4 7 -2] *
Ab[2 5 6 -3; 9] * ARd[8; -4 7 6] for (FR, ARu, ARd, At, Ab) in zip(FRi, ARui, ARdir, Ati, Abi)]

return circshift(FRm, -1)
end

"""
Rm = Rmap(FRi::Vector{<:AbstractTensorMap},
ARui::Vector{<:AbstractTensorMap},
ARdir::Vector{<:AbstractTensorMap},
)

```
── ARuᵢⱼ ──┐ ──┐
│ Rᵢⱼ = Rᵢⱼ₋₁
── ARdᵢᵣⱼ ──┘ ──┘
```
"""
function Rmap(Ri::Vector{<:AbstractTensorMap},
ARui::Vector{<:AbstractTensorMap},
ARdir::Vector{<:AbstractTensorMap})
Rm = [@tensoropt R[-1; -5] := ARu[-1 2 3; 4] * R[4; 6] * ARd[6; -5 2 3] for (R, ARu, ARd) in zip(Ri, ARui, ARdir)]

return circshift(Rm, -1)
end

"""
ACm = ACmap(ACj::Vector{<:AbstractTensorMap},
FLj::Vector{<:AbstractTensorMap},
FRj::Vector{<:AbstractTensorMap},
Atj::Vector{<:AbstractTensorMap},
Abj::Vector{<:AbstractTensorMap})

```
┌─────── ACᵢⱼ ─────┐
┌───── ACᵢ₊₁ⱼ ─────┐ │ │ │
│ │ │ = FLᵢⱼ ─── Oᵢⱼ ───── FRᵢⱼ
│ │ │

```
"""
function ACmap(ACj::Vector{<:AbstractTensorMap},
FLj::Vector{<:AbstractTensorMap},
FRj::Vector{<:AbstractTensorMap},
Atj::Vector{<:AbstractTensorMap},
Abj::Vector{<:AbstractTensorMap})
ACm = [@tensoropt AC[-1 -2 -3; -4] := AC[1 2 3; 4] * FL[-1 6 5; 1]* At[9; 2 7 -2 6] *
Ab[3 8 -3 5; 9] * FR[4 7 8; -4] for (AC, FL, FR, At, Ab) in zip(ACj, FLj, FRj, Atj, Abj)]

return circshift(ACm, 1)
end

"""
Cmap(Cij, FLjp, FRj, II)

```
┌────Cᵢⱼ ───┐
┌── Cᵢ₊₁ⱼ ──┐ │ │
│ │ = FLᵢⱼ₊₁ ──── FRᵢⱼ
│ │

```
"""
function Cmap(Cj::Vector{<:AbstractTensorMap},
FLjr::Vector{<:AbstractTensorMap},
FRj::Vector{<:AbstractTensorMap})
Cm = [@tensoropt C[-1; -2] := C[1; 2] * FL[-1 3 4; 1] * FR[2 3 4; -2] for (C, FL, FR) in zip(Cj, FLjr, FRj)]

return circshift(Cm, 1)
end

"""
nearest_neighbour_energy(ipeps::InfinitePEPS, Hh, Hv, env::VUMPSEnv)

Compute the energy of the nearest neighbour Hamiltonian for an infinite PEPS.

```
┌────── ACuᵢⱼ ────ARuᵢᵣ ──────┐
│ │ │ │
FLoᵢⱼ ── Oᵢⱼ ────── Oᵢᵣ ──── FRoᵢᵣ ir = Ni + 1 - i
│ │ │ │ jr = j + 1
└───── ACdᵢᵣⱼ ─────ARdᵢᵣᵣ─────┘

┌─────── ACuᵢⱼ ─────┐
│ │ │
FLuᵢⱼ ─── Oᵢⱼ ───── FRuᵢⱼ ir = i + 1
│ │ │ irr = Ni - i
FLoᵢᵣⱼ ── Oᵢᵣⱼ ─── FRoᵢᵣⱼ
│ │ │
└────── ACdᵢᵣᵣⱼ ───┘
```
"""
function nearest_neighbour_energy(ipeps::InfinitePEPS, Hh, Hv, env::VUMPSEnv)
@unpack ACu, ARu, ACd, ARd, FLu, FRu, FLo, FRo = env
Ni, Nj = size(ipeps)

energy_tol = 0
for j in 1:Nj, i in 1:Ni
# horizontal contraction
ir = Ni + 1 - i
jr = mod1(j + 1, Nj)
@tensoropt oph[-1 -2; -3 -4] := FLo[i,j][18 12 15; 5] * ACu[i,j][5 6 7; 8] * ipeps.A[i,j][-1; 6 13 19 12] *
conj(ipeps.A[i,j][-3; 7 16 20 15]) * conj(ACd[ir,j][18 19 20; 21]) *
ARu[i,jr][8 9 10; 11] * ipeps.A[i,jr][-2; 9 14 22 13] *
conj(ipeps.A[i,jr][-4; 10 17 23 16]) * conj(ARd[ir,jr][21 22 23; 24]) * FRo[i,jr][11 14 17; 24]

@tensor eh = oph[1 2; 3 4] * Hh[3 4; 1 2]
@tensor nh = oph[1 2; 1 2]
energy_tol += eh / nh
# @show eh / nh eh nh

# vertical contraction
ir = mod1(i + 1, Ni)
irr = mod1(Ni - i, Ni)
@tensoropt opv[-1 -2; -3 -4] := FLu[i,j][21 19 20; 18] * ACu[i,j][18 12 15 5] * ipeps.A[i,j][-1; 12 6 13 19] *
conj(ipeps.A[i,j][-3; 15 7 16 20]) * FRu[i,j][5 6 7; 8] * FLo[ir,j][24 22 23; 21] *
ipeps.A[ir,j][-2; 13 9 14 22] * conj(ipeps.A[ir,j][-4; 16 10 17 23]) *
FRo[ir,j][8 9 10; 11] * conj(ACd[irr,j][24 14 17; 11])

@tensor ev = opv[1 2; 3 4] * Hv[3 4; 1 2]
@tensor nv = opv[1 2; 1 2]
energy_tol += ev / nv
# @show ev / nv ev nv

# penalty term
# energy_tol += 0.1 * abs(eh / nh - eh / nh)
end

return energy_tol/Ni/Nj
end
Loading