diff --git a/docs/src/changelog.md b/docs/src/changelog.md index 6005693f5..64bb570b0 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -24,9 +24,16 @@ When releasing a new version, move the "Unreleased" changes to a new version sec - `ClusterExpansion` time-evolution algorithm for `make_time_mpo`, constructing evolution MPOs for finite and infinite nearest-neighbor Hamiltonians through a nonperturbative cluster expansion that matches exact exponentials on clusters of up to `N` sites. ([#252](https://github.com/QuantumKitHub/MPSKit.jl/pull/252)) +- `changebonds(::InfiniteMPOHamiltonian, ::SvdCut)` compresses infinite Hamiltonians with the + almost-Schmidt decomposition of [Parker et al.](https://arxiv.org/abs/1909.06341), and + `left_canonicalize!`/`right_canonicalize!` bring an `InfiniteMPOHamiltonian` into canonical form. ([#539](https://github.com/QuantumKitHub/MPSKit.jl/pull/539)) ### Changed +- `left_canonicalize!`/`right_canonicalize!(::FiniteMPOHamiltonian)` now sweep over the whole chain + instead of a single site, and return the gauge transformations alongside the operator, matching + the infinite versions. ([#539](https://github.com/QuantumKitHub/MPSKit.jl/pull/539)) + ### Deprecated ### Removed @@ -35,6 +42,8 @@ When releasing a new version, move the "Unreleased" changes to a new version sec - `convert(TensorMap, H)` for single-site `FiniteMPOHamiltonian`s no longer contracts the site tensor with itself. ([#252](https://github.com/QuantumKitHub/MPSKit.jl/pull/252)) - `show` methods for finite MPS/MPOs no longer check unit space condition through `oneunit`, which is not defined for multifusion spaces. ([#538](https://github.com/QuantumKitHub/MPSKit.jl/pull/538)) +- The canonical forms of `MPOHamiltonian`s, and therefore `SvdCut` compression, are now correct when + interactions contain identity components or redundant channels. ([#539](https://github.com/QuantumKitHub/MPSKit.jl/pull/539)) ### Performance diff --git a/src/algorithms/changebonds/svdcut.jl b/src/algorithms/changebonds/svdcut.jl index 79ad973ba..611ff87a6 100644 --- a/src/algorithms/changebonds/svdcut.jl +++ b/src/algorithms/changebonds/svdcut.jl @@ -153,22 +153,15 @@ end changebonds(mpo::FiniteMPOHamiltonian, alg::SvdCut) = changebonds!(copy(mpo), alg) function changebonds!(H::FiniteMPOHamiltonian, alg::SvdCut) # orthogonality center to the left - for i in length(H):-1:2 - H = right_canonicalize!(H, i) - end + right_canonicalize!(H) # a Jordan bond has a start and a finish level on top of its interaction channels channels_before = [dim(left_virtualspace(H, i)) - 2 for i in 2:length(H)] - # swipe right alg_trunc = MatrixAlgebraKit.TruncatedAlgorithm(alg.alg_svd, alg.trunc) - for i in 1:(length(H) - 1) - H = left_canonicalize!(H, i; alg = MatrixAlgebraKit.LeftOrthViaSVD(alg_trunc)) - end - # swipe left -- TODO: do we really need this double sweep? - for i in length(H):-1:2 - H = right_canonicalize!(H, i; alg = MatrixAlgebraKit.RightOrthViaSVD(alg_trunc)) - end + left_canonicalize!(H; alg = MatrixAlgebraKit.LeftOrthViaSVD(alg_trunc)) + # TODO: do we really need this double sweep? + right_canonicalize!(H; alg = MatrixAlgebraKit.RightOrthViaSVD(alg_trunc)) emptied = findall( i -> channels_before[i] > 0 && dim(left_virtualspace(H, i + 1)) - 2 == 0, @@ -178,3 +171,55 @@ function changebonds!(H::FiniteMPOHamiltonian, alg::SvdCut) return H end + +changebonds(H::InfiniteMPOHamiltonian, alg::SvdCut) = changebonds!(copy(H), alg) +function changebonds!(H::InfiniteMPOHamiltonian, alg::SvdCut) + channels_before = [dim(right_virtualspace(H[i])) - 2 for i in 1:length(H)] + + # the interaction channels of the gauge between the right- and left-canonical forms + right_canonicalize!(H) + _, Cs = left_canonicalize!(H) + + Us = map(C -> first(svd_trunc(C; trunc = alg.trunc, alg = alg.alg_svd)), Cs) + for i in 1:length(H) + H[i] = _project_channels(Us[i - 1], H[i], Us[i]) + end + + emptied = findall( + i -> channels_before[i] > 0 && dim(right_virtualspace(H[i])) - 2 == 0, + eachindex(channels_before) + ) + isempty(emptied) || @warn "`SvdCut` discarded every interaction crossing bond(s) $emptied; those terms are gone from the compressed operator. Loosen `trunc` or check the scale of the input Hamiltonian." + + return H +end + +# Compute `diag(1, Uₗ', 1) * W * diag(1, Uᵣ, 1)`. +function _project_channels(Uₗ, W::JordanMPOTensor, Uᵣ) + P = physicalspace(W) + Vₗ, Vᵣ = _channelspace(Uₗ), _channelspace(Uᵣ) + WA, WB, WC = W.A, W.B, W.C + Uₗᴴ = adjoint(Uₗ) + + if dim(Vₗ) == 0 || dim(Vᵣ) == 0 + A = similar(WA, Vₗ ⊗ P ← P ⊗ Vᵣ) + else + @plansor A[l p; p' r] := Uₗᴴ[l; l'] * WA[l' p; p' r'] * Uᵣ[r'; r] + end + if dim(Vₗ) == 0 + B = similar(WB, Vₗ ⊗ P ← P) + else + @plansor B[l p; p'] := Uₗᴴ[l; l'] * WB[l' p; p'] + end + if dim(Vᵣ) == 0 + C = similar(WC, P ← P ⊗ Vᵣ) + else + @plansor C[p; p' r] := WC[p; p' r'] * Uᵣ[r'; r] + end + + unit = rightunitspace(P) + V = ⊞(unit, Vₗ, unit) ⊗ P ← P ⊗ ⊞(unit, Vᵣ, unit) + return JordanMPOTensor(V, A, B, C, W.D) +end + +_channelspace(U) = dim(U) == 0 ? SumSpace{spacetype(U)}() : domain(U)[1] diff --git a/src/operators/ortho.jl b/src/operators/ortho.jl index 9ea1f64a5..5a4779a51 100644 --- a/src/operators/ortho.jl +++ b/src/operators/ortho.jl @@ -2,75 +2,205 @@ _dense(t::SparseBlockTensorMap) = BlockTensorMap(t) _dense(t) = t +# Finite MPOHamiltonians +# ---------------------- + +""" + left_canonicalize!(H::FiniteMPOHamiltonian; alg = Defaults.alg_orth()) -> H, L + left_canonicalize!(H::InfiniteMPOHamiltonian; tol, maxiter, verbosity) -> H, L + +Bring `H` into left-canonical form in place, where the upper-left block `[1 C; 0 A]` of every site +has orthonormal columns with respect to the trace inner product normalized by the physical +dimension. The identity components of `C` are gauged away, and the decompositions respect the block +structure. + +Also returns the interaction-channel blocks `L` of the gauge transformations on every bond, such +that `diag(1, L[i - 1], 1) * H₀[i] = H[i] * diag(1, L[i], 1)` up to the identity components that +are gauged away, with `L[0]` the identity for a `FiniteMPOHamiltonian`. + +For a `FiniteMPOHamiltonian`, this is a single sweep of QR decompositions from left to right. Any +other (possibly truncating) algorithm `alg` is applied to the interaction channels only, as used by +[`SvdCut`](@ref). + +For an `InfiniteMPOHamiltonian`, the sweeps over the unit cell are repeated until the gauge +transformation is unitary, removing redundant interaction channels along the way. + +# Keyword arguments + +- `alg = Defaults.alg_orth()`: (finite) algorithm for the decomposition of the interaction channels +- `tol = Defaults.tolgauge`: (infinite) convergence tolerance on the deviation of the gauge + transformation from a unitary one, which is also the relative tolerance for discarding redundant + channels +- `maxiter = Defaults.maxiter`: (infinite) maximal number of sweeps over the unit cell +- `verbosity = VERBOSE_WARN`: (infinite) setting for how much information is displayed + +# References + +* [Parker et al. Phys. Rev. B 102 (2020)](@cite parker2020) +""" +function left_canonicalize!( + H::FiniteMPOHamiltonian; alg::MatrixAlgebraKit.AbstractAlgorithm = Defaults.alg_orth() + ) + Ls = map(1:(length(H) - 1)) do i + H[i], t, R = _left_orth_jordan(H[i]; alg) + H[i + 1] = _absorb_left_gauge(H[i], t, R, H[i + 1]) + return R + end + return H, Ls +end + +""" + right_canonicalize!(H::FiniteMPOHamiltonian; alg = Defaults.alg_orth()) -> H, L + right_canonicalize!(H::InfiniteMPOHamiltonian; tol, maxiter, verbosity) -> H, L + +Bring `H` into right-canonical form in place, where the lower-right block `[A B; 0 1]` of every site +has orthonormal rows with respect to the trace inner product normalized by the physical dimension. +The identity components of `B` are gauged away, and the decompositions respect the block structure. + +Also returns the interaction-channel blocks `L` of the gauge transformations on every bond, such +that `H₀[i] * diag(1, L[i], 1) = diag(1, L[i - 1], 1) * H[i]` up to the identity components that +are gauged away, with `L[length(H)]` the identity for a `FiniteMPOHamiltonian`. + +For a `FiniteMPOHamiltonian`, this is a single sweep of LQ decompositions from right to left. For an +`InfiniteMPOHamiltonian`, the sweeps over the unit cell are repeated until the gauge transformation +is unitary, removing redundant interaction channels along the way. + +See [`left_canonicalize!`](@ref) for the keyword arguments. +""" +function right_canonicalize!( + H::FiniteMPOHamiltonian; alg::MatrixAlgebraKit.AbstractAlgorithm = Defaults.alg_orth() + ) + Ls = map(length(H):-1:2) do i + R, t, H[i] = _right_orth_jordan(H[i]; alg) + H[i - 1] = _absorb_right_gauge(H[i - 1], R, t, H[i]) + return R + end + return H, reverse!(Ls) +end + +# Infinite MPOHamiltonians +# ------------------------ + function left_canonicalize!( - H::FiniteMPOHamiltonian, i::Int; - alg::MatrixAlgebraKit.AbstractAlgorithm = Defaults.alg_orth() + H::InfiniteMPOHamiltonian; + tol::Real = Defaults.tolgauge, maxiter::Int = Defaults.maxiter, verbosity::Int = VERBOSE_WARN ) - 1 ≤ i < length(H) || throw(ArgumentError("Bounds error in canonicalize")) + alg_trunc = MatrixAlgebraKit.TruncatedAlgorithm(Defaults.alg_svd(), trunctol(; rtol = tol)) + alg = MatrixAlgebraKit.LeftOrthViaSVD(alg_trunc) + Ls = PeriodicVector(Vector{Any}(nothing, length(H))) + return with_verbosity(; verbosity) do + log = IterLog("MPO LC") + @log_initialization loginit!(log, Inf) + for iter in 1:maxiter + local t, R + for i in 1:length(H) + H[i], t, R = _left_orth_jordan(H[i]; alg) + H[i + 1] = _absorb_left_gauge(H[i], t, R, H[i + 1]) + Ls[i] = isnothing(Ls[i]) || dim(R) == 0 ? R : R * Ls[i] + end + ϵ = _gauge_distance(t, R) + if ϵ <= tol + @log_convergence logfinish!(log, iter, ϵ) + break + elseif iter == maxiter + @log_nonconvergence logcancel!(log, iter, ϵ) + else + @log_iteration logiter!(log, iter, ϵ) + end + end + return H, map(identity, Ls) + end +end - W = H[i] +function right_canonicalize!( + H::InfiniteMPOHamiltonian; + tol::Real = Defaults.tolgauge, maxiter::Int = Defaults.maxiter, verbosity::Int = VERBOSE_WARN + ) + alg_trunc = MatrixAlgebraKit.TruncatedAlgorithm(Defaults.alg_svd(), trunctol(; rtol = tol)) + alg = MatrixAlgebraKit.RightOrthViaSVD(alg_trunc) + Ls = PeriodicVector(Vector{Any}(nothing, length(H))) + return with_verbosity(; verbosity) do + log = IterLog("MPO RC") + @log_initialization loginit!(log, Inf) + for iter in 1:maxiter + local t, R + for i in length(H):-1:1 + R, t, H[i] = _right_orth_jordan(H[i]; alg) + H[i - 1] = _absorb_right_gauge(H[i - 1], R, t, H[i]) + Ls[i - 1] = isnothing(Ls[i - 1]) || dim(R) == 0 ? R : Ls[i - 1] * R + end + ϵ = _gauge_distance(t, R) + if ϵ <= tol + @log_convergence logfinish!(log, iter, ϵ) + break + elseif iter == maxiter + @log_nonconvergence logcancel!(log, iter, ϵ) + else + @log_iteration logiter!(log, iter, ϵ) + end + end + return H, map(identity, Ls) + end +end + +# distance of the gauge transformation `[1 t 0; 0 R 0; 0 0 1]` (or its mirror) to a unitary one +function _gauge_distance(t, R) + codomain(R) == domain(R) || return Inf + return norm(t) + norm(R' * R - id(storagetype(R), domain(R))) +end + +# Implementation +# -------------- + +# Decompose `W = Q * G` with `G = [1 t 0; 0 R 0; 0 0 1]` and `Q` left-canonical. +function _left_orth_jordan(W::JordanMPOTensor; alg = Defaults.alg_orth()) P = physicalspace(W) _rightunit = rightunitspace(P) S = spacetype(W) d = sqrt(dim(P)) + WA, WB, WD = W.A, W.B, W.D - # block accessors recompute a fresh `SparseBlockTensorMap` on every access, so bind - # them once and reuse the locals throughout - WA, WB, WC, WD = W.A, W.B, W.C, W.D - - # orthogonalize second column against first - Wi = W[1, 1, 1, 1] - WI = removeunit(Wi, 1) - @plansor t[l; r] := conj(WI[p; p' l]) * WC[p; p' r] - # TODO: the following is currently broken due to a TensorKit bug - # @plansor C′[p; p' r] := WC[p; p' r] - WI[p; p' l] * t[l; r] - @plansor C′[p; p' r] := -WI[p; p' l] * t[l; r] - add!(C′, _dense(WC)) - - # QR of second column - if size(W, 1) == 1 - tmp = transpose(C′, ((2, 1), (3,))) - Q, R = left_orth!(tmp; alg) - - if dim(R) == 0 # fully truncated - V = _rightunit ⊞ _rightunit - Q1 = typeof(WA)(undef, SumSpace{S}() ⊗ P ← P ⊗ SumSpace{S}()) - Q2 = typeof(WC)(undef, P ← P ⊗ SumSpace{S}()) - else - V = ⊞(_rightunit, space(R, 1), _rightunit) - scale!(Q, d) - scale!(R, inv(d)) - Q1 = typeof(WA)(undef, SumSpace{S}() ⊗ P ← P ⊗ space(R, 1)) - Q2 = transpose(Q, ((2,), (1, 3))) - end - H[i] = JordanMPOTensor(codomain(W) ← physicalspace(W) ⊗ V, Q1, WB, Q2, WD) + # QR of the upper-left block `[1 C; 0 A]` with the identity as first column, which keeps all other + # columns of `Q` orthogonal to it, even if the block is rank-deficient + rows = 1:max(size(W, 1) - 1, 1) + V = BlockTensorMap(SparseBlockTensorMap(W)[rows, 1, 1, 1:(end - 1)]) + Q, R = left_orth!(transpose(V, ((3, 1, 2), (4,))); alg = Defaults.alg_orth()) + Π = _split_unit(domain(Q)[1], _rightunit, storagetype(W)) + Q = (Q * Π')[:, :, :, 2:2] + R = Π * R + t = scale!(R[1:1, 2:end], inv(d)) + R = R[2:2, 2:end] + + # other (possibly truncating) algorithms only act on the interaction channels + if alg != Defaults.alg_orth() + U, R = left_orth!(R; alg) + Q = Q * U + end + scale!(R, inv(d)) + + if dim(R) == 0 # fully truncated + V = _rightunit ⊞ _rightunit + Q1 = typeof(WA)(undef, codomain(WA) ← P ⊗ SumSpace{S}()) + Q2 = typeof(W.C)(undef, P ← P ⊗ SumSpace{S}()) else - tmp = transpose(cat(insertleftunit(C′, 1), WA; dims = 1), ((3, 1, 2), (4,))) - Q, R = left_orth!(tmp; alg) - - if dim(R) == 0 # fully truncated - V = _rightunit ⊞ _rightunit - Q1 = typeof(WA)(undef, codomain(WA) ← P ⊗ SumSpace{S}()) - Q2 = typeof(WC)(undef, P ← P ⊗ SumSpace{S}()) - else - scale!(Q, d) - scale!(R, inv(d)) - Q′ = transpose(Q, ((2, 3), (1, 4))) - Q1 = Q′[2:end, 1, 1, 1] - Q2 = removeunit(SparseBlockTensorMap(Q′[1:1, 1, 1, 1]), 1) - V = ⊞(_rightunit, right_virtualspace(Q′), _rightunit) - end - H[i] = JordanMPOTensor(codomain(W) ← P ⊗ V, Q1, WB, Q2, WD) + Q = scale!(transpose(Q, ((2, 3), (1, 4))), d) + Q1 = size(W, 1) == 1 ? typeof(WA)(undef, SumSpace{S}() ⊗ P ← P ⊗ right_virtualspace(Q)) : + SparseBlockTensorMap(Q[2:end, 1, 1, 1]) + Q2 = removeunit(SparseBlockTensorMap(Q[1:1, 1, 1, 1]), 1) + V = ⊞(_rightunit, right_virtualspace(Q), _rightunit) end + return JordanMPOTensor(codomain(W) ← P ⊗ V, Q1, WB, Q2, WD), t, R +end - # absorb into next site - W′ = H[i + 1] +# Compute `G * W′` for the gauge `G = [1 t 0; 0 R 0; 0 0 1]` obtained from `_left_orth_jordan`, +# where `Q` is the tensor to the left of `W′`. +function _absorb_left_gauge(Q::JordanMPOTensor, t, R, W′::JordanMPOTensor) W′A, W′B, W′C, W′D = W′.A, W′.B, W′.C, W′.D if size(W′, 4) > 1 && dim(R) != 0 @plansor A′[l p; p' r] := R[l; r'] * W′A[r' p; p' r] else - A′ = similar(W′A, right_virtualspace(H[i].A) ⊗ physicalspace(W′) ← domain(W′A)) + A′ = similar(W′A, right_virtualspace(Q.A) ⊗ physicalspace(W′) ← domain(W′A)) end if size(W′, 4) > 1 @@ -83,88 +213,67 @@ function left_canonicalize!( if dim(R) != 0 @plansor B′[l p; p'] := R[l; r] * W′B[r p; p'] else - B′ = similar(W′B, right_virtualspace(H[i].A) ⊗ physicalspace(W′) ← domain(W′B)) + B′ = similar(W′B, right_virtualspace(Q.A) ⊗ physicalspace(W′) ← domain(W′B)) end @plansor D′[l p; p'] := t[l; r] * W′B[r p; p'] D′ = add!(removeunit(D′, 1), _dense(W′D)) - H[i + 1] = JordanMPOTensor( - right_virtualspace(H[i]) ⊗ physicalspace(W′) ← domain(W′), - A′, B′, C′, D′ + return JordanMPOTensor( + right_virtualspace(Q) ⊗ physicalspace(W′) ← domain(W′), A′, B′, C′, D′ ) - return H end -function right_canonicalize!( - H::FiniteMPOHamiltonian, i::Int; - alg::MatrixAlgebraKit.AbstractAlgorithm = Defaults.alg_orth() - ) - 1 < i ≤ length(H) || throw(ArgumentError("Bounds error in canonicalize")) - - W = H[i] +# Decompose `W = G * Q` with `G = [1 0 0; 0 R t; 0 0 1]` and `Q` right-canonical. +function _right_orth_jordan(W::JordanMPOTensor; alg = Defaults.alg_orth()) P = physicalspace(W) _rightunit = rightunitspace(P) S = spacetype(W) d = sqrt(dim(P)) + WA, WC, WD = W.A, W.C, W.D - # block accessors recompute a fresh `SparseBlockTensorMap` on every access, so bind - # them once and reuse the locals throughout - WA, WB, WC, WD = W.A, W.B, W.C, W.D - - # orthogonalize second row against last - Wi = W[end, 1, 1, end] - WI = removeunit(Wi, 4) - @plansor t[l; r] := conj(WI[r p; p']) * WB[l p; p'] - # TODO: the following is currently broken due to a TensorKit bug - # @plansor B′[l p; p'] := WB[l p; p'] - WI[r p; p'] * t[l; r] - @plansor B′[l p; p'] := -WI[r p; p'] * t[l; r] - add!(B′, _dense(WB)) - - # LQ of second row - if size(W, 4) == 1 - tmp = transpose(B′, ((1,), (3, 2))) - R, Q = right_orth!(tmp; alg) - - if dim(R) == 0 - V = _rightunit ⊞ _rightunit - Q1 = typeof(WA)(undef, SumSpace{S}() ⊗ P ← P ⊗ SumSpace{S}()) - Q2 = typeof(WB)(undef, SumSpace{S}() ⊗ P ← P) - else - V = ⊞(_rightunit, space(Q, 1), _rightunit) - scale!(Q, d) - scale!(R, inv(d)) - Q1 = typeof(WA)(undef, space(Q, 1) ⊗ P ← P ⊗ SumSpace{S}()) - Q2 = transpose(Q, ((1, 3), (2,))) - end - H[i] = JordanMPOTensor(V ⊗ P ← domain(W), Q1, Q2, WC, WD) + # LQ of the lower-right block `[A B; 0 1]` with the identity as first row, which keeps all other + # rows of `Q` orthogonal to it, even if the block is rank-deficient + cols = min(2, size(W, 4)):size(W, 4) + V = SparseBlockTensorMap(W)[2:end, 1, 1, cols] + V = BlockTensorMap(cat(V[end:end, 1, 1, :], V[1:(end - 1), 1, 1, :]; dims = 1)) + R, Q = right_orth!(transpose(V, ((1,), (3, 4, 2))); alg = Defaults.alg_orth()) + Π = _split_unit(codomain(Q)[1], _rightunit, storagetype(W)) + Q = (Π * Q)[2:2, :, :, :] + R = R * Π' + t = scale!(R[2:end, 1:1], inv(d)) + R = R[2:end, 2:2] + + # other (possibly truncating) algorithms only act on the interaction channels + if alg != Defaults.alg_orth() + R, U = right_orth!(R; alg) + Q = U * Q + end + scale!(R, inv(d)) + + if dim(R) == 0 # fully truncated + V = _rightunit ⊞ _rightunit + Q1 = typeof(WA)(undef, SumSpace{S}() ⊗ P ← domain(WA)) + Q2 = typeof(W.B)(undef, SumSpace{S}() ⊗ P ← P) else - B′′ = insertleftunit(B′, 4) - tmp = transpose(cat(B′′, WA; dims = 4), ((1,), (3, 4, 2))) - R, Q = right_orth!(tmp; alg) - if dim(R) == 0 - V = _rightunit ⊞ _rightunit - Q1 = typeof(WA)(undef, SumSpace{S}() ⊗ P ← domain(WA)) - Q2 = typeof(WB)(undef, SumSpace{S}() ⊗ P ← P) - else - scale!(Q, d) - scale!(R, inv(d)) - Q′ = transpose(Q, ((1, 4), (2, 3))) - Q1 = SparseBlockTensorMap(Q′[1, 1, 1, 2:end]) - Q2 = removeunit(SparseBlockTensorMap(Q′[1, 1, 1, 1:1]), 4) - V = ⊞(_rightunit, left_virtualspace(Q′), _rightunit) - end - H[i] = JordanMPOTensor(V ⊗ P ← domain(W), Q1, Q2, WC, WD) + Q = scale!(transpose(Q, ((1, 4), (2, 3))), d) + Q1 = size(W, 4) == 1 ? typeof(WA)(undef, left_virtualspace(Q) ⊗ P ← P ⊗ SumSpace{S}()) : + SparseBlockTensorMap(Q[1, 1, 1, 1:(end - 1)]) + Q2 = removeunit(SparseBlockTensorMap(Q[1, 1, 1, end:end]), 4) + V = ⊞(_rightunit, left_virtualspace(Q), _rightunit) end + return R, t, JordanMPOTensor(V ⊗ P ← domain(W), Q1, Q2, WC, WD) +end - # absorb into previous site - W′ = H[i - 1] +# Compute `W′ * G` for the gauge `G = [1 0 0; 0 R t; 0 0 1]` obtained from `_right_orth_jordan`, +# where `Q` is the tensor to the right of `W′`. +function _absorb_right_gauge(W′::JordanMPOTensor, R, t, Q::JordanMPOTensor) W′A, W′B, W′C, W′D = W′.A, W′.B, W′.C, W′.D if size(W′, 1) > 1 && dim(R) != 0 @plansor A′[l p; p' r] := W′A[l p; p' r'] * R[r'; r] else - A′ = similar(W′A, codomain(W′A) ← physicalspace(W′A) ⊗ left_virtualspace(H[i].A)) + A′ = similar(W′A, codomain(W′A) ← physicalspace(W′A) ⊗ left_virtualspace(Q.A)) end if size(W′, 1) > 1 @@ -177,11 +286,32 @@ function right_canonicalize!( if dim(R) != 0 @plansor C′[p; p' r] := W′C[p; p' r'] * R[r'; r] else - C′ = similar(W′C, codomain(W′C) ← physicalspace(W′) ⊗ left_virtualspace(H[i].A)) + C′ = similar(W′C, codomain(W′C) ← physicalspace(W′) ⊗ left_virtualspace(Q.A)) end @plansor D′[p; p' r] := W′C[p; p' r'] * t[r'; r] D′ = add!(removeunit(D′, 3), _dense(W′D)) - H[i - 1] = JordanMPOTensor(codomain(W′) ← physicalspace(W′) ⊗ V, A′, B′, C′, D′) - return H + return JordanMPOTensor( + codomain(W′) ← physicalspace(W′) ⊗ left_virtualspace(Q), A′, B′, C′, D′ + ) +end + +# Isometry that splits the first unit-sector basis vector of `V = ⊞(X)` off into its own block. +function _split_unit(V::SumSpace, unit, ::Type{A}) where {A} + X = only(V) + c₀ = only(sectors(unit)) + rest = X ⊖ unit + T = scalartype(A) + TT = tensormaptype(spacetype(X), 1, 1, A) + Π₁ = TT(undef, unit ← X) + copyto!(block(Π₁, c₀), [one(T) zeros(T, 1, blockdim(X, c₀) - 1)]) + Π₂ = TT(undef, rest ← X) + for (c, b) in blocks(Π₂) + Id = Matrix{T}(LinearAlgebra.I, size(b, 1), size(b, 1)) + copyto!(b, c == c₀ ? [zeros(T, size(b, 1)) Id] : Id) + end + Π = BlockTensorMap{TT}(undef, (unit ⊞ rest) ← V) + Π[1, 1] = Π₁ + Π[2, 1] = Π₂ + return Π end diff --git a/test/algorithms/changebonds.jl b/test/algorithms/changebonds.jl index ab451e990..93b894b84 100644 --- a/test/algorithms/changebonds.jl +++ b/test/algorithms/changebonds.jl @@ -249,6 +249,37 @@ end @test expectation_value(ψ, H′) ≈ expectation_value(ψ, H) rtol = 1.0e-10 end +@testset "SvdCut on InfiniteMPOHamiltonian $(spacetype(pspace))" for (pspace, Dspace) in spacelist + Random.seed!(2468) + nn = rand(ComplexF64, pspace * pspace, pspace * pspace) + nn += nn' + L = 3 + lattice = PeriodicVector(fill(pspace, L)) + H0 = InfiniteMPOHamiltonian(lattice, (i, i + 1) => nn for i in 1:L) + + InfiniteMPOHamiltonian(lattice, (i, i + 2) => 0.5^i * nn for i in 1:L) + # adding the terms twice duplicates every interaction channel + H = H0 + 2 * H0 + ψ = InfiniteMPS(rand, ComplexF64, fill(pspace, L), fill(Dspace, L)) + + H′ = changebonds(H, SvdCut(; trunc = trunctol(; atol = 1.0e-12))) + @test all(i -> dim(right_virtualspace(H′[i])) ≤ dim(right_virtualspace(H0[i])), 1:L) + @test expectation_value(ψ, H′) ≈ expectation_value(ψ, H) rtol = 1.0e-10 + + # exponentially decaying interactions, written out as explicit terms + Z = rand(ComplexF64, pspace, pspace) + Z += Z' + H_exp = InfiniteMPOHamiltonian(PeriodicVector([pspace]), (1, 1 + r) => 0.5^r * (Z ⊗ Z) for r in 1:20) + H_exp′ = changebonds(H_exp, SvdCut(; trunc = trunctol(; rtol = 1.0e-5))) + @test dim(right_virtualspace(H_exp′[1])) < dim(right_virtualspace(H_exp[1])) + ψ1 = InfiniteMPS(rand, ComplexF64, [pspace], [Dspace]) + @test expectation_value(ψ1, H_exp′) ≈ expectation_value(ψ1, H_exp) rtol = 1.0e-4 + + # purely on-site terms have no channels to compress + H_onsite = InfiniteMPOHamiltonian(PeriodicVector([pspace]), (1,) => Z) + @test expectation_value(ψ1, changebonds(H_onsite, SvdCut(; trunc = truncrank(0)))) ≈ + expectation_value(ψ1, H_onsite) +end + # Regression: RandExpand used to be able to introduce a NaN entanglement entropy. @testset "RandExpand does not introduce NaN entropy" begin ψ = InfiniteMPS([ℂ^2], [ℂ^5]) diff --git a/test/hamiltonian/finite.jl b/test/hamiltonian/finite.jl index 14f7c9794..3fe4cd3ef 100644 --- a/test/hamiltonian/finite.jl +++ b/test/hamiltonian/finite.jl @@ -113,8 +113,7 @@ using TensorKit: ℙ end end -@testset "Finite MPOHamiltonian canonicalization with decoupled sites" begin - # regression test for https://github.com/QuantumKitHub/MPSKit.jl/issues/522 +@testset "Finite MPOHamiltonian canonicalization" begin T = ComplexF64 for V in (ℂ^2, U1Space(-1 => 1, 0 => 1, 1 => 1), Rep[SU₂](1 // 2 => 1)) O₁ = randn(T, V, V) @@ -123,33 +122,56 @@ end O₂ += O₂' E = id(storagetype(O₁), domain(O₁)) + # regression test for https://github.com/QuantumKitHub/MPSKit.jl/issues/522 # bulk block of the middle site is empty on the left but not on the right H_right = FiniteMPOHamiltonian(fill(V, 3), (1,) => O₁, (2, 3) => O₂) H_right_tm = O₁ ⊗ E ⊗ E + E ⊗ O₂ - H = MPSKit.right_canonicalize!(copy(H_right), 2) + H, = MPSKit.right_canonicalize!(copy(H_right)) @test convert(TensorMap, H) ≈ H_right_tm # bulk block of the middle site is empty on the right but not on the left H_left = FiniteMPOHamiltonian(fill(V, 3), (1, 2) => O₂, (3,) => O₁) H_left_tm = O₂ ⊗ E + E ⊗ E ⊗ O₁ - H = MPSKit.left_canonicalize!(copy(H_left), 2) + H, = MPSKit.left_canonicalize!(copy(H_left)) @test convert(TensorMap, H) ≈ H_left_tm L = 5 H5 = FiniteMPOHamiltonian(fill(V, L), (1, 2) => O₂, (3,) => O₁, (4, 5) => O₂) H5_tm = O₂ ⊗ E ⊗ E ⊗ E + E ⊗ E ⊗ O₁ ⊗ E ⊗ E + E ⊗ E ⊗ E ⊗ O₂ - H = copy(H5) - for i in L:-1:2 - MPSKit.right_canonicalize!(H, i) - end + H, = MPSKit.right_canonicalize!(copy(H5)) @test convert(TensorMap, H) ≈ H5_tm - for i in 1:(L - 1) - MPSKit.left_canonicalize!(H, i) - end + H, = MPSKit.left_canonicalize!(H) @test convert(TensorMap, H) ≈ H5_tm H5_trunc = changebonds(H5, SvdCut(; trunc = truncrank(16))) @test convert(TensorMap, H5_trunc) ≈ H5_tm + + # canonical forms, also when `C`/`B` have identity components and the blocks are rank-deficient + L = 4 + H = FiniteMPOHamiltonian(fill(V, L), (i, i + 1) => O₂ for i in 1:(L - 1)) + H_tm = convert(TensorMap, H) + + HL, Ls = MPSKit.left_canonicalize!(copy(H)) + @test length(Ls) == L - 1 + @test all(i -> codomain(Ls[i])[1] == right_virtualspace(HL[i].A), 1:(L - 1)) + @test convert(TensorMap, HL) ≈ H_tm + for W in HL[2:(L - 1)] + Vb = MPSKit.SparseBlockTensorMap(W)[1:(end - 1), 1, 1, 1:(end - 1)] + @plansor G[r; r′] := conj(Vb[l p; p′ r]) * Vb[l p; p′ r′] + G = TensorMap(G) / dim(physicalspace(W)) + @test G ≈ id(domain(G)) + end + + HR, Rs = MPSKit.right_canonicalize!(copy(H)) + @test length(Rs) == L - 1 + @test all(i -> domain(Rs[i])[1] == left_virtualspace(HR[i + 1].A), 1:(L - 1)) + @test convert(TensorMap, HR) ≈ H_tm + for W in HR[2:(L - 1)] + Vb = MPSKit.SparseBlockTensorMap(W)[2:end, 1, 1, 2:end] + @plansor G[l; l′] := Vb[l p; p′ r] * conj(Vb[l′ p; p′ r]) + G = TensorMap(G) / dim(physicalspace(W)) + @test G ≈ id(codomain(G)) + end end end diff --git a/test/hamiltonian/infinite.jl b/test/hamiltonian/infinite.jl index f44ed42b4..4dc2f4c51 100644 --- a/test/hamiltonian/infinite.jl +++ b/test/hamiltonian/infinite.jl @@ -68,3 +68,36 @@ end @test expectation_value(ψ1, Hreal_long_range + Hcomplex_local) ≈ expectation_value(ψ1, Hreal_long_range) + expectation_value(ψ1, Hcomplex_local) atol = 1.0e-10 end + +@testset "InfiniteMPOHamiltonian canonical forms $(sectortype(pspace))" for (pspace, Dspace) in zip(pspaces, vspaces) + O₂ = rand(ComplexF64, pspace^2, pspace^2) + O₂ += O₂' + O₃ = rand(ComplexF64, pspace^3, pspace^3) + O₃ += O₃' + lattice = PeriodicVector(fill(pspace, 2)) + H = InfiniteMPOHamiltonian(lattice, (i, i + 1) => O₂ for i in 1:2) + + InfiniteMPOHamiltonian(lattice, (i, i + 1, i + 2) => O₃ for i in 1:2) + + ψ = InfiniteMPS(rand, ComplexF64, fill(pspace, 2), fill(Dspace, 2)) + E = expectation_value(ψ, H) + + # the upper-left block `[1 C; 0 A]` has orthonormal columns + HL, = MPSKit.left_canonicalize!(copy(H)) + @test expectation_value(ψ, HL) ≈ E + for W in parent(HL) + V = MPSKit.SparseBlockTensorMap(W)[1:(end - 1), 1, 1, 1:(end - 1)] + @plansor G[r; r′] := conj(V[l p; p′ r]) * V[l p; p′ r′] + G = TensorMap(G) / dim(physicalspace(W)) + @test G ≈ id(domain(G)) + end + + # the lower-right block `[A B; 0 1]` has orthonormal rows + HR, = MPSKit.right_canonicalize!(copy(H)) + @test expectation_value(ψ, HR) ≈ E + for W in parent(HR) + V = MPSKit.SparseBlockTensorMap(W)[2:end, 1, 1, 2:end] + @plansor G[l; l′] := V[l p; p′ r] * conj(V[l′ p; p′ r]) + G = TensorMap(G) / dim(physicalspace(W)) + @test G ≈ id(codomain(G)) + end +end