From c1478ed3d3e7d350a38e68d8c8c86814dc954b67 Mon Sep 17 00:00:00 2001 From: afeuerpfeil Date: Mon, 6 Oct 2025 19:31:19 -0400 Subject: [PATCH 01/19] precompute contraction GL*A and A*GR in JordanMPO_AC and _AC2 --- .../derivatives/hamiltonian_derivatives.jl | 33 ++++++++++++------- 1 file changed, 21 insertions(+), 12 deletions(-) diff --git a/src/algorithms/derivatives/hamiltonian_derivatives.jl b/src/algorithms/derivatives/hamiltonian_derivatives.jl index b922205f8..edf8d9bcb 100644 --- a/src/algorithms/derivatives/hamiltonian_derivatives.jl +++ b/src/algorithms/derivatives/hamiltonian_derivatives.jl @@ -138,9 +138,12 @@ function AC_hamiltonian( finished = removeunit(GL[end], 2) end - # continuing - A = W.A - continuing = (GL[2:(end - 1)], A, GR[2:(end - 1)]) + if nonzero_length(W.A) > 0 + @tensor GLW[-1 -2 -3; -4 -5] ≔ GL[2:(end - 1)][-1 1; -4] * W.A[1 -2; -5 -3] + continuing = (GLW, permute(GR[2:(end - 1)], (2, 1), (3,))) + else + continuing = missing + end return JordanMPO_AC_Hamiltonian( onsite, not_started, finished, starting, ending, continuing @@ -287,7 +290,15 @@ function AC2_hamiltonian( # continuing - continuing # TODO: MPODerivativeOperator code reuse + optimization - AA = (GL[2:(end - 1)], W1.A, W2.A, GR[2:(end - 1)]) + # TODO: One should only calculate these operators if necessary. One can do this by checking nonzero_length(A1) > 0 && nonzero_length(A2) > 0 and then setting AA=missing or as we discussed via bit matrix multiplication. + ## TODO: Think about how one could and whether one should store these objects and use them for (a) advancing environments in iDMRG, (b) reuse ind backwards-sweep in IDMRG, (c) subspace expansion + if nonzero_length(W1.A) > 0 && nonzero_length(W2.A) > 0 + @tensor GLW[-1 -2 -3; -4 -5] ≔ GL[2:(end -1)][-1 1; -4] * W1.A[1 -2; -5 -3] + @tensor GWR[-1 -2 -3; -4 -5] ≔ W2.A[-1 -5; -3 1] * GR[2:(end -1)][-2 1; -4] + AA = (GLW, GWR) + else + AA = missing + end return JordanMPO_AC2_Hamiltonian(II, IC, ID, CB, CA, AB, AA, BE, DE, EE) end @@ -303,9 +314,9 @@ function (H::JordanMPO_AC_Hamiltonian)(x::MPSTensor) ismissing(H.C) || @plansor y[-1 -2; -3] += x[-1 2; 1] * H.C[-2 -3; 2 1] ismissing(H.B) || @plansor y[-1 -2; -3] += H.B[-1 -2; 1 2] * x[1 2; -3] - GL, A, GR = H.A - if nonzero_length(A) > 0 - @plansor y[-1 -2; -3] += GL[-1 5; 4] * x[4 2; 1] * A[5 -2; 2 3] * GR[1 3; -3] + if !ismissing(H.A) + GLW, GR = H.A + @tensor y[-1 -2; -3] += GLW[-1 -2 3; 1 2] * x[1 2; 4] * GR[3 4; -3] end return y @@ -323,11 +334,9 @@ function (H::JordanMPO_AC2_Hamiltonian)(x::MPOTensor) ismissing(H.DE) || @plansor y[-1 -2; -3 -4] += x[-1 1; -3 -4] * H.DE[-2; 1] ismissing(H.EE) || @plansor y[-1 -2; -3 -4] += x[1 -2; -3 -4] * H.EE[-1; 1] - GL, A1, A2, GR = H.AA - if nonzero_length(A1) > 0 && nonzero_length(A2) > 0 - # TODO: there are too many entries here, this could be further optimized - @plansor y[-1 -2; -3 -4] += GL[-1 7; 6] * x[6 5; 1 3] * A1[7 -2; 5 4] * - A2[4 -4; 3 2] * GR[1 2; -3] + if !ismissing(H.AA) + GLW, GWR = H.AA + @tensor y[-1 -2; -3 -4] += GLW[-1 -2 3; 1 2] * x[1 2; 4 5] * GWR[3 4 5; -3 -4] end return y From 88cab3eb2c9eda00495f11d4d3e3452dc75613ec Mon Sep 17 00:00:00 2001 From: afeuerpfeil Date: Mon, 6 Oct 2025 19:42:42 -0400 Subject: [PATCH 02/19] fix formatting and make einstein labels for contracted indices consistent --- src/algorithms/derivatives/hamiltonian_derivatives.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/algorithms/derivatives/hamiltonian_derivatives.jl b/src/algorithms/derivatives/hamiltonian_derivatives.jl index edf8d9bcb..ab74a525b 100644 --- a/src/algorithms/derivatives/hamiltonian_derivatives.jl +++ b/src/algorithms/derivatives/hamiltonian_derivatives.jl @@ -293,8 +293,8 @@ function AC2_hamiltonian( # TODO: One should only calculate these operators if necessary. One can do this by checking nonzero_length(A1) > 0 && nonzero_length(A2) > 0 and then setting AA=missing or as we discussed via bit matrix multiplication. ## TODO: Think about how one could and whether one should store these objects and use them for (a) advancing environments in iDMRG, (b) reuse ind backwards-sweep in IDMRG, (c) subspace expansion if nonzero_length(W1.A) > 0 && nonzero_length(W2.A) > 0 - @tensor GLW[-1 -2 -3; -4 -5] ≔ GL[2:(end -1)][-1 1; -4] * W1.A[1 -2; -5 -3] - @tensor GWR[-1 -2 -3; -4 -5] ≔ W2.A[-1 -5; -3 1] * GR[2:(end -1)][-2 1; -4] + @tensor GLW[-1 -2 -3; -4 -5] ≔ GL[2:(end - 1)][-1 1; -4] * W1.A[1 -2; -5 -3] + @tensor GWR[-1 -2 -3; -4 -5] ≔ W2.A[-1 -5; -3 1] * GR[2:(end - 1)][-2 1; -4] AA = (GLW, GWR) else AA = missing @@ -316,7 +316,7 @@ function (H::JordanMPO_AC_Hamiltonian)(x::MPSTensor) if !ismissing(H.A) GLW, GR = H.A - @tensor y[-1 -2; -3] += GLW[-1 -2 3; 1 2] * x[1 2; 4] * GR[3 4; -3] + @tensor y[-1 -2; -3] += GLW[-1 -2 4; 3 2] * x[3 2; 1] * GR[4 1; -3] end return y @@ -336,7 +336,7 @@ function (H::JordanMPO_AC2_Hamiltonian)(x::MPOTensor) if !ismissing(H.AA) GLW, GWR = H.AA - @tensor y[-1 -2; -3 -4] += GLW[-1 -2 3; 1 2] * x[1 2; 4 5] * GWR[3 4 5; -3 -4] + @tensor y[-1 -2; -3 -4] += GLW[-1 -2 5; 4 3] * x[4 3; 2 1] * GWR[5 2 1; -3 -4] end return y From a975ba39bf32de5db2ace27c9d40dc9db08242bb Mon Sep 17 00:00:00 2001 From: afeuerpfeil Date: Mon, 6 Oct 2025 19:48:36 -0400 Subject: [PATCH 03/19] revert change for contraction order --- src/algorithms/derivatives/hamiltonian_derivatives.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/algorithms/derivatives/hamiltonian_derivatives.jl b/src/algorithms/derivatives/hamiltonian_derivatives.jl index ab74a525b..aeb934420 100644 --- a/src/algorithms/derivatives/hamiltonian_derivatives.jl +++ b/src/algorithms/derivatives/hamiltonian_derivatives.jl @@ -316,7 +316,7 @@ function (H::JordanMPO_AC_Hamiltonian)(x::MPSTensor) if !ismissing(H.A) GLW, GR = H.A - @tensor y[-1 -2; -3] += GLW[-1 -2 4; 3 2] * x[3 2; 1] * GR[4 1; -3] + @tensor y[-1 -2; -3] += GLW[-1 -2 3; 1 2] * x[1 2; 4] * GR[3 4; -3] end return y @@ -336,7 +336,7 @@ function (H::JordanMPO_AC2_Hamiltonian)(x::MPOTensor) if !ismissing(H.AA) GLW, GWR = H.AA - @tensor y[-1 -2; -3 -4] += GLW[-1 -2 5; 4 3] * x[4 3; 2 1] * GWR[5 2 1; -3 -4] + @tensor y[-1 -2; -3 -4] += GLW[-1 -2 3; 1 2] * x[1 2; 4 5] * GWR[3 4 5; -3 -4] end return y From 46a9dc729d8046c70ec1bff3d29122555b0c2d78 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Mon, 6 Oct 2025 21:59:40 -0400 Subject: [PATCH 04/19] update planar contractions --- .../derivatives/hamiltonian_derivatives.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/algorithms/derivatives/hamiltonian_derivatives.jl b/src/algorithms/derivatives/hamiltonian_derivatives.jl index aeb934420..523e377df 100644 --- a/src/algorithms/derivatives/hamiltonian_derivatives.jl +++ b/src/algorithms/derivatives/hamiltonian_derivatives.jl @@ -139,8 +139,8 @@ function AC_hamiltonian( end if nonzero_length(W.A) > 0 - @tensor GLW[-1 -2 -3; -4 -5] ≔ GL[2:(end - 1)][-1 1; -4] * W.A[1 -2; -5 -3] - continuing = (GLW, permute(GR[2:(end - 1)], (2, 1), (3,))) + @plansor GLW[-1 -2 -3; -4 -5] ≔ GL[2:(end - 1)][-1 1; -4] * W.A[1 -2; -5 -3] + continuing = (GLW, GR[2:(end - 1)]) else continuing = missing end @@ -293,8 +293,8 @@ function AC2_hamiltonian( # TODO: One should only calculate these operators if necessary. One can do this by checking nonzero_length(A1) > 0 && nonzero_length(A2) > 0 and then setting AA=missing or as we discussed via bit matrix multiplication. ## TODO: Think about how one could and whether one should store these objects and use them for (a) advancing environments in iDMRG, (b) reuse ind backwards-sweep in IDMRG, (c) subspace expansion if nonzero_length(W1.A) > 0 && nonzero_length(W2.A) > 0 - @tensor GLW[-1 -2 -3; -4 -5] ≔ GL[2:(end - 1)][-1 1; -4] * W1.A[1 -2; -5 -3] - @tensor GWR[-1 -2 -3; -4 -5] ≔ W2.A[-1 -5; -3 1] * GR[2:(end - 1)][-2 1; -4] + @plansor GLW[-1 -2 -3; -4 -5] ≔ GL[2:(end - 1)][-1 1; -4] * W1.A[1 -2; -5 -3] + @plansor GWR[-1 -2 -3; -4 -5] ≔ W2.A[-3 -5; -2 1] * GR[2:(end - 1)][-1 1; -4] AA = (GLW, GWR) else AA = missing @@ -316,7 +316,7 @@ function (H::JordanMPO_AC_Hamiltonian)(x::MPSTensor) if !ismissing(H.A) GLW, GR = H.A - @tensor y[-1 -2; -3] += GLW[-1 -2 3; 1 2] * x[1 2; 4] * GR[3 4; -3] + @tensor y[-1 -2; -3] += GLW[-1 -2 4; 1 2] * x[1 2; 3] * GR[3 4; -3] end return y @@ -336,7 +336,7 @@ function (H::JordanMPO_AC2_Hamiltonian)(x::MPOTensor) if !ismissing(H.AA) GLW, GWR = H.AA - @tensor y[-1 -2; -3 -4] += GLW[-1 -2 3; 1 2] * x[1 2; 4 5] * GWR[3 4 5; -3 -4] + @plansor y[-1 -2; -3 -4] += GLW[-1 -2 5; 1 2] * x[1 2; 3 4] * GWR[3 4 5; -3 -4] end return y From 0a8b6d0affbf4bdd02e1e327b79f6b615d123009 Mon Sep 17 00:00:00 2001 From: afeuerpfeil Date: Tue, 7 Oct 2025 07:58:11 -0400 Subject: [PATCH 05/19] add missing plansor and fix indexing into messing error --- .../derivatives/mpo_derivatives copy.jl | 122 ++++++++++++++++++ 1 file changed, 122 insertions(+) create mode 100644 src/algorithms/derivatives/mpo_derivatives copy.jl diff --git a/src/algorithms/derivatives/mpo_derivatives copy.jl b/src/algorithms/derivatives/mpo_derivatives copy.jl new file mode 100644 index 000000000..ac252803c --- /dev/null +++ b/src/algorithms/derivatives/mpo_derivatives copy.jl @@ -0,0 +1,122 @@ +""" + struct MPODerivativeOperator{L,O<:Tuple,R} + +Effective local operator obtained from taking the partial derivative of an MPS-MPO-MPS sandwich. +""" +struct MPODerivativeOperator{L, R, O} <: DerivativeOperator + leftenv::L + rightenv::R +end + +Base.length(H::MPODerivativeOperator) = length(H.operators) + +const MPO_C_Hamiltonian{L, R} = MPODerivativeOperator{L, R, 0} +MPO_C_Hamiltonian(GL::L, GR::R) where {L, R} = MPODerivativeOperator{L,R,0}(GL, GR) + +const MPO_AC_Hamiltonian{L, R, O} = MPODerivativeOperator{L, R, O} +MPO_AC_Hamiltonian(GL::L, ::O, GR::R) where {L, R, O} = MPODerivativeOperator{L,R,O}(GL, GR) + +const MPO_AC2_Hamiltonian{L, R} = MPODerivativeOperator{L, R, 2} +MPO_AC2_Hamiltonian(GL::L, ::O₁, ::O₂, GR::R) where {L, R, O₁, O₂} = MPODerivativeOperator{L,R,Tuple{O₁,O₂}}(GL, GR) + + +# Constructors +# ------------ +function C_hamiltonian(site::Int, below, operator, above, envs) + return MPO_C_Hamiltonian(leftenv(envs, site + 1, below), rightenv(envs, site, below)) +end +function AC_hamiltonian(site::Int, below, operator<:MPOTensor, above, envs) + O = operator[site] + L = leftenv(envs, site, below) + @plansor L[-1 -2 -3; -4 -5] ≔ L[-1 1; -4] * O[1 -2; -5 -3] + R = rightenv(envs, site, below) + R = permute(R, (2, 1), (3,)) + return MPO_AC_Hamiltonian(L, O, R) +end +function AC_hamiltonian(site::Int, below, operator, above, envs) + O = isnothing(operator) ? nothing : operator[site] + return MPO_AC_Hamiltonian(leftenv(envs, site, below), O, rightenv(envs, site, below)) +end +function AC2_hamiltonian(site::Int, below, operator<:MPOTensor, above, envs) + O1, O2 = operator[site], operator[site + 1] + L = leftenv(envs, site, below) + @tensor GLW[-1 -2 -3; -4 -5] ≔ L[-1 1; -4] * O1[1 -2; -5 -3] + @tensor GWR[-1 -2 -3; -4 -5] ≔ O2[-1 -5; -3 1] * R[-2 1; -4] + return MPO_AC2_Hamiltonian(L, O1, O2, R) +end +function AC2_hamiltonian(site::Int, below, operator, above, envs) + O1, O2 = isnothing(operator) ? (nothing, nothing) : (operator[site], operator[site + 1]) + return MPO_AC2_Hamiltonian( + leftenv(envs, site, below), O1, O2, rightenv(envs, site + 1, below) + ) +end + +# Properties +# ---------- +function TensorKit.domain(H::MPODerivativeOperator{L, R, O}) where {L, R, O} + V_l = right_virtualspace(H.leftenv) + V_r = left_virtualspace(H.rightenv) + V_o = prod(physicalspace, H.O; init = one(V_l)) + return V_l ⊗ V_o ⊗ V_r +end +function TensorKit.codomain(H::MPODerivativeOperator) + V_l = left_virtualspace(H.leftenv) + V_r = right_virtualspace(H.rightenv) + V_o = prod(physicalspace, H.O; init = one(V_l)) + return V_l ⊗ V_o ⊗ V_r +end + +# Actions +# ------- +function (h::MPO_C_Hamiltonian{<:MPSBondTensor, <:MPSBondTensor})(x::MPSBondTensor) + @plansor y[-1; -2] ≔ h.leftenv[-1; 1] * x[1; 2] * h.rightenv[2; -2] + return y isa AbstractBlockTensorMap ? only(y) : y +end +function (h::MPO_C_Hamiltonian{<:MPSTensor, <:MPSTensor})(x::MPSBondTensor) + @plansor y[-1; -2] ≔ h.leftenv[-1 3; 1] * x[1; 2] * h.rightenv[2 3; -2] + return y isa AbstractBlockTensorMap ? only(y) : y +end + +function (h::MPO_AC_Hamiltonian{<:MPSTensor, <:MPOTensor, <:MPSTensor})(x::MPSTensor) + @plansor y[-1 -2; -3] ≔ h.leftenv[-1 5; 4] * x[4 2; 1] * h.operators[1][5 -2; 2 3] * + h.rightenv[1 3; -3] + return y isa AbstractBlockTensorMap ? only(y) : y +end +function (h::MPO_AC_Hamiltonian{<:MPSTensor, <:Number, <:MPSTensor})(x::MPSTensor) + @plansor y[-1 -2; -3] ≔ ( + h.leftenv[-1 5; 4] * x[4 6; 1] * τ[6 5; 7 -2] * h.rightenv[1 7; -3] + ) * only(h.operators) + return y isa AbstractBlockTensorMap ? only(y) : y +end +function (h::MPO_AC_Hamiltonian{<:MPSBondTensor, Nothing, <:MPSBondTensor})(x::MPSTensor) + return @plansor y[-1 -2; -3] ≔ h.leftenv[-1; 2] * x[2 -2; 1] * h.rightenv[1; -3] +end +function (h::MPO_AC_Hamiltonian{<:MPSTensor, <:MPOTensor, <:MPSTensor})( + x::GenericMPSTensor{<:Any, 3} + ) + @plansor y[-1 -2 -3; -4] ≔ h.leftenv[-1 7; 6] * x[6 4 2; 1] * + h.operators[1][7 -2; 4 5] * τ[5 -3; 2 3] * h.rightenv[1 3; -4] + return y isa AbstractBlockTensorMap ? only(y) : y +end + +function (h::MPO_AC2_Hamiltonian{<:MPSBondTensor, Nothing, Nothing, <:MPSBondTensor})( + x::MPOTensor + ) + @plansor y[-1 -2; -3 -4] ≔ h.leftenv[-1; 1] * x[1 -2; 2 -4] * h.rightenv[2 -3] + return y isa AbstractBlockTensorMap ? only(y) : y +end +function (h::MPO_AC2_Hamiltonian{<:MPSTensor, <:MPOTensor, <:MPOTensor, <:MPSTensor})( + x::MPOTensor + ) + @plansor y[-1 -2; -3 -4] ≔ h.leftenv[-1 7; 6] * x[6 5; 1 3] * + h.operators[1][7 -2; 5 4] * h.operators[2][4 -4; 3 2] * h.rightenv[1 2; -3] + return y isa AbstractBlockTensorMap ? only(y) : y +end +function (h::MPO_AC2_Hamiltonian{<:MPSTensor, <:MPOTensor, <:MPOTensor, <:MPSTensor})( + x::AbstractTensorMap{<:Any, <:Any, 3, 3} + ) + @plansor y[-1 -2 -3; -4 -5 -6] ≔ h.leftenv[-1 11; 10] * x[10 8 6; 1 2 4] * + h.rightenv[1 3; -4] * h.operators[1][11 -2; 8 9] * τ[9 -3; 6 7] * + h.operators[2][7 -6; 4 5] * τ[5 -5; 2 3] + return y isa AbstractBlockTensorMap ? only(y) : y +end From 5ea9055d38f9b9e3509227efe9d00170db7d91ea Mon Sep 17 00:00:00 2001 From: afeuerpfeil Date: Tue, 7 Oct 2025 07:58:50 -0400 Subject: [PATCH 06/19] Commit correct file --- .../derivatives/hamiltonian_derivatives.jl | 6 +- .../derivatives/mpo_derivatives copy.jl | 122 ------------------ 2 files changed, 3 insertions(+), 125 deletions(-) delete mode 100644 src/algorithms/derivatives/mpo_derivatives copy.jl diff --git a/src/algorithms/derivatives/hamiltonian_derivatives.jl b/src/algorithms/derivatives/hamiltonian_derivatives.jl index 523e377df..194e1e000 100644 --- a/src/algorithms/derivatives/hamiltonian_derivatives.jl +++ b/src/algorithms/derivatives/hamiltonian_derivatives.jl @@ -16,7 +16,7 @@ struct JordanMPO_AC_Hamiltonian{O1, O2, O3} <: DerivativeOperator ) # obtaining storagetype of environments since these should have already mixed # the types of the operator and state - gl = continuing[1] + gl = @coalesce(onsite, not_started, finished, starting, ending) S = spacetype(gl) M = storagetype(gl) O1 = tensormaptype(S, 1, 1, M) @@ -52,7 +52,7 @@ struct JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4} <: DerivativeOperator function JordanMPO_AC2_Hamiltonian(II, IC, ID, CB, CA, AB, AA, BE, DE, EE) # obtaining storagetype of environments since these should have already mixed # the types of the operator and state - gl = AA[1] + gl = @coalesce(II, IC, ID, CB, CA, AB, BE, DE, EE) S = spacetype(gl) M = storagetype(gl) O1 = tensormaptype(S, 1, 1, M) @@ -316,7 +316,7 @@ function (H::JordanMPO_AC_Hamiltonian)(x::MPSTensor) if !ismissing(H.A) GLW, GR = H.A - @tensor y[-1 -2; -3] += GLW[-1 -2 4; 1 2] * x[1 2; 3] * GR[3 4; -3] + @plansor y[-1 -2; -3] += GLW[-1 -2 3; 1 2] * x[1 2; 4] * GR[4 3; -3] end return y diff --git a/src/algorithms/derivatives/mpo_derivatives copy.jl b/src/algorithms/derivatives/mpo_derivatives copy.jl deleted file mode 100644 index ac252803c..000000000 --- a/src/algorithms/derivatives/mpo_derivatives copy.jl +++ /dev/null @@ -1,122 +0,0 @@ -""" - struct MPODerivativeOperator{L,O<:Tuple,R} - -Effective local operator obtained from taking the partial derivative of an MPS-MPO-MPS sandwich. -""" -struct MPODerivativeOperator{L, R, O} <: DerivativeOperator - leftenv::L - rightenv::R -end - -Base.length(H::MPODerivativeOperator) = length(H.operators) - -const MPO_C_Hamiltonian{L, R} = MPODerivativeOperator{L, R, 0} -MPO_C_Hamiltonian(GL::L, GR::R) where {L, R} = MPODerivativeOperator{L,R,0}(GL, GR) - -const MPO_AC_Hamiltonian{L, R, O} = MPODerivativeOperator{L, R, O} -MPO_AC_Hamiltonian(GL::L, ::O, GR::R) where {L, R, O} = MPODerivativeOperator{L,R,O}(GL, GR) - -const MPO_AC2_Hamiltonian{L, R} = MPODerivativeOperator{L, R, 2} -MPO_AC2_Hamiltonian(GL::L, ::O₁, ::O₂, GR::R) where {L, R, O₁, O₂} = MPODerivativeOperator{L,R,Tuple{O₁,O₂}}(GL, GR) - - -# Constructors -# ------------ -function C_hamiltonian(site::Int, below, operator, above, envs) - return MPO_C_Hamiltonian(leftenv(envs, site + 1, below), rightenv(envs, site, below)) -end -function AC_hamiltonian(site::Int, below, operator<:MPOTensor, above, envs) - O = operator[site] - L = leftenv(envs, site, below) - @plansor L[-1 -2 -3; -4 -5] ≔ L[-1 1; -4] * O[1 -2; -5 -3] - R = rightenv(envs, site, below) - R = permute(R, (2, 1), (3,)) - return MPO_AC_Hamiltonian(L, O, R) -end -function AC_hamiltonian(site::Int, below, operator, above, envs) - O = isnothing(operator) ? nothing : operator[site] - return MPO_AC_Hamiltonian(leftenv(envs, site, below), O, rightenv(envs, site, below)) -end -function AC2_hamiltonian(site::Int, below, operator<:MPOTensor, above, envs) - O1, O2 = operator[site], operator[site + 1] - L = leftenv(envs, site, below) - @tensor GLW[-1 -2 -3; -4 -5] ≔ L[-1 1; -4] * O1[1 -2; -5 -3] - @tensor GWR[-1 -2 -3; -4 -5] ≔ O2[-1 -5; -3 1] * R[-2 1; -4] - return MPO_AC2_Hamiltonian(L, O1, O2, R) -end -function AC2_hamiltonian(site::Int, below, operator, above, envs) - O1, O2 = isnothing(operator) ? (nothing, nothing) : (operator[site], operator[site + 1]) - return MPO_AC2_Hamiltonian( - leftenv(envs, site, below), O1, O2, rightenv(envs, site + 1, below) - ) -end - -# Properties -# ---------- -function TensorKit.domain(H::MPODerivativeOperator{L, R, O}) where {L, R, O} - V_l = right_virtualspace(H.leftenv) - V_r = left_virtualspace(H.rightenv) - V_o = prod(physicalspace, H.O; init = one(V_l)) - return V_l ⊗ V_o ⊗ V_r -end -function TensorKit.codomain(H::MPODerivativeOperator) - V_l = left_virtualspace(H.leftenv) - V_r = right_virtualspace(H.rightenv) - V_o = prod(physicalspace, H.O; init = one(V_l)) - return V_l ⊗ V_o ⊗ V_r -end - -# Actions -# ------- -function (h::MPO_C_Hamiltonian{<:MPSBondTensor, <:MPSBondTensor})(x::MPSBondTensor) - @plansor y[-1; -2] ≔ h.leftenv[-1; 1] * x[1; 2] * h.rightenv[2; -2] - return y isa AbstractBlockTensorMap ? only(y) : y -end -function (h::MPO_C_Hamiltonian{<:MPSTensor, <:MPSTensor})(x::MPSBondTensor) - @plansor y[-1; -2] ≔ h.leftenv[-1 3; 1] * x[1; 2] * h.rightenv[2 3; -2] - return y isa AbstractBlockTensorMap ? only(y) : y -end - -function (h::MPO_AC_Hamiltonian{<:MPSTensor, <:MPOTensor, <:MPSTensor})(x::MPSTensor) - @plansor y[-1 -2; -3] ≔ h.leftenv[-1 5; 4] * x[4 2; 1] * h.operators[1][5 -2; 2 3] * - h.rightenv[1 3; -3] - return y isa AbstractBlockTensorMap ? only(y) : y -end -function (h::MPO_AC_Hamiltonian{<:MPSTensor, <:Number, <:MPSTensor})(x::MPSTensor) - @plansor y[-1 -2; -3] ≔ ( - h.leftenv[-1 5; 4] * x[4 6; 1] * τ[6 5; 7 -2] * h.rightenv[1 7; -3] - ) * only(h.operators) - return y isa AbstractBlockTensorMap ? only(y) : y -end -function (h::MPO_AC_Hamiltonian{<:MPSBondTensor, Nothing, <:MPSBondTensor})(x::MPSTensor) - return @plansor y[-1 -2; -3] ≔ h.leftenv[-1; 2] * x[2 -2; 1] * h.rightenv[1; -3] -end -function (h::MPO_AC_Hamiltonian{<:MPSTensor, <:MPOTensor, <:MPSTensor})( - x::GenericMPSTensor{<:Any, 3} - ) - @plansor y[-1 -2 -3; -4] ≔ h.leftenv[-1 7; 6] * x[6 4 2; 1] * - h.operators[1][7 -2; 4 5] * τ[5 -3; 2 3] * h.rightenv[1 3; -4] - return y isa AbstractBlockTensorMap ? only(y) : y -end - -function (h::MPO_AC2_Hamiltonian{<:MPSBondTensor, Nothing, Nothing, <:MPSBondTensor})( - x::MPOTensor - ) - @plansor y[-1 -2; -3 -4] ≔ h.leftenv[-1; 1] * x[1 -2; 2 -4] * h.rightenv[2 -3] - return y isa AbstractBlockTensorMap ? only(y) : y -end -function (h::MPO_AC2_Hamiltonian{<:MPSTensor, <:MPOTensor, <:MPOTensor, <:MPSTensor})( - x::MPOTensor - ) - @plansor y[-1 -2; -3 -4] ≔ h.leftenv[-1 7; 6] * x[6 5; 1 3] * - h.operators[1][7 -2; 5 4] * h.operators[2][4 -4; 3 2] * h.rightenv[1 2; -3] - return y isa AbstractBlockTensorMap ? only(y) : y -end -function (h::MPO_AC2_Hamiltonian{<:MPSTensor, <:MPOTensor, <:MPOTensor, <:MPSTensor})( - x::AbstractTensorMap{<:Any, <:Any, 3, 3} - ) - @plansor y[-1 -2 -3; -4 -5 -6] ≔ h.leftenv[-1 11; 10] * x[10 8 6; 1 2 4] * - h.rightenv[1 3; -4] * h.operators[1][11 -2; 8 9] * τ[9 -3; 6 7] * - h.operators[2][7 -6; 4 5] * τ[5 -5; 2 3] - return y isa AbstractBlockTensorMap ? only(y) : y -end From f3290c3a708f2bed69d0d07f0070f10b50610a09 Mon Sep 17 00:00:00 2001 From: afeuerpfeil Date: Tue, 7 Oct 2025 09:24:22 -0400 Subject: [PATCH 07/19] Port changes from hamiltonian_derivatives to mpo_derivatives and reuse code from mpo_derivatives --- .../derivatives/hamiltonian_derivatives.jl | 42 ++++------ src/algorithms/derivatives/mpo_derivatives.jl | 78 ++++++++++++++++++- 2 files changed, 93 insertions(+), 27 deletions(-) diff --git a/src/algorithms/derivatives/hamiltonian_derivatives.jl b/src/algorithms/derivatives/hamiltonian_derivatives.jl index 194e1e000..2358c5aa5 100644 --- a/src/algorithms/derivatives/hamiltonian_derivatives.jl +++ b/src/algorithms/derivatives/hamiltonian_derivatives.jl @@ -69,12 +69,6 @@ end # Constructors # ------------ -const _HAM_MPS_TYPES = Union{ - FiniteMPS{<:MPSTensor}, - WindowMPS{<:MPSTensor}, - InfiniteMPS{<:MPSTensor}, -} - function AC_hamiltonian( site::Int, below::_HAM_MPS_TYPES, operator::MPOHamiltonian{<:JordanMPOTensor}, above::_HAM_MPS_TYPES, envs @@ -83,9 +77,11 @@ function AC_hamiltonian( GR = rightenv(envs, site, below) W = operator[site] + GR_2 = GR[2:(end - 1)] + GL_2 = GL[2:(end - 1)] + # starting if nonzero_length(W.C) > 0 - GR_2 = GR[2:(end - 1)] @plansor starting_[-1 -2; -3 -4] ≔ W.C[-1; -3 1] * GR_2[-4 1; -2] starting = only(starting_) else @@ -94,7 +90,6 @@ function AC_hamiltonian( # ending if nonzero_length(W.B) > 0 - GL_2 = GL[2:(end - 1)] @plansor ending_[-1 -2; -3 -4] ≔ GL_2[-1 1; -3] * W.B[1 -2; -4] ending = only(ending_) else @@ -139,8 +134,7 @@ function AC_hamiltonian( end if nonzero_length(W.A) > 0 - @plansor GLW[-1 -2 -3; -4 -5] ≔ GL[2:(end - 1)][-1 1; -4] * W.A[1 -2; -5 -3] - continuing = (GLW, GR[2:(end - 1)]) + continuing = AC_hamiltonian(GL_2, W.A, GR_2) else continuing = missing end @@ -159,10 +153,13 @@ function AC2_hamiltonian( W1 = operator[site] W2 = operator[site + 1] + GR_2 = GR[2:(end - 1)] + GL_2 = GL[2:(end - 1)] + # starting left - continuing right if nonzero_length(W1.C) > 0 && nonzero_length(W2.A) > 0 @plansor CA_[-1 -2 -3; -4 -5 -6] ≔ W1.C[-1; -4 2] * W2.A[2 -2; -5 1] * - GR[2:(end - 1)][-6 1; -3] + GR_2[-6 1; -3] CA = only(CA_) else CA = missing @@ -170,7 +167,7 @@ function AC2_hamiltonian( # continuing left - ending right if nonzero_length(W1.A) > 0 && nonzero_length(W2.B) > 0 - @plansor AB_[-1 -2 -3; -4 -5 -6] ≔ GL[2:(end - 1)][-1 2; -4] * W1.A[2 -2; -5 1] * + @plansor AB_[-1 -2 -3; -4 -5 -6] ≔ GL_2[-1 2; -4] * W1.A[2 -2; -5 1] * W2.B[1 -3; -6] AB = only(AB_) else @@ -200,10 +197,10 @@ function AC2_hamiltonian( if !ismissing(CA) I = id(storagetype(GR[1]), physicalspace(W1)) @plansor CA[-1 -2 -3; -4 -5 -6] += (I[-1; -4] * W2.C[-2; -5 1]) * - GR[2:(end - 1)][-6 1; -3] + GR_2[-6 1; -3] IC = missing else - @plansor IC[-1 -2; -3 -4] ≔ W2.C[-1; -3 1] * GR[2:(end - 1)][-4 1; -2] + @plansor IC[-1 -2; -3 -4] ≔ W2.C[-1; -3 1] * GR_2[-4 1; -2] end else IC = missing @@ -213,11 +210,11 @@ function AC2_hamiltonian( if nonzero_length(W1.B) > 0 if !ismissing(AB) I = id(storagetype(GL[end]), physicalspace(W2)) - @plansor AB[-1 -2 -3; -4 -5 -6] += GL[2:(end - 1)][-1 1; -4] * + @plansor AB[-1 -2 -3; -4 -5 -6] += GL_2[-1 1; -4] * (W1.B[1 -2; -5] * I[-3; -6]) BE = missing else - @plansor BE[-1 -2; -3 -4] ≔ GL[2:(end - 1)][-1 2; -3] * W1.B[2 -2; -4] + @plansor BE[-1 -2; -3 -4] ≔ GL_2[-1 2; -3] * W1.B[2 -2; -4] end else BE = missing @@ -290,12 +287,9 @@ function AC2_hamiltonian( # continuing - continuing # TODO: MPODerivativeOperator code reuse + optimization - # TODO: One should only calculate these operators if necessary. One can do this by checking nonzero_length(A1) > 0 && nonzero_length(A2) > 0 and then setting AA=missing or as we discussed via bit matrix multiplication. ## TODO: Think about how one could and whether one should store these objects and use them for (a) advancing environments in iDMRG, (b) reuse ind backwards-sweep in IDMRG, (c) subspace expansion if nonzero_length(W1.A) > 0 && nonzero_length(W2.A) > 0 - @plansor GLW[-1 -2 -3; -4 -5] ≔ GL[2:(end - 1)][-1 1; -4] * W1.A[1 -2; -5 -3] - @plansor GWR[-1 -2 -3; -4 -5] ≔ W2.A[-3 -5; -2 1] * GR[2:(end - 1)][-1 1; -4] - AA = (GLW, GWR) + AA = AC2_hamiltonian(GL_2, W1.A, W2.A, GR_2) else AA = missing end @@ -313,10 +307,8 @@ function (H::JordanMPO_AC_Hamiltonian)(x::MPSTensor) ismissing(H.I) || @plansor y[-1 -2; -3] += x[-1 -2; 1] * H.I[1; -3] ismissing(H.C) || @plansor y[-1 -2; -3] += x[-1 2; 1] * H.C[-2 -3; 2 1] ismissing(H.B) || @plansor y[-1 -2; -3] += H.B[-1 -2; 1 2] * x[1 2; -3] - if !ismissing(H.A) - GLW, GR = H.A - @plansor y[-1 -2; -3] += GLW[-1 -2 3; 1 2] * x[1 2; 4] * GR[4 3; -3] + y += H.A(x) end return y @@ -333,10 +325,8 @@ function (H::JordanMPO_AC2_Hamiltonian)(x::MPOTensor) ismissing(H.BE) || @plansor y[-1 -2; -3 -4] += x[1 2; -3 -4] * H.BE[-1 -2; 1 2] ismissing(H.DE) || @plansor y[-1 -2; -3 -4] += x[-1 1; -3 -4] * H.DE[-2; 1] ismissing(H.EE) || @plansor y[-1 -2; -3 -4] += x[1 -2; -3 -4] * H.EE[-1; 1] - if !ismissing(H.AA) - GLW, GWR = H.AA - @plansor y[-1 -2; -3 -4] += GLW[-1 -2 5; 1 2] * x[1 2; 3 4] * GWR[3 4 5; -3 -4] + y += H.AA(x) end return y diff --git a/src/algorithms/derivatives/mpo_derivatives.jl b/src/algorithms/derivatives/mpo_derivatives.jl index 12b030b1e..8bd6dd4f5 100644 --- a/src/algorithms/derivatives/mpo_derivatives.jl +++ b/src/algorithms/derivatives/mpo_derivatives.jl @@ -9,9 +9,15 @@ struct MPODerivativeOperator{L, O <: Tuple, R} <: DerivativeOperator rightenv::R end +struct MPOContractedDerivativeOperator{L, R, N} <: DerivativeOperator + leftenv::L + rightenv::R +end + Base.length(H::MPODerivativeOperator) = length(H.operators) +Base.length(::MPOContractedDerivativeOperator{L, R, N}) where {L, R, N} = N -const MPO_C_Hamiltonian{L, R} = MPODerivativeOperator{L, Tuple{}, R} +const MPO_C_Hamiltonian{L, R} = MPOContractedDerivativeOperator{L, Tuple{}, R} MPO_C_Hamiltonian(GL, GR) = MPODerivativeOperator(GL, (), GR) const MPO_AC_Hamiltonian{L, O, R} = MPODerivativeOperator{L, Tuple{O}, R} @@ -20,8 +26,23 @@ MPO_AC_Hamiltonian(GL, O, GR) = MPODerivativeOperator(GL, (O,), GR) const MPO_AC2_Hamiltonian{L, O₁, O₂, R} = MPODerivativeOperator{L, Tuple{O₁, O₂}, R} MPO_AC2_Hamiltonian(GL, O1, O2, GR) = MPODerivativeOperator(GL, (O1, O2), GR) +const MPO_Contracted_C_Hamiltonian{L, R} = MPOContractedDerivativeOperator{L, R, 0} +MPO_Contracted_C_Hamiltonian(GL::L, GR::R) where {L, R} = MPOContractedDerivativeOperator{L,R,0}(GL, GR) + +const MPO_Contracted_AC_Hamiltonian{L, R} = MPOContractedDerivativeOperator{L, R, 1} +MPO_Contracted_AC_Hamiltonian(GL::L, GR::R) where {L, R} = MPOContractedDerivativeOperator{L,R,1}(GL, GR) + +const MPO_Contracted_AC2_Hamiltonian{L, R} = MPOContractedDerivativeOperator{L, R, 2} +MPO_Contracted_AC2_Hamiltonian(GL::L, GR::R) where {L, R} = MPOContractedDerivativeOperator{L,R,2}(GL, GR) + # Constructors # ------------ +const _HAM_MPS_TYPES = Union{ + FiniteMPS{<:MPSTensor}, + WindowMPS{<:MPSTensor}, + InfiniteMPS{<:MPSTensor}, +} + function C_hamiltonian(site::Int, below, operator, above, envs) return MPO_C_Hamiltonian(leftenv(envs, site + 1, below), rightenv(envs, site, below)) end @@ -35,6 +56,30 @@ function AC2_hamiltonian(site::Int, below, operator, above, envs) leftenv(envs, site, below), O1, O2, rightenv(envs, site + 1, below) ) end +function C_hamiltonian(site::Int, below::_HAM_MPS_TYPES, operator::MPOHamiltonian{<:MPOTensor}, above::_HAM_MPS_TYPES, envs) + return MPO_Contracted_C_Hamiltonian(leftenv(envs, site + 1, below), rightenv(envs, site, below)) +end +function AC_hamiltonian(site::Int, below::_HAM_MPS_TYPES, operator::MPOHamiltonian{<:MPOTensor}, above::_HAM_MPS_TYPES, envs) + O = operator[site] + GL = leftenv(envs, site, below) + return AC_hamiltonian(GL, O, rightenv(envs, site, below)) +end +function AC_hamiltonian(GL::MPSTensor, O::MPOTensor, GR::MPSTensor) + @plansor GLW[-1 -2 -3; -4 -5] ≔ GL[-1 1; -4] * O[1 -2; -5 -3] + return MPO_Contracted_AC_Hamiltonian(GLW, O, GR) +end +function AC2_hamiltonian(site::Int, below::_HAM_MPS_TYPES, operator::MPOHamiltonian{<:MPOTensor}, above::_HAM_MPS_TYPES, envs) + O1 = operator[site] + O2 = operator[site + 1] + GL = leftenv(envs, site, below) + GR = rightenv(envs, site + 1, below) + return AC2_hamiltonian(GL, O1, O2, GR) +end +function AC2_hamiltonian(GL::MPSTensor, O1::MPOTensor, O2::MPOTensor, GR::MPSTensor) + @plansor GLW[-1 -2 -3; -4 -5] ≔ GL[-1 1; -4] * O1[1 -2; -5 -3] + @plansor GWR[-1 -2 -3; -4 -5] ≔ O2[-3 -5; -2 1] * GR[-1 1; -4] + return MPO_Contracted_AC2_Hamiltonian(GLW, GWR) +end # Properties # ---------- @@ -50,6 +95,20 @@ function TensorKit.codomain(H::MPODerivativeOperator) V_o = prod(physicalspace, H.O; init = one(V_l)) return V_l ⊗ V_o ⊗ V_r end +function TensorKit.domain(H::MPOContractedDerivativeOperator) + V_l = right_virtualspace(H.leftenv) + V_r = left_virtualspace(H.rightenv) + ## TODO: How to deal with the H.O here? + V_o = prod(physicalspace, H.O; init = one(V_l)) + return V_l ⊗ V_o ⊗ V_r +end +function TensorKit.codomain(H::MPOContractedDerivativeOperator) + V_l = left_virtualspace(H.leftenv) + V_r = right_virtualspace(H.rightenv) + ## TODO: How to deal with the H.O here? + V_o = prod(physicalspace, H.O; init = one(V_l)) + return V_l ⊗ V_o ⊗ V_r +end # Actions # ------- @@ -105,3 +164,20 @@ function (h::MPO_AC2_Hamiltonian{<:MPSTensor, <:MPOTensor, <:MPOTensor, <:MPSTen h.operators[2][7 -6; 4 5] * τ[5 -5; 2 3] return y isa AbstractBlockTensorMap ? only(y) : y end +function (h::MPO_Contracted_C_Hamiltonian)(x::MPSBondTensor) + @plansor y[-1; -2] ≔ h.leftenv[-1 3; 1] * x[1; 2] * h.rightenv[2 3; -2] + return y isa AbstractBlockTensorMap ? only(y) : y +end +function (h::MPO_Contracted_AC_Hamiltonian)( + x::GenericMPSTensor{<:Any, 3} + ) + @plansor y[-1 -2; -3] ≔ h.leftenv[-1 -2 3; 1 2] * x[1 2; 4] * h.rightenv[4 3; -3] + return y isa AbstractBlockTensorMap ? only(y) : y +end +function (h::MPO_Contracted_AC2_Hamiltonian)( + x::MPOTensor + ) + @plansor y[-1 -2; -3 -4] ≔ h.leftenv[-1 -2 5; 1 2] * x[1 2; 3 4] * h.rightenv[3 4 5; -3 -4] + return y isa AbstractBlockTensorMap ? only(y) : y +end +## TODO: Which interface should we use for an inplace += version to be used in hamiltonian_derivatives.jl? \ No newline at end of file From ade24dee89b9f511033fe5fb448c8ffd3e957565 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Tue, 7 Oct 2025 09:42:05 -0400 Subject: [PATCH 08/19] avoid additional allocation --- src/algorithms/derivatives/hamiltonian_derivatives.jl | 11 ++--------- 1 file changed, 2 insertions(+), 9 deletions(-) diff --git a/src/algorithms/derivatives/hamiltonian_derivatives.jl b/src/algorithms/derivatives/hamiltonian_derivatives.jl index 2358c5aa5..21ff79c6b 100644 --- a/src/algorithms/derivatives/hamiltonian_derivatives.jl +++ b/src/algorithms/derivatives/hamiltonian_derivatives.jl @@ -300,22 +300,18 @@ end # Actions # ------- function (H::JordanMPO_AC_Hamiltonian)(x::MPSTensor) - y = zerovector(x) - + y = ismissing(H.A) ? zerovector(x) : H.A(x) ismissing(H.D) || @plansor y[-1 -2; -3] += x[-1 1; -3] * H.D[-2; 1] ismissing(H.E) || @plansor y[-1 -2; -3] += H.E[-1; 1] * x[1 -2; -3] ismissing(H.I) || @plansor y[-1 -2; -3] += x[-1 -2; 1] * H.I[1; -3] ismissing(H.C) || @plansor y[-1 -2; -3] += x[-1 2; 1] * H.C[-2 -3; 2 1] ismissing(H.B) || @plansor y[-1 -2; -3] += H.B[-1 -2; 1 2] * x[1 2; -3] - if !ismissing(H.A) - y += H.A(x) - end return y end function (H::JordanMPO_AC2_Hamiltonian)(x::MPOTensor) - y = zerovector(x) + y = ismissing(H.AA) ? zerovector(x) : H.AA(x) ismissing(H.II) || @plansor y[-1 -2; -3 -4] += x[-1 -2; 1 -4] * H.II[-3; 1] ismissing(H.IC) || @plansor y[-1 -2; -3 -4] += x[-1 -2; 1 2] * H.IC[-4 -3; 2 1] ismissing(H.ID) || @plansor y[-1 -2; -3 -4] += x[-1 -2; -3 1] * H.ID[-4; 1] @@ -325,9 +321,6 @@ function (H::JordanMPO_AC2_Hamiltonian)(x::MPOTensor) ismissing(H.BE) || @plansor y[-1 -2; -3 -4] += x[1 2; -3 -4] * H.BE[-1 -2; 1 2] ismissing(H.DE) || @plansor y[-1 -2; -3 -4] += x[-1 1; -3 -4] * H.DE[-2; 1] ismissing(H.EE) || @plansor y[-1 -2; -3 -4] += x[1 -2; -3 -4] * H.EE[-1; 1] - if !ismissing(H.AA) - y += H.AA(x) - end return y end From 6fc9e7acb87c14b72cbb64dcc151a147450e4bf5 Mon Sep 17 00:00:00 2001 From: afeuerpfeil Date: Tue, 7 Oct 2025 14:05:40 -0400 Subject: [PATCH 09/19] remove redundanct contracted_C operator and dispatch to the contracted derivates for MPO{<:MPOTensor} --- src/algorithms/derivatives/mpo_derivatives.jl | 23 +++++-------------- 1 file changed, 6 insertions(+), 17 deletions(-) diff --git a/src/algorithms/derivatives/mpo_derivatives.jl b/src/algorithms/derivatives/mpo_derivatives.jl index 8bd6dd4f5..02ec00b34 100644 --- a/src/algorithms/derivatives/mpo_derivatives.jl +++ b/src/algorithms/derivatives/mpo_derivatives.jl @@ -17,7 +17,7 @@ end Base.length(H::MPODerivativeOperator) = length(H.operators) Base.length(::MPOContractedDerivativeOperator{L, R, N}) where {L, R, N} = N -const MPO_C_Hamiltonian{L, R} = MPOContractedDerivativeOperator{L, Tuple{}, R} +const MPO_C_Hamiltonian{L, R} = MPODerivativeOperator{L, Tuple{}, R} MPO_C_Hamiltonian(GL, GR) = MPODerivativeOperator(GL, (), GR) const MPO_AC_Hamiltonian{L, O, R} = MPODerivativeOperator{L, Tuple{O}, R} @@ -26,9 +26,6 @@ MPO_AC_Hamiltonian(GL, O, GR) = MPODerivativeOperator(GL, (O,), GR) const MPO_AC2_Hamiltonian{L, O₁, O₂, R} = MPODerivativeOperator{L, Tuple{O₁, O₂}, R} MPO_AC2_Hamiltonian(GL, O1, O2, GR) = MPODerivativeOperator(GL, (O1, O2), GR) -const MPO_Contracted_C_Hamiltonian{L, R} = MPOContractedDerivativeOperator{L, R, 0} -MPO_Contracted_C_Hamiltonian(GL::L, GR::R) where {L, R} = MPOContractedDerivativeOperator{L,R,0}(GL, GR) - const MPO_Contracted_AC_Hamiltonian{L, R} = MPOContractedDerivativeOperator{L, R, 1} MPO_Contracted_AC_Hamiltonian(GL::L, GR::R) where {L, R} = MPOContractedDerivativeOperator{L,R,1}(GL, GR) @@ -56,19 +53,16 @@ function AC2_hamiltonian(site::Int, below, operator, above, envs) leftenv(envs, site, below), O1, O2, rightenv(envs, site + 1, below) ) end -function C_hamiltonian(site::Int, below::_HAM_MPS_TYPES, operator::MPOHamiltonian{<:MPOTensor}, above::_HAM_MPS_TYPES, envs) - return MPO_Contracted_C_Hamiltonian(leftenv(envs, site + 1, below), rightenv(envs, site, below)) -end -function AC_hamiltonian(site::Int, below::_HAM_MPS_TYPES, operator::MPOHamiltonian{<:MPOTensor}, above::_HAM_MPS_TYPES, envs) +function AC_hamiltonian(site::Int, below::_HAM_MPS_TYPES, operator::MPO{<:MPOTensor}, above::_HAM_MPS_TYPES, envs) O = operator[site] GL = leftenv(envs, site, below) return AC_hamiltonian(GL, O, rightenv(envs, site, below)) end function AC_hamiltonian(GL::MPSTensor, O::MPOTensor, GR::MPSTensor) @plansor GLW[-1 -2 -3; -4 -5] ≔ GL[-1 1; -4] * O[1 -2; -5 -3] - return MPO_Contracted_AC_Hamiltonian(GLW, O, GR) + return MPO_Contracted_AC_Hamiltonian(GLW, GR) end -function AC2_hamiltonian(site::Int, below::_HAM_MPS_TYPES, operator::MPOHamiltonian{<:MPOTensor}, above::_HAM_MPS_TYPES, envs) +function AC2_hamiltonian(site::Int, below::_HAM_MPS_TYPES, operator::MPO{<:MPOTensor}, above::_HAM_MPS_TYPES, envs) O1 = operator[site] O2 = operator[site + 1] GL = leftenv(envs, site, below) @@ -164,12 +158,8 @@ function (h::MPO_AC2_Hamiltonian{<:MPSTensor, <:MPOTensor, <:MPOTensor, <:MPSTen h.operators[2][7 -6; 4 5] * τ[5 -5; 2 3] return y isa AbstractBlockTensorMap ? only(y) : y end -function (h::MPO_Contracted_C_Hamiltonian)(x::MPSBondTensor) - @plansor y[-1; -2] ≔ h.leftenv[-1 3; 1] * x[1; 2] * h.rightenv[2 3; -2] - return y isa AbstractBlockTensorMap ? only(y) : y -end function (h::MPO_Contracted_AC_Hamiltonian)( - x::GenericMPSTensor{<:Any, 3} + x::AbstractTensorMap{<:Any, <:Any, 2, 1} ) @plansor y[-1 -2; -3] ≔ h.leftenv[-1 -2 3; 1 2] * x[1 2; 4] * h.rightenv[4 3; -3] return y isa AbstractBlockTensorMap ? only(y) : y @@ -179,5 +169,4 @@ function (h::MPO_Contracted_AC2_Hamiltonian)( ) @plansor y[-1 -2; -3 -4] ≔ h.leftenv[-1 -2 5; 1 2] * x[1 2; 3 4] * h.rightenv[3 4 5; -3 -4] return y isa AbstractBlockTensorMap ? only(y) : y -end -## TODO: Which interface should we use for an inplace += version to be used in hamiltonian_derivatives.jl? \ No newline at end of file +end \ No newline at end of file From 1f3ffc828c2d990169c6b3b07f929e858407192d Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Tue, 7 Oct 2025 18:54:39 -0400 Subject: [PATCH 10/19] restore type determination --- .../derivatives/hamiltonian_derivatives.jl | 57 ++++++------------- 1 file changed, 17 insertions(+), 40 deletions(-) diff --git a/src/algorithms/derivatives/hamiltonian_derivatives.jl b/src/algorithms/derivatives/hamiltonian_derivatives.jl index 21ff79c6b..8af9eae69 100644 --- a/src/algorithms/derivatives/hamiltonian_derivatives.jl +++ b/src/algorithms/derivatives/hamiltonian_derivatives.jl @@ -10,26 +10,7 @@ struct JordanMPO_AC_Hamiltonian{O1, O2, O3} <: DerivativeOperator E::Union{O1, Missing} # finished C::Union{O2, Missing} # starting B::Union{O2, Missing} # ending - A::O3 # continuing - function JordanMPO_AC_Hamiltonian( - onsite, not_started, finished, starting, ending, continuing - ) - # obtaining storagetype of environments since these should have already mixed - # the types of the operator and state - gl = @coalesce(onsite, not_started, finished, starting, ending) - S = spacetype(gl) - M = storagetype(gl) - O1 = tensormaptype(S, 1, 1, M) - O2 = tensormaptype(S, 2, 2, M) - return new{O1, O2, typeof(continuing)}( - onsite, not_started, finished, starting, ending, continuing - ) - end - function JordanMPO_AC_Hamiltonian{O1, O2, O3}( - onsite, not_started, finished, starting, ending, continuing - ) where {O1, O2, O3} - return new{O1, O2, O3}(onsite, not_started, finished, starting, ending, continuing) - end + A::Union{O3, Missing} # continuing end """ @@ -45,26 +26,10 @@ struct JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4} <: DerivativeOperator CB::Union{O2, Missing} # starting left - ending right CA::Union{O3, Missing} # starting left - continuing right AB::Union{O3, Missing} # continuing left - ending right - AA::O4 # continuing left - continuing right + AA::Union{O4, Missing} # continuing left - continuing right BE::Union{O2, Missing} # ending left DE::Union{O1, Missing} # onsite left EE::Union{O1, Missing} # finished - function JordanMPO_AC2_Hamiltonian(II, IC, ID, CB, CA, AB, AA, BE, DE, EE) - # obtaining storagetype of environments since these should have already mixed - # the types of the operator and state - gl = @coalesce(II, IC, ID, CB, CA, AB, BE, DE, EE) - S = spacetype(gl) - M = storagetype(gl) - O1 = tensormaptype(S, 1, 1, M) - O2 = tensormaptype(S, 2, 2, M) - O3 = tensormaptype(S, 3, 3, M) - return new{O1, O2, O3, typeof(AA)}(II, IC, ID, CB, CA, AB, AA, BE, DE, EE) - end - function JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4}( - II, IC, ID, CB, CA, AB, AA, BE, DE, EE - ) where {O1, O2, O3, O4} - return new{O1, O2, O3, O4}(II, IC, ID, CB, CA, AB, AA, BE, DE, EE) - end end # Constructors @@ -139,7 +104,13 @@ function AC_hamiltonian( continuing = missing end - return JordanMPO_AC_Hamiltonian( + S = spacetype(GL) + M = storagetype(GL) + O1 = tensormaptype(S, 1, 1, M) + O2 = tensormaptype(S, 2, 2, M) + O3 = Core.Compiler.return_type(AC_hamiltonian, typeof((GL_2, W.A, GR_2))) + + return JordanMPO_AC_Hamiltonian{O1, O2, O3}( onsite, not_started, finished, starting, ending, continuing ) end @@ -286,7 +257,6 @@ function AC2_hamiltonian( end # continuing - continuing - # TODO: MPODerivativeOperator code reuse + optimization ## TODO: Think about how one could and whether one should store these objects and use them for (a) advancing environments in iDMRG, (b) reuse ind backwards-sweep in IDMRG, (c) subspace expansion if nonzero_length(W1.A) > 0 && nonzero_length(W2.A) > 0 AA = AC2_hamiltonian(GL_2, W1.A, W2.A, GR_2) @@ -294,7 +264,14 @@ function AC2_hamiltonian( AA = missing end - return JordanMPO_AC2_Hamiltonian(II, IC, ID, CB, CA, AB, AA, BE, DE, EE) + S = spacetype(GL) + M = storagetype(GL) + O1 = tensormaptype(S, 1, 1, M) + O2 = tensormaptype(S, 2, 2, M) + O3 = tensormaptype(S, 3, 3, M) + O4 = Core.Compiler.return_type(AC2_hamiltonian, typeof((GL_2, W1.A, W2.A, GR_2))) + + return JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4}(II, IC, ID, CB, CA, AB, AA, BE, DE, EE) end # Actions From bca4f46e7dc34d86a3fac953aae51fb0c0fd59be Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Tue, 7 Oct 2025 19:19:35 -0400 Subject: [PATCH 11/19] fix Aqua complaints --- .../derivatives/hamiltonian_derivatives.jl | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/src/algorithms/derivatives/hamiltonian_derivatives.jl b/src/algorithms/derivatives/hamiltonian_derivatives.jl index 8af9eae69..1a611d1b9 100644 --- a/src/algorithms/derivatives/hamiltonian_derivatives.jl +++ b/src/algorithms/derivatives/hamiltonian_derivatives.jl @@ -11,6 +11,14 @@ struct JordanMPO_AC_Hamiltonian{O1, O2, O3} <: DerivativeOperator C::Union{O2, Missing} # starting B::Union{O2, Missing} # ending A::Union{O3, Missing} # continuing + + # need inner constructor to prohibit no-type-param constructor with unbound vars + function JordanMPO_AC_Hamiltonian{O1, O2, O3}( + D::Union{O1, Missing}, I::Union{O1, Missing}, E::Union{O1, Missing}, + C::Union{O2, Missing}, B::Union{O2, Missing}, A::Union{O3, Missing}, + ) where {O1, O2, O3} + return new{O1, O2, O3}(D, I, E, C, B, A) + end end """ @@ -30,6 +38,16 @@ struct JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4} <: DerivativeOperator BE::Union{O2, Missing} # ending left DE::Union{O1, Missing} # onsite left EE::Union{O1, Missing} # finished + + # need inner constructor to prohibit no-type-param constructor with unbound vars + function JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4}( + II::Union{O1, Missing}, IC::Union{O2, Missing}, ID::Union{O1, Missing}, + CB::Union{O2, Missing}, CA::Union{O3, Missing}, AB::Union{O3, Missing}, + AA::Union{O4, Missing}, BE::Union{O2, Missing}, DE::Union{O1, Missing}, + EE::Union{O1, Missing} + ) where {O1, O2, O3, O4} + return new{O1, O2, O3, O4}(II, IC, ID, CB, CA, AB, AA, BE, DE, EE) + end end # Constructors From e5216fa369b28384f4e6137c320f25daf7210f88 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Tue, 7 Oct 2025 19:19:35 -0400 Subject: [PATCH 12/19] fix Aqua complaints --- .../derivatives/hamiltonian_derivatives.jl | 20 +++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/src/algorithms/derivatives/hamiltonian_derivatives.jl b/src/algorithms/derivatives/hamiltonian_derivatives.jl index 8af9eae69..c1db0eb13 100644 --- a/src/algorithms/derivatives/hamiltonian_derivatives.jl +++ b/src/algorithms/derivatives/hamiltonian_derivatives.jl @@ -11,6 +11,14 @@ struct JordanMPO_AC_Hamiltonian{O1, O2, O3} <: DerivativeOperator C::Union{O2, Missing} # starting B::Union{O2, Missing} # ending A::Union{O3, Missing} # continuing + + # need inner constructor to prohibit no-type-param constructor with unbound vars + function JordanMPO_AC_Hamiltonian{O1, O2, O3}( + D::Union{O1, Missing}, I::Union{O1, Missing}, E::Union{O1, Missing}, + C::Union{O2, Missing}, B::Union{O2, Missing}, A::Union{O3, Missing}, + ) where {O1, O2, O3} + return new{O1, O2, O3}(D, I, E, C, B, A) + end end """ @@ -30,6 +38,16 @@ struct JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4} <: DerivativeOperator BE::Union{O2, Missing} # ending left DE::Union{O1, Missing} # onsite left EE::Union{O1, Missing} # finished + + # need inner constructor to prohibit no-type-param constructor with unbound vars + function JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4}( + II::Union{O1, Missing}, IC::Union{O2, Missing}, ID::Union{O1, Missing}, + CB::Union{O2, Missing}, CA::Union{O3, Missing}, AB::Union{O3, Missing}, + AA::Union{O4, Missing}, BE::Union{O2, Missing}, DE::Union{O1, Missing}, + EE::Union{O1, Missing} + ) where {O1, O2, O3, O4} + return new{O1, O2, O3, O4}(II, IC, ID, CB, CA, AB, AA, BE, DE, EE) + end end # Constructors @@ -172,6 +190,7 @@ function AC2_hamiltonian( IC = missing else @plansor IC[-1 -2; -3 -4] ≔ W2.C[-1; -3 1] * GR_2[-4 1; -2] + IC = only(IC) end else IC = missing @@ -186,6 +205,7 @@ function AC2_hamiltonian( BE = missing else @plansor BE[-1 -2; -3 -4] ≔ GL_2[-1 2; -3] * W1.B[2 -2; -4] + BE = only(BE) end else BE = missing From bab244a1d6d5aa26508e04c99ce768f094de4c56 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Fri, 10 Oct 2025 13:21:01 -0400 Subject: [PATCH 13/19] attempt to fix type stability --- .../derivatives/hamiltonian_derivatives.jl | 24 +++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/src/algorithms/derivatives/hamiltonian_derivatives.jl b/src/algorithms/derivatives/hamiltonian_derivatives.jl index c1db0eb13..801f9b6eb 100644 --- a/src/algorithms/derivatives/hamiltonian_derivatives.jl +++ b/src/algorithms/derivatives/hamiltonian_derivatives.jl @@ -20,6 +20,15 @@ struct JordanMPO_AC_Hamiltonian{O1, O2, O3} <: DerivativeOperator return new{O1, O2, O3}(D, I, E, C, B, A) end end +function JordanMPO_AC_Hamiltonian{O1, O2, O3}( + D, I, E, C, B, A, + ) where {O1, O2, O3} + return JordanMPO_AC_Hamiltonian{O1, O2, O3}( + ismissing(D) ? D : convert(O1, D), ismissing(I) ? I : convert(O1, I), + ismissing(E) ? E : convert(O1, E), ismissing(C) ? C : convert(O2, C), + ismissing(B) ? B : convert(O2, B), ismissing(A) ? A : convert(O3, A) + ) +end """ JordanMPO_AC2_Hamiltonian{O1,O2,O3,O4} @@ -50,6 +59,21 @@ struct JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4} <: DerivativeOperator end end +function JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4}( + II::Union{O1, Missing}, IC::Union{O2, Missing}, ID::Union{O1, Missing}, + CB::Union{O2, Missing}, CA::Union{O3, Missing}, AB::Union{O3, Missing}, + AA::Union{O4, Missing}, BE::Union{O2, Missing}, DE::Union{O1, Missing}, + EE::Union{O1, Missing} + ) where {O1, O2, O3, O4} + return JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4}( + ismissing(II) ? II : convert(O1, II), ismissing(IC) ? IC : convert(O2, IC), + ismissing(ID) ? ID : convert(O1, II), ismissing(CB) ? CB : convert(O2, CB), + ismissing(CA) ? CA : convert(O3, CA), ismissing(AB) ? AB : convert(O3, AB), + ismissing(AA) ? AA : convert(P4, AA), ismissing(BE) ? BE : convert(O2, BE), + ismissing(DE) ? DE : convert(O1, DE), ismissing(EE) ? EE : convert(O1, EE) + ) +end + # Constructors # ------------ function AC_hamiltonian( From edc66455904d057218b85183cd31f6619e73d7ae Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Fri, 10 Oct 2025 14:34:06 -0400 Subject: [PATCH 14/19] small fix --- src/algorithms/derivatives/hamiltonian_derivatives.jl | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/algorithms/derivatives/hamiltonian_derivatives.jl b/src/algorithms/derivatives/hamiltonian_derivatives.jl index 801f9b6eb..42a92eee6 100644 --- a/src/algorithms/derivatives/hamiltonian_derivatives.jl +++ b/src/algorithms/derivatives/hamiltonian_derivatives.jl @@ -60,10 +60,7 @@ struct JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4} <: DerivativeOperator end function JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4}( - II::Union{O1, Missing}, IC::Union{O2, Missing}, ID::Union{O1, Missing}, - CB::Union{O2, Missing}, CA::Union{O3, Missing}, AB::Union{O3, Missing}, - AA::Union{O4, Missing}, BE::Union{O2, Missing}, DE::Union{O1, Missing}, - EE::Union{O1, Missing} + II, IC, ID, CB, CA, AB, AA, BE, DE, EE ) where {O1, O2, O3, O4} return JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4}( ismissing(II) ? II : convert(O1, II), ismissing(IC) ? IC : convert(O2, IC), From 1688b7bb708c3d7686f16432d7da1b0f7bf95c69 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Fri, 10 Oct 2025 16:06:26 -0400 Subject: [PATCH 15/19] type stability improvements --- .../derivatives/hamiltonian_derivatives.jl | 29 ++----------------- 1 file changed, 2 insertions(+), 27 deletions(-) diff --git a/src/algorithms/derivatives/hamiltonian_derivatives.jl b/src/algorithms/derivatives/hamiltonian_derivatives.jl index 42a92eee6..b35bf85d2 100644 --- a/src/algorithms/derivatives/hamiltonian_derivatives.jl +++ b/src/algorithms/derivatives/hamiltonian_derivatives.jl @@ -14,21 +14,11 @@ struct JordanMPO_AC_Hamiltonian{O1, O2, O3} <: DerivativeOperator # need inner constructor to prohibit no-type-param constructor with unbound vars function JordanMPO_AC_Hamiltonian{O1, O2, O3}( - D::Union{O1, Missing}, I::Union{O1, Missing}, E::Union{O1, Missing}, - C::Union{O2, Missing}, B::Union{O2, Missing}, A::Union{O3, Missing}, + D, I, E, C, B, A, ) where {O1, O2, O3} return new{O1, O2, O3}(D, I, E, C, B, A) end end -function JordanMPO_AC_Hamiltonian{O1, O2, O3}( - D, I, E, C, B, A, - ) where {O1, O2, O3} - return JordanMPO_AC_Hamiltonian{O1, O2, O3}( - ismissing(D) ? D : convert(O1, D), ismissing(I) ? I : convert(O1, I), - ismissing(E) ? E : convert(O1, E), ismissing(C) ? C : convert(O2, C), - ismissing(B) ? B : convert(O2, B), ismissing(A) ? A : convert(O3, A) - ) -end """ JordanMPO_AC2_Hamiltonian{O1,O2,O3,O4} @@ -50,27 +40,12 @@ struct JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4} <: DerivativeOperator # need inner constructor to prohibit no-type-param constructor with unbound vars function JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4}( - II::Union{O1, Missing}, IC::Union{O2, Missing}, ID::Union{O1, Missing}, - CB::Union{O2, Missing}, CA::Union{O3, Missing}, AB::Union{O3, Missing}, - AA::Union{O4, Missing}, BE::Union{O2, Missing}, DE::Union{O1, Missing}, - EE::Union{O1, Missing} + II, IC, ID, CB, CA, AB, AA, BE, DE, EE ) where {O1, O2, O3, O4} return new{O1, O2, O3, O4}(II, IC, ID, CB, CA, AB, AA, BE, DE, EE) end end -function JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4}( - II, IC, ID, CB, CA, AB, AA, BE, DE, EE - ) where {O1, O2, O3, O4} - return JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4}( - ismissing(II) ? II : convert(O1, II), ismissing(IC) ? IC : convert(O2, IC), - ismissing(ID) ? ID : convert(O1, II), ismissing(CB) ? CB : convert(O2, CB), - ismissing(CA) ? CA : convert(O3, CA), ismissing(AB) ? AB : convert(O3, AB), - ismissing(AA) ? AA : convert(P4, AA), ismissing(BE) ? BE : convert(O2, BE), - ismissing(DE) ? DE : convert(O1, DE), ismissing(EE) ? EE : convert(O1, EE) - ) -end - # Constructors # ------------ function AC_hamiltonian( From 03695358ea0b0be0b9236f3be60bd4846753adc9 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Fri, 10 Oct 2025 16:06:32 -0400 Subject: [PATCH 16/19] fix deprecation warnings --- src/algorithms/statmech/idmrg.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/algorithms/statmech/idmrg.jl b/src/algorithms/statmech/idmrg.jl index 5aa7b0b74..4f4204733 100644 --- a/src/algorithms/statmech/idmrg.jl +++ b/src/algorithms/statmech/idmrg.jl @@ -19,7 +19,7 @@ function leading_boundary( for row in 1:size(ψ, 1) ac = ψ.AC[row, col] (col == size(ψ, 2)) && (ac = copy(ac)) # needed in next sweep - ψ.AL[row, col], ψ.C[row, col] = leftorth!(ac) + ψ.AL[row, col], ψ.C[row, col] = qr_compact!(ac) end transfer_leftenv!(envs, ψ, operator, ψ, col + 1) @@ -31,7 +31,7 @@ function leading_boundary( _, ψ.AC[:, col] = fixedpoint(Hac, ψ.AC[:, col], :LM, alg_eigsolve) for row in 1:size(ψ, 1) - ψ.C[row, col - 1], temp = rightorth!(_transpose_tail(ψ.AC[row, col])) + ψ.C[row, col - 1], temp = lq_compact!(_transpose_tail(ψ.AC[row, col])) ψ.AR[row, col] = _transpose_front(temp) end From 4fc23f5342073a1f2c18e462cf6eb38a9d965326 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Fri, 10 Oct 2025 16:06:36 -0400 Subject: [PATCH 17/19] stabilize test --- test/algorithms.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/algorithms.jl b/test/algorithms.jl index c1d0ad9e9..c91be0656 100644 --- a/test/algorithms.jl +++ b/test/algorithms.jl @@ -190,7 +190,7 @@ module TestAlgorithms @testset "LazySum FiniteMPS groundstate" verbose = true begin tol = 1.0e-8 - D = 15 + D = 32 atol = 1.0e-2 L = 10 From 69502678ffab8b6ec4992119d8851d690e4d2fc0 Mon Sep 17 00:00:00 2001 From: afeuerpfeil Date: Thu, 30 Oct 2025 17:58:18 -0400 Subject: [PATCH 18/19] sync --- Project.toml | 4 +- src/MPSKit.jl | 4 +- src/algorithms/changebonds/changebonds.jl | 2 +- src/algorithms/changebonds/randexpand.jl | 134 +++++++- .../hamiltonian_derivatives copy.jl | 319 ++++++++++++++++++ .../derivatives/hamiltonian_derivatives.jl | 299 +--------------- src/algorithms/derivatives/mpo_derivatives.jl | 6 +- src/algorithms/groundstate/idmrg.jl | 4 +- src/algorithms/groundstate/vumps.jl | 2 +- src/environments/abstract_envs.jl | 5 +- src/states/abstractmps.jl | 38 ++- src/states/infinitemps.jl | 5 +- src/states/ortho.jl | 2 +- src/utility/dynamictols.jl | 5 +- 14 files changed, 490 insertions(+), 339 deletions(-) create mode 100644 src/algorithms/derivatives/hamiltonian_derivatives copy.jl diff --git a/Project.toml b/Project.toml index 23fc5e15f..f97025cba 100644 --- a/Project.toml +++ b/Project.toml @@ -26,7 +26,7 @@ VectorInterface = "409d34a3-91d5-4945-b6ec-7529ddf182d8" [compat] Accessors = "0.1" Aqua = "0.8.9" -BlockTensorKit = "0.3" +BlockTensorKit = "0.2,0.3" Combinatorics = "1" Compat = "3.47, 4.10" DocStringExtensions = "0.9.3" @@ -42,7 +42,7 @@ Plots = "1.40" Printf = "1" Random = "1" RecipesBase = "1.1" -TensorKit = "0.15.1" +TensorKit = "0.15.1,0.14" TensorKitManifolds = "0.7" TensorOperations = "5" Test = "1" diff --git a/src/MPSKit.jl b/src/MPSKit.jl index a6c586144..1d09467fa 100644 --- a/src/MPSKit.jl +++ b/src/MPSKit.jl @@ -36,7 +36,7 @@ export FiniteExcited, QuasiparticleAnsatz, ChepigaAnsatz, ChepigaAnsatz2 export time_evolve, timestep, timestep!, make_time_mpo export TDVP, TDVP2, WI, WII, TaylorCluster export changebonds, changebonds! -export VUMPSSvdCut, OptimalExpand, SvdCut, RandExpand +export VUMPSSvdCut, OptimalExpand, SvdCut, RandExpand, RandPerturbedExpand export propagator export DynamicalDMRG, NaiveInvert, Jeckelmann export exact_diagonalization, fidelity_susceptibility @@ -194,4 +194,4 @@ function __init__() return nothing end -end +end \ No newline at end of file diff --git a/src/algorithms/changebonds/changebonds.jl b/src/algorithms/changebonds/changebonds.jl index 700b5ab0b..d20e7399c 100644 --- a/src/algorithms/changebonds/changebonds.jl +++ b/src/algorithms/changebonds/changebonds.jl @@ -51,4 +51,4 @@ function _expand!(ψ::MultilineMPS, AL′::PeriodicMatrix, AR′::PeriodicMatrix _expand!(ψ[i], AL′[i, :], AR′[i, :]) end return ψ -end +end \ No newline at end of file diff --git a/src/algorithms/changebonds/randexpand.jl b/src/algorithms/changebonds/randexpand.jl index 28e87b3d6..73ed83f2a 100644 --- a/src/algorithms/changebonds/randexpand.jl +++ b/src/algorithms/changebonds/randexpand.jl @@ -2,8 +2,9 @@ $(TYPEDEF) An algorithm that expands the bond dimension by adding random unitary vectors that are -orthogonal to the existing state. This is achieved by performing a truncated SVD on a random -two-site MPS tensor, which is made orthogonal to the existing state. +orthogonal to the existing state. This means that additional directions are added to +`AL` and `AR` that are contained in the nullspace of both. Note that this is happens in +parallel, and therefore the expansion will never go beyond the local two-site subspace. ## Fields @@ -17,39 +18,42 @@ $(TYPEDFIELDS) trscheme::TruncationStrategy end -function changebonds(ψ::InfiniteMPS, alg::RandExpand) +function changebonds!(ψ::InfiniteMPS, alg::RandExpand) T = eltype(ψ.AL) AL′ = similar(ψ.AL) AR′ = similar(ψ.AR, tensormaptype(spacetype(T), 1, numind(T) - 1, storagetype(T))) - for i in 1:length(ψ) - # determine optimal expansion spaces around bond i - AC2 = randomize!(_transpose_front(ψ.AC[i]) * _transpose_tail(ψ.AR[i + 1])) - # Use the nullspaces and SVD decomposition to determine the optimal expansion space + for i in 1:length(ψ) + # obtain spaces by sampling the support of both the left and right nullspace VL = left_null(ψ.AL[i]) - VR = right_null!(_transpose_tail(ψ.AR[i + 1])) - intermediate = normalize!(adjoint(VL) * AC2 * adjoint(VR)) - U, _, Vᴴ = svd_trunc!(intermediate; trunc = alg.trscheme, alg = alg.alg_svd) - - AL′[i] = VL * U - AR′[i + 1] = Vᴴ * VR + VR = right_null!(_transpose_tail(ψ.AR[i + 1]; copy = true)) + V = sample_space(infimum(right_virtualspace(VL), space(VR, 1)), alg.trscheme) + + # obtain (orthogonal) directions as isometries in that direction + XL = randisometry(scalartype(VL), right_virtualspace(VL) ← V) + AL′[i] = VL * XL + XR = randisometry(storagetype(VL), space(VR, 1) ← V) + AR′[i + 1] = XR * VR end - return _expand(ψ, AL′, AR′) + return _expand!(ψ, AL′, AR′) end -function changebonds(ψ::MultilineMPS, alg::RandExpand) - return Multiline(map(x -> changebonds(x, alg), ψ.data)) +function changebonds!(ψ::MultilineMPS, alg::RandExpand) + foreach(Base.Fix2(changebonds!, alg), ψ.data) + return ψ end -changebonds(ψ::AbstractFiniteMPS, alg::RandExpand) = changebonds!(copy(ψ), alg) +changebonds(ψ::AbstractMPS, alg::RandExpand) = changebonds!(copy(ψ), alg) +changebonds(ψ::MultilineMPS, alg::RandExpand) = changebonds!(copy(ψ), alg) + function changebonds!(ψ::AbstractFiniteMPS, alg::RandExpand) for i in 1:(length(ψ) - 1) AC2 = randomize!(_transpose_front(ψ.AC[i]) * _transpose_tail(ψ.AR[i + 1])) #Calculate nullspaces for left and right NL = left_null(ψ.AC[i]) - NR = right_null!(_transpose_tail(ψ.AR[i + 1])) + NR = right_null!(_transpose_tail(ψ.AR[i + 1]; copy = true)) #Use this nullspaces and SVD decomposition to determine the optimal expansion space intermediate = normalize!(adjoint(NL) * AC2 * adjoint(NR)) @@ -67,3 +71,97 @@ function changebonds!(ψ::AbstractFiniteMPS, alg::RandExpand) return normalize!(ψ) end + +function sample_space(V, strategy) + S = TensorKit.SectorDict(c => sort!(Random.rand(dim(V, c)); rev = true) for c in sectors(V)) + ind = MatrixAlgebraKit.findtruncated_svd(S, strategy) + return TensorKit.Factorizations.truncate_space(V, ind) +end + + +""" +$(TYPEDEF) + +An algorithm that expands the bond dimension by adding random unitary vectors that are +orthogonal to the existing state, in a sweeping fashion. Additionally, some random noise +is added to the state in order for it to remain gauge-able. + +## Fields + +$(TYPEDFIELDS) +""" +@kwdef struct RandPerturbedExpand{S} <: Algorithm + "algorithm used for the singular value decomposition" + alg_svd::S = Defaults.alg_svd() + + "algorithm used for [truncation](@extref MatrixAlgebraKit.TruncationStrategy] the expanded space" + trscheme::TruncationStrategy + + "amount of noise that is added to the current state" + noisefactor::Float64 = eps()^(3 / 4) + + "algorithm used for gauging the state" + alg_gauge = Defaults.alg_gauge(; dynamic_tols = false) +end + +function changebonds!(ψ::InfiniteMPS, alg::RandPerturbedExpand) + for i in 1:length(ψ) + # obtain space by sampling the support of left nullspace + # add (orthogonal) directions as isometries in that direction + VL = left_null(ψ.AL[i]) + V = sample_space(right_virtualspace(VL), alg.trscheme) + XL = randisometry(scalartype(VL), right_virtualspace(VL) ← V) + ψ.AL[i] = catdomain(ψ.AL[i], VL * XL) + + # make sure the next site fits, by "absorbing" into a larger tensor + # with some random noise to ensure state is still gauge-able + AL = ψ.AL[i + 1] + AL′ = similar(AL, right_virtualspace(ψ.AL[i]) ⊗ physicalspace(AL) ← right_virtualspace(AL)) + scale!(randn!(AL′), alg.noisefactor) + ψ.AL[i + 1] = TensorKit.absorb!(AL′, AL) + end + + # properly regauge the state: + makefullrank!(ψ.AL) + ψ.AR .= similar.(ψ.AL) + # ψ.AC .= similar.(ψ.AL) + + # initial guess for gauge is embedded original C + C₀ = similar(ψ.C[0], right_virtualspace(ψ.AL[end]) ← left_virtualspace(ψ.AL[1])) + absorb!(id!(C₀), ψ.C[0]) + + gaugefix!(ψ, ψ.AL, C₀; order = :R, alg.alg_gauge.maxiter, alg.alg_gauge.tol) + + for i in reverse(1:length(ψ)) + # obtain space by sampling the support of left nullspace + # add (orthogonal) directions as isometries in that direction + AR_tail = _transpose_tail(ψ.AR[i]) + VR = right_null(AR_tail) + V = sample_space(space(VR, 1), alg.trscheme) + XR = randisometry(scalartype(VR), space(VR, 1) ← V) + ψ.AR[i] = _transpose_front(catcodomain(AR_tail, XR' * VR)) + + # make sure the next site fits, by "absorbing" into a larger tensor + # with some random noise to ensure state is still gauge-able + AR = ψ.AR[i - 1] + AR′ = similar(AR, left_virtualspace(AR) ⊗ physicalspace(AR) ← left_virtualspace(ψ.AR[i])) + scale!(randn!(AR′), alg.noisefactor) + ψ.AR[i - 1] = TensorKit.absorb!(AR′, AR) + end + + # properly regauge the state: + makefullrank!(ψ.AR) + ψ.AL .= similar.(ψ.AR) + ψ.AC .= similar.(ψ.AR) + + # initial guess for gauge is embedded original C + C₀ = similar(ψ.C[0], right_virtualspace(ψ.AR[end]) ← left_virtualspace(ψ.AR[1])) + absorb!(id!(C₀), ψ.C[0]) + + gaugefix!(ψ, ψ.AR, C₀; order = :LR, alg.alg_gauge.maxiter, alg.alg_gauge.tol) + mul!.(ψ.AC, ψ.AL, ψ.C) + + return ψ +end + +changebonds(ψ::InfiniteMPS, alg::RandPerturbedExpand) = changebonds!(copy(ψ), alg) \ No newline at end of file diff --git a/src/algorithms/derivatives/hamiltonian_derivatives copy.jl b/src/algorithms/derivatives/hamiltonian_derivatives copy.jl new file mode 100644 index 000000000..f5290ebff --- /dev/null +++ b/src/algorithms/derivatives/hamiltonian_derivatives copy.jl @@ -0,0 +1,319 @@ +""" + JordanMPO_AC_Hamiltonian{O1,O2,O3} + +Efficient operator for representing the single-site derivative of a `MPOHamiltonian` sandwiched between two MPSs. +In particular, this operator aims to make maximal use of the structure of the `MPOHamiltonian` to reduce the number of operations required to apply the operator to a tensor. +""" +struct JordanMPO_AC_Hamiltonian{O1, O2, O3} <: DerivativeOperator + D::Union{O1, Missing} # onsite + I::Union{O1, Missing} # not started + E::Union{O1, Missing} # finished + C::Union{O2, Missing} # starting + B::Union{O2, Missing} # ending + A::Union{O3, Missing} # continuing + + # need inner constructor to prohibit no-type-param constructor with unbound vars + function JordanMPO_AC_Hamiltonian{O1, O2, O3}( + D, I, E, C, B, A, + ) where {O1, O2, O3} + return new{O1, O2, O3}(D, I, E, C, B, A) + end +end + +""" + JordanMPO_AC2_Hamiltonian{O1,O2,O3,O4} + +Efficient operator for representing the single-site derivative of a `MPOHamiltonian` sandwiched between two MPSs. +In particular, this operator aims to make maximal use of the structure of the `MPOHamiltonian` to reduce the number of operations required to apply the operator to a tensor. +""" +struct JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4} <: DerivativeOperator + II::Union{O1, Missing} # not_started + IC::Union{O2, Missing} # starting right + ID::Union{O1, Missing} # onsite right + CB::Union{O2, Missing} # starting left - ending right + CA::Union{O3, Missing} # starting left - continuing right + AB::Union{O3, Missing} # continuing left - ending right + AA::Union{O4, Missing} # continuing left - continuing right + BE::Union{O2, Missing} # ending left + DE::Union{O1, Missing} # onsite left + EE::Union{O1, Missing} # finished + + # need inner constructor to prohibit no-type-param constructor with unbound vars + function JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4}( + II, IC, ID, CB, CA, AB, AA, BE, DE, EE + ) where {O1, O2, O3, O4} + return new{O1, O2, O3, O4}(II, IC, ID, CB, CA, AB, AA, BE, DE, EE) + end +end + +# Constructors +# ------------ +function AC_hamiltonian( + site::Int, below::_HAM_MPS_TYPES, operator::MPOHamiltonian{<:JordanMPOTensor}, + above::_HAM_MPS_TYPES, envs + ) + GL = leftenv(envs, site, below) + GR = rightenv(envs, site, below) + W = operator[site] + + GR_2 = GR[2:(end - 1)] + GL_2 = GL[2:(end - 1)] + + # starting + if nonzero_length(W.C) > 0 + @plansor starting_[-1 -2; -3 -4] ≔ W.C[-1; -3 1] * GR_2[-4 1; -2] + starting = only(starting_) + else + starting = missing + end + + # ending + if nonzero_length(W.B) > 0 + @plansor ending_[-1 -2; -3 -4] ≔ GL_2[-1 1; -3] * W.B[1 -2; -4] + ending = only(ending_) + else + ending = missing + end + + # onsite + if nonzero_length(W.D) > 0 + if !ismissing(starting) + @plansor starting[-1 -2; -3 -4] += W.D[-1; -3] * removeunit(GR[end], 2)[-4; -2] + onsite = missing + elseif !ismissing(ending) + @plansor ending[-1 -2; -3 -4] += removeunit(GL[1], 2)[-1; -3] * W.D[-2; -4] + onsite = missing + else + onsite = W.D + end + else + onsite = missing + end + + # not_started + if isfinite(operator) && site == length(operator) + not_started = missing + elseif !ismissing(starting) + I = id(storagetype(GR[1]), physicalspace(W)) + @plansor starting[-1 -2; -3 -4] += I[-1; -3] * removeunit(GR[1], 2)[-4; -2] + not_started = missing + else + not_started = removeunit(GR[1], 2) + end + + # finished + if isfinite(operator) && site == 1 + finished = missing + elseif !ismissing(ending) + I = id(storagetype(GL[end]), physicalspace(W)) + @plansor ending[-1 -2; -3 -4] += removeunit(GL[end], 2)[-1; -3] * I[-2; -4] + finished = missing + else + finished = removeunit(GL[end], 2) + end + + if nonzero_length(W.A) > 0 + continuing = AC_hamiltonian(GL_2, W.A, GR_2) + else + continuing = missing + end + + S = spacetype(GL) + M = storagetype(GL) + O1 = tensormaptype(S, 1, 1, M) + O2 = tensormaptype(S, 2, 2, M) + O3 = Core.Compiler.return_type(AC_hamiltonian, typeof((GL_2, W.A, GR_2))) + + return JordanMPO_AC_Hamiltonian{O1, O2, O3}( + onsite, not_started, finished, starting, ending, continuing + ) +end + +function AC2_hamiltonian( + site::Int, below::_HAM_MPS_TYPES, operator::MPOHamiltonian{<:JordanMPOTensor}, + above::_HAM_MPS_TYPES, envs + ) + GL = leftenv(envs, site, below) + GR = rightenv(envs, site + 1, below) + W1 = operator[site] + W2 = operator[site + 1] + + GR_2 = GR[2:(end - 1)] + GL_2 = GL[2:(end - 1)] + + # starting left - continuing right + if nonzero_length(W1.C) > 0 && nonzero_length(W2.A) > 0 + @plansor CA_[-1 -2 -3; -4 -5 -6] ≔ W1.C[-1; -4 2] * W2.A[2 -2; -5 1] * + GR_2[-6 1; -3] + CA = only(CA_) + else + CA = missing + end + + # continuing left - ending right + if nonzero_length(W1.A) > 0 && nonzero_length(W2.B) > 0 + @plansor AB_[-1 -2 -3; -4 -5 -6] ≔ GL_2[-1 2; -4] * W1.A[2 -2; -5 1] * + W2.B[1 -3; -6] + AB = only(AB_) + else + AB = missing + end + + # middle + if nonzero_length(W1.C) > 0 && nonzero_length(W2.B) > 0 + if !ismissing(CA) + @plansor CA[-1 -2 -3; -4 -5 -6] += W1.C[-1; -4 1] * W2.B[1 -2; -5] * + removeunit(GR[end], 2)[-6; -3] + CB = missing + elseif !ismissing(AB) + @plansor AB[-1 -2 -3; -4 -5 -6] += removeunit(GL[1], 2)[-1; -4] * + W1.C[-2; -5 1] * W2.B[1 -3; -6] + CB = missing + else + @plansor CB_[-1 -2; -3 -4] ≔ W1.C[-1; -3 1] * W2.B[1 -2; -4] + CB = only(CB_) + end + else + CB = missing + end + + # starting right + if nonzero_length(W2.C) > 0 + if !ismissing(CA) + I = id(storagetype(GR[1]), physicalspace(W1)) + @plansor CA[-1 -2 -3; -4 -5 -6] += (I[-1; -4] * W2.C[-2; -5 1]) * + GR_2[-6 1; -3] + IC = missing + else + @plansor IC[-1 -2; -3 -4] ≔ W2.C[-1; -3 1] * GR_2[-4 1; -2] + IC = only(IC) + end + else + IC = missing + end + + # ending left + if nonzero_length(W1.B) > 0 + if !ismissing(AB) + I = id(storagetype(GL[end]), physicalspace(W2)) + @plansor AB[-1 -2 -3; -4 -5 -6] += GL_2[-1 1; -4] * + (W1.B[1 -2; -5] * I[-3; -6]) + BE = missing + else + @plansor BE[-1 -2; -3 -4] ≔ GL_2[-1 2; -3] * W1.B[2 -2; -4] + BE = only(BE) + end + else + BE = missing + end + + # onsite left + if nonzero_length(W1.D) > 0 + if !ismissing(BE) + @plansor BE[-1 -2; -3 -4] += removeunit(GL[1], 2)[-1; -3] * W1.D[-2; -4] + DE = missing + elseif !ismissing(AB) + I = id(storagetype(GL[end]), physicalspace(W2)) + @plansor AB[-1 -2 -3; -4 -5 -6] += removeunit(GL[1], 2)[-1; -4] * + (W1.D[-2; -5] * I[-3; -6]) + DE = missing + # TODO: could also try in CA? + else + DE = only(W1.D) + end + else + DE = missing + end + + # onsite right + if nonzero_length(W2.D) > 0 + if !ismissing(IC) + @plansor IC[-1 -2; -3 -4] += W2.D[-1; -3] * removeunit(GR[end], 2)[-4; -2] + ID = missing + elseif !ismissing(CA) + I = id(storagetype(GR[1]), physicalspace(W1)) + @plansor CA[-1 -2 -3; -4 -5 -6] += (I[-1; -4] * W2.D[-2; -5]) * + removeunit(GR[end], 2)[-6; -3] + ID = missing + else + ID = only(W2.D) + end + else + ID = missing + end + + # finished + if isfinite(operator) && site + 1 == length(operator) + II = missing + elseif !ismissing(IC) + I = id(storagetype(GR[1]), physicalspace(W2)) + @plansor IC[-1 -2; -3 -4] += I[-1; -3] * removeunit(GR[1], 2)[-4; -2] + II = missing + elseif !ismissing(CA) + I = id(storagetype(GR[1]), physicalspace(W1) ⊗ physicalspace(W2)) + @plansor CA[-1 -2 -3; -4 -5 -6] += I[-1 -2; -4 -5] * removeunit(GR[1], 2)[-6; -3] + II = missing + else + II = transpose(removeunit(GR[1], 2)) + end + + # unstarted + if isfinite(operator) && site == 1 + EE = missing + elseif !ismissing(BE) + I = id(storagetype(GL[end]), physicalspace(W1)) + @plansor BE[-1 -2; -3 -4] += removeunit(GL[end], 2)[-1; -3] * I[-2; -4] + EE = missing + elseif !ismissing(AB) + I = id(storagetype(GL[end]), physicalspace(W1) ⊗ physicalspace(W2)) + @plansor AB[-1 -2 -3; -4 -5 -6] += removeunit(GL[end], 2)[-1; -4] * I[-2 -3; -5 -6] + EE = missing + else + EE = removeunit(GL[end], 2) + end + + # continuing - continuing + ## TODO: Think about how one could and whether one should store these objects and use them for (a) advancing environments in iDMRG, (b) reuse ind backwards-sweep in IDMRG, (c) subspace expansion + if nonzero_length(W1.A) > 0 && nonzero_length(W2.A) > 0 + AA = AC2_hamiltonian(GL_2, W1.A, W2.A, GR_2) + else + AA = missing + end + + S = spacetype(GL) + M = storagetype(GL) + O1 = tensormaptype(S, 1, 1, M) + O2 = tensormaptype(S, 2, 2, M) + O3 = tensormaptype(S, 3, 3, M) + O4 = Core.Compiler.return_type(AC2_hamiltonian, typeof((GL_2, W1.A, W2.A, GR_2))) + + return JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4}(II, IC, ID, CB, CA, AB, AA, BE, DE, EE) +end + +# Actions +# ------- +function (H::JordanMPO_AC_Hamiltonian)(x::MPSTensor) + y = ismissing(H.A) ? zerovector(x) : H.A(x) + ismissing(H.D) || @plansor y[-1 -2; -3] += x[-1 1; -3] * H.D[-2; 1] + ismissing(H.E) || @plansor y[-1 -2; -3] += H.E[-1; 1] * x[1 -2; -3] + ismissing(H.I) || @plansor y[-1 -2; -3] += x[-1 -2; 1] * H.I[1; -3] + ismissing(H.C) || @plansor y[-1 -2; -3] += x[-1 2; 1] * H.C[-2 -3; 2 1] + ismissing(H.B) || @plansor y[-1 -2; -3] += H.B[-1 -2; 1 2] * x[1 2; -3] + + return y +end + +function (H::JordanMPO_AC2_Hamiltonian)(x::MPOTensor) + y = ismissing(H.AA) ? zerovector(x) : H.AA(x) + ismissing(H.II) || @plansor y[-1 -2; -3 -4] += x[-1 -2; 1 -4] * H.II[-3; 1] + ismissing(H.IC) || @plansor y[-1 -2; -3 -4] += x[-1 -2; 1 2] * H.IC[-4 -3; 2 1] + ismissing(H.ID) || @plansor y[-1 -2; -3 -4] += x[-1 -2; -3 1] * H.ID[-4; 1] + ismissing(H.CB) || @plansor y[-1 -2; -3 -4] += x[-1 1; -3 2] * H.CB[-2 -4; 1 2] + ismissing(H.CA) || @plansor y[-1 -2; -3 -4] += x[-1 1; 3 2] * H.CA[-2 -4 -3; 1 2 3] + ismissing(H.AB) || @plansor y[-1 -2; -3 -4] += x[1 2; -3 3] * H.AB[-1 -2 -4; 1 2 3] + ismissing(H.BE) || @plansor y[-1 -2; -3 -4] += x[1 2; -3 -4] * H.BE[-1 -2; 1 2] + ismissing(H.DE) || @plansor y[-1 -2; -3 -4] += x[-1 1; -3 -4] * H.DE[-2; 1] + ismissing(H.EE) || @plansor y[-1 -2; -3 -4] += x[1 -2; -3 -4] * H.EE[-1; 1] + + return y +end \ No newline at end of file diff --git a/src/algorithms/derivatives/hamiltonian_derivatives.jl b/src/algorithms/derivatives/hamiltonian_derivatives.jl index b35bf85d2..009bdd069 100644 --- a/src/algorithms/derivatives/hamiltonian_derivatives.jl +++ b/src/algorithms/derivatives/hamiltonian_derivatives.jl @@ -1,50 +1,3 @@ -""" - JordanMPO_AC_Hamiltonian{O1,O2,O3} - -Efficient operator for representing the single-site derivative of a `MPOHamiltonian` sandwiched between two MPSs. -In particular, this operator aims to make maximal use of the structure of the `MPOHamiltonian` to reduce the number of operations required to apply the operator to a tensor. -""" -struct JordanMPO_AC_Hamiltonian{O1, O2, O3} <: DerivativeOperator - D::Union{O1, Missing} # onsite - I::Union{O1, Missing} # not started - E::Union{O1, Missing} # finished - C::Union{O2, Missing} # starting - B::Union{O2, Missing} # ending - A::Union{O3, Missing} # continuing - - # need inner constructor to prohibit no-type-param constructor with unbound vars - function JordanMPO_AC_Hamiltonian{O1, O2, O3}( - D, I, E, C, B, A, - ) where {O1, O2, O3} - return new{O1, O2, O3}(D, I, E, C, B, A) - end -end - -""" - JordanMPO_AC2_Hamiltonian{O1,O2,O3,O4} - -Efficient operator for representing the single-site derivative of a `MPOHamiltonian` sandwiched between two MPSs. -In particular, this operator aims to make maximal use of the structure of the `MPOHamiltonian` to reduce the number of operations required to apply the operator to a tensor. -""" -struct JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4} <: DerivativeOperator - II::Union{O1, Missing} # not_started - IC::Union{O2, Missing} # starting right - ID::Union{O1, Missing} # onsite right - CB::Union{O2, Missing} # starting left - ending right - CA::Union{O3, Missing} # starting left - continuing right - AB::Union{O3, Missing} # continuing left - ending right - AA::Union{O4, Missing} # continuing left - continuing right - BE::Union{O2, Missing} # ending left - DE::Union{O1, Missing} # onsite left - EE::Union{O1, Missing} # finished - - # need inner constructor to prohibit no-type-param constructor with unbound vars - function JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4}( - II, IC, ID, CB, CA, AB, AA, BE, DE, EE - ) where {O1, O2, O3, O4} - return new{O1, O2, O3, O4}(II, IC, ID, CB, CA, AB, AA, BE, DE, EE) - end -end # Constructors # ------------ @@ -56,77 +9,7 @@ function AC_hamiltonian( GR = rightenv(envs, site, below) W = operator[site] - GR_2 = GR[2:(end - 1)] - GL_2 = GL[2:(end - 1)] - - # starting - if nonzero_length(W.C) > 0 - @plansor starting_[-1 -2; -3 -4] ≔ W.C[-1; -3 1] * GR_2[-4 1; -2] - starting = only(starting_) - else - starting = missing - end - - # ending - if nonzero_length(W.B) > 0 - @plansor ending_[-1 -2; -3 -4] ≔ GL_2[-1 1; -3] * W.B[1 -2; -4] - ending = only(ending_) - else - ending = missing - end - - # onsite - if nonzero_length(W.D) > 0 - if !ismissing(starting) - @plansor starting[-1 -2; -3 -4] += W.D[-1; -3] * removeunit(GR[end], 2)[-4; -2] - onsite = missing - elseif !ismissing(ending) - @plansor ending[-1 -2; -3 -4] += removeunit(GL[1], 2)[-1; -3] * W.D[-2; -4] - onsite = missing - else - onsite = W.D - end - else - onsite = missing - end - - # not_started - if isfinite(operator) && site == length(operator) - not_started = missing - elseif !ismissing(starting) - I = id(storagetype(GR[1]), physicalspace(W)) - @plansor starting[-1 -2; -3 -4] += I[-1; -3] * removeunit(GR[1], 2)[-4; -2] - not_started = missing - else - not_started = removeunit(GR[1], 2) - end - - # finished - if isfinite(operator) && site == 1 - finished = missing - elseif !ismissing(ending) - I = id(storagetype(GL[end]), physicalspace(W)) - @plansor ending[-1 -2; -3 -4] += removeunit(GL[end], 2)[-1; -3] * I[-2; -4] - finished = missing - else - finished = removeunit(GL[end], 2) - end - - if nonzero_length(W.A) > 0 - continuing = AC_hamiltonian(GL_2, W.A, GR_2) - else - continuing = missing - end - - S = spacetype(GL) - M = storagetype(GL) - O1 = tensormaptype(S, 1, 1, M) - O2 = tensormaptype(S, 2, 2, M) - O3 = Core.Compiler.return_type(AC_hamiltonian, typeof((GL_2, W.A, GR_2))) - - return JordanMPO_AC_Hamiltonian{O1, O2, O3}( - onsite, not_started, finished, starting, ending, continuing - ) + return AC_hamiltonian(TensorMap(GL), TensorMap(W), TensorMap(GR)) end function AC2_hamiltonian( @@ -137,183 +20,5 @@ function AC2_hamiltonian( GR = rightenv(envs, site + 1, below) W1 = operator[site] W2 = operator[site + 1] - - GR_2 = GR[2:(end - 1)] - GL_2 = GL[2:(end - 1)] - - # starting left - continuing right - if nonzero_length(W1.C) > 0 && nonzero_length(W2.A) > 0 - @plansor CA_[-1 -2 -3; -4 -5 -6] ≔ W1.C[-1; -4 2] * W2.A[2 -2; -5 1] * - GR_2[-6 1; -3] - CA = only(CA_) - else - CA = missing - end - - # continuing left - ending right - if nonzero_length(W1.A) > 0 && nonzero_length(W2.B) > 0 - @plansor AB_[-1 -2 -3; -4 -5 -6] ≔ GL_2[-1 2; -4] * W1.A[2 -2; -5 1] * - W2.B[1 -3; -6] - AB = only(AB_) - else - AB = missing - end - - # middle - if nonzero_length(W1.C) > 0 && nonzero_length(W2.B) > 0 - if !ismissing(CA) - @plansor CA[-1 -2 -3; -4 -5 -6] += W1.C[-1; -4 1] * W2.B[1 -2; -5] * - removeunit(GR[end], 2)[-6; -3] - CB = missing - elseif !ismissing(AB) - @plansor AB[-1 -2 -3; -4 -5 -6] += removeunit(GL[1], 2)[-1; -4] * - W1.C[-2; -5 1] * W2.B[1 -3; -6] - CB = missing - else - @plansor CB_[-1 -2; -3 -4] ≔ W1.C[-1; -3 1] * W2.B[1 -2; -4] - CB = only(CB_) - end - else - CB = missing - end - - # starting right - if nonzero_length(W2.C) > 0 - if !ismissing(CA) - I = id(storagetype(GR[1]), physicalspace(W1)) - @plansor CA[-1 -2 -3; -4 -5 -6] += (I[-1; -4] * W2.C[-2; -5 1]) * - GR_2[-6 1; -3] - IC = missing - else - @plansor IC[-1 -2; -3 -4] ≔ W2.C[-1; -3 1] * GR_2[-4 1; -2] - IC = only(IC) - end - else - IC = missing - end - - # ending left - if nonzero_length(W1.B) > 0 - if !ismissing(AB) - I = id(storagetype(GL[end]), physicalspace(W2)) - @plansor AB[-1 -2 -3; -4 -5 -6] += GL_2[-1 1; -4] * - (W1.B[1 -2; -5] * I[-3; -6]) - BE = missing - else - @plansor BE[-1 -2; -3 -4] ≔ GL_2[-1 2; -3] * W1.B[2 -2; -4] - BE = only(BE) - end - else - BE = missing - end - - # onsite left - if nonzero_length(W1.D) > 0 - if !ismissing(BE) - @plansor BE[-1 -2; -3 -4] += removeunit(GL[1], 2)[-1; -3] * W1.D[-2; -4] - DE = missing - elseif !ismissing(AB) - I = id(storagetype(GL[end]), physicalspace(W2)) - @plansor AB[-1 -2 -3; -4 -5 -6] += removeunit(GL[1], 2)[-1; -4] * - (W1.D[-2; -5] * I[-3; -6]) - DE = missing - # TODO: could also try in CA? - else - DE = only(W1.D) - end - else - DE = missing - end - - # onsite right - if nonzero_length(W2.D) > 0 - if !ismissing(IC) - @plansor IC[-1 -2; -3 -4] += W2.D[-1; -3] * removeunit(GR[end], 2)[-4; -2] - ID = missing - elseif !ismissing(CA) - I = id(storagetype(GR[1]), physicalspace(W1)) - @plansor CA[-1 -2 -3; -4 -5 -6] += (I[-1; -4] * W2.D[-2; -5]) * - removeunit(GR[end], 2)[-6; -3] - ID = missing - else - ID = only(W2.D) - end - else - ID = missing - end - - # finished - if isfinite(operator) && site + 1 == length(operator) - II = missing - elseif !ismissing(IC) - I = id(storagetype(GR[1]), physicalspace(W2)) - @plansor IC[-1 -2; -3 -4] += I[-1; -3] * removeunit(GR[1], 2)[-4; -2] - II = missing - elseif !ismissing(CA) - I = id(storagetype(GR[1]), physicalspace(W1) ⊗ physicalspace(W2)) - @plansor CA[-1 -2 -3; -4 -5 -6] += I[-1 -2; -4 -5] * removeunit(GR[1], 2)[-6; -3] - II = missing - else - II = transpose(removeunit(GR[1], 2)) - end - - # unstarted - if isfinite(operator) && site == 1 - EE = missing - elseif !ismissing(BE) - I = id(storagetype(GL[end]), physicalspace(W1)) - @plansor BE[-1 -2; -3 -4] += removeunit(GL[end], 2)[-1; -3] * I[-2; -4] - EE = missing - elseif !ismissing(AB) - I = id(storagetype(GL[end]), physicalspace(W1) ⊗ physicalspace(W2)) - @plansor AB[-1 -2 -3; -4 -5 -6] += removeunit(GL[end], 2)[-1; -4] * I[-2 -3; -5 -6] - EE = missing - else - EE = removeunit(GL[end], 2) - end - - # continuing - continuing - ## TODO: Think about how one could and whether one should store these objects and use them for (a) advancing environments in iDMRG, (b) reuse ind backwards-sweep in IDMRG, (c) subspace expansion - if nonzero_length(W1.A) > 0 && nonzero_length(W2.A) > 0 - AA = AC2_hamiltonian(GL_2, W1.A, W2.A, GR_2) - else - AA = missing - end - - S = spacetype(GL) - M = storagetype(GL) - O1 = tensormaptype(S, 1, 1, M) - O2 = tensormaptype(S, 2, 2, M) - O3 = tensormaptype(S, 3, 3, M) - O4 = Core.Compiler.return_type(AC2_hamiltonian, typeof((GL_2, W1.A, W2.A, GR_2))) - - return JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4}(II, IC, ID, CB, CA, AB, AA, BE, DE, EE) -end - -# Actions -# ------- -function (H::JordanMPO_AC_Hamiltonian)(x::MPSTensor) - y = ismissing(H.A) ? zerovector(x) : H.A(x) - ismissing(H.D) || @plansor y[-1 -2; -3] += x[-1 1; -3] * H.D[-2; 1] - ismissing(H.E) || @plansor y[-1 -2; -3] += H.E[-1; 1] * x[1 -2; -3] - ismissing(H.I) || @plansor y[-1 -2; -3] += x[-1 -2; 1] * H.I[1; -3] - ismissing(H.C) || @plansor y[-1 -2; -3] += x[-1 2; 1] * H.C[-2 -3; 2 1] - ismissing(H.B) || @plansor y[-1 -2; -3] += H.B[-1 -2; 1 2] * x[1 2; -3] - - return y -end - -function (H::JordanMPO_AC2_Hamiltonian)(x::MPOTensor) - y = ismissing(H.AA) ? zerovector(x) : H.AA(x) - ismissing(H.II) || @plansor y[-1 -2; -3 -4] += x[-1 -2; 1 -4] * H.II[-3; 1] - ismissing(H.IC) || @plansor y[-1 -2; -3 -4] += x[-1 -2; 1 2] * H.IC[-4 -3; 2 1] - ismissing(H.ID) || @plansor y[-1 -2; -3 -4] += x[-1 -2; -3 1] * H.ID[-4; 1] - ismissing(H.CB) || @plansor y[-1 -2; -3 -4] += x[-1 1; -3 2] * H.CB[-2 -4; 1 2] - ismissing(H.CA) || @plansor y[-1 -2; -3 -4] += x[-1 1; 3 2] * H.CA[-2 -4 -3; 1 2 3] - ismissing(H.AB) || @plansor y[-1 -2; -3 -4] += x[1 2; -3 3] * H.AB[-1 -2 -4; 1 2 3] - ismissing(H.BE) || @plansor y[-1 -2; -3 -4] += x[1 2; -3 -4] * H.BE[-1 -2; 1 2] - ismissing(H.DE) || @plansor y[-1 -2; -3 -4] += x[-1 1; -3 -4] * H.DE[-2; 1] - ismissing(H.EE) || @plansor y[-1 -2; -3 -4] += x[1 -2; -3 -4] * H.EE[-1; 1] - - return y + return AC2_hamiltonian(TensorMap(GL), TensorMap(W1), TensorMap(W2), TensorMap(GR)) end diff --git a/src/algorithms/derivatives/mpo_derivatives.jl b/src/algorithms/derivatives/mpo_derivatives.jl index 02ec00b34..2232d221e 100644 --- a/src/algorithms/derivatives/mpo_derivatives.jl +++ b/src/algorithms/derivatives/mpo_derivatives.jl @@ -71,7 +71,8 @@ function AC2_hamiltonian(site::Int, below::_HAM_MPS_TYPES, operator::MPO{<:MPOTe end function AC2_hamiltonian(GL::MPSTensor, O1::MPOTensor, O2::MPOTensor, GR::MPSTensor) @plansor GLW[-1 -2 -3; -4 -5] ≔ GL[-1 1; -4] * O1[1 -2; -5 -3] - @plansor GWR[-1 -2 -3; -4 -5] ≔ O2[-3 -5; -2 1] * GR[-1 1; -4] + # @plansor GWR[-1 -2 -3; -4 -5] ≔ O2[-3 -5; -2 1] * GR[-1 1; -4] + @tensor GWR[-1 -2 -3; -4 -5] ≔ O2[-1 -5; -3 1] * GR[-2 1; -4] return MPO_Contracted_AC2_Hamiltonian(GLW, GWR) end @@ -167,6 +168,7 @@ end function (h::MPO_Contracted_AC2_Hamiltonian)( x::MPOTensor ) - @plansor y[-1 -2; -3 -4] ≔ h.leftenv[-1 -2 5; 1 2] * x[1 2; 3 4] * h.rightenv[3 4 5; -3 -4] + # @plansor y[-1 -2; -3 -4] ≔ h.leftenv[-1 -2 5; 1 2] * x[1 2; 3 4] * h.rightenv[3 4 5; -3 -4] + @tensor y[-1 -2; -3 -4] ≔ h.leftenv[-1 -2 3; 1 2] * x[1 2; 4 5] * h.rightenv[3 4 5; -3 -4] return y isa AbstractBlockTensorMap ? only(y) : y end \ No newline at end of file diff --git a/src/algorithms/groundstate/idmrg.jl b/src/algorithms/groundstate/idmrg.jl index 0b1e848fd..0ad55c23b 100644 --- a/src/algorithms/groundstate/idmrg.jl +++ b/src/algorithms/groundstate/idmrg.jl @@ -24,7 +24,7 @@ $(TYPEDFIELDS) alg_eigsolve::A = Defaults.alg_eigsolve() end -function find_groundstate(ost::InfiniteMPS, H, alg::IDMRG, envs = environments(ost, H)) +function find_groundstate(ost::AbstractInfiniteMPS, H, alg::IDMRG, envs = environments(ost, H)) ϵ::Float64 = calc_galerkin(ost, H, ost, envs) ψ = copy(ost) log = IterLog("IDMRG") @@ -131,7 +131,7 @@ $(TYPEDFIELDS) trscheme::TruncationStrategy end -function find_groundstate(ost::InfiniteMPS, H, alg::IDMRG2, envs = environments(ost, H)) +function find_groundstate(ost::AbstractInfiniteMPS, H, alg::IDMRG2, envs = environments(ost, H)) length(ost) < 2 && throw(ArgumentError("unit cell should be >= 2")) ϵ::Float64 = calc_galerkin(ost, H, ost, envs) diff --git a/src/algorithms/groundstate/vumps.jl b/src/algorithms/groundstate/vumps.jl index 99da8234c..418337047 100644 --- a/src/algorithms/groundstate/vumps.jl +++ b/src/algorithms/groundstate/vumps.jl @@ -46,7 +46,7 @@ struct VUMPSState{S, O, E} end function find_groundstate( - mps::InfiniteMPS, operator, alg::VUMPS, envs = environments(mps, operator) + mps::AbstractInfiniteMPS, operator, alg::VUMPS, envs = environments(mps, operator) ) return dominant_eigsolve(operator, mps, alg, envs; which = :SR) end diff --git a/src/environments/abstract_envs.jl b/src/environments/abstract_envs.jl index cbdcafee1..97caefe03 100644 --- a/src/environments/abstract_envs.jl +++ b/src/environments/abstract_envs.jl @@ -81,7 +81,8 @@ function environment_alg( tol = Defaults.tol, maxiter = Defaults.maxiter, krylovdim = Defaults.krylovdim, verbosity = Defaults.VERBOSE_NONE ) - return GMRES(; tol, maxiter, krylovdim, verbosity) + max_krylovdim = dim(left_virtualspace(above, 1)) * dim(left_virtualspace(below, 1)) + return GMRES(; tol, maxiter, krylovdim = min(max_krylovdim, krylovdim), verbosity) end function environment_alg( ::Union{InfiniteQP, MultilineQP}, ::Union{InfiniteMPO, MultilineMPO}, @@ -90,4 +91,4 @@ function environment_alg( verbosity = Defaults.VERBOSE_NONE ) return GMRES(; tol, maxiter, krylovdim, verbosity) -end +end \ No newline at end of file diff --git a/src/states/abstractmps.jl b/src/states/abstractmps.jl index dd756b2f4..2ac1cdeb8 100644 --- a/src/states/abstractmps.jl +++ b/src/states/abstractmps.jl @@ -15,7 +15,7 @@ const MPSTensor{S} = GenericMPSTensor{S, 2} # the usual mps tensors on which we """ MPSTensor([f, eltype], d::Int, left_D::Int, [right_D]::Int]) - MPSTensor([f, eltype], physicalspace::Union{S,CompositeSpace{S}}, + MPSTensor([f, eltype], physicalspace::Union{S,CompositeSpace{S}}, left_virtualspace::S, [right_virtualspace]::S) where {S<:ElementarySpace} Construct an `MPSTensor` with given physical and virtual spaces. @@ -93,10 +93,11 @@ Determine whether the given tensor is full rank, i.e. whether both the map from virtual space and the physical space to the right virtual space, and the map from the right virtual space and the physical space to the left virtual space are injective. """ -function isfullrank(A::GenericMPSTensor; side = :both) - Vₗ = _firstspace(A) - Vᵣ = _lastspace(A) - P = ⊗(space.(Ref(A), 2:(numind(A) - 1))...) +isfullrank(A::GenericMPSTensor; kwargs...) = isfullrank(space(A); kwargs...) +function isfullrank(V::TensorKit.TensorMapSpace; side = :both) + Vₗ = V[1] + Vᵣ = V[numind(V)] + P = ⊗(getindex.(Ref(V), 2:(numind(V) - 1))...) return if side === :both Vₗ ⊗ P ≿ Vᵣ' && Vₗ' ≾ P ⊗ Vᵣ elseif side === :right @@ -129,6 +130,29 @@ function makefullrank!(A::PeriodicVector{<:GenericMPSTensor}; alg = Defaults.alg return A end +function makefullrank!(virtualspaces::PeriodicVector{S}, physicalspaces::PeriodicVector{S}) where {S <: ElementarySpace} + haschanged = true + while haschanged + haschanged = false + for i in 1:length(virtualspaces) + Vmax = fuse(virtualspaces[i - 1], physicalspaces[i - 1]) + if !(virtualspaces[i] ≾ Vmax) + virtualspaces[i] = infimum(virtualspaces[i], Vmax) + haschanged = true + end + end + for i in reverse(1:length(virtualspaces)) + Vmax = fuse(dual(physicalspaces[i - 1]), virtualspaces[i]) + if !(virtualspaces[i - 1] ≾ Vmax) + virtualspaces[i - 1] = infimum(virtualspaces[i - 1], Vmax) + haschanged = true + end + end + end + + return virtualspaces +end + # Tensor accessors # ---------------- @doc """ @@ -177,7 +201,7 @@ TensorKit.sectortype(ψtype::Type{<:AbstractMPS}) = sectortype(site_type(ψtype) """ left_virtualspace(ψ::AbstractMPS, [pos=1:length(ψ)]) - + Return the virtual space of the bond to the left of sites `pos`. !!! warning @@ -220,4 +244,4 @@ physicalspace(ψ::AbstractMPS) = map(Base.Fix1(physicalspace, ψ), eachsite(ψ)) Return an iterator over the sites of the MPS `state`. """ -eachsite(ψ::AbstractMPS) = eachindex(ψ) +eachsite(ψ::AbstractMPS) = eachindex(ψ) \ No newline at end of file diff --git a/src/states/infinitemps.jl b/src/states/infinitemps.jl index c9a4cf906..79620ebac 100644 --- a/src/states/infinitemps.jl +++ b/src/states/infinitemps.jl @@ -1,5 +1,6 @@ +abstract type AbstractInfiniteMPS <: AbstractMPS end """ - InfiniteMPS{A<:GenericMPSTensor,B<:MPSBondTensor} <: AbtractMPS + InfiniteMPS{A<:GenericMPSTensor,B<:MPSBondTensor} <: AbstractInfiniteMPS Type that represents an infinite Matrix Product State. @@ -43,7 +44,7 @@ tensors `As`, or a list of left-gauged tensors `ALs`. - `tol`: gauge fixing tolerance - `maxiter`: gauge fixing maximum iterations """ -struct InfiniteMPS{A <: GenericMPSTensor, B <: MPSBondTensor} <: AbstractMPS +struct InfiniteMPS{A <: GenericMPSTensor, B <: MPSBondTensor} <: AbstractInfiniteMPS AL::PeriodicVector{A} AR::PeriodicVector{A} C::PeriodicVector{B} diff --git a/src/states/ortho.jl b/src/states/ortho.jl index a29942339..3e935f456 100644 --- a/src/states/ortho.jl +++ b/src/states/ortho.jl @@ -151,7 +151,7 @@ end Bring updated `AC` and `C` tensors back into a consistent set of left or right canonical tensors. This minimizes `∥AC_i - AL_i * C_i∥` or `∥AC_i - C_{i-1} * AR_i∥`. -The `alg` is passed on to `left_orth!` and `right_orth!`, and can be used to control the kind of +The `alg` is passed on to `left_orth!` and `right_orth!`, and can be used to control the kind of factorization used. By default, this is set to a (positive) QR/LQ, even though the optimal algorithm would use a polar decompositions instead, sacrificing a bit of performance for accuracy. diff --git a/src/utility/dynamictols.jl b/src/utility/dynamictols.jl index dc2a15781..61aae1f99 100644 --- a/src/utility/dynamictols.jl +++ b/src/utility/dynamictols.jl @@ -58,10 +58,11 @@ end Update the tolerance of the algorithm `alg` based on the current iteration `iter` and the current error `ϵ`, where the new tolerance is given by - + new_tol = clamp(ϵ * alg.tol_factor / sqrt(iter), alg.tol_min, alg.tol_max) """ function updatetol(alg::DynamicTol, iter::Integer, ϵ::Real) + iter = max(iter, one(iter)) new_tol = clamp(ϵ * alg.tol_factor / sqrt(iter), alg.tol_min, alg.tol_max) return _updatetol(alg.alg, new_tol) end @@ -71,4 +72,4 @@ function _updatetol(alg, tol::Real) return Accessors.@set alg.tol = tol end -end +end \ No newline at end of file From 23f9c232bc4d4615f677fb333155cb814d41f86a Mon Sep 17 00:00:00 2001 From: afeuerpfeil Date: Thu, 30 Oct 2025 18:04:33 -0400 Subject: [PATCH 19/19] Revert "sync" This reverts commit 69502678ffab8b6ec4992119d8851d690e4d2fc0. --- Project.toml | 4 +- src/MPSKit.jl | 4 +- src/algorithms/changebonds/changebonds.jl | 2 +- src/algorithms/changebonds/randexpand.jl | 134 +------- .../hamiltonian_derivatives copy.jl | 319 ------------------ .../derivatives/hamiltonian_derivatives.jl | 299 +++++++++++++++- src/algorithms/derivatives/mpo_derivatives.jl | 6 +- src/algorithms/groundstate/idmrg.jl | 4 +- src/algorithms/groundstate/vumps.jl | 2 +- src/environments/abstract_envs.jl | 5 +- src/states/abstractmps.jl | 38 +-- src/states/infinitemps.jl | 5 +- src/states/ortho.jl | 2 +- src/utility/dynamictols.jl | 5 +- 14 files changed, 339 insertions(+), 490 deletions(-) delete mode 100644 src/algorithms/derivatives/hamiltonian_derivatives copy.jl diff --git a/Project.toml b/Project.toml index f97025cba..23fc5e15f 100644 --- a/Project.toml +++ b/Project.toml @@ -26,7 +26,7 @@ VectorInterface = "409d34a3-91d5-4945-b6ec-7529ddf182d8" [compat] Accessors = "0.1" Aqua = "0.8.9" -BlockTensorKit = "0.2,0.3" +BlockTensorKit = "0.3" Combinatorics = "1" Compat = "3.47, 4.10" DocStringExtensions = "0.9.3" @@ -42,7 +42,7 @@ Plots = "1.40" Printf = "1" Random = "1" RecipesBase = "1.1" -TensorKit = "0.15.1,0.14" +TensorKit = "0.15.1" TensorKitManifolds = "0.7" TensorOperations = "5" Test = "1" diff --git a/src/MPSKit.jl b/src/MPSKit.jl index 1d09467fa..a6c586144 100644 --- a/src/MPSKit.jl +++ b/src/MPSKit.jl @@ -36,7 +36,7 @@ export FiniteExcited, QuasiparticleAnsatz, ChepigaAnsatz, ChepigaAnsatz2 export time_evolve, timestep, timestep!, make_time_mpo export TDVP, TDVP2, WI, WII, TaylorCluster export changebonds, changebonds! -export VUMPSSvdCut, OptimalExpand, SvdCut, RandExpand, RandPerturbedExpand +export VUMPSSvdCut, OptimalExpand, SvdCut, RandExpand export propagator export DynamicalDMRG, NaiveInvert, Jeckelmann export exact_diagonalization, fidelity_susceptibility @@ -194,4 +194,4 @@ function __init__() return nothing end -end \ No newline at end of file +end diff --git a/src/algorithms/changebonds/changebonds.jl b/src/algorithms/changebonds/changebonds.jl index d20e7399c..700b5ab0b 100644 --- a/src/algorithms/changebonds/changebonds.jl +++ b/src/algorithms/changebonds/changebonds.jl @@ -51,4 +51,4 @@ function _expand!(ψ::MultilineMPS, AL′::PeriodicMatrix, AR′::PeriodicMatrix _expand!(ψ[i], AL′[i, :], AR′[i, :]) end return ψ -end \ No newline at end of file +end diff --git a/src/algorithms/changebonds/randexpand.jl b/src/algorithms/changebonds/randexpand.jl index 73ed83f2a..28e87b3d6 100644 --- a/src/algorithms/changebonds/randexpand.jl +++ b/src/algorithms/changebonds/randexpand.jl @@ -2,9 +2,8 @@ $(TYPEDEF) An algorithm that expands the bond dimension by adding random unitary vectors that are -orthogonal to the existing state. This means that additional directions are added to -`AL` and `AR` that are contained in the nullspace of both. Note that this is happens in -parallel, and therefore the expansion will never go beyond the local two-site subspace. +orthogonal to the existing state. This is achieved by performing a truncated SVD on a random +two-site MPS tensor, which is made orthogonal to the existing state. ## Fields @@ -18,42 +17,39 @@ $(TYPEDFIELDS) trscheme::TruncationStrategy end -function changebonds!(ψ::InfiniteMPS, alg::RandExpand) +function changebonds(ψ::InfiniteMPS, alg::RandExpand) T = eltype(ψ.AL) AL′ = similar(ψ.AL) AR′ = similar(ψ.AR, tensormaptype(spacetype(T), 1, numind(T) - 1, storagetype(T))) - for i in 1:length(ψ) - # obtain spaces by sampling the support of both the left and right nullspace + # determine optimal expansion spaces around bond i + AC2 = randomize!(_transpose_front(ψ.AC[i]) * _transpose_tail(ψ.AR[i + 1])) + + # Use the nullspaces and SVD decomposition to determine the optimal expansion space VL = left_null(ψ.AL[i]) - VR = right_null!(_transpose_tail(ψ.AR[i + 1]; copy = true)) - V = sample_space(infimum(right_virtualspace(VL), space(VR, 1)), alg.trscheme) - - # obtain (orthogonal) directions as isometries in that direction - XL = randisometry(scalartype(VL), right_virtualspace(VL) ← V) - AL′[i] = VL * XL - XR = randisometry(storagetype(VL), space(VR, 1) ← V) - AR′[i + 1] = XR * VR + VR = right_null!(_transpose_tail(ψ.AR[i + 1])) + intermediate = normalize!(adjoint(VL) * AC2 * adjoint(VR)) + U, _, Vᴴ = svd_trunc!(intermediate; trunc = alg.trscheme, alg = alg.alg_svd) + + AL′[i] = VL * U + AR′[i + 1] = Vᴴ * VR end - return _expand!(ψ, AL′, AR′) + return _expand(ψ, AL′, AR′) end -function changebonds!(ψ::MultilineMPS, alg::RandExpand) - foreach(Base.Fix2(changebonds!, alg), ψ.data) - return ψ +function changebonds(ψ::MultilineMPS, alg::RandExpand) + return Multiline(map(x -> changebonds(x, alg), ψ.data)) end -changebonds(ψ::AbstractMPS, alg::RandExpand) = changebonds!(copy(ψ), alg) -changebonds(ψ::MultilineMPS, alg::RandExpand) = changebonds!(copy(ψ), alg) - +changebonds(ψ::AbstractFiniteMPS, alg::RandExpand) = changebonds!(copy(ψ), alg) function changebonds!(ψ::AbstractFiniteMPS, alg::RandExpand) for i in 1:(length(ψ) - 1) AC2 = randomize!(_transpose_front(ψ.AC[i]) * _transpose_tail(ψ.AR[i + 1])) #Calculate nullspaces for left and right NL = left_null(ψ.AC[i]) - NR = right_null!(_transpose_tail(ψ.AR[i + 1]; copy = true)) + NR = right_null!(_transpose_tail(ψ.AR[i + 1])) #Use this nullspaces and SVD decomposition to determine the optimal expansion space intermediate = normalize!(adjoint(NL) * AC2 * adjoint(NR)) @@ -71,97 +67,3 @@ function changebonds!(ψ::AbstractFiniteMPS, alg::RandExpand) return normalize!(ψ) end - -function sample_space(V, strategy) - S = TensorKit.SectorDict(c => sort!(Random.rand(dim(V, c)); rev = true) for c in sectors(V)) - ind = MatrixAlgebraKit.findtruncated_svd(S, strategy) - return TensorKit.Factorizations.truncate_space(V, ind) -end - - -""" -$(TYPEDEF) - -An algorithm that expands the bond dimension by adding random unitary vectors that are -orthogonal to the existing state, in a sweeping fashion. Additionally, some random noise -is added to the state in order for it to remain gauge-able. - -## Fields - -$(TYPEDFIELDS) -""" -@kwdef struct RandPerturbedExpand{S} <: Algorithm - "algorithm used for the singular value decomposition" - alg_svd::S = Defaults.alg_svd() - - "algorithm used for [truncation](@extref MatrixAlgebraKit.TruncationStrategy] the expanded space" - trscheme::TruncationStrategy - - "amount of noise that is added to the current state" - noisefactor::Float64 = eps()^(3 / 4) - - "algorithm used for gauging the state" - alg_gauge = Defaults.alg_gauge(; dynamic_tols = false) -end - -function changebonds!(ψ::InfiniteMPS, alg::RandPerturbedExpand) - for i in 1:length(ψ) - # obtain space by sampling the support of left nullspace - # add (orthogonal) directions as isometries in that direction - VL = left_null(ψ.AL[i]) - V = sample_space(right_virtualspace(VL), alg.trscheme) - XL = randisometry(scalartype(VL), right_virtualspace(VL) ← V) - ψ.AL[i] = catdomain(ψ.AL[i], VL * XL) - - # make sure the next site fits, by "absorbing" into a larger tensor - # with some random noise to ensure state is still gauge-able - AL = ψ.AL[i + 1] - AL′ = similar(AL, right_virtualspace(ψ.AL[i]) ⊗ physicalspace(AL) ← right_virtualspace(AL)) - scale!(randn!(AL′), alg.noisefactor) - ψ.AL[i + 1] = TensorKit.absorb!(AL′, AL) - end - - # properly regauge the state: - makefullrank!(ψ.AL) - ψ.AR .= similar.(ψ.AL) - # ψ.AC .= similar.(ψ.AL) - - # initial guess for gauge is embedded original C - C₀ = similar(ψ.C[0], right_virtualspace(ψ.AL[end]) ← left_virtualspace(ψ.AL[1])) - absorb!(id!(C₀), ψ.C[0]) - - gaugefix!(ψ, ψ.AL, C₀; order = :R, alg.alg_gauge.maxiter, alg.alg_gauge.tol) - - for i in reverse(1:length(ψ)) - # obtain space by sampling the support of left nullspace - # add (orthogonal) directions as isometries in that direction - AR_tail = _transpose_tail(ψ.AR[i]) - VR = right_null(AR_tail) - V = sample_space(space(VR, 1), alg.trscheme) - XR = randisometry(scalartype(VR), space(VR, 1) ← V) - ψ.AR[i] = _transpose_front(catcodomain(AR_tail, XR' * VR)) - - # make sure the next site fits, by "absorbing" into a larger tensor - # with some random noise to ensure state is still gauge-able - AR = ψ.AR[i - 1] - AR′ = similar(AR, left_virtualspace(AR) ⊗ physicalspace(AR) ← left_virtualspace(ψ.AR[i])) - scale!(randn!(AR′), alg.noisefactor) - ψ.AR[i - 1] = TensorKit.absorb!(AR′, AR) - end - - # properly regauge the state: - makefullrank!(ψ.AR) - ψ.AL .= similar.(ψ.AR) - ψ.AC .= similar.(ψ.AR) - - # initial guess for gauge is embedded original C - C₀ = similar(ψ.C[0], right_virtualspace(ψ.AR[end]) ← left_virtualspace(ψ.AR[1])) - absorb!(id!(C₀), ψ.C[0]) - - gaugefix!(ψ, ψ.AR, C₀; order = :LR, alg.alg_gauge.maxiter, alg.alg_gauge.tol) - mul!.(ψ.AC, ψ.AL, ψ.C) - - return ψ -end - -changebonds(ψ::InfiniteMPS, alg::RandPerturbedExpand) = changebonds!(copy(ψ), alg) \ No newline at end of file diff --git a/src/algorithms/derivatives/hamiltonian_derivatives copy.jl b/src/algorithms/derivatives/hamiltonian_derivatives copy.jl deleted file mode 100644 index f5290ebff..000000000 --- a/src/algorithms/derivatives/hamiltonian_derivatives copy.jl +++ /dev/null @@ -1,319 +0,0 @@ -""" - JordanMPO_AC_Hamiltonian{O1,O2,O3} - -Efficient operator for representing the single-site derivative of a `MPOHamiltonian` sandwiched between two MPSs. -In particular, this operator aims to make maximal use of the structure of the `MPOHamiltonian` to reduce the number of operations required to apply the operator to a tensor. -""" -struct JordanMPO_AC_Hamiltonian{O1, O2, O3} <: DerivativeOperator - D::Union{O1, Missing} # onsite - I::Union{O1, Missing} # not started - E::Union{O1, Missing} # finished - C::Union{O2, Missing} # starting - B::Union{O2, Missing} # ending - A::Union{O3, Missing} # continuing - - # need inner constructor to prohibit no-type-param constructor with unbound vars - function JordanMPO_AC_Hamiltonian{O1, O2, O3}( - D, I, E, C, B, A, - ) where {O1, O2, O3} - return new{O1, O2, O3}(D, I, E, C, B, A) - end -end - -""" - JordanMPO_AC2_Hamiltonian{O1,O2,O3,O4} - -Efficient operator for representing the single-site derivative of a `MPOHamiltonian` sandwiched between two MPSs. -In particular, this operator aims to make maximal use of the structure of the `MPOHamiltonian` to reduce the number of operations required to apply the operator to a tensor. -""" -struct JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4} <: DerivativeOperator - II::Union{O1, Missing} # not_started - IC::Union{O2, Missing} # starting right - ID::Union{O1, Missing} # onsite right - CB::Union{O2, Missing} # starting left - ending right - CA::Union{O3, Missing} # starting left - continuing right - AB::Union{O3, Missing} # continuing left - ending right - AA::Union{O4, Missing} # continuing left - continuing right - BE::Union{O2, Missing} # ending left - DE::Union{O1, Missing} # onsite left - EE::Union{O1, Missing} # finished - - # need inner constructor to prohibit no-type-param constructor with unbound vars - function JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4}( - II, IC, ID, CB, CA, AB, AA, BE, DE, EE - ) where {O1, O2, O3, O4} - return new{O1, O2, O3, O4}(II, IC, ID, CB, CA, AB, AA, BE, DE, EE) - end -end - -# Constructors -# ------------ -function AC_hamiltonian( - site::Int, below::_HAM_MPS_TYPES, operator::MPOHamiltonian{<:JordanMPOTensor}, - above::_HAM_MPS_TYPES, envs - ) - GL = leftenv(envs, site, below) - GR = rightenv(envs, site, below) - W = operator[site] - - GR_2 = GR[2:(end - 1)] - GL_2 = GL[2:(end - 1)] - - # starting - if nonzero_length(W.C) > 0 - @plansor starting_[-1 -2; -3 -4] ≔ W.C[-1; -3 1] * GR_2[-4 1; -2] - starting = only(starting_) - else - starting = missing - end - - # ending - if nonzero_length(W.B) > 0 - @plansor ending_[-1 -2; -3 -4] ≔ GL_2[-1 1; -3] * W.B[1 -2; -4] - ending = only(ending_) - else - ending = missing - end - - # onsite - if nonzero_length(W.D) > 0 - if !ismissing(starting) - @plansor starting[-1 -2; -3 -4] += W.D[-1; -3] * removeunit(GR[end], 2)[-4; -2] - onsite = missing - elseif !ismissing(ending) - @plansor ending[-1 -2; -3 -4] += removeunit(GL[1], 2)[-1; -3] * W.D[-2; -4] - onsite = missing - else - onsite = W.D - end - else - onsite = missing - end - - # not_started - if isfinite(operator) && site == length(operator) - not_started = missing - elseif !ismissing(starting) - I = id(storagetype(GR[1]), physicalspace(W)) - @plansor starting[-1 -2; -3 -4] += I[-1; -3] * removeunit(GR[1], 2)[-4; -2] - not_started = missing - else - not_started = removeunit(GR[1], 2) - end - - # finished - if isfinite(operator) && site == 1 - finished = missing - elseif !ismissing(ending) - I = id(storagetype(GL[end]), physicalspace(W)) - @plansor ending[-1 -2; -3 -4] += removeunit(GL[end], 2)[-1; -3] * I[-2; -4] - finished = missing - else - finished = removeunit(GL[end], 2) - end - - if nonzero_length(W.A) > 0 - continuing = AC_hamiltonian(GL_2, W.A, GR_2) - else - continuing = missing - end - - S = spacetype(GL) - M = storagetype(GL) - O1 = tensormaptype(S, 1, 1, M) - O2 = tensormaptype(S, 2, 2, M) - O3 = Core.Compiler.return_type(AC_hamiltonian, typeof((GL_2, W.A, GR_2))) - - return JordanMPO_AC_Hamiltonian{O1, O2, O3}( - onsite, not_started, finished, starting, ending, continuing - ) -end - -function AC2_hamiltonian( - site::Int, below::_HAM_MPS_TYPES, operator::MPOHamiltonian{<:JordanMPOTensor}, - above::_HAM_MPS_TYPES, envs - ) - GL = leftenv(envs, site, below) - GR = rightenv(envs, site + 1, below) - W1 = operator[site] - W2 = operator[site + 1] - - GR_2 = GR[2:(end - 1)] - GL_2 = GL[2:(end - 1)] - - # starting left - continuing right - if nonzero_length(W1.C) > 0 && nonzero_length(W2.A) > 0 - @plansor CA_[-1 -2 -3; -4 -5 -6] ≔ W1.C[-1; -4 2] * W2.A[2 -2; -5 1] * - GR_2[-6 1; -3] - CA = only(CA_) - else - CA = missing - end - - # continuing left - ending right - if nonzero_length(W1.A) > 0 && nonzero_length(W2.B) > 0 - @plansor AB_[-1 -2 -3; -4 -5 -6] ≔ GL_2[-1 2; -4] * W1.A[2 -2; -5 1] * - W2.B[1 -3; -6] - AB = only(AB_) - else - AB = missing - end - - # middle - if nonzero_length(W1.C) > 0 && nonzero_length(W2.B) > 0 - if !ismissing(CA) - @plansor CA[-1 -2 -3; -4 -5 -6] += W1.C[-1; -4 1] * W2.B[1 -2; -5] * - removeunit(GR[end], 2)[-6; -3] - CB = missing - elseif !ismissing(AB) - @plansor AB[-1 -2 -3; -4 -5 -6] += removeunit(GL[1], 2)[-1; -4] * - W1.C[-2; -5 1] * W2.B[1 -3; -6] - CB = missing - else - @plansor CB_[-1 -2; -3 -4] ≔ W1.C[-1; -3 1] * W2.B[1 -2; -4] - CB = only(CB_) - end - else - CB = missing - end - - # starting right - if nonzero_length(W2.C) > 0 - if !ismissing(CA) - I = id(storagetype(GR[1]), physicalspace(W1)) - @plansor CA[-1 -2 -3; -4 -5 -6] += (I[-1; -4] * W2.C[-2; -5 1]) * - GR_2[-6 1; -3] - IC = missing - else - @plansor IC[-1 -2; -3 -4] ≔ W2.C[-1; -3 1] * GR_2[-4 1; -2] - IC = only(IC) - end - else - IC = missing - end - - # ending left - if nonzero_length(W1.B) > 0 - if !ismissing(AB) - I = id(storagetype(GL[end]), physicalspace(W2)) - @plansor AB[-1 -2 -3; -4 -5 -6] += GL_2[-1 1; -4] * - (W1.B[1 -2; -5] * I[-3; -6]) - BE = missing - else - @plansor BE[-1 -2; -3 -4] ≔ GL_2[-1 2; -3] * W1.B[2 -2; -4] - BE = only(BE) - end - else - BE = missing - end - - # onsite left - if nonzero_length(W1.D) > 0 - if !ismissing(BE) - @plansor BE[-1 -2; -3 -4] += removeunit(GL[1], 2)[-1; -3] * W1.D[-2; -4] - DE = missing - elseif !ismissing(AB) - I = id(storagetype(GL[end]), physicalspace(W2)) - @plansor AB[-1 -2 -3; -4 -5 -6] += removeunit(GL[1], 2)[-1; -4] * - (W1.D[-2; -5] * I[-3; -6]) - DE = missing - # TODO: could also try in CA? - else - DE = only(W1.D) - end - else - DE = missing - end - - # onsite right - if nonzero_length(W2.D) > 0 - if !ismissing(IC) - @plansor IC[-1 -2; -3 -4] += W2.D[-1; -3] * removeunit(GR[end], 2)[-4; -2] - ID = missing - elseif !ismissing(CA) - I = id(storagetype(GR[1]), physicalspace(W1)) - @plansor CA[-1 -2 -3; -4 -5 -6] += (I[-1; -4] * W2.D[-2; -5]) * - removeunit(GR[end], 2)[-6; -3] - ID = missing - else - ID = only(W2.D) - end - else - ID = missing - end - - # finished - if isfinite(operator) && site + 1 == length(operator) - II = missing - elseif !ismissing(IC) - I = id(storagetype(GR[1]), physicalspace(W2)) - @plansor IC[-1 -2; -3 -4] += I[-1; -3] * removeunit(GR[1], 2)[-4; -2] - II = missing - elseif !ismissing(CA) - I = id(storagetype(GR[1]), physicalspace(W1) ⊗ physicalspace(W2)) - @plansor CA[-1 -2 -3; -4 -5 -6] += I[-1 -2; -4 -5] * removeunit(GR[1], 2)[-6; -3] - II = missing - else - II = transpose(removeunit(GR[1], 2)) - end - - # unstarted - if isfinite(operator) && site == 1 - EE = missing - elseif !ismissing(BE) - I = id(storagetype(GL[end]), physicalspace(W1)) - @plansor BE[-1 -2; -3 -4] += removeunit(GL[end], 2)[-1; -3] * I[-2; -4] - EE = missing - elseif !ismissing(AB) - I = id(storagetype(GL[end]), physicalspace(W1) ⊗ physicalspace(W2)) - @plansor AB[-1 -2 -3; -4 -5 -6] += removeunit(GL[end], 2)[-1; -4] * I[-2 -3; -5 -6] - EE = missing - else - EE = removeunit(GL[end], 2) - end - - # continuing - continuing - ## TODO: Think about how one could and whether one should store these objects and use them for (a) advancing environments in iDMRG, (b) reuse ind backwards-sweep in IDMRG, (c) subspace expansion - if nonzero_length(W1.A) > 0 && nonzero_length(W2.A) > 0 - AA = AC2_hamiltonian(GL_2, W1.A, W2.A, GR_2) - else - AA = missing - end - - S = spacetype(GL) - M = storagetype(GL) - O1 = tensormaptype(S, 1, 1, M) - O2 = tensormaptype(S, 2, 2, M) - O3 = tensormaptype(S, 3, 3, M) - O4 = Core.Compiler.return_type(AC2_hamiltonian, typeof((GL_2, W1.A, W2.A, GR_2))) - - return JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4}(II, IC, ID, CB, CA, AB, AA, BE, DE, EE) -end - -# Actions -# ------- -function (H::JordanMPO_AC_Hamiltonian)(x::MPSTensor) - y = ismissing(H.A) ? zerovector(x) : H.A(x) - ismissing(H.D) || @plansor y[-1 -2; -3] += x[-1 1; -3] * H.D[-2; 1] - ismissing(H.E) || @plansor y[-1 -2; -3] += H.E[-1; 1] * x[1 -2; -3] - ismissing(H.I) || @plansor y[-1 -2; -3] += x[-1 -2; 1] * H.I[1; -3] - ismissing(H.C) || @plansor y[-1 -2; -3] += x[-1 2; 1] * H.C[-2 -3; 2 1] - ismissing(H.B) || @plansor y[-1 -2; -3] += H.B[-1 -2; 1 2] * x[1 2; -3] - - return y -end - -function (H::JordanMPO_AC2_Hamiltonian)(x::MPOTensor) - y = ismissing(H.AA) ? zerovector(x) : H.AA(x) - ismissing(H.II) || @plansor y[-1 -2; -3 -4] += x[-1 -2; 1 -4] * H.II[-3; 1] - ismissing(H.IC) || @plansor y[-1 -2; -3 -4] += x[-1 -2; 1 2] * H.IC[-4 -3; 2 1] - ismissing(H.ID) || @plansor y[-1 -2; -3 -4] += x[-1 -2; -3 1] * H.ID[-4; 1] - ismissing(H.CB) || @plansor y[-1 -2; -3 -4] += x[-1 1; -3 2] * H.CB[-2 -4; 1 2] - ismissing(H.CA) || @plansor y[-1 -2; -3 -4] += x[-1 1; 3 2] * H.CA[-2 -4 -3; 1 2 3] - ismissing(H.AB) || @plansor y[-1 -2; -3 -4] += x[1 2; -3 3] * H.AB[-1 -2 -4; 1 2 3] - ismissing(H.BE) || @plansor y[-1 -2; -3 -4] += x[1 2; -3 -4] * H.BE[-1 -2; 1 2] - ismissing(H.DE) || @plansor y[-1 -2; -3 -4] += x[-1 1; -3 -4] * H.DE[-2; 1] - ismissing(H.EE) || @plansor y[-1 -2; -3 -4] += x[1 -2; -3 -4] * H.EE[-1; 1] - - return y -end \ No newline at end of file diff --git a/src/algorithms/derivatives/hamiltonian_derivatives.jl b/src/algorithms/derivatives/hamiltonian_derivatives.jl index 009bdd069..b35bf85d2 100644 --- a/src/algorithms/derivatives/hamiltonian_derivatives.jl +++ b/src/algorithms/derivatives/hamiltonian_derivatives.jl @@ -1,3 +1,50 @@ +""" + JordanMPO_AC_Hamiltonian{O1,O2,O3} + +Efficient operator for representing the single-site derivative of a `MPOHamiltonian` sandwiched between two MPSs. +In particular, this operator aims to make maximal use of the structure of the `MPOHamiltonian` to reduce the number of operations required to apply the operator to a tensor. +""" +struct JordanMPO_AC_Hamiltonian{O1, O2, O3} <: DerivativeOperator + D::Union{O1, Missing} # onsite + I::Union{O1, Missing} # not started + E::Union{O1, Missing} # finished + C::Union{O2, Missing} # starting + B::Union{O2, Missing} # ending + A::Union{O3, Missing} # continuing + + # need inner constructor to prohibit no-type-param constructor with unbound vars + function JordanMPO_AC_Hamiltonian{O1, O2, O3}( + D, I, E, C, B, A, + ) where {O1, O2, O3} + return new{O1, O2, O3}(D, I, E, C, B, A) + end +end + +""" + JordanMPO_AC2_Hamiltonian{O1,O2,O3,O4} + +Efficient operator for representing the single-site derivative of a `MPOHamiltonian` sandwiched between two MPSs. +In particular, this operator aims to make maximal use of the structure of the `MPOHamiltonian` to reduce the number of operations required to apply the operator to a tensor. +""" +struct JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4} <: DerivativeOperator + II::Union{O1, Missing} # not_started + IC::Union{O2, Missing} # starting right + ID::Union{O1, Missing} # onsite right + CB::Union{O2, Missing} # starting left - ending right + CA::Union{O3, Missing} # starting left - continuing right + AB::Union{O3, Missing} # continuing left - ending right + AA::Union{O4, Missing} # continuing left - continuing right + BE::Union{O2, Missing} # ending left + DE::Union{O1, Missing} # onsite left + EE::Union{O1, Missing} # finished + + # need inner constructor to prohibit no-type-param constructor with unbound vars + function JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4}( + II, IC, ID, CB, CA, AB, AA, BE, DE, EE + ) where {O1, O2, O3, O4} + return new{O1, O2, O3, O4}(II, IC, ID, CB, CA, AB, AA, BE, DE, EE) + end +end # Constructors # ------------ @@ -9,7 +56,77 @@ function AC_hamiltonian( GR = rightenv(envs, site, below) W = operator[site] - return AC_hamiltonian(TensorMap(GL), TensorMap(W), TensorMap(GR)) + GR_2 = GR[2:(end - 1)] + GL_2 = GL[2:(end - 1)] + + # starting + if nonzero_length(W.C) > 0 + @plansor starting_[-1 -2; -3 -4] ≔ W.C[-1; -3 1] * GR_2[-4 1; -2] + starting = only(starting_) + else + starting = missing + end + + # ending + if nonzero_length(W.B) > 0 + @plansor ending_[-1 -2; -3 -4] ≔ GL_2[-1 1; -3] * W.B[1 -2; -4] + ending = only(ending_) + else + ending = missing + end + + # onsite + if nonzero_length(W.D) > 0 + if !ismissing(starting) + @plansor starting[-1 -2; -3 -4] += W.D[-1; -3] * removeunit(GR[end], 2)[-4; -2] + onsite = missing + elseif !ismissing(ending) + @plansor ending[-1 -2; -3 -4] += removeunit(GL[1], 2)[-1; -3] * W.D[-2; -4] + onsite = missing + else + onsite = W.D + end + else + onsite = missing + end + + # not_started + if isfinite(operator) && site == length(operator) + not_started = missing + elseif !ismissing(starting) + I = id(storagetype(GR[1]), physicalspace(W)) + @plansor starting[-1 -2; -3 -4] += I[-1; -3] * removeunit(GR[1], 2)[-4; -2] + not_started = missing + else + not_started = removeunit(GR[1], 2) + end + + # finished + if isfinite(operator) && site == 1 + finished = missing + elseif !ismissing(ending) + I = id(storagetype(GL[end]), physicalspace(W)) + @plansor ending[-1 -2; -3 -4] += removeunit(GL[end], 2)[-1; -3] * I[-2; -4] + finished = missing + else + finished = removeunit(GL[end], 2) + end + + if nonzero_length(W.A) > 0 + continuing = AC_hamiltonian(GL_2, W.A, GR_2) + else + continuing = missing + end + + S = spacetype(GL) + M = storagetype(GL) + O1 = tensormaptype(S, 1, 1, M) + O2 = tensormaptype(S, 2, 2, M) + O3 = Core.Compiler.return_type(AC_hamiltonian, typeof((GL_2, W.A, GR_2))) + + return JordanMPO_AC_Hamiltonian{O1, O2, O3}( + onsite, not_started, finished, starting, ending, continuing + ) end function AC2_hamiltonian( @@ -20,5 +137,183 @@ function AC2_hamiltonian( GR = rightenv(envs, site + 1, below) W1 = operator[site] W2 = operator[site + 1] - return AC2_hamiltonian(TensorMap(GL), TensorMap(W1), TensorMap(W2), TensorMap(GR)) + + GR_2 = GR[2:(end - 1)] + GL_2 = GL[2:(end - 1)] + + # starting left - continuing right + if nonzero_length(W1.C) > 0 && nonzero_length(W2.A) > 0 + @plansor CA_[-1 -2 -3; -4 -5 -6] ≔ W1.C[-1; -4 2] * W2.A[2 -2; -5 1] * + GR_2[-6 1; -3] + CA = only(CA_) + else + CA = missing + end + + # continuing left - ending right + if nonzero_length(W1.A) > 0 && nonzero_length(W2.B) > 0 + @plansor AB_[-1 -2 -3; -4 -5 -6] ≔ GL_2[-1 2; -4] * W1.A[2 -2; -5 1] * + W2.B[1 -3; -6] + AB = only(AB_) + else + AB = missing + end + + # middle + if nonzero_length(W1.C) > 0 && nonzero_length(W2.B) > 0 + if !ismissing(CA) + @plansor CA[-1 -2 -3; -4 -5 -6] += W1.C[-1; -4 1] * W2.B[1 -2; -5] * + removeunit(GR[end], 2)[-6; -3] + CB = missing + elseif !ismissing(AB) + @plansor AB[-1 -2 -3; -4 -5 -6] += removeunit(GL[1], 2)[-1; -4] * + W1.C[-2; -5 1] * W2.B[1 -3; -6] + CB = missing + else + @plansor CB_[-1 -2; -3 -4] ≔ W1.C[-1; -3 1] * W2.B[1 -2; -4] + CB = only(CB_) + end + else + CB = missing + end + + # starting right + if nonzero_length(W2.C) > 0 + if !ismissing(CA) + I = id(storagetype(GR[1]), physicalspace(W1)) + @plansor CA[-1 -2 -3; -4 -5 -6] += (I[-1; -4] * W2.C[-2; -5 1]) * + GR_2[-6 1; -3] + IC = missing + else + @plansor IC[-1 -2; -3 -4] ≔ W2.C[-1; -3 1] * GR_2[-4 1; -2] + IC = only(IC) + end + else + IC = missing + end + + # ending left + if nonzero_length(W1.B) > 0 + if !ismissing(AB) + I = id(storagetype(GL[end]), physicalspace(W2)) + @plansor AB[-1 -2 -3; -4 -5 -6] += GL_2[-1 1; -4] * + (W1.B[1 -2; -5] * I[-3; -6]) + BE = missing + else + @plansor BE[-1 -2; -3 -4] ≔ GL_2[-1 2; -3] * W1.B[2 -2; -4] + BE = only(BE) + end + else + BE = missing + end + + # onsite left + if nonzero_length(W1.D) > 0 + if !ismissing(BE) + @plansor BE[-1 -2; -3 -4] += removeunit(GL[1], 2)[-1; -3] * W1.D[-2; -4] + DE = missing + elseif !ismissing(AB) + I = id(storagetype(GL[end]), physicalspace(W2)) + @plansor AB[-1 -2 -3; -4 -5 -6] += removeunit(GL[1], 2)[-1; -4] * + (W1.D[-2; -5] * I[-3; -6]) + DE = missing + # TODO: could also try in CA? + else + DE = only(W1.D) + end + else + DE = missing + end + + # onsite right + if nonzero_length(W2.D) > 0 + if !ismissing(IC) + @plansor IC[-1 -2; -3 -4] += W2.D[-1; -3] * removeunit(GR[end], 2)[-4; -2] + ID = missing + elseif !ismissing(CA) + I = id(storagetype(GR[1]), physicalspace(W1)) + @plansor CA[-1 -2 -3; -4 -5 -6] += (I[-1; -4] * W2.D[-2; -5]) * + removeunit(GR[end], 2)[-6; -3] + ID = missing + else + ID = only(W2.D) + end + else + ID = missing + end + + # finished + if isfinite(operator) && site + 1 == length(operator) + II = missing + elseif !ismissing(IC) + I = id(storagetype(GR[1]), physicalspace(W2)) + @plansor IC[-1 -2; -3 -4] += I[-1; -3] * removeunit(GR[1], 2)[-4; -2] + II = missing + elseif !ismissing(CA) + I = id(storagetype(GR[1]), physicalspace(W1) ⊗ physicalspace(W2)) + @plansor CA[-1 -2 -3; -4 -5 -6] += I[-1 -2; -4 -5] * removeunit(GR[1], 2)[-6; -3] + II = missing + else + II = transpose(removeunit(GR[1], 2)) + end + + # unstarted + if isfinite(operator) && site == 1 + EE = missing + elseif !ismissing(BE) + I = id(storagetype(GL[end]), physicalspace(W1)) + @plansor BE[-1 -2; -3 -4] += removeunit(GL[end], 2)[-1; -3] * I[-2; -4] + EE = missing + elseif !ismissing(AB) + I = id(storagetype(GL[end]), physicalspace(W1) ⊗ physicalspace(W2)) + @plansor AB[-1 -2 -3; -4 -5 -6] += removeunit(GL[end], 2)[-1; -4] * I[-2 -3; -5 -6] + EE = missing + else + EE = removeunit(GL[end], 2) + end + + # continuing - continuing + ## TODO: Think about how one could and whether one should store these objects and use them for (a) advancing environments in iDMRG, (b) reuse ind backwards-sweep in IDMRG, (c) subspace expansion + if nonzero_length(W1.A) > 0 && nonzero_length(W2.A) > 0 + AA = AC2_hamiltonian(GL_2, W1.A, W2.A, GR_2) + else + AA = missing + end + + S = spacetype(GL) + M = storagetype(GL) + O1 = tensormaptype(S, 1, 1, M) + O2 = tensormaptype(S, 2, 2, M) + O3 = tensormaptype(S, 3, 3, M) + O4 = Core.Compiler.return_type(AC2_hamiltonian, typeof((GL_2, W1.A, W2.A, GR_2))) + + return JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4}(II, IC, ID, CB, CA, AB, AA, BE, DE, EE) +end + +# Actions +# ------- +function (H::JordanMPO_AC_Hamiltonian)(x::MPSTensor) + y = ismissing(H.A) ? zerovector(x) : H.A(x) + ismissing(H.D) || @plansor y[-1 -2; -3] += x[-1 1; -3] * H.D[-2; 1] + ismissing(H.E) || @plansor y[-1 -2; -3] += H.E[-1; 1] * x[1 -2; -3] + ismissing(H.I) || @plansor y[-1 -2; -3] += x[-1 -2; 1] * H.I[1; -3] + ismissing(H.C) || @plansor y[-1 -2; -3] += x[-1 2; 1] * H.C[-2 -3; 2 1] + ismissing(H.B) || @plansor y[-1 -2; -3] += H.B[-1 -2; 1 2] * x[1 2; -3] + + return y +end + +function (H::JordanMPO_AC2_Hamiltonian)(x::MPOTensor) + y = ismissing(H.AA) ? zerovector(x) : H.AA(x) + ismissing(H.II) || @plansor y[-1 -2; -3 -4] += x[-1 -2; 1 -4] * H.II[-3; 1] + ismissing(H.IC) || @plansor y[-1 -2; -3 -4] += x[-1 -2; 1 2] * H.IC[-4 -3; 2 1] + ismissing(H.ID) || @plansor y[-1 -2; -3 -4] += x[-1 -2; -3 1] * H.ID[-4; 1] + ismissing(H.CB) || @plansor y[-1 -2; -3 -4] += x[-1 1; -3 2] * H.CB[-2 -4; 1 2] + ismissing(H.CA) || @plansor y[-1 -2; -3 -4] += x[-1 1; 3 2] * H.CA[-2 -4 -3; 1 2 3] + ismissing(H.AB) || @plansor y[-1 -2; -3 -4] += x[1 2; -3 3] * H.AB[-1 -2 -4; 1 2 3] + ismissing(H.BE) || @plansor y[-1 -2; -3 -4] += x[1 2; -3 -4] * H.BE[-1 -2; 1 2] + ismissing(H.DE) || @plansor y[-1 -2; -3 -4] += x[-1 1; -3 -4] * H.DE[-2; 1] + ismissing(H.EE) || @plansor y[-1 -2; -3 -4] += x[1 -2; -3 -4] * H.EE[-1; 1] + + return y end diff --git a/src/algorithms/derivatives/mpo_derivatives.jl b/src/algorithms/derivatives/mpo_derivatives.jl index 2232d221e..02ec00b34 100644 --- a/src/algorithms/derivatives/mpo_derivatives.jl +++ b/src/algorithms/derivatives/mpo_derivatives.jl @@ -71,8 +71,7 @@ function AC2_hamiltonian(site::Int, below::_HAM_MPS_TYPES, operator::MPO{<:MPOTe end function AC2_hamiltonian(GL::MPSTensor, O1::MPOTensor, O2::MPOTensor, GR::MPSTensor) @plansor GLW[-1 -2 -3; -4 -5] ≔ GL[-1 1; -4] * O1[1 -2; -5 -3] - # @plansor GWR[-1 -2 -3; -4 -5] ≔ O2[-3 -5; -2 1] * GR[-1 1; -4] - @tensor GWR[-1 -2 -3; -4 -5] ≔ O2[-1 -5; -3 1] * GR[-2 1; -4] + @plansor GWR[-1 -2 -3; -4 -5] ≔ O2[-3 -5; -2 1] * GR[-1 1; -4] return MPO_Contracted_AC2_Hamiltonian(GLW, GWR) end @@ -168,7 +167,6 @@ end function (h::MPO_Contracted_AC2_Hamiltonian)( x::MPOTensor ) - # @plansor y[-1 -2; -3 -4] ≔ h.leftenv[-1 -2 5; 1 2] * x[1 2; 3 4] * h.rightenv[3 4 5; -3 -4] - @tensor y[-1 -2; -3 -4] ≔ h.leftenv[-1 -2 3; 1 2] * x[1 2; 4 5] * h.rightenv[3 4 5; -3 -4] + @plansor y[-1 -2; -3 -4] ≔ h.leftenv[-1 -2 5; 1 2] * x[1 2; 3 4] * h.rightenv[3 4 5; -3 -4] return y isa AbstractBlockTensorMap ? only(y) : y end \ No newline at end of file diff --git a/src/algorithms/groundstate/idmrg.jl b/src/algorithms/groundstate/idmrg.jl index 0ad55c23b..0b1e848fd 100644 --- a/src/algorithms/groundstate/idmrg.jl +++ b/src/algorithms/groundstate/idmrg.jl @@ -24,7 +24,7 @@ $(TYPEDFIELDS) alg_eigsolve::A = Defaults.alg_eigsolve() end -function find_groundstate(ost::AbstractInfiniteMPS, H, alg::IDMRG, envs = environments(ost, H)) +function find_groundstate(ost::InfiniteMPS, H, alg::IDMRG, envs = environments(ost, H)) ϵ::Float64 = calc_galerkin(ost, H, ost, envs) ψ = copy(ost) log = IterLog("IDMRG") @@ -131,7 +131,7 @@ $(TYPEDFIELDS) trscheme::TruncationStrategy end -function find_groundstate(ost::AbstractInfiniteMPS, H, alg::IDMRG2, envs = environments(ost, H)) +function find_groundstate(ost::InfiniteMPS, H, alg::IDMRG2, envs = environments(ost, H)) length(ost) < 2 && throw(ArgumentError("unit cell should be >= 2")) ϵ::Float64 = calc_galerkin(ost, H, ost, envs) diff --git a/src/algorithms/groundstate/vumps.jl b/src/algorithms/groundstate/vumps.jl index 418337047..99da8234c 100644 --- a/src/algorithms/groundstate/vumps.jl +++ b/src/algorithms/groundstate/vumps.jl @@ -46,7 +46,7 @@ struct VUMPSState{S, O, E} end function find_groundstate( - mps::AbstractInfiniteMPS, operator, alg::VUMPS, envs = environments(mps, operator) + mps::InfiniteMPS, operator, alg::VUMPS, envs = environments(mps, operator) ) return dominant_eigsolve(operator, mps, alg, envs; which = :SR) end diff --git a/src/environments/abstract_envs.jl b/src/environments/abstract_envs.jl index 97caefe03..cbdcafee1 100644 --- a/src/environments/abstract_envs.jl +++ b/src/environments/abstract_envs.jl @@ -81,8 +81,7 @@ function environment_alg( tol = Defaults.tol, maxiter = Defaults.maxiter, krylovdim = Defaults.krylovdim, verbosity = Defaults.VERBOSE_NONE ) - max_krylovdim = dim(left_virtualspace(above, 1)) * dim(left_virtualspace(below, 1)) - return GMRES(; tol, maxiter, krylovdim = min(max_krylovdim, krylovdim), verbosity) + return GMRES(; tol, maxiter, krylovdim, verbosity) end function environment_alg( ::Union{InfiniteQP, MultilineQP}, ::Union{InfiniteMPO, MultilineMPO}, @@ -91,4 +90,4 @@ function environment_alg( verbosity = Defaults.VERBOSE_NONE ) return GMRES(; tol, maxiter, krylovdim, verbosity) -end \ No newline at end of file +end diff --git a/src/states/abstractmps.jl b/src/states/abstractmps.jl index 2ac1cdeb8..dd756b2f4 100644 --- a/src/states/abstractmps.jl +++ b/src/states/abstractmps.jl @@ -15,7 +15,7 @@ const MPSTensor{S} = GenericMPSTensor{S, 2} # the usual mps tensors on which we """ MPSTensor([f, eltype], d::Int, left_D::Int, [right_D]::Int]) - MPSTensor([f, eltype], physicalspace::Union{S,CompositeSpace{S}}, + MPSTensor([f, eltype], physicalspace::Union{S,CompositeSpace{S}}, left_virtualspace::S, [right_virtualspace]::S) where {S<:ElementarySpace} Construct an `MPSTensor` with given physical and virtual spaces. @@ -93,11 +93,10 @@ Determine whether the given tensor is full rank, i.e. whether both the map from virtual space and the physical space to the right virtual space, and the map from the right virtual space and the physical space to the left virtual space are injective. """ -isfullrank(A::GenericMPSTensor; kwargs...) = isfullrank(space(A); kwargs...) -function isfullrank(V::TensorKit.TensorMapSpace; side = :both) - Vₗ = V[1] - Vᵣ = V[numind(V)] - P = ⊗(getindex.(Ref(V), 2:(numind(V) - 1))...) +function isfullrank(A::GenericMPSTensor; side = :both) + Vₗ = _firstspace(A) + Vᵣ = _lastspace(A) + P = ⊗(space.(Ref(A), 2:(numind(A) - 1))...) return if side === :both Vₗ ⊗ P ≿ Vᵣ' && Vₗ' ≾ P ⊗ Vᵣ elseif side === :right @@ -130,29 +129,6 @@ function makefullrank!(A::PeriodicVector{<:GenericMPSTensor}; alg = Defaults.alg return A end -function makefullrank!(virtualspaces::PeriodicVector{S}, physicalspaces::PeriodicVector{S}) where {S <: ElementarySpace} - haschanged = true - while haschanged - haschanged = false - for i in 1:length(virtualspaces) - Vmax = fuse(virtualspaces[i - 1], physicalspaces[i - 1]) - if !(virtualspaces[i] ≾ Vmax) - virtualspaces[i] = infimum(virtualspaces[i], Vmax) - haschanged = true - end - end - for i in reverse(1:length(virtualspaces)) - Vmax = fuse(dual(physicalspaces[i - 1]), virtualspaces[i]) - if !(virtualspaces[i - 1] ≾ Vmax) - virtualspaces[i - 1] = infimum(virtualspaces[i - 1], Vmax) - haschanged = true - end - end - end - - return virtualspaces -end - # Tensor accessors # ---------------- @doc """ @@ -201,7 +177,7 @@ TensorKit.sectortype(ψtype::Type{<:AbstractMPS}) = sectortype(site_type(ψtype) """ left_virtualspace(ψ::AbstractMPS, [pos=1:length(ψ)]) - + Return the virtual space of the bond to the left of sites `pos`. !!! warning @@ -244,4 +220,4 @@ physicalspace(ψ::AbstractMPS) = map(Base.Fix1(physicalspace, ψ), eachsite(ψ)) Return an iterator over the sites of the MPS `state`. """ -eachsite(ψ::AbstractMPS) = eachindex(ψ) \ No newline at end of file +eachsite(ψ::AbstractMPS) = eachindex(ψ) diff --git a/src/states/infinitemps.jl b/src/states/infinitemps.jl index 79620ebac..c9a4cf906 100644 --- a/src/states/infinitemps.jl +++ b/src/states/infinitemps.jl @@ -1,6 +1,5 @@ -abstract type AbstractInfiniteMPS <: AbstractMPS end """ - InfiniteMPS{A<:GenericMPSTensor,B<:MPSBondTensor} <: AbstractInfiniteMPS + InfiniteMPS{A<:GenericMPSTensor,B<:MPSBondTensor} <: AbtractMPS Type that represents an infinite Matrix Product State. @@ -44,7 +43,7 @@ tensors `As`, or a list of left-gauged tensors `ALs`. - `tol`: gauge fixing tolerance - `maxiter`: gauge fixing maximum iterations """ -struct InfiniteMPS{A <: GenericMPSTensor, B <: MPSBondTensor} <: AbstractInfiniteMPS +struct InfiniteMPS{A <: GenericMPSTensor, B <: MPSBondTensor} <: AbstractMPS AL::PeriodicVector{A} AR::PeriodicVector{A} C::PeriodicVector{B} diff --git a/src/states/ortho.jl b/src/states/ortho.jl index 3e935f456..a29942339 100644 --- a/src/states/ortho.jl +++ b/src/states/ortho.jl @@ -151,7 +151,7 @@ end Bring updated `AC` and `C` tensors back into a consistent set of left or right canonical tensors. This minimizes `∥AC_i - AL_i * C_i∥` or `∥AC_i - C_{i-1} * AR_i∥`. -The `alg` is passed on to `left_orth!` and `right_orth!`, and can be used to control the kind of +The `alg` is passed on to `left_orth!` and `right_orth!`, and can be used to control the kind of factorization used. By default, this is set to a (positive) QR/LQ, even though the optimal algorithm would use a polar decompositions instead, sacrificing a bit of performance for accuracy. diff --git a/src/utility/dynamictols.jl b/src/utility/dynamictols.jl index 61aae1f99..dc2a15781 100644 --- a/src/utility/dynamictols.jl +++ b/src/utility/dynamictols.jl @@ -58,11 +58,10 @@ end Update the tolerance of the algorithm `alg` based on the current iteration `iter` and the current error `ϵ`, where the new tolerance is given by - + new_tol = clamp(ϵ * alg.tol_factor / sqrt(iter), alg.tol_min, alg.tol_max) """ function updatetol(alg::DynamicTol, iter::Integer, ϵ::Real) - iter = max(iter, one(iter)) new_tol = clamp(ϵ * alg.tol_factor / sqrt(iter), alg.tol_min, alg.tol_max) return _updatetol(alg.alg, new_tol) end @@ -72,4 +71,4 @@ function _updatetol(alg, tol::Real) return Accessors.@set alg.tol = tol end -end \ No newline at end of file +end