From 8bd563db10021ab9df67f564918d25c750f612ba Mon Sep 17 00:00:00 2001 From: Atsushi Ueda <101638759+dartsushi@users.noreply.github.com> Date: Wed, 19 Aug 2026 12:06:30 +0200 Subject: [PATCH 1/4] Faster optimization for SloopTNR Modify optimize_S to avoid calculating a constant over and over again --- src/schemes/symmetric_looptnr.jl | 18 +++++++++++++++++- 1 file changed, 17 insertions(+), 1 deletion(-) diff --git a/src/schemes/symmetric_looptnr.jl b/src/schemes/symmetric_looptnr.jl index 2481cbae..32b577d9 100644 --- a/src/schemes/symmetric_looptnr.jl +++ b/src/schemes/symmetric_looptnr.jl @@ -74,6 +74,21 @@ function TTtoNorm(TT) return tr(T8 * b) end +function TtoNorm(T) + @tensor TT[-1 -2; -3 -4] := T[1 -1 2 -3] * conj(T[1 -2 2 -4]) + return TTtoNorm(TT) +end + +function cost_looptnr(S, T, n_TT) + @assert eltype(S) == Float64 "Modification is needed for complex numbers!" + SS = StoSS(S) + + @tensor TSS[-1 -2; -3 -4] := T[1 -1 2 -3] * conj(SS[1 -2 2 -4]) + @tensor S4[-1 -2; -3 -4] := SS[1 -1 2 -3] * conj(SS[1 -2 2 -4]) + # T + return n_TT + TTtoNorm(S4) - 2 * TTtoNorm(TSS) +end + function cost_looptnr(S, T) @assert eltype(S) == Float64 "Modification is needed for complex numbers!" SS = StoSS(S) @@ -92,7 +107,8 @@ function fg(f, A) end function optimize_S(scheme, S) - opt_fun(x) = cost_looptnr(x, scheme.T) + n_TT = TtoNorm(scheme.T) + opt_fun(x) = cost_looptnr(x, scheme.T, n_TT) opt_fg(x) = fg(opt_fun, x) Sopt, fx, gx, numfg, normgradhistory = optimize( opt_fg, S, From 2d6a4f42884342f53aa3ba68ffdca290ec1f447e Mon Sep 17 00:00:00 2001 From: Atsushi Ueda <101638759+dartsushi@users.noreply.github.com> Date: Wed, 19 Aug 2026 12:08:32 +0200 Subject: [PATCH 2/4] Update gradient tolerance and max iterations in SLoopTNR --- src/schemes/symmetric_looptnr.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/schemes/symmetric_looptnr.jl b/src/schemes/symmetric_looptnr.jl index 32b577d9..6d201fba 100644 --- a/src/schemes/symmetric_looptnr.jl +++ b/src/schemes/symmetric_looptnr.jl @@ -27,7 +27,7 @@ mutable struct SLoopTNR{E, S, TT <: AbstractTensorMap{E, S, 4, 0}} <: TNRScheme{ "Gradient optimization algorithm" gradalg::OptimKit.LBFGS - function SLoopTNR(T::TT; gradalg = LBFGS(10; verbosity = 0, gradtol = 6.0e-7, maxiter = 40000)) where {E, S, TT <: AbstractTensorMap{E, S, 4, 0}} + function SLoopTNR(T::TT; gradalg = LBFGS(10; verbosity = 0, gradtol = 1.0e-6, maxiter = 200)) where {E, S, TT <: AbstractTensorMap{E, S, 4, 0}} return new{E, S, TT}(T, gradalg) end end From f24c0a729e0d90893011571b1ac4e5cbeccb872a Mon Sep 17 00:00:00 2001 From: Atsushi Ueda <101638759+dartsushi@users.noreply.github.com> Date: Wed, 19 Aug 2026 14:59:35 +0200 Subject: [PATCH 3/4] Got rid of Zygote I manually put the gradient and removed AD. It is now two times faster --- src/schemes/symmetric_looptnr.jl | 79 +++++++++++++++++++++----------- 1 file changed, 52 insertions(+), 27 deletions(-) diff --git a/src/schemes/symmetric_looptnr.jl b/src/schemes/symmetric_looptnr.jl index 6d201fba..499fba0a 100644 --- a/src/schemes/symmetric_looptnr.jl +++ b/src/schemes/symmetric_looptnr.jl @@ -48,7 +48,7 @@ classical_ising_inv() = classical_ising_inv(ising_βc) ########## utility functions ########## function trnorm_2x2(T) - @tensor TT[-1 -2; -3 -4] := T[1 -1 2 -3] * conj(T[1 -2 2 -4]) + @tensoropt TT[-1 -2; -3 -4] := T[1 -1 2 -3] * conj(T[1 -2 2 -4]) return sqrt(TTtoNorm(TT)) end @@ -56,7 +56,7 @@ end function StoSS(S) V = domain(S)[1] b = isomorphism(V, V') - @tensor SS[-1 -2 -3 -4] := S[-1 -2; 1] * S[-3 -4; 2] * b[1 2] + @tensoropt SS[-1 -2 -3 -4] := S[-1 -2; 1] * S[-3 -4; 2] * b[1 2] return SS end @@ -64,52 +64,77 @@ function TTtoNorm(TT) V = domain(TT) b = isomorphism(V[1] ⊗ V[2], V[1]' ⊗ V[2]') TTb = TT * b - @tensor T4[-1 -2; -3 -4] := TT[-1 -2; 1 2] * TTb[-3 -4; 1 2] + @tensoropt T4[-1 -2; -3 -4] := TT[-1 -2; 1 2] * TTb[-3 -4; 1 2] V = domain(T4) b = isomorphism(V[1] ⊗ V[2], V[1]' ⊗ V[2]') T4b = T4 * b - @tensor T8[-1 -2; -3 -4] := T4[-1 -2; 1 2] * T4b[-3 -4; 1 2] + @tensoropt T8[-1 -2; -3 -4] := T4[-1 -2; 1 2] * T4b[-3 -4; 1 2] V = domain(T8) b = isomorphism(V[1] ⊗ V[2], V[1]' ⊗ V[2]') return tr(T8 * b) end function TtoNorm(T) - @tensor TT[-1 -2; -3 -4] := T[1 -1 2 -3] * conj(T[1 -2 2 -4]) + @tensoropt TT[-1 -2; -3 -4] := T[1 -1 2 -3] * conj(T[1 -2 2 -4]) return TTtoNorm(TT) end - + function cost_looptnr(S, T, n_TT) @assert eltype(S) == Float64 "Modification is needed for complex numbers!" SS = StoSS(S) - @tensor TSS[-1 -2; -3 -4] := T[1 -1 2 -3] * conj(SS[1 -2 2 -4]) - @tensor S4[-1 -2; -3 -4] := SS[1 -1 2 -3] * conj(SS[1 -2 2 -4]) - # T + @tensoropt TSS[-1 -2; -3 -4] := T[1 -1 2 -3] * conj(SS[1 -2 2 -4]) + @tensoropt S4[-1 -2; -3 -4] := SS[1 -1 2 -3] * conj(SS[1 -2 2 -4]) return n_TT + TTtoNorm(S4) - 2 * TTtoNorm(TSS) end -function cost_looptnr(S, T) - @assert eltype(S) == Float64 "Modification is needed for complex numbers!" - SS = StoSS(S) - @tensor TT[-1 -2; -3 -4] := T[1 -1 2 -3] * conj(T[1 -2 2 -4]) - @tensor TSS[-1 -2; -3 -4] := T[1 -1 2 -3] * conj(SS[1 -2 2 -4]) - @tensor S4[-1 -2; -3 -4] := SS[1 -1 2 -3] * conj(SS[1 -2 2 -4]) - # T - return TTtoNorm(TT) + TTtoNorm(S4) - 2 * TTtoNorm(TSS) +########## Gradient Optimization ########## +function loop_environment(TT) + V = domain(TT) + b = isomorphism(V[1] ⊗ V[2], V[1]' ⊗ V[2]') + TTb = TT * b + @tensoropt T4[-1 -2; -3 -4] := TT[-1 -2; 1 2] * TTb[-3 -4; 1 2] + return T4 * b * TT end -########## Gradient Optimization ########## -function fg(f, A) - f_out, g = Zygote.withgradient(f, A) +function StoSS_pullback(S, gSS) + V = domain(S)[1] + b = isomorphism(V, V') + Sb = S * b + @tensoropt gS[-1 -2; -3] := gSS[-1 -2 1 2] * conj(Sb[1 2; -3]) + return 2 * gS +end + +function gradient_looptnr(S, T, SS, TSS, S4) + # One of the eight equal contributions from the S-S overlap. + env_SS = loop_environment(S4) + @tensoropt gSS[-1 -2 -3 -4] := SS[-1 1 -3 2] * env_SS[-2 1; -4 2] + gS_SS = StoSS_pullback(S, gSS) + + # One of the four equal contributions from the T-S overlap. + env_TS = loop_environment(TSS) + env_TS_transposed = permute(env_TS, ((2, 1), (4, 3))) + @tensoropt gTS_dual[-1 -2 -3 -4] := conj(T[-1 1 -3 2]) * env_TS_transposed[-2 1; -4 2] + gTS = flip(gTS_dual, (1, 2, 3, 4)) + gS_TS = StoSS_pullback(S, gTS) + + return 8 * gS_SS - 2 * 4 * gS_TS +end + +function cost_looptnr_fg(S, T, n_TT) + @assert eltype(S) == Float64 "Modification is needed for complex numbers!" + SS = StoSS(S) - return f_out, g[1] + @tensoropt TSS[-1 -2; -3 -4] := T[1 -1 2 -3] * conj(SS[1 -2 2 -4]) + @tensoropt S4[-1 -2; -3 -4] := SS[1 -1 2 -3] * conj(SS[1 -2 2 -4]) + cost = n_TT + TTtoNorm(S4) - 2 * TTtoNorm(TSS) + grad = gradient_looptnr(S, T, SS, TSS, S4) + return cost, grad end function optimize_S(scheme, S) n_TT = TtoNorm(scheme.T) - opt_fun(x) = cost_looptnr(x, scheme.T, n_TT) - opt_fg(x) = fg(opt_fun, x) + opt_fg(x) = cost_looptnr_fg(x, scheme.T, n_TT) Sopt, fx, gx, numfg, normgradhistory = optimize( opt_fg, S, scheme.gradalg @@ -142,7 +167,7 @@ end function entanglement_filtering(T; trunc = trunctol(atol = 1.0e-12)) entanglement_function(steps, data) = abs(data[end]) - entanglement_criterion = maxiter(100) & convcrit(1.0e-12, entanglement_function) + entanglement_criterion = maxiter(200) & convcrit(1.0e-12, entanglement_function) psi_center = Ψ_center(T) psi_corner = Ψ_corner(T) @@ -162,7 +187,7 @@ function entanglement_filtering(T; trunc = trunctol(atol = 1.0e-12)) P_top = PL_list[3] P_left = PL_list[3] - @tensor T_new[-1 -2 -3 -4] := T[1 2 3 4] * P_left[-1; 1] * P_bottom[-2; 2] * + @tensoropt T_new[-1 -2 -3 -4] := T[1 2 3 4] * P_left[-1; 1] * P_bottom[-2; 2] * P_top[-3; 3] * P_right[-4; 4] return T_new end @@ -184,14 +209,14 @@ function ef_oneloop(T, trunc::TruncationStrategy) criterion, trunc ) i = 1 - @tensor S[-2 -1; -3] := ΨB[i][-1; -2 2] * PRs[mod(i, 8) + 1][2; -3] + @tensoropt S[-2 -1; -3] := ΨB[i][-1; -2 2] * PRs[mod(i, 8) + 1][2; -3] return S end ########## Updating the tensor ########## function combine_4S(S) Sflip = flip(S, (1, 2)) - @tensor Tnew[-1 -2 -3 -4] := S[1 2; -4] * Sflip[1 4; -3] * S[3 4; -1] * Sflip[3 2; -2] + @tensoropt Tnew[-1 -2 -3 -4] := S[1 2; -4] * Sflip[1 4; -3] * S[3 4; -1] * Sflip[3 2; -2] return Tnew end From ffd49f4e1ebb1539709717c6dd911ef8c28df19c Mon Sep 17 00:00:00 2001 From: dartsushi Date: Wed, 19 Aug 2026 19:20:09 +0200 Subject: [PATCH 4/4] fptensor --- src/TNRKit.jl | 5 +- src/utility/fixed_point_tensor.jl | 322 ++++++++++++++++++++++++++++++ test/schemes/schemes.jl | 79 ++++++++ 3 files changed, 405 insertions(+), 1 deletion(-) create mode 100644 src/utility/fixed_point_tensor.jl diff --git a/src/TNRKit.jl b/src/TNRKit.jl index 78f47824..8b7cadd2 100644 --- a/src/TNRKit.jl +++ b/src/TNRKit.jl @@ -5,7 +5,7 @@ using MatrixAlgebraKit using MatrixAlgebraKit: TruncationStrategy using LoggingExtras, Printf using KrylovKit -using OptimKit, Zygote +using OptimKit using DocStringExtensions using SpecialFunctions using FastGaussQuadrature @@ -138,6 +138,9 @@ include("utility/transfer_matrix.jl") include("utility/cft.jl") export CFTData, extract_tau_and_c +include("utility/fixed_point_tensor.jl") +export fixed_point_tensor, fixed_point_tensor_4x4 + include("utility/gs_degeneracy.jl") export ground_state_degeneracy, gu_wen_ratio diff --git a/src/utility/fixed_point_tensor.jl b/src/utility/fixed_point_tensor.jl new file mode 100644 index 00000000..3a735707 --- /dev/null +++ b/src/utility/fixed_point_tensor.jl @@ -0,0 +1,322 @@ +""" + fixed_point_tensor(T; nstates = 3, eig_tol = 1.0e-12, eig_krylovdim = 40, + return_basis = false) + fixed_point_tensor(scheme::SLoopTNR; kwargs...) + +Compute the normalized fixed-point tensor elements of a four-leg tensor `T` in +the transfer-matrix eigenbasis. The construction follows Eq. (2) of +[Ueda and Yamazaki (2023)](https://arxiv.org/abs/2307.02523): the four-tensor +`2 × 2` transfer matrices in the horizontal and vertical directions are +diagonalized, their leading `nstates` eigenvectors are used as boundary +projectors for a mirrored `2 × 2` tensor patch, and the result is normalized by +its `(1, 1, 1, 1)` element. The returned leg order is left-bottom-top-right, +matching the TNRKit convention. + +For the critical Ising model, the three leading states are ordered as +`(1, σ, ε)`, so `fixed_point_tensor(scheme)[2, 2, 2, 2]` is the `σσσσ` +element. + +Set `return_basis = true` to return a named tuple containing `elements`, the +horizontal and vertical bases, and their corresponding transfer-matrix +eigenvalues. +""" +function fixed_point_tensor( + T::AbstractTensorMap{E, S, 4, 0}; nstates::Int = 3, + eig_tol::Real = 1.0e-12, eig_krylovdim::Int = 40, + return_basis::Bool = false + ) where {E, S} + nstates > 0 || throw(ArgumentError("nstates must be positive")) + eig_tol > 0 || throw(ArgumentError("eig_tol must be positive")) + eig_krylovdim > nstates || throw(ArgumentError("eig_krylovdim must exceed nstates")) + + A = convert(Array, T) + allequal(size(A)) || throw(DimensionMismatch("all four tensor legs must have equal dimension")) + nstates <= size(A, 1)^2 || throw(DimensionMismatch( + "cannot retain $nstates states from a transfer matrix of dimension $(size(A, 1)^2)" + )) + + horizontal_basis, horizontal_eigenvalues = _fixed_point_basis( + T, true, nstates, eig_tol, eig_krylovdim + ) + vertical_basis, vertical_eigenvalues = _fixed_point_basis( + T, false, nstates, eig_tol, eig_krylovdim + ) + horizontal_projector = conj.(horizontal_basis) + vertical_projector = conj.(vertical_basis) + Aflip = convert(Array, flip(T, (1, 2, 3, 4))) + + # TNRKit orders the legs as left-bottom-top-right. The four tensors are + # related by mirror symmetry and arranged as + # + # top + # u -------- v + # / \ + # left a--T--------Tf--b right + # | | + # c--Tf-------T--d + # \ / + # w -------- r + # bottom + # + # and each pair of boundary indices is projected onto the corresponding + # transfer-matrix eigenbasis. Keeping this as one optimized contraction + # avoids materializing the eight-index boundary tensor. + @tensoropt elements[left, bottom, top, right] := + A[a, x, u, y] * Aflip[b, z, v, y] * Aflip[c, x, w, q] * A[d, z, r, q] * + horizontal_projector[a, c, left] * vertical_projector[w, r, bottom] * + vertical_projector[u, v, top] * horizontal_projector[b, d, right] + + normalization = elements[1, 1, 1, 1] + iszero(normalization) && throw(ArgumentError("the fixed-point identity element is zero")) + elements ./= normalization + + if return_basis + return (; + elements, horizontal_basis, vertical_basis, + horizontal_eigenvalues, vertical_eigenvalues, + ) + end + return elements +end + +fixed_point_tensor(scheme::SLoopTNR; kwargs...) = fixed_point_tensor(scheme.T; kwargs...) + +""" + fixed_point_tensor_4x4(T; nstates = 3, eig_tol = 1.0e-12, + eig_krylovdim = 40, return_basis = false) + fixed_point_tensor_4x4(scheme::SLoopTNR; kwargs...) + +Compute fixed-point tensor elements from a mirrored `4 × 4` patch. Each CFT +basis state is an eigenvector on four boundary bonds. The horizontal and +vertical transfer matrices are applied matrix-free as four successive column +or row tensor-network contractions, so a dense matrix of size `D^4 × D^4` is +never constructed. + +The result has left-bottom-top-right leg order and is normalized by its +`(1, 1, 1, 1)` element. With `return_basis = true`, the return value has the +same named-tuple layout as [`fixed_point_tensor`](@ref), but each basis has +shape `D × D × D × D × nstates`. +""" +function fixed_point_tensor_4x4( + T::AbstractTensorMap{E, S, 4, 0}; nstates::Int = 3, + eig_tol::Real = 1.0e-12, eig_krylovdim::Int = 40, + return_basis::Bool = false + ) where {E, S} + nstates > 0 || throw(ArgumentError("nstates must be positive")) + eig_tol > 0 || throw(ArgumentError("eig_tol must be positive")) + eig_krylovdim > nstates || throw(ArgumentError("eig_krylovdim must exceed nstates")) + + A = convert(Array, T) + allequal(size(A)) || throw(DimensionMismatch("all four tensor legs must have equal dimension")) + d = size(A, 1) + nstates <= d^4 || throw(DimensionMismatch( + "cannot retain $nstates states from a transfer matrix of dimension $(d^4)" + )) + + horizontal_basis, horizontal_eigenvalues = _fixed_point_basis_4x4( + T, true, nstates, eig_tol, eig_krylovdim + ) + vertical_basis, vertical_eigenvalues = _fixed_point_basis_4x4( + T, false, nstates, eig_tol, eig_krylovdim + ) + elements = _project_fixed_point_patch_4x4( + A, convert(Array, flip(T, (1, 2, 3, 4))), + conj.(horizontal_basis), conj.(vertical_basis) + ) + + normalization = elements[1, 1, 1, 1] + iszero(normalization) && throw(ArgumentError("the fixed-point identity element is zero")) + elements ./= normalization + + if return_basis + return (; + elements, horizontal_basis, vertical_basis, + horizontal_eigenvalues, vertical_eigenvalues, + ) + end + return elements +end + +fixed_point_tensor_4x4(scheme::SLoopTNR; kwargs...) = + fixed_point_tensor_4x4(scheme.T; kwargs...) + +function _fixed_point_transfer_matrix( + T::AbstractTensorMap{E, S, 4, 0}, horizontal::Bool + ) where {E, S} + Tflip = flip(T, (1, 2, 3, 4)) + if horizontal + @tensoropt transfer[-1 -2; -3 -4] := + T[-1 1; 3 2] * Tflip[-3 4; 5 2] * + Tflip[-2 1; 3 6] * T[-4 4; 5 6] + else + @tensoropt transfer[-1 -2; -3 -4] := + T[1 3; -1 2] * Tflip[1 4; -2 2] * + Tflip[5 3; -3 6] * T[5 4; -4 6] + end + dense_transfer = convert(Array, transfer) + return reshape( + dense_transfer, size(dense_transfer, 1) * size(dense_transfer, 2), : + ) +end + +function _fixed_point_basis( + T::AbstractTensorMap{E, S, 4, 0}, horizontal::Bool, nstates::Int, + eig_tol::Real, eig_krylovdim::Int + ) where {E, S} + transfer = _fixed_point_transfer_matrix(T, horizontal) + hermitian_transfer = Hermitian((transfer + transfer') / 2) + x0 = convert.(eltype(transfer), sin.(eachindex(axes(transfer, 1)))) + eigenvalues, eigenvectors, info = eigsolve( + hermitian_transfer, x0, nstates, :LM; + krylovdim = min(size(transfer, 1), eig_krylovdim), maxiter = 300, + tol = eig_tol, verbosity = 0 + ) + info.converged < nstates && @warn "Fixed-point transfer-matrix eigensolver did not converge" horizontal info + + order = sortperm(real.(eigenvalues); rev = true)[1:nstates] + eigenvalues = eigenvalues[order] + eigenvectors = reduce(hcat, eigenvectors[order]) + + # Fix the otherwise arbitrary phase of every state. This makes tensor + # elements with an odd number of a given state reproducible as well. + for state in axes(eigenvectors, 2) + vector = @view eigenvectors[:, state] + pivot = vector[argmax(abs.(vector))] + iszero(pivot) || (vector .*= conj(pivot) / abs(pivot)) + end + + d = dim(codomain(T)[1]) + return reshape(eigenvectors, d, d, nstates), eigenvalues +end + +function _fixed_point_transfer_action_4x4( + T::AbstractTensorMap{E, S, 4, 0}, horizontal::Bool + ) where {E, S} + A = convert(Array, T) + Aflip = convert(Array, flip(T, (1, 2, 3, 4))) + d = size(A, 1) + + if horizontal + return function (vector) + boundary = reshape(vector, d, d, d, d) + for column in 4:-1:1 + boundary = _apply_fixed_point_column_4x4(boundary, A, Aflip, column) + end + return vec(boundary) + end + end + return function (vector) + boundary = reshape(vector, d, d, d, d) + for row in 4:-1:1 + boundary = _apply_fixed_point_row_4x4(boundary, A, Aflip, row) + end + return vec(boundary) + end +end + +function _apply_fixed_point_row_4x4(boundary, A, Aflip, row::Int) + tensors = [iseven(row + column) ? A : Aflip for column in 1:4] + if isodd(row) + indices = [ + [4, 5, -1, 1], [2, 6, -2, 1], + [2, 7, -3, 3], [4, 8, -4, 3], [5, 6, 7, 8], + ] + else + indices = [ + [4, -1, 5, 1], [2, -2, 6, 1], + [2, -3, 7, 3], [4, -4, 8, 3], [5, 6, 7, 8], + ] + end + return ncon([tensors..., boundary], indices) +end + +function _apply_fixed_point_column_4x4(boundary, A, Aflip, column::Int) + tensors = [iseven(row + column) ? A : Aflip for row in 1:4] + if iseven(column) + indices = [ + [5, 1, 4, -1], [6, 1, 2, -2], + [7, 3, 2, -3], [8, 3, 4, -4], [5, 6, 7, 8], + ] + else + indices = [ + [-1, 1, 4, 5], [-2, 1, 2, 6], + [-3, 3, 2, 7], [-4, 3, 4, 8], [5, 6, 7, 8], + ] + end + return ncon([tensors..., boundary], indices) +end + +function _fixed_point_basis_4x4( + T::AbstractTensorMap{E, S, 4, 0}, horizontal::Bool, nstates::Int, + eig_tol::Real, eig_krylovdim::Int + ) where {E, S} + d = dim(codomain(T)[1]) + transfer = _fixed_point_transfer_action_4x4(T, horizontal) + x0 = convert.(E, sin.(1:(d^4))) + eigenvalues, eigenvectors, info = eigsolve( + transfer, x0, nstates, :LM; + krylovdim = min(d^4, eig_krylovdim), maxiter = 300, + tol = eig_tol, verbosity = 0 + ) + info.converged < nstates && @warn "4 × 4 fixed-point transfer-matrix eigensolver did not converge" horizontal info + + order = sortperm(real.(eigenvalues); rev = true)[1:nstates] + eigenvalues = eigenvalues[order] + eigenvectors = reduce(hcat, eigenvectors[order]) + for state in axes(eigenvectors, 2) + vector = @view eigenvectors[:, state] + pivot = vector[argmax(abs.(vector))] + iszero(pivot) || (vector .*= conj(pivot) / abs(pivot)) + end + return reshape(eigenvectors, d, d, d, d, nstates), eigenvalues +end + +function _project_fixed_point_patch_4x4(A, Aflip, horizontal_projector, vertical_projector) + site_indices = [zeros(Int, 4) for _ in 1:4, _ in 1:4] + label = 1 + + for row in 1:4, column in 1:3 + leg = isodd(column) ? 4 : 1 + site_indices[row, column][leg] = label + site_indices[row, column + 1][leg] = label + label += 1 + end + for row in 1:3, column in 1:4 + leg = isodd(row) ? 2 : 3 + site_indices[row, column][leg] = label + site_indices[row + 1, column][leg] = label + label += 1 + end + + left = Int[] + right = Int[] + for row in 1:4 + push!(left, label) + site_indices[row, 1][1] = label + label += 1 + push!(right, label) + site_indices[row, 4][1] = label + label += 1 + end + bottom = Int[] + top = Int[] + for column in 1:4 + push!(bottom, label) + site_indices[4, column][3] = label + label += 1 + push!(top, label) + site_indices[1, column][3] = label + label += 1 + end + + tensors = Any[ + iseven(row + column) ? A : Aflip for row in 1:4 for column in 1:4 + ] + indices = [site_indices[row, column] for row in 1:4 for column in 1:4] + append!( + tensors, + (horizontal_projector, vertical_projector, vertical_projector, horizontal_projector) + ) + append!(indices, ([left; -1], [bottom; -2], [top; -3], [right; -4])) + return ncon(tensors, indices) +end diff --git a/test/schemes/schemes.jl b/test/schemes/schemes.jl index b428a0a6..7cdc161d 100644 --- a/test/schemes/schemes.jl +++ b/test/schemes/schemes.jl @@ -1,6 +1,8 @@ using Test using TNRKit using TensorKit +using LinearAlgebra +using OptimKit # This tests every scheme in the library on the Z2 symmetric Ising model. println("---------------------") @@ -369,6 +371,71 @@ end end # SLoopTNR +@testset "Fixed-point tensor basis" begin + T_inv = classical_ising_inv() + Tflip = flip(T_inv, (1, 2, 3, 4)) + result = fixed_point_tensor(T_inv; return_basis = true) + horizontal_transfer = TNRKit._fixed_point_transfer_matrix( + T_inv, true + ) + horizontal_basis = reshape( + result.horizontal_basis, size(horizontal_transfer, 1), : + ) + vertical_transfer = TNRKit._fixed_point_transfer_matrix( + T_inv, false + ) + vertical_basis = reshape(result.vertical_basis, size(vertical_transfer, 1), :) + @tensoropt vertical_tensor[-1 -2; -3 -4] := + T_inv[1 3; -1 2] * Tflip[1 4; -2 2] * + Tflip[5 3; -3 6] * T_inv[5 4; -4 6] + expected_vertical_transfer = reshape(convert(Array, vertical_tensor), size(vertical_transfer)) + + @test result.elements[1, 1, 1, 1] ≈ 1 + @test vertical_transfer ≈ expected_vertical_transfer + @test horizontal_transfer * horizontal_basis ≈ + horizontal_basis * Diagonal(result.horizontal_eigenvalues) + @test vertical_transfer * vertical_basis ≈ + vertical_basis * Diagonal(result.vertical_eigenvalues) + @test result.horizontal_eigenvalues ≈ result.vertical_eigenvalues +end + +@testset "4 × 4 fixed-point tensor basis" begin + T_inv = classical_ising_inv() + result = fixed_point_tensor_4x4(T_inv; return_basis = true, eig_tol = 1.0e-10) + horizontal_transfer = TNRKit._fixed_point_transfer_action_4x4(T_inv, true) + vertical_transfer = TNRKit._fixed_point_transfer_action_4x4(T_inv, false) + horizontal_basis = reshape(result.horizontal_basis, :, 3) + vertical_basis = reshape(result.vertical_basis, :, 3) + + @test result.elements[1, 1, 1, 1] ≈ 1 + @test real(result.elements[2, 2, 2, 2]) ≈ 0.3964205 atol = 1.0e-6 + for state in 1:3 + @test horizontal_transfer(horizontal_basis[:, state]) ≈ + result.horizontal_eigenvalues[state] * horizontal_basis[:, state] + @test vertical_transfer(vertical_basis[:, state]) ≈ + result.vertical_eigenvalues[state] * vertical_basis[:, state] + end + @test result.horizontal_eigenvalues ≈ result.vertical_eigenvalues +end + +@testset "SLoopTNR - Manual gradient" begin + V = ℝ^2 + T_inv = ones(Float64, V ⊗ V ⊗ V ⊗ V ← one(V)) + S = ones(Float64, V ⊗ V ← V) + dS = ones(Float64, space(S)) + n_TT = TNRKit.TtoNorm(T_inv) + + cost, grad = TNRKit.cost_looptnr_fg(S, T_inv, n_TT) + ϵ = 1.0e-6 + directional_derivative = ( + TNRKit.cost_looptnr(S + ϵ * dS, T_inv, n_TT) - + TNRKit.cost_looptnr(S - ϵ * dS, T_inv, n_TT) + ) / (2ϵ) + + @test cost ≈ TNRKit.cost_looptnr(S, T_inv, n_TT) + @test real(dot(grad, dS)) ≈ directional_derivative rtol = 1.0e-9 +end + @testset "SLoopTNR - Ising Model" begin @info "SLoopTNR ising free energy" T_inv = classical_ising_inv() @@ -377,6 +444,18 @@ end data = run!(scheme, truncrank(4), maxiter(25)) @test free_energy(data, ising_βc) ≈ f_onsager rtol = 1.0e-5 + + @info "SLoopTNR Ising fixed-point tensor" + gradalg = LBFGS(10; verbosity = 0, gradtol = 6.0e-7, maxiter = 2000) + scheme = SLoopTNR(classical_ising_inv(); gradalg) + run!(scheme, truncrank(16), maxiter(16)) + fp = fixed_point_tensor(scheme) + + # The finite-χ value approaches the exact 0.645 from below (the paper + # reports 0.610 at D = 96 and finite L). + σ4 = real(fp[2, 2, 2, 2]) + @test σ4 ≈ 0.5967 atol = 5.0e-3 + @test σ4 ≈ 0.645 atol = 6.0e-2 end # ctm