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/7] 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/7] 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/7] 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/7] 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 From 7ed42758a65a2a5700bf38169d9e4b75653d82cc Mon Sep 17 00:00:00 2001 From: Atsushi Ueda <101638759+dartsushi@users.noreply.github.com> Date: Thu, 20 Aug 2026 12:09:37 +0200 Subject: [PATCH 5/7] Remove fixed_point_tensor inclusion and export Removed fixed_point_tensor and its export from TNRKit. --- src/TNRKit.jl | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/TNRKit.jl b/src/TNRKit.jl index 8b7cadd2..1bccfd87 100644 --- a/src/TNRKit.jl +++ b/src/TNRKit.jl @@ -138,9 +138,6 @@ 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 From dde4720f41e3b408fda97ab638beb5b734774386 Mon Sep 17 00:00:00 2001 From: Atsushi Ueda <101638759+dartsushi@users.noreply.github.com> Date: Thu, 20 Aug 2026 12:11:26 +0200 Subject: [PATCH 6/7] Remove fixed-point tensor basis tests Removed tests for fixed-point tensor basis and 4x4 fixed-point tensor basis. --- test/schemes/schemes.jl | 48 ----------------------------------------- 1 file changed, 48 deletions(-) diff --git a/test/schemes/schemes.jl b/test/schemes/schemes.jl index 7cdc161d..130b7ae4 100644 --- a/test/schemes/schemes.jl +++ b/test/schemes/schemes.jl @@ -370,54 +370,6 @@ end @test free_energy(data, ising_βc; initial_size = 2) ≈ f_onsager rtol = 1.0e-6 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)) From e4782531ce9f441f376b9b13019cac2671180155 Mon Sep 17 00:00:00 2001 From: Atsushi Ueda <101638759+dartsushi@users.noreply.github.com> Date: Thu, 20 Aug 2026 12:12:01 +0200 Subject: [PATCH 7/7] Delete src/utility/fixed_point_tensor.jl --- src/utility/fixed_point_tensor.jl | 322 ------------------------------ 1 file changed, 322 deletions(-) delete mode 100644 src/utility/fixed_point_tensor.jl diff --git a/src/utility/fixed_point_tensor.jl b/src/utility/fixed_point_tensor.jl deleted file mode 100644 index 3a735707..00000000 --- a/src/utility/fixed_point_tensor.jl +++ /dev/null @@ -1,322 +0,0 @@ -""" - 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