From 8c152d51545279d97e8d9dd8035ebbf2c31ad034 Mon Sep 17 00:00:00 2001 From: leburgel Date: Mon, 10 Mar 2025 12:05:23 +0100 Subject: [PATCH 01/10] Implement norm-preserving retractions --- .../optimization/peps_optimization.jl | 90 +++++++++++++++++-- 1 file changed, 84 insertions(+), 6 deletions(-) diff --git a/src/algorithms/optimization/peps_optimization.jl b/src/algorithms/optimization/peps_optimization.jl index aecb2132b..3fc60dc25 100644 --- a/src/algorithms/optimization/peps_optimization.jl +++ b/src/algorithms/optimization/peps_optimization.jl @@ -105,9 +105,17 @@ function fixedpoint( gradnorms_unitcell = Vector{Matrix{T}}() times = Vector{Float64}() + # normalize the initial guess + peps₀ = peps_normalize(peps₀) + # optimize operator cost function (peps_final, env_final), cost, ∂cost, numfg, convergence_history = optimize( - (peps₀, env₀), alg.optimizer; retract, inner=real_inner, finalize! + (peps₀, env₀), + alg.optimizer; + retract, + inner=real_inner, + finalize!, + (transport!)=(peps_transport!), ) do (peps, env) start_time = time_ns() E, gs = withgradient(peps) do ψ @@ -144,13 +152,83 @@ function fixedpoint( return peps_final, env_final, cost, info end -# Update PEPS unit cell in non-mutating way -# Note: Both x and η are InfinitePEPS during optimization +""" + peps_normalize(A::InfinitePEPS) + +Normalize the individual tensors in the unit cell of an `InfinitePEPS` such that they each +have unit Euclidean norm. +""" +function peps_normalize(A::InfinitePEPS) + normalized_tensors = map(unitcell(A)) do t + return t / norm(t) + end + return InfinitePEPS(normalized_tensors) +end + +""" + peps_retract(x, η, α) + +Performs a norm-preserving retraction of an infinite PEPS `A = x[1]` along `η` with step +size `α`, giving a new PEPS `A´`, +```math +A' \\leftarrow \\cos(α ||η||) A + \\sin(α ||η||) \\frac{η}{||η||}, +``` +and corresponding directional derivative `ξ`, +```math +ξ = \\cos(α ||η||) η - \\sin(α ||η||) ||η|| A, +``` +such that ``\\langle A', ξ \\rangle = 0``. If the initial state was properly normalized, we +have that +```math +||A'|| = ||A|| = 1. +``` +""" function peps_retract(x, η, α) - peps = deepcopy(x[1]) - peps.A .+= η.A .* α + peps = x[1] + norms_η = norm.(η.A) + + peps´ = similar(x[1]) + peps´.A .= cos.(α .* norms_η) .* peps.A .+ sin.(α .* norms_η) .* η.A ./ norms_η + env = deepcopy(x[2]) - return (peps, env), η + + ξ = similar(η) + ξ.A .= cos.(α .* norms_η) .* η.A .- sin.(α .* norms_η) .* norms_η .* peps.A + + return (peps´, env), ξ +end + +""" + peps_transport!(ξ, x, η, α, x′) + +Transports a direction at `A = x[1]` to a valid direction at `A´ = x´[1]` corresponding to +the norm-preserving retraction of `A` along `η` with step size `α`. In particular, starting +from a direction `η` of the form +```math +ξ = \\left\\langle \\frac{η}{||η||}, ξ \\right\\rangle \\frac{η}{||η||} + Δξ +``` +where ``\\langle Δξ, A \\rangle = \\langle Δξ, η \\rangle = 0``, it returns +```math +ξ(α) = \\left\\langle \\frac{η}{||η||}, ξ \\right\\rangle \\left( \\cos(α ||η||) \\frac{η}{||η||} - \\sin(α ||η||) A \\right) + Δξ +``` +such that ``||ξ(α)|| = ||ξ||, \\langle A', ξ(α) \\rangle = 0``. +""" +function peps_transport!(ξ, x, η, α, x´) + peps = x[1] + + norms_η = norm.(η.A) + normalized_η = η.A ./ norms_η + overlaps_η_ξ = inner.(normalized_η, ξ.A) + + # isolate the orthogonal component + Δξ = ξ.A .- overlaps_η_ξ .* normalized_η + + # keep orthogonal component fixed, modify the rest by the proper directional derivative + ξ.A .= + overlaps_η_ξ .* + (cos.(α .* norms_η) .* normalized_η .- sin.(α .* norms_η) .* peps.A) .+ Δξ + + return ξ end # Take real valued part of dot product From 5b4123a21c07d3e4c997df687c295ca97ca99bd0 Mon Sep 17 00:00:00 2001 From: leburgel Date: Mon, 10 Mar 2025 13:59:29 +0100 Subject: [PATCH 02/10] Drop implicit assumption that state has unit norm --- .../optimization/peps_optimization.jl | 28 +++++++++++-------- 1 file changed, 16 insertions(+), 12 deletions(-) diff --git a/src/algorithms/optimization/peps_optimization.jl b/src/algorithms/optimization/peps_optimization.jl index 3fc60dc25..00d2fca55 100644 --- a/src/algorithms/optimization/peps_optimization.jl +++ b/src/algorithms/optimization/peps_optimization.jl @@ -171,29 +171,30 @@ end Performs a norm-preserving retraction of an infinite PEPS `A = x[1]` along `η` with step size `α`, giving a new PEPS `A´`, ```math -A' \\leftarrow \\cos(α ||η||) A + \\sin(α ||η||) \\frac{η}{||η||}, +A' \\leftarrow \\cos \\left( α \\frac{||η||}{||A||} \\right) A + \\sin \\left( α \\frac{||η||}{||A||} \\right) ||A|| \\frac{η}{||η||}, ``` and corresponding directional derivative `ξ`, ```math -ξ = \\cos(α ||η||) η - \\sin(α ||η||) ||η|| A, -``` -such that ``\\langle A', ξ \\rangle = 0``. If the initial state was properly normalized, we -have that -```math -||A'|| = ||A|| = 1. +ξ = \\cos \\left( α \\frac{||η||}{||A||} \\right) η - \\sin \\left( α \\frac{||η||}{||A||} \\right) ||η|| \\frac{A}{||A||}, ``` +such that ``\\langle A', ξ \\rangle = 0`` and ``||A'|| = ||A||``. """ function peps_retract(x, η, α) peps = x[1] + norms_peps = norm.(peps.A) norms_η = norm.(η.A) peps´ = similar(x[1]) - peps´.A .= cos.(α .* norms_η) .* peps.A .+ sin.(α .* norms_η) .* η.A ./ norms_η + peps´.A .= + cos.(α .* norms_η ./ norms_peps) .* peps.A .+ + sin.(α .* norms_η ./ norms_peps) .* norms_peps .* η.A ./ norms_η env = deepcopy(x[2]) ξ = similar(η) - ξ.A .= cos.(α .* norms_η) .* η.A .- sin.(α .* norms_η) .* norms_η .* peps.A + ξ.A .= + cos.(α .* norms_η ./ norms_peps) .* η.A .- + sin.(α .* norms_η ./ norms_peps) .* norms_η .* peps.A ./ norms_peps return (peps´, env), ξ end @@ -209,12 +210,13 @@ from a direction `η` of the form ``` where ``\\langle Δξ, A \\rangle = \\langle Δξ, η \\rangle = 0``, it returns ```math -ξ(α) = \\left\\langle \\frac{η}{||η||}, ξ \\right\\rangle \\left( \\cos(α ||η||) \\frac{η}{||η||} - \\sin(α ||η||) A \\right) + Δξ +ξ(α) = \\left\\langle \\frac{η}{||η||}, ξ \\right \\rangle \\left( \\cos \\left( α \\frac{||η||}{||A||} \\right) \\frac{η}{||η||} - \\sin( \\left( α \\frac{||η||}{||A||} \\right) \\frac{A}{||A||} \\right) + Δξ ``` such that ``||ξ(α)|| = ||ξ||, \\langle A', ξ(α) \\rangle = 0``. """ function peps_transport!(ξ, x, η, α, x´) peps = x[1] + norms_peps = norm.(peps.A) norms_η = norm.(η.A) normalized_η = η.A ./ norms_η @@ -225,8 +227,10 @@ function peps_transport!(ξ, x, η, α, x´) # keep orthogonal component fixed, modify the rest by the proper directional derivative ξ.A .= - overlaps_η_ξ .* - (cos.(α .* norms_η) .* normalized_η .- sin.(α .* norms_η) .* peps.A) .+ Δξ + overlaps_η_ξ .* ( + cos.(α .* norms_η ./ norms_peps) .* normalized_η .- + sin.(α .* norms_η ./ norms_peps) .* peps.A ./ norms_peps + ) .+ Δξ return ξ end From 391e7fbbe988fec49ff7cb5b2ba2c4dbb6d94164 Mon Sep 17 00:00:00 2001 From: Lander Burgelman <39218680+leburgel@users.noreply.github.com> Date: Wed, 12 Mar 2025 21:16:35 +0100 Subject: [PATCH 03/10] Update src/algorithms/optimization/peps_optimization.jl Co-authored-by: Lukas Devos --- src/algorithms/optimization/peps_optimization.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/algorithms/optimization/peps_optimization.jl b/src/algorithms/optimization/peps_optimization.jl index 03e6193e4..a5b5e3d82 100644 --- a/src/algorithms/optimization/peps_optimization.jl +++ b/src/algorithms/optimization/peps_optimization.jl @@ -272,9 +272,7 @@ Normalize the individual tensors in the unit cell of an `InfinitePEPS` such that have unit Euclidean norm. """ function peps_normalize(A::InfinitePEPS) - normalized_tensors = map(unitcell(A)) do t - return t / norm(t) - end + normalized_tensors = normalize.(unitcell(A)) return InfinitePEPS(normalized_tensors) end From 4a0641130d3987b66bd7076fe23eb36d3bb6b157 Mon Sep 17 00:00:00 2001 From: leburgel Date: Wed, 12 Mar 2025 21:46:20 +0100 Subject: [PATCH 04/10] Increase tol for cotangent gauge sensitivity warnings --- src/utility/svd.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/utility/svd.jl b/src/utility/svd.jl index 5193ca0f5..1c86bf500 100644 --- a/src/utility/svd.jl +++ b/src/utility/svd.jl @@ -290,8 +290,8 @@ function ChainRulesCore.rrule( Ũ, S̃, Ṽ⁺, info = tsvd(t, alg; trunc, p) U, S, V⁺ = info.U_full, info.S_full, info.V_full # untruncated SVD decomposition - smallest_sval = minimum(minimum(abs.(diag(b))) for (_, b) in blocks(S̃)) - pullback_tol = max(1e-14, 1e-2 * smallest_sval) + smallest_sval = minimum(((_, b),) -> minimum(diag(b)), blocks(S)) + pullback_tol = max(eps(scalartype(S))^(3 / 4), smallest_sval) function tsvd!_nothing_pullback(ΔUSVi) ΔU, ΔS, ΔV⁺, = unthunk.(ΔUSVi) From b8c62fc4d1c0ebfb26aeb95de900729f198581ad Mon Sep 17 00:00:00 2001 From: leburgel Date: Wed, 12 Mar 2025 21:47:52 +0100 Subject: [PATCH 05/10] Fix typo --- src/utility/svd.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/utility/svd.jl b/src/utility/svd.jl index 1c86bf500..a179641ff 100644 --- a/src/utility/svd.jl +++ b/src/utility/svd.jl @@ -290,8 +290,8 @@ function ChainRulesCore.rrule( Ũ, S̃, Ṽ⁺, info = tsvd(t, alg; trunc, p) U, S, V⁺ = info.U_full, info.S_full, info.V_full # untruncated SVD decomposition - smallest_sval = minimum(((_, b),) -> minimum(diag(b)), blocks(S)) - pullback_tol = max(eps(scalartype(S))^(3 / 4), smallest_sval) + smallest_sval = minimum(((_, b),) -> minimum(diag(b)), blocks(S̃)) + pullback_tol = max(eps(scalartype(S̃))^(3 / 4), smallest_sval) function tsvd!_nothing_pullback(ΔUSVi) ΔU, ΔS, ΔV⁺, = unthunk.(ΔUSVi) From c7e804da0c76af09650613ab48b0b0469b494021 Mon Sep 17 00:00:00 2001 From: leburgel Date: Wed, 12 Mar 2025 22:06:55 +0100 Subject: [PATCH 06/10] Not too high, not too low? --- src/utility/svd.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/utility/svd.jl b/src/utility/svd.jl index a179641ff..db39c3a5c 100644 --- a/src/utility/svd.jl +++ b/src/utility/svd.jl @@ -291,7 +291,11 @@ function ChainRulesCore.rrule( U, S, V⁺ = info.U_full, info.S_full, info.V_full # untruncated SVD decomposition smallest_sval = minimum(((_, b),) -> minimum(diag(b)), blocks(S̃)) - pullback_tol = max(eps(scalartype(S̃))^(3 / 4), smallest_sval) + pullback_tol = clamp( + smallest_sval, + eps(scalartype(S̃))^(3 / 4), + TensorKitCRCExt.default_pullback_gaugetol(S), + ) function tsvd!_nothing_pullback(ΔUSVi) ΔU, ΔS, ΔV⁺, = unthunk.(ΔUSVi) From d761218a1221d58789d38b5611dcd9e17a9a553d Mon Sep 17 00:00:00 2001 From: leburgel Date: Wed, 12 Mar 2025 22:20:23 +0100 Subject: [PATCH 07/10] Definitely too low --- src/utility/svd.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/utility/svd.jl b/src/utility/svd.jl index db39c3a5c..8511465d1 100644 --- a/src/utility/svd.jl +++ b/src/utility/svd.jl @@ -292,9 +292,7 @@ function ChainRulesCore.rrule( smallest_sval = minimum(((_, b),) -> minimum(diag(b)), blocks(S̃)) pullback_tol = clamp( - smallest_sval, - eps(scalartype(S̃))^(3 / 4), - TensorKitCRCExt.default_pullback_gaugetol(S), + smallest_sval, eps(scalartype(S̃))^(3 / 4), eps(scalartype(S̃))^(1 / 2) ) function tsvd!_nothing_pullback(ΔUSVi) From 7c490cdeb4d5f757244366af100df7525497c542 Mon Sep 17 00:00:00 2001 From: leburgel Date: Thu, 13 Mar 2025 07:28:11 +0100 Subject: [PATCH 08/10] Very arbitrary --- src/utility/svd.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utility/svd.jl b/src/utility/svd.jl index 8511465d1..f214d4273 100644 --- a/src/utility/svd.jl +++ b/src/utility/svd.jl @@ -292,7 +292,7 @@ function ChainRulesCore.rrule( smallest_sval = minimum(((_, b),) -> minimum(diag(b)), blocks(S̃)) pullback_tol = clamp( - smallest_sval, eps(scalartype(S̃))^(3 / 4), eps(scalartype(S̃))^(1 / 2) + 10 * smallest_sval, eps(scalartype(S̃))^(3 / 4), eps(scalartype(S̃))^(1 / 2) ) function tsvd!_nothing_pullback(ΔUSVi) From e357dde2f882f37327c5c3788c6ebb2ab9f92bda Mon Sep 17 00:00:00 2001 From: leburgel Date: Thu, 13 Mar 2025 08:02:46 +0100 Subject: [PATCH 09/10] Restore tolerance, reduce bond dimension instead? --- src/utility/svd.jl | 2 +- test/heisenberg.jl | 2 +- test/tf_ising.jl | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/utility/svd.jl b/src/utility/svd.jl index f214d4273..8511465d1 100644 --- a/src/utility/svd.jl +++ b/src/utility/svd.jl @@ -292,7 +292,7 @@ function ChainRulesCore.rrule( smallest_sval = minimum(((_, b),) -> minimum(diag(b)), blocks(S̃)) pullback_tol = clamp( - 10 * smallest_sval, eps(scalartype(S̃))^(3 / 4), eps(scalartype(S̃))^(1 / 2) + smallest_sval, eps(scalartype(S̃))^(3 / 4), eps(scalartype(S̃))^(1 / 2) ) function tsvd!_nothing_pullback(ΔUSVi) diff --git a/test/heisenberg.jl b/test/heisenberg.jl index 01b344d59..fa2fdc571 100644 --- a/test/heisenberg.jl +++ b/test/heisenberg.jl @@ -8,7 +8,7 @@ using OptimKit # initialize parameters Dbond = 2 -χenv = 16 +χenv = 14 gradtol = 1e-3 # compare against Juraj Hasik's data: # https://github.com/jurajHasik/j1j2_ipeps_states/blob/main/single-site_pg-C4v-A1/j20.0/state_1s_A1_j20.0_D2_chi_opt48.dat diff --git a/test/tf_ising.jl b/test/tf_ising.jl index fa3b730ca..ed6a4d84b 100644 --- a/test/tf_ising.jl +++ b/test/tf_ising.jl @@ -18,7 +18,7 @@ mˣ = 0.91 # initialize parameters χbond = 2 -χenv = 16 +χenv = 14 gradtol = 1e-3 # initialize states From 8b37be790663971915d9cc04f298138284a3ae53 Mon Sep 17 00:00:00 2001 From: leburgel Date: Thu, 13 Mar 2025 10:20:30 +0100 Subject: [PATCH 10/10] No real difference, so restore bond dimension --- test/heisenberg.jl | 2 +- test/tf_ising.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/test/heisenberg.jl b/test/heisenberg.jl index fa2fdc571..01b344d59 100644 --- a/test/heisenberg.jl +++ b/test/heisenberg.jl @@ -8,7 +8,7 @@ using OptimKit # initialize parameters Dbond = 2 -χenv = 14 +χenv = 16 gradtol = 1e-3 # compare against Juraj Hasik's data: # https://github.com/jurajHasik/j1j2_ipeps_states/blob/main/single-site_pg-C4v-A1/j20.0/state_1s_A1_j20.0_D2_chi_opt48.dat diff --git a/test/tf_ising.jl b/test/tf_ising.jl index ed6a4d84b..fa3b730ca 100644 --- a/test/tf_ising.jl +++ b/test/tf_ising.jl @@ -18,7 +18,7 @@ mˣ = 0.91 # initialize parameters χbond = 2 -χenv = 14 +χenv = 16 gradtol = 1e-3 # initialize states