diff --git a/README.md b/README.md index b4b94c64c..377b5e7b5 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,4 @@ - - - PEPSKit.jl logo - + # PEPSKit.jl @@ -20,7 +17,7 @@ **Tools for working with projected entangled-pair states** -It contracts, it optimizes, it may break. +*It contracts, it optimizes, it evolves.* ## Installation diff --git a/docs/make.jl b/docs/make.jl index a437d7edc..c432c28b8 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -56,7 +56,7 @@ makedocs(; modules=[PEPSKit, MPSKitModels], sitename="PEPSKit.jl", format=Documenter.HTML(; - prettyurls=get(ENV, "CI", nothing) == "true", mathengine, size_threshold=1024000 + prettyurls=true, mathengine, assets=["assets/custom.css"], size_threshold=1024000 ), pages=[ "Home" => "index.md", diff --git a/docs/src/assets/custom.css b/docs/src/assets/custom.css new file mode 100644 index 000000000..ef9f18328 --- /dev/null +++ b/docs/src/assets/custom.css @@ -0,0 +1,42 @@ +/* set fixed non-trivial inversion and hue rotation */ +:root { + --inversionFraction: 100%; + --hueRotation: -180deg; +} + +/* color-invert images */ +.color-invertible { + filter: invert(var(--inversionFraction)) hue-rotate(var(--hueRotation)) !important; + -ms-filter: invert(var(--inversionFraction)) !important; + -webkit-filter: invert(var(--inversionFraction)) hue-rotate(var(--hueRotation)) !important; +} +/* including the logo */ +.docs-logo { + filter: invert(var(--inversionFraction)) hue-rotate(var(--hueRotation)) !important; + -ms-filter: invert(var(--inversionFraction)) !important; + -webkit-filter: invert(var(--inversionFraction)) hue-rotate(var(--hueRotation)) !important; +} + + +/* but disable for the two light themes */ +.theme--documenter-light .color-invertible { + filter: invert(0%) hue-rotate(0deg) !important; + -ms-filter: invert(var(0%)) !important; + -webkit-filter: invert(var(0%)) hue-rotate(0deg) !important; +} +.theme--catppuccin-latte .color-invertible { + filter: invert(0%) hue-rotate(0deg) !important; + -ms-filter: invert(var(0%)) !important; + -webkit-filter: invert(var(0%)) hue-rotate(0deg) !important; +} +/* for the logo as well */ +.theme--documenter-light .docs-logo { + filter: invert(0%) hue-rotate(0deg) !important; + -ms-filter: invert(var(0%)) !important; + -webkit-filter: invert(var(0%)) hue-rotate(0deg) !important; +} +.theme--catppuccin-latte .docs-logo { + filter: invert(0%) hue-rotate(0deg) !important; + -ms-filter: invert(var(0%)) !important; + -webkit-filter: invert(var(0%)) hue-rotate(0deg) !important; +} diff --git a/docs/src/assets/figures/pepo.svg b/docs/src/assets/figures/pepo.svg new file mode 100644 index 000000000..bfcd30c65 --- /dev/null +++ b/docs/src/assets/figures/pepo.svg @@ -0,0 +1,1637 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/src/assets/figures/pepo_fixedpoint_equation.svg b/docs/src/assets/figures/pepo_fixedpoint_equation.svg new file mode 100644 index 000000000..352a29713 --- /dev/null +++ b/docs/src/assets/figures/pepo_fixedpoint_equation.svg @@ -0,0 +1,4055 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/src/assets/figures/peps_norm_network.svg b/docs/src/assets/figures/peps_norm_network.svg new file mode 100644 index 000000000..2f1ef66f0 --- /dev/null +++ b/docs/src/assets/figures/peps_norm_network.svg @@ -0,0 +1,2494 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/src/assets/figures/peps_transfer_fixedpoint_equation.svg b/docs/src/assets/figures/peps_transfer_fixedpoint_equation.svg new file mode 100644 index 000000000..c4c0c1e95 --- /dev/null +++ b/docs/src/assets/figures/peps_transfer_fixedpoint_equation.svg @@ -0,0 +1,1925 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/src/assets/figures/peps_transfer_operator.svg b/docs/src/assets/figures/peps_transfer_operator.svg new file mode 100644 index 000000000..191ca5237 --- /dev/null +++ b/docs/src/assets/figures/peps_transfer_operator.svg @@ -0,0 +1,1126 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/src/assets/logo-dark.svg b/docs/src/assets/logo-dark.svg deleted file mode 100644 index 9d8f3844f..000000000 --- a/docs/src/assets/logo-dark.svg +++ /dev/null @@ -1,111 +0,0 @@ - - - - - - image/svg+xml - - - - - - - - - - - - - - - - diff --git a/docs/src/assets/logo_readme.svg b/docs/src/assets/logo_readme.svg new file mode 100644 index 000000000..6d6025309 --- /dev/null +++ b/docs/src/assets/logo_readme.svg @@ -0,0 +1,18 @@ + + +
+
+ + PEPSKit Logo +
+
+ + \ No newline at end of file diff --git a/docs/src/examples/3d_ising_partition_function/index.md b/docs/src/examples/3d_ising_partition_function/index.md index ee45eb022..9fab4c32c 100644 --- a/docs/src/examples/3d_ising_partition_function/index.md +++ b/docs/src/examples/3d_ising_partition_function/index.md @@ -106,26 +106,40 @@ true ## Contracting the partition function To contract our infinite 3D partition function, we first reinterpret it as an infinite power -of a slice-to-slice transfer operator ``T``, where ``T`` can be seen as an infinite 2D +of a slice-to-slice transfer operator ``\mathbb{T}``, where ``\mathbb{T}`` can be seen as an infinite 2D projected entangled-pair operator (PEPO) which consists of the rank-6 tensor `O` at each -site of an infinite 2D square lattice. In the same spirit as the boundary MPS approach, all -we need to contract the whole partition function is to find the leading eigenvector of this -PEPO. The fixed point of such a PEPO can be parametrized as a PEPS, and for the case of a -Hermitian transfer operator we can find this PEPS through [variational optimization](@cite -vanderstraeten_residual_2018). +site of an infinite 2D square lattice, -Indeed, for a Hermitian transfer operator ``T`` we can characterize the fixed point PEPS +```@raw html +
+pepo +
+``` + +To contract the original infinite network, all we need to do is to find the leading +eigenvector of the PEPO ``\mathbb{T}``, The fixed point of such a PEPO can be parametrized +as a PEPS ``\psi``, which should then satisfy the eigenvalue equation +``\mathbb{T} |\psi\rangle = \Lambda |\psi\rangle``, or diagrammatically: + +```@raw html +
+pepo fixedpoint equation +
+``` + +For a Hermitian transfer operator ``\mathbb{T}`` we can characterize the fixed point PEPS ``|\psi\rangle`` which satisfies the eigenvalue equation -``T |\psi\rangle = \Lambda |\psi\rangle`` corresponding to the largest magnitude eigenvalue -``\Lambda`` as the solution of a variational problem +``\mathbb{T} |\psi\rangle = \Lambda |\psi\rangle`` corresponding to the largest magnitude +eigenvalue ``\Lambda`` as the solution of a variational problem ```math -|\psi\rangle = \text{argmin}_{|\psi\rangle} \left ( \lim_{N \to ∞} - \frac{1}{N} \log \left( \frac{\langle \psi | T | \psi \rangle}{\langle \psi | \psi \rangle} \right) \right ) , +|\psi\rangle = \text{argmin}_{|\psi\rangle} \left ( \lim_{N \to ∞} - \frac{1}{N} \log \left( \frac{\langle \psi | \mathbb{T} | \psi \rangle}{\langle \psi | \psi \rangle} \right) \right ) , ``` -where ``N`` is the diverging number of sites of the 2D transfer operator ``T``. The function -minimized in this expression is exactly the free energy per site of the partition function, -so we essentially find the fixed-point PEPS by variationally minimizing the free energy. +where ``N`` is the diverging number of sites of the 2D transfer operator ``\mathbb{T}``. The function +minimized in this expression is exactly the free energy per site of the partition function. +This means we can directly find the fixed-point PEPS by +[variationally minimizing the free energy](@cite vanderstraeten_residual_2018). ### Defining the cost function @@ -134,7 +148,7 @@ use [OptimKit.jl](https://github.com/Jutho/OptimKit.jl) to actually optimize it. immediately recognize the denominator ``\langle \psi | \psi \rangle`` as the familiar PEPS norm, where we can compute the norm per site as the [`network_value`](@ref) of the corresponding [`InfiniteSquareNetwork`](@ref) by contracting it with the CTMRG algorithm. -Similarly, the numerator ``\langle \psi | T | \psi \rangle`` is nothing more than an +Similarly, the numerator ``\langle \psi | \mathbb{T} | \psi \rangle`` is nothing more than an `InfiniteSquareNetwork` consisting of three layers corresponding to the ket, transfer operator and bra objects. This object can also be constructed and contracted in a straightforward way, so we can again compute its `network_value`. @@ -279,83 +293,85 @@ optimizer_alg = LBFGS(32; maxiter=100, gradtol=1e-5, verbosity=3) ```` [ Info: LBFGS: initializing with f = -0.554073395182, ‖∇f‖ = 7.7844e-01 -┌ Warning: CTMRG cancel 150: obj = +1.702942228759e+01 +1.443123032606e-07im err = 2.4386740957e-05 time = 1.12 sec -└ @ PEPSKit ~/repos/PEPSKit.jl/src/algorithms/ctmrg/ctmrg.jl:155 -[ Info: LBFGS: iter 1, time 117.39 s: f = -0.777080930369, ‖∇f‖ = 3.1305e-02, α = 7.10e+02, m = 0, nfg = 7 -[ Info: LBFGS: iter 2, time 118.47 s: f = -0.784111515961, ‖∇f‖ = 2.0103e-02, α = 1.00e+00, m = 1, nfg = 1 -[ Info: LBFGS: iter 3, time 118.64 s: f = -0.792705733484, ‖∇f‖ = 2.3327e-02, α = 1.00e+00, m = 2, nfg = 1 -[ Info: LBFGS: iter 4, time 118.78 s: f = -0.796289732476, ‖∇f‖ = 2.2475e-02, α = 1.00e+00, m = 3, nfg = 1 -[ Info: LBFGS: iter 5, time 118.89 s: f = -0.799674902374, ‖∇f‖ = 7.0288e-03, α = 1.00e+00, m = 4, nfg = 1 -[ Info: LBFGS: iter 6, time 118.98 s: f = -0.800082100121, ‖∇f‖ = 1.2717e-03, α = 1.00e+00, m = 5, nfg = 1 -[ Info: LBFGS: iter 7, time 119.08 s: f = -0.800110603125, ‖∇f‖ = 1.3384e-03, α = 1.00e+00, m = 6, nfg = 1 -[ Info: LBFGS: iter 8, time 119.18 s: f = -0.800262201996, ‖∇f‖ = 2.4945e-03, α = 1.00e+00, m = 7, nfg = 1 -[ Info: LBFGS: iter 9, time 119.27 s: f = -0.800450505448, ‖∇f‖ = 2.9259e-03, α = 1.00e+00, m = 8, nfg = 1 -[ Info: LBFGS: iter 10, time 119.36 s: f = -0.800764917087, ‖∇f‖ = 1.7221e-03, α = 1.00e+00, m = 9, nfg = 1 -[ Info: LBFGS: iter 11, time 119.45 s: f = -0.800876048838, ‖∇f‖ = 2.2475e-03, α = 1.00e+00, m = 10, nfg = 1 -[ Info: LBFGS: iter 12, time 119.53 s: f = -0.801100867467, ‖∇f‖ = 1.5561e-03, α = 1.00e+00, m = 11, nfg = 1 -[ Info: LBFGS: iter 13, time 119.63 s: f = -0.801317048856, ‖∇f‖ = 1.1561e-03, α = 1.00e+00, m = 12, nfg = 1 -[ Info: LBFGS: iter 14, time 119.73 s: f = -0.801373050545, ‖∇f‖ = 7.1300e-04, α = 1.00e+00, m = 13, nfg = 1 -[ Info: LBFGS: iter 15, time 119.82 s: f = -0.801388615264, ‖∇f‖ = 2.8462e-04, α = 1.00e+00, m = 14, nfg = 1 -[ Info: LBFGS: iter 16, time 119.91 s: f = -0.801394633333, ‖∇f‖ = 2.7607e-04, α = 1.00e+00, m = 15, nfg = 1 -[ Info: LBFGS: iter 17, time 119.99 s: f = -0.801408061564, ‖∇f‖ = 3.6096e-04, α = 1.00e+00, m = 16, nfg = 1 -[ Info: LBFGS: iter 18, time 120.09 s: f = -0.801509542169, ‖∇f‖ = 1.9822e-03, α = 1.00e+00, m = 17, nfg = 1 -[ Info: LBFGS: iter 19, time 120.20 s: f = -0.801578405251, ‖∇f‖ = 1.8040e-03, α = 1.00e+00, m = 18, nfg = 1 -[ Info: LBFGS: iter 20, time 120.59 s: f = -0.801694524424, ‖∇f‖ = 2.9356e-03, α = 5.48e-01, m = 19, nfg = 3 -[ Info: LBFGS: iter 21, time 121.07 s: f = -0.801761920683, ‖∇f‖ = 1.1993e-03, α = 3.82e-01, m = 20, nfg = 2 -[ Info: LBFGS: iter 22, time 121.19 s: f = -0.801797785494, ‖∇f‖ = 6.0337e-04, α = 1.00e+00, m = 21, nfg = 1 -[ Info: LBFGS: iter 23, time 121.45 s: f = -0.801808747834, ‖∇f‖ = 3.7053e-04, α = 5.24e-01, m = 22, nfg = 2 -[ Info: LBFGS: iter 24, time 121.58 s: f = -0.801812729173, ‖∇f‖ = 3.0781e-04, α = 1.00e+00, m = 23, nfg = 1 -[ Info: LBFGS: iter 25, time 121.71 s: f = -0.801816445211, ‖∇f‖ = 2.9994e-04, α = 1.00e+00, m = 24, nfg = 1 -[ Info: LBFGS: iter 26, time 121.84 s: f = -0.801824713130, ‖∇f‖ = 3.6496e-04, α = 1.00e+00, m = 25, nfg = 1 -[ Info: LBFGS: iter 27, time 121.98 s: f = -0.801839673823, ‖∇f‖ = 5.4222e-04, α = 1.00e+00, m = 26, nfg = 1 -[ Info: LBFGS: iter 28, time 122.12 s: f = -0.801857478904, ‖∇f‖ = 2.7917e-04, α = 1.00e+00, m = 27, nfg = 1 -[ Info: LBFGS: iter 29, time 122.29 s: f = -0.801864555224, ‖∇f‖ = 1.2319e-04, α = 1.00e+00, m = 28, nfg = 1 -[ Info: LBFGS: iter 30, time 122.48 s: f = -0.801865598736, ‖∇f‖ = 8.6048e-05, α = 1.00e+00, m = 29, nfg = 1 -[ Info: LBFGS: iter 31, time 122.63 s: f = -0.801867571755, ‖∇f‖ = 8.8636e-05, α = 1.00e+00, m = 30, nfg = 1 -[ Info: LBFGS: iter 32, time 122.78 s: f = -0.801870393528, ‖∇f‖ = 2.6554e-04, α = 1.00e+00, m = 31, nfg = 1 -[ Info: LBFGS: iter 33, time 122.93 s: f = -0.801874797039, ‖∇f‖ = 2.7841e-04, α = 1.00e+00, m = 32, nfg = 1 -[ Info: LBFGS: iter 34, time 123.13 s: f = -0.801877566644, ‖∇f‖ = 1.8523e-04, α = 1.00e+00, m = 32, nfg = 1 -[ Info: LBFGS: iter 35, time 123.33 s: f = -0.801878506245, ‖∇f‖ = 2.0638e-04, α = 1.00e+00, m = 32, nfg = 1 -[ Info: LBFGS: iter 36, time 123.51 s: f = -0.801878995097, ‖∇f‖ = 5.6081e-05, α = 1.00e+00, m = 32, nfg = 1 -[ Info: LBFGS: iter 37, time 123.67 s: f = -0.801879153573, ‖∇f‖ = 6.2356e-05, α = 1.00e+00, m = 32, nfg = 1 -[ Info: LBFGS: iter 38, time 123.82 s: f = -0.801879355075, ‖∇f‖ = 6.0528e-05, α = 1.00e+00, m = 32, nfg = 1 -[ Info: LBFGS: iter 39, time 123.97 s: f = -0.801880115100, ‖∇f‖ = 6.2768e-05, α = 1.00e+00, m = 32, nfg = 1 -[ Info: LBFGS: iter 40, time 124.12 s: f = -0.801881475065, ‖∇f‖ = 6.2301e-05, α = 1.00e+00, m = 32, nfg = 1 -[ Info: LBFGS: iter 41, time 124.29 s: f = -0.801882272425, ‖∇f‖ = 9.5267e-05, α = 1.00e+00, m = 32, nfg = 1 -[ Info: LBFGS: iter 42, time 124.44 s: f = -0.801882600033, ‖∇f‖ = 5.1283e-05, α = 1.00e+00, m = 32, nfg = 1 -[ Info: LBFGS: iter 43, time 124.59 s: f = -0.801882711875, ‖∇f‖ = 2.6091e-05, α = 1.00e+00, m = 32, nfg = 1 -[ Info: LBFGS: iter 44, time 124.73 s: f = -0.801882805828, ‖∇f‖ = 2.9316e-05, α = 1.00e+00, m = 32, nfg = 1 -[ Info: LBFGS: iter 45, time 124.85 s: f = -0.801883027060, ‖∇f‖ = 2.7982e-05, α = 1.00e+00, m = 32, nfg = 1 -[ Info: LBFGS: iter 46, time 124.98 s: f = -0.801883402178, ‖∇f‖ = 3.8102e-05, α = 1.00e+00, m = 32, nfg = 1 -[ Info: LBFGS: iter 47, time 125.13 s: f = -0.801883718321, ‖∇f‖ = 5.3658e-05, α = 1.00e+00, m = 32, nfg = 1 -[ Info: LBFGS: iter 48, time 125.26 s: f = -0.801883962887, ‖∇f‖ = 2.8728e-05, α = 1.00e+00, m = 32, nfg = 1 -[ Info: LBFGS: iter 49, time 125.40 s: f = -0.801884158085, ‖∇f‖ = 3.0680e-05, α = 1.00e+00, m = 32, nfg = 1 -[ Info: LBFGS: iter 50, time 125.55 s: f = -0.801884385940, ‖∇f‖ = 4.1973e-05, α = 1.00e+00, m = 32, nfg = 1 -[ Info: LBFGS: iter 51, time 125.70 s: f = -0.801884810459, ‖∇f‖ = 6.8881e-05, α = 1.00e+00, m = 32, nfg = 1 -[ Info: LBFGS: iter 52, time 125.85 s: f = -0.801885011014, ‖∇f‖ = 3.8651e-05, α = 1.00e+00, m = 32, nfg = 1 -[ Info: LBFGS: iter 53, time 126.00 s: f = -0.801885126625, ‖∇f‖ = 1.9013e-05, α = 1.00e+00, m = 32, nfg = 1 -[ Info: LBFGS: iter 54, time 126.14 s: f = -0.801885186489, ‖∇f‖ = 3.2919e-05, α = 1.00e+00, m = 32, nfg = 1 -[ Info: LBFGS: iter 55, time 126.29 s: f = -0.801885309713, ‖∇f‖ = 4.8521e-05, α = 1.00e+00, m = 32, nfg = 1 -[ Info: LBFGS: iter 56, time 126.45 s: f = -0.801885491631, ‖∇f‖ = 1.1478e-04, α = 1.00e+00, m = 32, nfg = 1 -[ Info: LBFGS: iter 57, time 126.61 s: f = -0.801885912857, ‖∇f‖ = 7.7221e-05, α = 1.00e+00, m = 32, nfg = 1 -[ Info: LBFGS: iter 58, time 126.77 s: f = -0.801886451980, ‖∇f‖ = 6.5316e-05, α = 1.00e+00, m = 32, nfg = 1 -[ Info: LBFGS: iter 59, time 126.91 s: f = -0.801886639804, ‖∇f‖ = 5.1567e-05, α = 1.00e+00, m = 32, nfg = 1 -[ Info: LBFGS: iter 60, time 127.25 s: f = -0.801886699372, ‖∇f‖ = 4.5540e-05, α = 3.68e-01, m = 32, nfg = 2 -[ Info: LBFGS: iter 61, time 127.43 s: f = -0.801886723992, ‖∇f‖ = 2.1992e-05, α = 1.00e+00, m = 32, nfg = 1 -[ Info: LBFGS: iter 62, time 127.61 s: f = -0.801886735202, ‖∇f‖ = 1.8064e-05, α = 1.00e+00, m = 32, nfg = 1 -[ Info: LBFGS: iter 63, time 127.80 s: f = -0.801886771395, ‖∇f‖ = 3.8651e-05, α = 1.00e+00, m = 32, nfg = 1 -[ Info: LBFGS: iter 64, time 127.97 s: f = -0.801886801952, ‖∇f‖ = 4.2630e-05, α = 1.00e+00, m = 32, nfg = 1 -[ Info: LBFGS: iter 65, time 128.14 s: f = -0.801886837856, ‖∇f‖ = 3.9318e-05, α = 1.00e+00, m = 32, nfg = 1 -[ Info: LBFGS: iter 66, time 128.30 s: f = -0.801886916784, ‖∇f‖ = 3.8747e-05, α = 1.00e+00, m = 32, nfg = 1 -[ Info: LBFGS: iter 67, time 128.45 s: f = -0.801887030055, ‖∇f‖ = 3.7139e-05, α = 1.00e+00, m = 32, nfg = 1 -[ Info: LBFGS: iter 68, time 128.62 s: f = -0.801887141198, ‖∇f‖ = 5.7017e-05, α = 1.00e+00, m = 32, nfg = 1 -[ Info: LBFGS: iter 69, time 128.95 s: f = -0.801887199205, ‖∇f‖ = 3.0700e-05, α = 5.24e-01, m = 32, nfg = 2 -[ Info: LBFGS: iter 70, time 129.10 s: f = -0.801887246613, ‖∇f‖ = 1.3885e-05, α = 1.00e+00, m = 32, nfg = 1 -[ Info: LBFGS: iter 71, time 129.25 s: f = -0.801887263716, ‖∇f‖ = 1.5769e-05, α = 1.00e+00, m = 32, nfg = 1 -[ Info: LBFGS: iter 72, time 129.40 s: f = -0.801887319464, ‖∇f‖ = 2.1424e-05, α = 1.00e+00, m = 32, nfg = 1 -[ Info: LBFGS: iter 73, time 129.56 s: f = -0.801887406143, ‖∇f‖ = 1.9896e-05, α = 1.00e+00, m = 32, nfg = 1 -[ Info: LBFGS: iter 74, time 129.88 s: f = -0.801887467460, ‖∇f‖ = 1.9800e-05, α = 3.61e-01, m = 32, nfg = 2 -[ Info: LBFGS: converged after 75 iterations and time 130.03 s: f = -0.801887535670, ‖∇f‖ = 9.9339e-06 +┌ Warning: CTMRG cancel 150: obj = +1.702942227203e+01 +1.438609518622e-07im err = 2.4390784946e-05 time = 25.66 sec +└ @ PEPSKit ~/git/PEPSKit.jl/src/algorithms/ctmrg/ctmrg.jl:153 +[ Info: LBFGS: iter 1, time 72.51 s: f = -0.777080930424, ‖∇f‖ = 3.1305e-02, α = 7.10e+02, m = 0, nfg = 7 +[ Info: LBFGS: iter 2, time 76.22 s: f = -0.784111515920, ‖∇f‖ = 2.0103e-02, α = 1.00e+00, m = 1, nfg = 1 +[ Info: LBFGS: iter 3, time 79.91 s: f = -0.792705733264, ‖∇f‖ = 2.3327e-02, α = 1.00e+00, m = 2, nfg = 1 +[ Info: LBFGS: iter 4, time 82.31 s: f = -0.796289732161, ‖∇f‖ = 2.2475e-02, α = 1.00e+00, m = 3, nfg = 1 +[ Info: LBFGS: iter 5, time 84.94 s: f = -0.799674902231, ‖∇f‖ = 7.0288e-03, α = 1.00e+00, m = 4, nfg = 1 +[ Info: LBFGS: iter 6, time 87.13 s: f = -0.800082100115, ‖∇f‖ = 1.2717e-03, α = 1.00e+00, m = 5, nfg = 1 +[ Info: LBFGS: iter 7, time 89.94 s: f = -0.800110603120, ‖∇f‖ = 1.3384e-03, α = 1.00e+00, m = 6, nfg = 1 +[ Info: LBFGS: iter 8, time 92.01 s: f = -0.800262202003, ‖∇f‖ = 2.4945e-03, α = 1.00e+00, m = 7, nfg = 1 +[ Info: LBFGS: iter 9, time 93.65 s: f = -0.800450505439, ‖∇f‖ = 2.9259e-03, α = 1.00e+00, m = 8, nfg = 1 +[ Info: LBFGS: iter 10, time 95.79 s: f = -0.800764917063, ‖∇f‖ = 1.7221e-03, α = 1.00e+00, m = 9, nfg = 1 +[ Info: LBFGS: iter 11, time 98.09 s: f = -0.800876048822, ‖∇f‖ = 2.2475e-03, α = 1.00e+00, m = 10, nfg = 1 +[ Info: LBFGS: iter 12, time 99.35 s: f = -0.801100867388, ‖∇f‖ = 1.5561e-03, α = 1.00e+00, m = 11, nfg = 1 +[ Info: LBFGS: iter 13, time 100.88 s: f = -0.801317048785, ‖∇f‖ = 1.1561e-03, α = 1.00e+00, m = 12, nfg = 1 +[ Info: LBFGS: iter 14, time 102.48 s: f = -0.801373050522, ‖∇f‖ = 7.1300e-04, α = 1.00e+00, m = 13, nfg = 1 +[ Info: LBFGS: iter 15, time 103.88 s: f = -0.801388615258, ‖∇f‖ = 2.8462e-04, α = 1.00e+00, m = 14, nfg = 1 +[ Info: LBFGS: iter 16, time 105.37 s: f = -0.801394633326, ‖∇f‖ = 2.7607e-04, α = 1.00e+00, m = 15, nfg = 1 +[ Info: LBFGS: iter 17, time 107.89 s: f = -0.801408061547, ‖∇f‖ = 3.6096e-04, α = 1.00e+00, m = 16, nfg = 1 +[ Info: LBFGS: iter 18, time 109.37 s: f = -0.801509542237, ‖∇f‖ = 1.9822e-03, α = 1.00e+00, m = 17, nfg = 1 +[ Info: LBFGS: iter 19, time 111.39 s: f = -0.801578405298, ‖∇f‖ = 1.8040e-03, α = 1.00e+00, m = 18, nfg = 1 +[ Info: LBFGS: iter 20, time 117.86 s: f = -0.801694524832, ‖∇f‖ = 2.9356e-03, α = 5.48e-01, m = 19, nfg = 3 +[ Info: LBFGS: iter 21, time 123.02 s: f = -0.801761920585, ‖∇f‖ = 1.1993e-03, α = 3.82e-01, m = 20, nfg = 2 +[ Info: LBFGS: iter 22, time 124.85 s: f = -0.801797785496, ‖∇f‖ = 6.0337e-04, α = 1.00e+00, m = 21, nfg = 1 +[ Info: LBFGS: iter 23, time 129.85 s: f = -0.801808747834, ‖∇f‖ = 3.7053e-04, α = 5.24e-01, m = 22, nfg = 2 +[ Info: LBFGS: iter 24, time 131.58 s: f = -0.801812729196, ‖∇f‖ = 3.0781e-04, α = 1.00e+00, m = 23, nfg = 1 +[ Info: LBFGS: iter 25, time 133.29 s: f = -0.801816445181, ‖∇f‖ = 2.9994e-04, α = 1.00e+00, m = 24, nfg = 1 +[ Info: LBFGS: iter 26, time 136.36 s: f = -0.801824712580, ‖∇f‖ = 3.6497e-04, α = 1.00e+00, m = 25, nfg = 1 +[ Info: LBFGS: iter 27, time 139.94 s: f = -0.801839673918, ‖∇f‖ = 5.4217e-04, α = 1.00e+00, m = 26, nfg = 1 +[ Info: LBFGS: iter 28, time 143.71 s: f = -0.801857480042, ‖∇f‖ = 2.7918e-04, α = 1.00e+00, m = 27, nfg = 1 +[ Info: LBFGS: iter 29, time 149.41 s: f = -0.801864556175, ‖∇f‖ = 1.2323e-04, α = 1.00e+00, m = 28, nfg = 1 +[ Info: LBFGS: iter 30, time 154.13 s: f = -0.801865599394, ‖∇f‖ = 8.6049e-05, α = 1.00e+00, m = 29, nfg = 1 +[ Info: LBFGS: iter 31, time 157.32 s: f = -0.801867571015, ‖∇f‖ = 8.8673e-05, α = 1.00e+00, m = 30, nfg = 1 +[ Info: LBFGS: iter 32, time 163.14 s: f = -0.801870391456, ‖∇f‖ = 2.6510e-04, α = 1.00e+00, m = 31, nfg = 1 +[ Info: LBFGS: iter 33, time 165.78 s: f = -0.801874799763, ‖∇f‖ = 2.7836e-04, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 34, time 167.98 s: f = -0.801877566434, ‖∇f‖ = 1.8512e-04, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 35, time 170.53 s: f = -0.801878506140, ‖∇f‖ = 2.0603e-04, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 36, time 172.48 s: f = -0.801878994665, ‖∇f‖ = 5.6086e-05, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 37, time 174.33 s: f = -0.801879153394, ‖∇f‖ = 6.2420e-05, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 38, time 176.63 s: f = -0.801879354902, ‖∇f‖ = 6.0532e-05, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 39, time 179.06 s: f = -0.801880114487, ‖∇f‖ = 6.2680e-05, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 40, time 181.77 s: f = -0.801881471358, ‖∇f‖ = 6.2392e-05, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 41, time 184.73 s: f = -0.801882270540, ‖∇f‖ = 9.4942e-05, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 42, time 187.02 s: f = -0.801882598674, ‖∇f‖ = 5.1910e-05, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 43, time 189.45 s: f = -0.801882711254, ‖∇f‖ = 2.6275e-05, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 44, time 191.07 s: f = -0.801882805203, ‖∇f‖ = 2.9425e-05, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 45, time 192.75 s: f = -0.801883026107, ‖∇f‖ = 2.7910e-05, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 46, time 194.86 s: f = -0.801883400172, ‖∇f‖ = 3.7426e-05, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 47, time 197.16 s: f = -0.801883717581, ‖∇f‖ = 5.4717e-05, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 48, time 199.08 s: f = -0.801883966703, ‖∇f‖ = 2.9045e-05, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 49, time 201.86 s: f = -0.801884163647, ‖∇f‖ = 3.0661e-05, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 50, time 205.72 s: f = -0.801884391105, ‖∇f‖ = 4.1905e-05, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 51, time 209.31 s: f = -0.801884815983, ‖∇f‖ = 6.9018e-05, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 52, time 212.00 s: f = -0.801885013426, ‖∇f‖ = 3.8025e-05, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 53, time 215.29 s: f = -0.801885126302, ‖∇f‖ = 1.9306e-05, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 54, time 219.02 s: f = -0.801885184513, ‖∇f‖ = 3.3083e-05, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 55, time 221.19 s: f = -0.801885308658, ‖∇f‖ = 4.9014e-05, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 56, time 223.60 s: f = -0.801885502272, ‖∇f‖ = 1.1303e-04, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 57, time 226.58 s: f = -0.801885922461, ‖∇f‖ = 7.5880e-05, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 58, time 230.06 s: f = -0.801886457901, ‖∇f‖ = 7.2957e-05, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 59, time 232.86 s: f = -0.801886614664, ‖∇f‖ = 6.8816e-05, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 60, time 239.07 s: f = -0.801886696733, ‖∇f‖ = 3.0687e-05, α = 4.26e-01, m = 32, nfg = 2 +[ Info: LBFGS: iter 61, time 243.89 s: f = -0.801886716271, ‖∇f‖ = 2.1581e-05, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 62, time 247.56 s: f = -0.801886732686, ‖∇f‖ = 1.7659e-05, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 63, time 251.23 s: f = -0.801886790357, ‖∇f‖ = 4.1045e-05, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 64, time 255.93 s: f = -0.801886827021, ‖∇f‖ = 4.0831e-05, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 65, time 259.48 s: f = -0.801886871471, ‖∇f‖ = 4.1034e-05, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 66, time 261.08 s: f = -0.801886949561, ‖∇f‖ = 5.1171e-05, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 67, time 262.97 s: f = -0.801887066611, ‖∇f‖ = 4.5902e-05, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 68, time 265.75 s: f = -0.801887172302, ‖∇f‖ = 7.4809e-05, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 69, time 269.35 s: f = -0.801887249758, ‖∇f‖ = 3.9621e-05, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 70, time 271.31 s: f = -0.801887292123, ‖∇f‖ = 1.3999e-05, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 71, time 273.21 s: f = -0.801887312573, ‖∇f‖ = 1.0813e-05, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 72, time 275.65 s: f = -0.801887349746, ‖∇f‖ = 1.1335e-05, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 73, time 278.92 s: f = -0.801887427063, ‖∇f‖ = 1.9028e-05, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 74, time 281.71 s: f = -0.801887495619, ‖∇f‖ = 1.9287e-05, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: iter 75, time 286.27 s: f = -0.801887521348, ‖∇f‖ = 2.5285e-05, α = 5.43e-01, m = 32, nfg = 2 +[ Info: LBFGS: iter 76, time 289.89 s: f = -0.801887560128, ‖∇f‖ = 2.7097e-05, α = 1.00e+00, m = 32, nfg = 1 +[ Info: LBFGS: converged after 77 iterations and time 292.11 s: f = -0.801887571683, ‖∇f‖ = 7.8821e-06 ```` @@ -370,7 +386,7 @@ the final value of the cost function we have just optimized. ```` ```` --0.8018875356699146 +-0.8018875716825512 ```` As another check, we can compute the magnetization per site and compare it to a [reference @@ -388,7 +404,7 @@ m_ref = 0.667162 ```` ```` -0.00011315233182807027 +0.00011017887044251218 ```` --- diff --git a/docs/src/examples/3d_ising_partition_function/main.ipynb b/docs/src/examples/3d_ising_partition_function/main.ipynb index 508a4b750..5f171655f 100644 --- a/docs/src/examples/3d_ising_partition_function/main.ipynb +++ b/docs/src/examples/3d_ising_partition_function/main.ipynb @@ -137,26 +137,40 @@ "## Contracting the partition function\n", "\n", "To contract our infinite 3D partition function, we first reinterpret it as an infinite power\n", - "of a slice-to-slice transfer operator $T$, where $T$ can be seen as an infinite 2D\n", + "of a slice-to-slice transfer operator $\\mathbb{T}$, where $\\mathbb{T}$ can be seen as an infinite 2D\n", "projected entangled-pair operator (PEPO) which consists of the rank-6 tensor `O` at each\n", - "site of an infinite 2D square lattice. In the same spirit as the boundary MPS approach, all\n", - "we need to contract the whole partition function is to find the leading eigenvector of this\n", - "PEPO. The fixed point of such a PEPO can be parametrized as a PEPS, and for the case of a\n", - "Hermitian transfer operator we can find this PEPS through [variational optimization](@cite\n", - "vanderstraeten_residual_2018).\n", + "site of an infinite 2D square lattice,\n", "\n", - "Indeed, for a Hermitian transfer operator $T$ we can characterize the fixed point PEPS\n", + "\n", + "
\n", + "\"pepo\"\n", + "
\n", + "\n", + "\n", + "To contract the original infinite network, all we need to do is to find the leading\n", + "eigenvector of the PEPO $\\mathbb{T}$, The fixed point of such a PEPO can be parametrized\n", + "as a PEPS $\\psi$, which should then satisfy the eigenvalue equation\n", + "$\\mathbb{T} |\\psi\\rangle = \\Lambda |\\psi\\rangle$, or diagrammatically:\n", + "\n", + "\n", + "
\n", + "\"pepo\n", + "
\n", + "\n", + "\n", + "For a Hermitian transfer operator $\\mathbb{T}$ we can characterize the fixed point PEPS\n", "$|\\psi\\rangle$ which satisfies the eigenvalue equation\n", - "$T |\\psi\\rangle = \\Lambda |\\psi\\rangle$ corresponding to the largest magnitude eigenvalue\n", - "$\\Lambda$ as the solution of a variational problem\n", + "$\\mathbb{T} |\\psi\\rangle = \\Lambda |\\psi\\rangle$ corresponding to the largest magnitude\n", + "eigenvalue $\\Lambda$ as the solution of a variational problem\n", "\n", "$$\n", - "|\\psi\\rangle = \\text{argmin}_{|\\psi\\rangle} \\left ( \\lim_{N \\to ∞} - \\frac{1}{N} \\log \\left( \\frac{\\langle \\psi | T | \\psi \\rangle}{\\langle \\psi | \\psi \\rangle} \\right) \\right ) ,\n", + "|\\psi\\rangle = \\text{argmin}_{|\\psi\\rangle} \\left ( \\lim_{N \\to ∞} - \\frac{1}{N} \\log \\left( \\frac{\\langle \\psi | \\mathbb{T} | \\psi \\rangle}{\\langle \\psi | \\psi \\rangle} \\right) \\right ) ,\n", "$$\n", "\n", - "where $N$ is the diverging number of sites of the 2D transfer operator $T$. The function\n", - "minimized in this expression is exactly the free energy per site of the partition function,\n", - "so we essentially find the fixed-point PEPS by variationally minimizing the free energy.\n", + "where $N$ is the diverging number of sites of the 2D transfer operator $\\mathbb{T}$. The function\n", + "minimized in this expression is exactly the free energy per site of the partition function.\n", + "This means we can directly find the fixed-point PEPS by\n", + "[variationally minimizing the free energy](@cite vanderstraeten_residual_2018).\n", "\n", "### Defining the cost function\n", "\n", @@ -165,7 +179,7 @@ "immediately recognize the denominator $\\langle \\psi | \\psi \\rangle$ as the familiar PEPS\n", "norm, where we can compute the norm per site as the `network_value` of the\n", "corresponding `InfiniteSquareNetwork` by contracting it with the CTMRG algorithm.\n", - "Similarly, the numerator $\\langle \\psi | T | \\psi \\rangle$ is nothing more than an\n", + "Similarly, the numerator $\\langle \\psi | \\mathbb{T} | \\psi \\rangle$ is nothing more than an\n", "`InfiniteSquareNetwork` consisting of three layers corresponding to the ket, transfer\n", "operator and bra objects. This object can also be constructed and contracted in a\n", "straightforward way, so we can again compute its `network_value`.\n", @@ -395,11 +409,11 @@ "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", - "version": "1.11.5" + "version": "1.10.4" }, "kernelspec": { - "name": "julia-1.11", - "display_name": "Julia 1.11.5", + "name": "julia-1.10", + "display_name": "Julia 1.10.4", "language": "julia" } }, diff --git a/docs/src/examples/boundary_mps/index.md b/docs/src/examples/boundary_mps/index.md index 0b15b13b6..8ddeeb1ad 100644 --- a/docs/src/examples/boundary_mps/index.md +++ b/docs/src/examples/boundary_mps/index.md @@ -32,70 +32,125 @@ using TensorKit, PEPSKit, MPSKit ## Computing a PEPS norm -We start by initializing a random infinite PEPS. Let us use uniformly distributed complex -entries using `rand` (which sometimes lead to better convergence than Gaussian distributed -`randn` elements): +We start by initializing a random infinite PEPS. Let us use normally distributed complex +entries using `randn`: ````julia -peps₀ = InfinitePEPS(rand, ComplexF64, ComplexSpace(2), ComplexSpace(2)) +ψ = InfinitePEPS(randn, ComplexF64, ComplexSpace(2), ComplexSpace(2)) ```` ```` InfinitePEPS{TensorKit.TensorMap{ComplexF64, TensorKit.ComplexSpace, 1, 4, Vector{ComplexF64}}}(TensorKit.TensorMap{ComplexF64, TensorKit.ComplexSpace, 1, 4, Vector{ComplexF64}}[TensorMap(ℂ^2 ← (ℂ^2 ⊗ ℂ^2 ⊗ (ℂ^2)' ⊗ (ℂ^2)')): [:, :, 1, 1, 1] = - 0.8343040072662887 + 0.15425705836788395im 0.4612746978522435 + 0.7411151918989216im - 0.6640771294125087 + 0.4428356798799721im 0.9163597170532635 + 0.24145695415210522im + -0.5524390176345264 - 0.07357188568178248im 0.34014501646081047 - 0.7552574870030472im + -0.5455245317233405 + 0.8946618856309984im 1.249282911658007 + 0.45352274131986825im [:, :, 2, 1, 1] = - 0.44289651954161835 + 0.5968081052313008im 0.5473659268881094 + 0.37528062658773764im - 0.00644367423621961 + 0.9414462569909486im 0.36006028879229457 + 0.6157267258321241im + 0.33621043661988675 + 0.4400876608299719im -0.9866664087107284 - 0.28688827761325675im + -0.0077250067072679235 + 1.7380910495900947im -0.19071062901939098 - 1.1367500834118434im [:, :, 1, 2, 1] = - 0.04956065285909117 + 0.26119820734171617im 0.9153298540884296 + 0.3990244910357601im - 0.17944112964295234 + 0.4233545106724528im 0.020358359069476473 + 0.6501897922267199im + -0.09149850722392933 + 0.3560942836258964im 1.6255618447281441 - 0.5689426732891244im + -0.19309251474097275 - 0.32363899914302613im -0.025356816648697236 + 0.5632279168368712im [:, :, 2, 2, 1] = - 0.040493161136161526 + 0.03501665486055905im 0.2591040734810338 + 0.8830094105726012im - 0.781658280511654 + 0.9662812119384394im 0.8169988652653896 + 0.674481616952991im + 0.07675114584269166 - 0.011479824536308164im -0.17779977372973318 + 1.1379201927122535im + -1.0116302866282385 - 0.9253070687198848im 1.1649047337212566 + 0.9936369101208083im [:, :, 1, 1, 2] = - 0.2242833355717867 + 0.14929928451790686im 0.6883051212688887 + 0.588769359105893im - 0.046322385671192734 + 0.8543796191082029im 0.6437874016748227 + 0.257253015722232im + 0.2510676919806213 - 0.182052326055189im -0.5792402993550532 - 0.4309109406268341im + 0.04501645227038913 - 0.8140971172854408im -0.5608346802110794 + 0.21262550530307248im [:, :, 2, 1, 2] = - 0.8719996187768273 + 0.25052026742300637im 0.5714417314833022 + 0.9944321644519715im - 0.4273547968422168 + 0.6068478826937488im 0.4946426302106661 + 0.8353867377249198im + 1.5061767210554262 + 0.17190948125245623im -0.8001234458239143 + 0.6764943808639017im + -0.8176938467062373 - 0.40919675695722396im -0.6692181340575689 + 0.6923370271564298im [:, :, 1, 2, 2] = - 0.6857354516279699 + 0.0952105576480895im 0.14591923452838773 + 0.0853564870114002im - 0.6779060054394721 + 0.4947495895268207im 0.9280821860668365 + 0.931316796924268im + -0.16556382071485704 + 0.2540132491548349im 0.05546115732751907 + 0.3723175507964387im + -0.29883021417599165 - 0.07229462525164528im -1.200173153698329 - 0.45509299328832953im [:, :, 2, 2, 2] = - 0.3716373366637086 + 0.2556099109043021im 0.7954831107819061 + 0.016909518973250215im - 0.9376161032047406 + 0.6320864548041844im 0.7900851372111909 + 0.5457560526661245im + 0.289873563752043 + 0.44718981087960125im 0.018357838612906643 + 0.9634127683557584im + 0.5128282969211142 - 0.2865462937979091im -0.44278618042821827 + 0.2612084385439659im ;;]) ```` -To compute its norm, usually we would construct a double-layer `InfiniteSquareNetwork` which -encodes the bra-ket PEPS overlap and then contract this infinite square network, for example -using CTMRG. Here however, we will use another approach. If we take out a single row of this -infinite norm network, we can interpret it as a 2D row-to-row transfer operator ``T``. Here, -this transfer operator consists of an effective local rank-4 tensor at every site of a 2D -square lattice, where the local effective tensor is given by the contraction of a bra and -ket [`PEPSKit.PEPSTensor`](@ref) across their physical leg. Since the network we want to -contract can be interpreted as the infinite power of ``T``, we can contract it by finding -its leading eigenvector as a 1D MPS, which we call the boundary MPS. - -In PEPSKit.jl, we can directly construct the transfer operator corresponding to a PEPS norm -network from a given infinite PEPS as an [`InfiniteTransferPEPS`](@ref) object. -Additionally, we need to specify which direction should be facing north (`dir=1` -corresponding to north, counting clockwise) and which row is selected from the north - but -since we have a trivial unit cell there is only one row: +To compute its norm, we have to contract a double-layer network which encodes the bra-ket +PEPS overlap ``\langle ψ | ψ \rangle``: + +```@raw html +
+peps norm network +
+``` + +In PEPSKit.jl, this structure is represented as an [`InfiniteSquareNetwork`](@ref) object, +whose effective local rank-4 constituent tensor is given by the contraction of a pair of bra +and ket [`PEPSKit.PEPSTensor`](@ref)s across their physical legs. Until now, we have always +contracted such a network using the CTMRG algorithm. Here however, we will use another +approach. + +If we take out a single row of this infinite norm network, we can interpret it as a 1D +row-to-row transfer operator ``\mathbb{T}``, + +```@raw html +
+peps transfer operator +
+``` + +This transfer operator can be seen as an infinite chain of the effective local rank-4 +tensors that make up the PEPS norm network. Since the network we want to contract can be +interpreted as the infinite power of ``\mathbb{T}``, we can contract it by finding its +leading eigenvector as a 1D MPS ``| \psi_{\text{MPS}} \rangle``, which we call the boundary +MPS. This boundary MPS should satisfy the eigenvalue equation +``\mathbb{T} | \psi_{\text{MPS}} \rangle \approx \Lambda | \psi_{\text{MPS}} \rangle``, or +diagrammatically: + +```@raw html +
+peps transfer fixedpoint equation +
+``` + +Note that if ``\mathbb{T}`` is Hermitian, we can formulate this eigenvalue equation in terms of a +variational problem for the free energy, + +```math +\begin{align} +f &= \lim_{N \to ∞} - \frac{1}{N} \log \left( \frac{\langle \psi_{\text{MPS}} | \mathbb{T} | \psi_{\text{MPS}} \rangle}{\langle \psi_{\text{MPS}} | \psi_{\text{MPS}} \rangle} \right), +\\ +&= -\log(\lambda) +\end{align} +``` + +where ``\lambda = \Lambda^{1/N}`` is the 'eigenvalue per site' of ``\mathbb{T}``, giving +``f`` the meaning of a free energy density. + +Since the contraction of a PEPS norm network is in essence exactly the same problem as the +contraction of a 2D classical partition function, we can directly use boundary MPS +algorithms designed for 2D statistical mechanics models in this context. In particular, +we'll use the [the VUMPS algorithm](@cite vanderstraeten_tangentspace_2019) to perform the +boundary MPS contraction, and we'll call it through the [`leading_boundary`](@ref) method +from MPSKit.jl. This method precisely finds the MPS fixed point of a 1D transfer operator. + +## Boundary MPS contractions with PEPSKit.jl + +To use [`leading_boundary`](@ref), we first need to contruct the transfer operator +``\mathbb{T}`` as an [`MPSKit.InfiniteMPO`](@extref) object. In PEPSKit.jl, we can directly +construct the transfer operator corresponding to a PEPS norm network from a given infinite +PEPS as an [`InfiniteTransferPEPS`](@ref) object, which is a specific kind of +[`MPSKit.InfiniteMPO`](@extref). + +To construct a 1D transfer operator from a 2D PEPS state, we need to specify which direction +should be facing north (`dir=1` corresponding to north, counting clockwise) and which row of +the network is selected from the north - but since we have a trivial unit cell there is only +one row here: ````julia dir = 1 ## does not rotate the partition function row = 1 -T = InfiniteTransferPEPS(peps₀, dir, row) +T = InfiniteTransferPEPS(ψ, dir, row) ```` ```` @@ -106,8 +161,8 @@ single site MPSKit.InfiniteMPO{Tuple{TensorKit.TensorMap{ComplexF64, TensorKit.C ```` -Since we'll find the leading eigenvector of ``T`` as a boundary MPS, we first need to -construct an initial guess to supply to our algorithm. We can do this using the +Since we'll find the leading eigenvector of ``\mathbb{T}`` as a boundary MPS, we first need +to construct an initial guess to supply to our algorithm. We can do this using the [`initialize_mps`](@ref) function, which constructs a random MPS with a specific virtual space for a given transfer operator. Here, we'll build an initial guess for the boundary MPS with a bond dimension of 20: @@ -129,43 +184,42 @@ Note that this will just construct a MPS with random Gaussian entries based on t spaces of the supplied transfer operator. Of course, one might come up with a better initial guess (leading to better convergence) depending on the application. To find the leading boundary MPS fixed point, we call [`leading_boundary`](@ref) using the -[`MPSKit.VUMPS`](@extref) algorithm from MPSKit: +[`MPSKit.VUMPS`](@extref) algorithm: ````julia mps, env, ϵ = leading_boundary(mps₀, T, VUMPS(; tol=1e-6, verbosity=2)); ```` ```` -[ Info: VUMPS init: obj = +5.052950412844e+00 +1.493192627823e-02im err = 8.4684e-01 -[ Info: VUMPS conv 4: obj = +1.744071149764e+01 +2.162482270345e-08im err = 2.0274771848e-07 time = 4.39 sec +[ Info: VUMPS init: obj = +1.674563752306e+00 +3.035692829590e+00im err = 7.5576e-01 +[ Info: VUMPS conv 120: obj = +6.831610878310e+00 -9.694386191727e-09im err = 9.5145748780e-07 time = 25.25 sec ```` The norm of the state per unit cell is then given by the expectation value -$\langle \psi_\text{MPS} | \mathbb{T} | \psi_\text{MPS} \rangle$: +$\langle \psi_\text{MPS} | \mathbb{T} | \psi_\text{MPS} \rangle$ per site: ````julia norm_vumps = abs(prod(expectation_value(mps, T))) ```` ```` -17.440711497641633 +6.831610878309706 ```` -This can be compared to the result obtained using CTMRG, where we see that the results match: +This can be compared to the result obtained using CTMRG, where we see that the results +match: ````julia -env_ctmrg, = leading_boundary( - CTMRGEnv(peps₀, ComplexSpace(20)), peps₀; tol=1e-6, verbosity=2 -) -norm_ctmrg = abs(norm(peps₀, env_ctmrg)) +env_ctmrg, = leading_boundary(CTMRGEnv(ψ, ComplexSpace(20)), ψ; tol=1e-6, verbosity=2) +norm_ctmrg = abs(norm(ψ, env_ctmrg)) @show abs(norm_vumps - norm_ctmrg) / norm_vumps; ```` ```` -[ Info: CTMRG init: obj = -5.556349490423e-01 +1.605938670370e+00im err = 1.0000e+00 -[ Info: CTMRG conv 37: obj = +1.744071151099e+01 err = 3.2046524013e-07 time = 0.19 sec -abs(norm_vumps - norm_ctmrg) / norm_vumps = 7.653151774414538e-10 +[ Info: CTMRG init: obj = -1.530193898578e+01 +2.533005049756e-01im err = 1.0000e+00 +[ Info: CTMRG conv 30: obj = +6.831603585666e+00 err = 4.7443353447e-07 time = 3.78 sec +abs(norm_vumps - norm_ctmrg) / norm_vumps = 1.0674852619316525e-6 ```` @@ -174,19 +228,18 @@ abs(norm_vumps - norm_ctmrg) / norm_vumps = 7.653151774414538e-10 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 rows or 'lines', each of which can be represented by an [`InfiniteTransferPEPS`](@ref) object. Such -a multi-line transfer operator is represented by a `MultilineTransferPEPS` object. In this -case, the boundary MPS is an [`MultilineMPS`](@extref) object, which should be initialized -by specifying a virtual space for each site in the partition function unit cell. +a multi-line transfer operator is represented by a [`PEPSKit.MultilineTransferPEPS`](@ref) +object. In this case, the boundary MPS is an [`MultilineMPS`](@extref) object, which should +be initialized by specifying a virtual space for each site in the partition function unit +cell. First, we construct a PEPS with a $2 \times 2$ unit cell using the `unitcell` keyword argument and then define the corresponding transfer operator, where we again specify the direction which will be facing north: ````julia -peps₀_2x2 = InfinitePEPS( - rand, ComplexF64, ComplexSpace(2), ComplexSpace(2); unitcell=(2, 2) -) -T_2x2 = PEPSKit.MultilineTransferPEPS(peps₀_2x2, dir); +ψ_2x2 = InfinitePEPS(rand, ComplexF64, ComplexSpace(2), ComplexSpace(2); unitcell=(2, 2)) +T_2x2 = PEPSKit.MultilineTransferPEPS(ψ_2x2, dir); ```` Now, the procedure is the same as before: We compute the norm once using VUMPS, once using CTMRG and then compare. @@ -197,20 +250,20 @@ mps_2x2, = leading_boundary(mps₀_2x2, T_2x2, VUMPS(; tol=1e-6, verbosity=2)) norm_2x2_vumps = abs(prod(expectation_value(mps_2x2, T_2x2))) env_ctmrg_2x2, = leading_boundary( - CTMRGEnv(peps₀_2x2, ComplexSpace(20)), peps₀_2x2; tol=1e-6, verbosity=2 + CTMRGEnv(ψ_2x2, ComplexSpace(20)), ψ_2x2; tol=1e-6, verbosity=2 ) -norm_2x2_ctmrg = abs(norm(peps₀_2x2, env_ctmrg_2x2)) +norm_2x2_ctmrg = abs(norm(ψ_2x2, env_ctmrg_2x2)) @show abs(norm_2x2_vumps - norm_2x2_ctmrg) / norm_2x2_vumps; ```` ```` -[ Info: VUMPS init: obj = +6.668046237341e+02 -1.267878277078e+01im err = 8.7901e-01 -┌ Warning: VUMPS cancel 200: obj = +9.713547617328e+04 -3.270767578613e+02im err = 4.7663743999e-02 time = 8.25 sec -└ @ MPSKit ~/.julia/packages/MPSKit/XpTWn/src/algorithms/groundstate/vumps.jl:73 -[ Info: CTMRG init: obj = +1.074898090007e+03 -2.096255594496e+02im err = 1.0000e+00 -[ Info: CTMRG conv 41: obj = +9.723959008610e+04 err = 6.0700367621e-07 time = 0.66 sec -abs(norm_2x2_vumps - norm_2x2_ctmrg) / norm_2x2_vumps = 0.0010661671847507165 +[ Info: VUMPS init: obj = +9.819004336879e+02 +3.170023011746e+01im err = 8.6013e-01 +┌ Warning: VUMPS cancel 200: obj = +1.050428151028e+05 -1.899053741715e+02im err = 2.8823293066e-02 time = 2.71 min +└ @ MPSKit ~/.julia/packages/MPSKit/cKwp2/src/algorithms/groundstate/vumps.jl:76 +[ Info: CTMRG init: obj = -5.668964998506e+01 -2.489382614044e+01im err = 1.0000e+00 +[ Info: CTMRG conv 82: obj = +1.046633714846e+05 err = 7.6507920704e-07 time = 36.93 sec +abs(norm_2x2_vumps - norm_2x2_ctmrg) / norm_2x2_vumps = 0.0036139041100238267 ```` @@ -223,8 +276,8 @@ case the CTMRG algorithm is recommended. Using exactly the same machinery, we can contract 2D networks 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](@ref e_3d_ising) -with our random PEPS from before and evaluate the overlap $\langle \psi_\text{PEPS} | -O_\text{PEPO} | \psi_\text{PEPS} \rangle$. +with our random PEPS from before and evaluate the overlap $\langle \psi | +T | \psi \rangle$. The classical Ising PEPO is defined as follows: @@ -248,8 +301,8 @@ To evaluate the overlap, we instantiate the PEPO and the corresponding [`Infinit in the right direction, on the right row of the partition function (trivial here): ````julia -pepo = ising_pepo(1) -transfer_pepo = InfiniteTransferPEPO(peps₀, pepo, 1, 1) +T = ising_pepo(1) +transfer_pepo = InfiniteTransferPEPO(ψ, T, 1, 1) ```` ```` @@ -270,9 +323,10 @@ norm_pepo = abs(prod(expectation_value(mps_pepo, transfer_pepo))); ```` ```` -[ Info: VUMPS init: obj = +3.309203535702e+01 -4.227375981212e-01im err = 9.3280e-01 -[ Info: VUMPS conv 5: obj = +2.483696258643e+02 +2.405968608492e-07im err = 5.0410779107e-08 time = 3.38 sec -norm_pepo = 248.36962586426358 +[ Info: VUMPS init: obj = +2.552585816795e+01 -2.455978462565e-01im err = 8.9526e-01 +┌ Warning: VUMPS cancel 200: obj = +6.211569029510e+01 +1.206390478837e+00im err = 6.9644987097e-01 time = 3.12 min +└ @ MPSKit ~/.julia/packages/MPSKit/cKwp2/src/algorithms/groundstate/vumps.jl:76 +norm_pepo = 62.12740424984486 ```` diff --git a/docs/src/examples/boundary_mps/main.ipynb b/docs/src/examples/boundary_mps/main.ipynb index b78452c21..db4ba69cc 100644 --- a/docs/src/examples/boundary_mps/main.ipynb +++ b/docs/src/examples/boundary_mps/main.ipynb @@ -57,9 +57,8 @@ "source": [ "## Computing a PEPS norm\n", "\n", - "We start by initializing a random infinite PEPS. Let us use uniformly distributed complex\n", - "entries using `rand` (which sometimes lead to better convergence than Gaussian distributed\n", - "`randn` elements):" + "We start by initializing a random infinite PEPS. Let us use normally distributed complex\n", + "entries using `randn`:" ], "metadata": {} }, @@ -67,7 +66,7 @@ "outputs": [], "cell_type": "code", "source": [ - "peps₀ = InfinitePEPS(rand, ComplexF64, ComplexSpace(2), ComplexSpace(2))" + "ψ = InfinitePEPS(randn, ComplexF64, ComplexSpace(2), ComplexSpace(2))" ], "metadata": {}, "execution_count": null @@ -75,21 +74,77 @@ { "cell_type": "markdown", "source": [ - "To compute its norm, usually we would construct a double-layer `InfiniteSquareNetwork` which\n", - "encodes the bra-ket PEPS overlap and then contract this infinite square network, for example\n", - "using CTMRG. Here however, we will use another approach. If we take out a single row of this\n", - "infinite norm network, we can interpret it as a 2D row-to-row transfer operator $T$. Here,\n", - "this transfer operator consists of an effective local rank-4 tensor at every site of a 2D\n", - "square lattice, where the local effective tensor is given by the contraction of a bra and\n", - "ket `PEPSKit.PEPSTensor` across their physical leg. Since the network we want to\n", - "contract can be interpreted as the infinite power of $T$, we can contract it by finding\n", - "its leading eigenvector as a 1D MPS, which we call the boundary MPS.\n", - "\n", - "In PEPSKit.jl, we can directly construct the transfer operator corresponding to a PEPS norm\n", - "network from a given infinite PEPS as an `InfiniteTransferPEPS` object.\n", - "Additionally, we need to specify which direction should be facing north (`dir=1`\n", - "corresponding to north, counting clockwise) and which row is selected from the north - but\n", - "since we have a trivial unit cell there is only one row:" + "To compute its norm, we have to contract a double-layer network which encodes the bra-ket\n", + "PEPS overlap $\\langle ψ | ψ \\rangle$:\n", + "\n", + "\n", + "
\n", + "\"peps\n", + "
\n", + "\n", + "\n", + "In PEPSKit.jl, this structure is represented as an `InfiniteSquareNetwork` object,\n", + "whose effective local rank-4 constituent tensor is given by the contraction of a pair of bra\n", + "and ket `PEPSKit.PEPSTensor`s across their physical legs. Until now, we have always\n", + "contracted such a network using the CTMRG algorithm. Here however, we will use another\n", + "approach.\n", + "\n", + "If we take out a single row of this infinite norm network, we can interpret it as a 1D\n", + "row-to-row transfer operator $\\mathbb{T}$,\n", + "\n", + "\n", + "
\n", + "\"peps\n", + "
\n", + "\n", + "\n", + "This transfer operator can be seen as an infinite chain of the effective local rank-4\n", + "tensors that make up the PEPS norm network. Since the network we want to contract can be\n", + "interpreted as the infinite power of $\\mathbb{T}$, we can contract it by finding its\n", + "leading eigenvector as a 1D MPS $| \\psi_{\\text{MPS}} \\rangle$, which we call the boundary\n", + "MPS. This boundary MPS should satisfy the eigenvalue equation\n", + "$\\mathbb{T} | \\psi_{\\text{MPS}} \\rangle \\approx \\Lambda | \\psi_{\\text{MPS}} \\rangle$, or\n", + "diagrammatically:\n", + "\n", + "\n", + "
\n", + "\"peps\n", + "
\n", + "\n", + "\n", + "Note that if $\\mathbb{T}$ is Hermitian, we can formulate this eigenvalue equation in terms of a\n", + "variational problem for the free energy,\n", + "\n", + "$$\n", + "\\begin{align}\n", + "f &= \\lim_{N \\to ∞} - \\frac{1}{N} \\log \\left( \\frac{\\langle \\psi_{\\text{MPS}} | \\mathbb{T} | \\psi_{\\text{MPS}} \\rangle}{\\langle \\psi_{\\text{MPS}} | \\psi_{\\text{MPS}} \\rangle} \\right),\n", + "\\\\\n", + "&= -\\log(\\lambda)\n", + "\\end{align}\n", + "$$\n", + "\n", + "where $\\lambda = \\Lambda^{1/N}$ is the 'eigenvalue per site' of $\\mathbb{T}$, giving\n", + "$f$ the meaning of a free energy density.\n", + "\n", + "Since the contraction of a PEPS norm network is in essence exactly the same problem as the\n", + "contraction of a 2D classical partition function, we can directly use boundary MPS\n", + "algorithms designed for 2D statistical mechanics models in this context. In particular,\n", + "we'll use the [the VUMPS algorithm](@cite vanderstraeten_tangentspace_2019) to perform the\n", + "boundary MPS contraction, and we'll call it through the `leading_boundary` method\n", + "from MPSKit.jl. This method precisely finds the MPS fixed point of a 1D transfer operator.\n", + "\n", + "## Boundary MPS contractions with PEPSKit.jl\n", + "\n", + "To use `leading_boundary`, we first need to contruct the transfer operator\n", + "$\\mathbb{T}$ as an `MPSKit.InfiniteMPO` object. In PEPSKit.jl, we can directly\n", + "construct the transfer operator corresponding to a PEPS norm network from a given infinite\n", + "PEPS as an `InfiniteTransferPEPS` object, which is a specific kind of\n", + "`MPSKit.InfiniteMPO`.\n", + "\n", + "To construct a 1D transfer operator from a 2D PEPS state, we need to specify which direction\n", + "should be facing north (`dir=1` corresponding to north, counting clockwise) and which row of\n", + "the network is selected from the north - but since we have a trivial unit cell there is only\n", + "one row here:" ], "metadata": {} }, @@ -99,7 +154,7 @@ "source": [ "dir = 1 ## does not rotate the partition function\n", "row = 1\n", - "T = InfiniteTransferPEPS(peps₀, dir, row)" + "T = InfiniteTransferPEPS(ψ, dir, row)" ], "metadata": {}, "execution_count": null @@ -107,8 +162,8 @@ { "cell_type": "markdown", "source": [ - "Since we'll find the leading eigenvector of $T$ as a boundary MPS, we first need to\n", - "construct an initial guess to supply to our algorithm. We can do this using the\n", + "Since we'll find the leading eigenvector of $\\mathbb{T}$ as a boundary MPS, we first need\n", + "to construct an initial guess to supply to our algorithm. We can do this using the\n", "`initialize_mps` function, which constructs a random MPS with a specific virtual\n", "space for a given transfer operator. Here, we'll build an initial guess for the boundary MPS\n", "with a bond dimension of 20:" @@ -131,7 +186,7 @@ "spaces of the supplied transfer operator. Of course, one might come up with a better initial\n", "guess (leading to better convergence) depending on the application. To find the leading\n", "boundary MPS fixed point, we call `leading_boundary` using the\n", - "`MPSKit.VUMPS` algorithm from MPSKit:" + "`MPSKit.VUMPS` algorithm:" ], "metadata": {} }, @@ -148,7 +203,7 @@ "cell_type": "markdown", "source": [ "The norm of the state per unit cell is then given by the expectation value\n", - "$\\langle \\psi_\\text{MPS} | \\mathbb{T} | \\psi_\\text{MPS} \\rangle$:" + "$\\langle \\psi_\\text{MPS} | \\mathbb{T} | \\psi_\\text{MPS} \\rangle$ per site:" ], "metadata": {} }, @@ -164,7 +219,8 @@ { "cell_type": "markdown", "source": [ - "This can be compared to the result obtained using CTMRG, where we see that the results match:" + "This can be compared to the result obtained using CTMRG, where we see that the results\n", + "match:" ], "metadata": {} }, @@ -172,10 +228,8 @@ "outputs": [], "cell_type": "code", "source": [ - "env_ctmrg, = leading_boundary(\n", - " CTMRGEnv(peps₀, ComplexSpace(20)), peps₀; tol=1e-6, verbosity=2\n", - ")\n", - "norm_ctmrg = abs(norm(peps₀, env_ctmrg))\n", + "env_ctmrg, = leading_boundary(CTMRGEnv(ψ, ComplexSpace(20)), ψ; tol=1e-6, verbosity=2)\n", + "norm_ctmrg = abs(norm(ψ, env_ctmrg))\n", "@show abs(norm_vumps - norm_ctmrg) / norm_vumps;" ], "metadata": {}, @@ -189,9 +243,10 @@ "For PEPS with non-trivial unit cells, the principle is exactly the same. The only difference\n", "is that now the transfer operator of the PEPS norm partition function has multiple rows or\n", "'lines', each of which can be represented by an `InfiniteTransferPEPS` object. Such\n", - "a multi-line transfer operator is represented by a `MultilineTransferPEPS` object. In this\n", - "case, the boundary MPS is an `MultilineMPS` object, which should be initialized\n", - "by specifying a virtual space for each site in the partition function unit cell.\n", + "a multi-line transfer operator is represented by a `PEPSKit.MultilineTransferPEPS`\n", + "object. In this case, the boundary MPS is an `MultilineMPS` object, which should\n", + "be initialized by specifying a virtual space for each site in the partition function unit\n", + "cell.\n", "\n", "First, we construct a PEPS with a $2 \\times 2$ unit cell using the `unitcell` keyword\n", "argument and then define the corresponding transfer operator, where we again specify the\n", @@ -203,10 +258,8 @@ "outputs": [], "cell_type": "code", "source": [ - "peps₀_2x2 = InfinitePEPS(\n", - " rand, ComplexF64, ComplexSpace(2), ComplexSpace(2); unitcell=(2, 2)\n", - ")\n", - "T_2x2 = PEPSKit.MultilineTransferPEPS(peps₀_2x2, dir);" + "ψ_2x2 = InfinitePEPS(rand, ComplexF64, ComplexSpace(2), ComplexSpace(2); unitcell=(2, 2))\n", + "T_2x2 = PEPSKit.MultilineTransferPEPS(ψ_2x2, dir);" ], "metadata": {}, "execution_count": null @@ -227,9 +280,9 @@ "norm_2x2_vumps = abs(prod(expectation_value(mps_2x2, T_2x2)))\n", "\n", "env_ctmrg_2x2, = leading_boundary(\n", - " CTMRGEnv(peps₀_2x2, ComplexSpace(20)), peps₀_2x2; tol=1e-6, verbosity=2\n", + " CTMRGEnv(ψ_2x2, ComplexSpace(20)), ψ_2x2; tol=1e-6, verbosity=2\n", ")\n", - "norm_2x2_ctmrg = abs(norm(peps₀_2x2, env_ctmrg_2x2))\n", + "norm_2x2_ctmrg = abs(norm(ψ_2x2, env_ctmrg_2x2))\n", "\n", "@show abs(norm_2x2_vumps - norm_2x2_ctmrg) / norm_2x2_vumps;" ], @@ -248,8 +301,8 @@ "Using exactly the same machinery, we can contract 2D networks which encode the expectation\n", "value of a PEPO for a given PEPS state. As an example, we can consider the overlap of the\n", "PEPO correponding to the partition function of 3D classical Ising model\n", - "with our random PEPS from before and evaluate the overlap $\\langle \\psi_\\text{PEPS} |\n", - "O_\\text{PEPO} | \\psi_\\text{PEPS} \\rangle$.\n", + "with our random PEPS from before and evaluate the overlap $\\langle \\psi |\n", + "T | \\psi \\rangle$.\n", "\n", "The classical Ising PEPO is defined as follows:" ], @@ -288,8 +341,8 @@ "outputs": [], "cell_type": "code", "source": [ - "pepo = ising_pepo(1)\n", - "transfer_pepo = InfiniteTransferPEPO(peps₀, pepo, 1, 1)" + "T = ising_pepo(1)\n", + "transfer_pepo = InfiniteTransferPEPO(ψ, T, 1, 1)" ], "metadata": {}, "execution_count": null @@ -337,11 +390,11 @@ "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", - "version": "1.11.5" + "version": "1.10.4" }, "kernelspec": { - "name": "julia-1.11", - "display_name": "Julia 1.11.5", + "name": "julia-1.10", + "display_name": "Julia 1.10.4", "language": "julia" } }, diff --git a/examples/3d_ising_partition_function/main.jl b/examples/3d_ising_partition_function/main.jl index 98c3e9b23..45693172c 100644 --- a/examples/3d_ising_partition_function/main.jl +++ b/examples/3d_ising_partition_function/main.jl @@ -95,26 +95,40 @@ md""" ## Contracting the partition function To contract our infinite 3D partition function, we first reinterpret it as an infinite power -of a slice-to-slice transfer operator ``T``, where ``T`` can be seen as an infinite 2D +of a slice-to-slice transfer operator ``\mathbb{T}``, where ``\mathbb{T}`` can be seen as an infinite 2D projected entangled-pair operator (PEPO) which consists of the rank-6 tensor `O` at each -site of an infinite 2D square lattice. In the same spirit as the boundary MPS approach, all -we need to contract the whole partition function is to find the leading eigenvector of this -PEPO. The fixed point of such a PEPO can be parametrized as a PEPS, and for the case of a -Hermitian transfer operator we can find this PEPS through [variational optimization](@cite -vanderstraeten_residual_2018). +site of an infinite 2D square lattice, -Indeed, for a Hermitian transfer operator ``T`` we can characterize the fixed point PEPS +```@raw html +
+pepo +
+``` + +To contract the original infinite network, all we need to do is to find the leading +eigenvector of the PEPO ``\mathbb{T}``, The fixed point of such a PEPO can be parametrized +as a PEPS ``\psi``, which should then satisfy the eigenvalue equation +``\mathbb{T} |\psi\rangle = \Lambda |\psi\rangle``, or diagrammatically: + +```@raw html +
+pepo fixedpoint equation +
+``` + +For a Hermitian transfer operator ``\mathbb{T}`` we can characterize the fixed point PEPS ``|\psi\rangle`` which satisfies the eigenvalue equation -``T |\psi\rangle = \Lambda |\psi\rangle`` corresponding to the largest magnitude eigenvalue -``\Lambda`` as the solution of a variational problem +``\mathbb{T} |\psi\rangle = \Lambda |\psi\rangle`` corresponding to the largest magnitude +eigenvalue ``\Lambda`` as the solution of a variational problem ```math -|\psi\rangle = \text{argmin}_{|\psi\rangle} \left ( \lim_{N \to ∞} - \frac{1}{N} \log \left( \frac{\langle \psi | T | \psi \rangle}{\langle \psi | \psi \rangle} \right) \right ) , +|\psi\rangle = \text{argmin}_{|\psi\rangle} \left ( \lim_{N \to ∞} - \frac{1}{N} \log \left( \frac{\langle \psi | \mathbb{T} | \psi \rangle}{\langle \psi | \psi \rangle} \right) \right ) , ``` -where ``N`` is the diverging number of sites of the 2D transfer operator ``T``. The function -minimized in this expression is exactly the free energy per site of the partition function, -so we essentially find the fixed-point PEPS by variationally minimizing the free energy. +where ``N`` is the diverging number of sites of the 2D transfer operator ``\mathbb{T}``. The function +minimized in this expression is exactly the free energy per site of the partition function. +This means we can directly find the fixed-point PEPS by +[variationally minimizing the free energy](@cite vanderstraeten_residual_2018). ### Defining the cost function @@ -123,7 +137,7 @@ use [OptimKit.jl](https://github.com/Jutho/OptimKit.jl) to actually optimize it. immediately recognize the denominator ``\langle \psi | \psi \rangle`` as the familiar PEPS norm, where we can compute the norm per site as the [`network_value`](@ref) of the corresponding [`InfiniteSquareNetwork`](@ref) by contracting it with the CTMRG algorithm. -Similarly, the numerator ``\langle \psi | T | \psi \rangle`` is nothing more than an +Similarly, the numerator ``\langle \psi | \mathbb{T} | \psi \rangle`` is nothing more than an `InfiniteSquareNetwork` consisting of three layers corresponding to the ket, transfer operator and bra objects. This object can also be constructed and contracted in a straightforward way, so we can again compute its `network_value`. diff --git a/examples/Cache.toml b/examples/Cache.toml index 7d57923f6..c761fe45a 100644 --- a/examples/Cache.toml +++ b/examples/Cache.toml @@ -1,8 +1,8 @@ bose_hubbard = "918635ca8e256252637bbdbc50ef909f93b2e7eec2275ba3bf662cb529f3988d" hubbard_su = "d0e20c521ba1965d1db0e38d828512890e0409bc18ea2a92b0abe21767ea0baf" 2d_ising_partition_function = "86bb41ebdc876c1c10b01dfde934c6888fd7025ccafdee921c02039c3b601d3b" -3d_ising_partition_function = "8e6e77412deae6490a7b63f3b7e0eff45363b2c884bca9672c1db87084f87a92" -boundary_mps = "b50d94ad1124dd647097d3ad8483345dc51ee707798d6ff5636ba7763631161f" +3d_ising_partition_function = "0555e05299dbe1f7abd403d845b84141684fd5dfc483d560de53850d94ad11a7" +boundary_mps = "7b944fd76d58dad280302d3693fee61bbef3fdf1768ebe21525bb2742d932efe" heisenberg_su = "01257e88355dd4822bfb380ee1b3ce79b5fea9862f4bee130e89a0d9805e1f3d" xxz = "80e537be77aabf93835dfe721dacc36b96a64142346c26771ebbb65a1c8f2bac" fermi_hubbard = "5836dda79371470b21347ad6952a875b4d48195f28ded9af4a1fa0e5b8f04c82" diff --git a/examples/boundary_mps/main.jl b/examples/boundary_mps/main.jl index 51a4cef73..2fc80221a 100644 --- a/examples/boundary_mps/main.jl +++ b/examples/boundary_mps/main.jl @@ -25,38 +25,94 @@ using TensorKit, PEPSKit, MPSKit md""" ## Computing a PEPS norm -We start by initializing a random infinite PEPS. Let us use uniformly distributed complex -entries using `rand` (which sometimes lead to better convergence than Gaussian distributed -`randn` elements): +We start by initializing a random infinite PEPS. Let us use normally distributed complex +entries using `randn`: """ -peps₀ = InfinitePEPS(rand, ComplexF64, ComplexSpace(2), ComplexSpace(2)) +ψ = InfinitePEPS(randn, ComplexF64, ComplexSpace(2), ComplexSpace(2)) md""" -To compute its norm, usually we would construct a double-layer `InfiniteSquareNetwork` which -encodes the bra-ket PEPS overlap and then contract this infinite square network, for example -using CTMRG. Here however, we will use another approach. If we take out a single row of this -infinite norm network, we can interpret it as a 2D row-to-row transfer operator ``T``. Here, -this transfer operator consists of an effective local rank-4 tensor at every site of a 2D -square lattice, where the local effective tensor is given by the contraction of a bra and -ket [`PEPSKit.PEPSTensor`](@ref) across their physical leg. Since the network we want to -contract can be interpreted as the infinite power of ``T``, we can contract it by finding -its leading eigenvector as a 1D MPS, which we call the boundary MPS. - -In PEPSKit.jl, we can directly construct the transfer operator corresponding to a PEPS norm -network from a given infinite PEPS as an [`InfiniteTransferPEPS`](@ref) object. -Additionally, we need to specify which direction should be facing north (`dir=1` -corresponding to north, counting clockwise) and which row is selected from the north - but -since we have a trivial unit cell there is only one row: + +To compute its norm, we have to contract a double-layer network which encodes the bra-ket +PEPS overlap ``\langle ψ | ψ \rangle``: + +```@raw html +
+peps norm network +
+``` + +In PEPSKit.jl, this structure is represented as an [`InfiniteSquareNetwork`](@ref) object, +whose effective local rank-4 constituent tensor is given by the contraction of a pair of bra +and ket [`PEPSKit.PEPSTensor`](@ref)s across their physical legs. Until now, we have always +contracted such a network using the CTMRG algorithm. Here however, we will use another +approach. + +If we take out a single row of this infinite norm network, we can interpret it as a 1D +row-to-row transfer operator ``\mathbb{T}``, + +```@raw html +
+peps transfer operator +
+``` + +This transfer operator can be seen as an infinite chain of the effective local rank-4 +tensors that make up the PEPS norm network. Since the network we want to contract can be +interpreted as the infinite power of ``\mathbb{T}``, we can contract it by finding its +leading eigenvector as a 1D MPS ``| \psi_{\text{MPS}} \rangle``, which we call the boundary +MPS. This boundary MPS should satisfy the eigenvalue equation +``\mathbb{T} | \psi_{\text{MPS}} \rangle \approx \Lambda | \psi_{\text{MPS}} \rangle``, or +diagrammatically: + +```@raw html +
+peps transfer fixedpoint equation +
+``` + +Note that if ``\mathbb{T}`` is Hermitian, we can formulate this eigenvalue equation in terms of a +variational problem for the free energy, + +```math +\begin{align} +f &= \lim_{N \to ∞} - \frac{1}{N} \log \left( \frac{\langle \psi_{\text{MPS}} | \mathbb{T} | \psi_{\text{MPS}} \rangle}{\langle \psi_{\text{MPS}} | \psi_{\text{MPS}} \rangle} \right), +\\ +&= -\log(\lambda) +\end{align} +``` + +where ``\lambda = \Lambda^{1/N}`` is the 'eigenvalue per site' of ``\mathbb{T}``, giving +``f`` the meaning of a free energy density. + +Since the contraction of a PEPS norm network is in essence exactly the same problem as the +contraction of a 2D classical partition function, we can directly use boundary MPS +algorithms designed for 2D statistical mechanics models in this context. In particular, +we'll use the [the VUMPS algorithm](@cite vanderstraeten_tangentspace_2019) to perform the +boundary MPS contraction, and we'll call it through the [`leading_boundary`](@ref) method +from MPSKit.jl. This method precisely finds the MPS fixed point of a 1D transfer operator. + +## Boundary MPS contractions with PEPSKit.jl + +To use [`leading_boundary`](@ref), we first need to contruct the transfer operator +``\mathbb{T}`` as an [`MPSKit.InfiniteMPO`](@extref) object. In PEPSKit.jl, we can directly +construct the transfer operator corresponding to a PEPS norm network from a given infinite +PEPS as an [`InfiniteTransferPEPS`](@ref) object, which is a specific kind of +[`MPSKit.InfiniteMPO`](@extref). + +To construct a 1D transfer operator from a 2D PEPS state, we need to specify which direction +should be facing north (`dir=1` corresponding to north, counting clockwise) and which row of +the network is selected from the north - but since we have a trivial unit cell there is only +one row here: """ dir = 1 ## does not rotate the partition function row = 1 -T = InfiniteTransferPEPS(peps₀, dir, row) +T = InfiniteTransferPEPS(ψ, dir, row) md""" -Since we'll find the leading eigenvector of ``T`` as a boundary MPS, we first need to -construct an initial guess to supply to our algorithm. We can do this using the +Since we'll find the leading eigenvector of ``\mathbb{T}`` as a boundary MPS, we first need +to construct an initial guess to supply to our algorithm. We can do this using the [`initialize_mps`](@ref) function, which constructs a random MPS with a specific virtual space for a given transfer operator. Here, we'll build an initial guess for the boundary MPS with a bond dimension of 20: @@ -65,31 +121,29 @@ with a bond dimension of 20: mps₀ = initialize_mps(T, [ComplexSpace(20)]) md""" - Note that this will just construct a MPS with random Gaussian entries based on the physical spaces of the supplied transfer operator. Of course, one might come up with a better initial guess (leading to better convergence) depending on the application. To find the leading boundary MPS fixed point, we call [`leading_boundary`](@ref) using the -[`MPSKit.VUMPS`](@extref) algorithm from MPSKit: +[`MPSKit.VUMPS`](@extref) algorithm: """ mps, env, ϵ = leading_boundary(mps₀, T, VUMPS(; tol=1e-6, verbosity=2)); md""" The norm of the state per unit cell is then given by the expectation value -$\langle \psi_\text{MPS} | \mathbb{T} | \psi_\text{MPS} \rangle$: +$\langle \psi_\text{MPS} | \mathbb{T} | \psi_\text{MPS} \rangle$ per site: """ norm_vumps = abs(prod(expectation_value(mps, T))) md""" -This can be compared to the result obtained using CTMRG, where we see that the results match: +This can be compared to the result obtained using CTMRG, where we see that the results +match: """ -env_ctmrg, = leading_boundary( - CTMRGEnv(peps₀, ComplexSpace(20)), peps₀; tol=1e-6, verbosity=2 -) -norm_ctmrg = abs(norm(peps₀, env_ctmrg)) +env_ctmrg, = leading_boundary(CTMRGEnv(ψ, ComplexSpace(20)), ψ; tol=1e-6, verbosity=2) +norm_ctmrg = abs(norm(ψ, env_ctmrg)) @show abs(norm_vumps - norm_ctmrg) / norm_vumps; md""" @@ -98,19 +152,18 @@ md""" 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 rows or 'lines', each of which can be represented by an [`InfiniteTransferPEPS`](@ref) object. Such -a multi-line transfer operator is represented by a `MultilineTransferPEPS` object. In this -case, the boundary MPS is an [`MultilineMPS`](@extref) object, which should be initialized -by specifying a virtual space for each site in the partition function unit cell. +a multi-line transfer operator is represented by a [`PEPSKit.MultilineTransferPEPS`](@ref) +object. In this case, the boundary MPS is an [`MultilineMPS`](@extref) object, which should +be initialized by specifying a virtual space for each site in the partition function unit +cell. First, we construct a PEPS with a $2 \times 2$ unit cell using the `unitcell` keyword argument and then define the corresponding transfer operator, where we again specify the direction which will be facing north: """ -peps₀_2x2 = InfinitePEPS( - rand, ComplexF64, ComplexSpace(2), ComplexSpace(2); unitcell=(2, 2) -) -T_2x2 = PEPSKit.MultilineTransferPEPS(peps₀_2x2, dir); +ψ_2x2 = InfinitePEPS(rand, ComplexF64, ComplexSpace(2), ComplexSpace(2); unitcell=(2, 2)) +T_2x2 = PEPSKit.MultilineTransferPEPS(ψ_2x2, dir); md""" Now, the procedure is the same as before: We compute the norm once using VUMPS, once using CTMRG and then compare. @@ -121,9 +174,9 @@ mps_2x2, = leading_boundary(mps₀_2x2, T_2x2, VUMPS(; tol=1e-6, verbosity=2)) norm_2x2_vumps = abs(prod(expectation_value(mps_2x2, T_2x2))) env_ctmrg_2x2, = leading_boundary( - CTMRGEnv(peps₀_2x2, ComplexSpace(20)), peps₀_2x2; tol=1e-6, verbosity=2 + CTMRGEnv(ψ_2x2, ComplexSpace(20)), ψ_2x2; tol=1e-6, verbosity=2 ) -norm_2x2_ctmrg = abs(norm(peps₀_2x2, env_ctmrg_2x2)) +norm_2x2_ctmrg = abs(norm(ψ_2x2, env_ctmrg_2x2)) @show abs(norm_2x2_vumps - norm_2x2_ctmrg) / norm_2x2_vumps; @@ -137,8 +190,8 @@ case the CTMRG algorithm is recommended. Using exactly the same machinery, we can contract 2D networks 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](@ref e_3d_ising) -with our random PEPS from before and evaluate the overlap $\langle \psi_\text{PEPS} | -O_\text{PEPO} | \psi_\text{PEPS} \rangle$. +with our random PEPS from before and evaluate the overlap $\langle \psi | +T | \psi \rangle$. The classical Ising PEPO is defined as follows: """ @@ -162,8 +215,8 @@ To evaluate the overlap, we instantiate the PEPO and the corresponding [`Infinit in the right direction, on the right row of the partition function (trivial here): """ -pepo = ising_pepo(1) -transfer_pepo = InfiniteTransferPEPO(peps₀, pepo, 1, 1) +T = ising_pepo(1) +transfer_pepo = InfiniteTransferPEPO(ψ, T, 1, 1) md""" As before, we converge the boundary MPS using VUMPS and then compute the expectation value: