diff --git a/docs/overrides.jl b/docs/overrides.jl index 36e0afcce..49e0d3967 100644 --- a/docs/overrides.jl +++ b/docs/overrides.jl @@ -6,7 +6,7 @@ function DocumenterVitepress.vitepress_config_transform(::VitepressOverrides, co occursin(marker, config) || error("DocumenterVitepress config has no default export to extend") return replace(config, marker => "const upstreamConfig = "; count = 1) * """ - import { withOverrides } from './overrides.mts' - export default withOverrides(upstreamConfig) - """ + import { withOverrides } from './overrides.mts' + export default withOverrides(upstreamConfig) + """ end diff --git a/docs/src/assets/mpskit.bib b/docs/src/assets/mpskit.bib index fd1352050..0ac909bac 100644 --- a/docs/src/assets/mpskit.bib +++ b/docs/src/assets/mpskit.bib @@ -773,6 +773,20 @@ @article{vanderstraeten2019 abstract = {In these lecture notes we give a technical overview of tangent-space methods for matrix product states in the thermodynamic limit. We introduce the manifold of uniform matrix product states, show how to compute different types of observables, and discuss the concept of a tangent space. We explain how to variationally optimize ground-state approximations, implement real-time evolution and describe elementary excitations for a given model Hamiltonian. Also, we explain how matrix product states approximate fixed points of one-dimensional transfer matrices. We show how all these methods can be translated to the language of continuous matrix product states for one-dimensional field theories. We conclude with some extensions of the tangent-space formalism and with an outlook to new applications.} } +@article{vanhecke2021cluster, + title = {Symmetric Cluster Expansions with Tensor Networks}, + author = {Vanhecke, Bram and Vanderstraeten, Laurens and Verstraete, Frank}, + year = {2021}, + month = feb, + journal = {Physical Review A}, + volume = {103}, + number = {2}, + pages = {L020402}, + publisher = {American Physical Society}, + doi = {10.1103/PhysRevA.103.L020402}, + url = {https://link.aps.org/doi/10.1103/PhysRevA.103.L020402} +} + @article{vanhecke2021, title = {Tangent-Space Methods for Truncating Uniform {{MPS}}}, author = {Vanhecke, Bram and Van Damme, Maarten and Haegeman, Jutho and Vanderstraeten, Laurens and Verstraete, Frank}, diff --git a/docs/src/changelog.md b/docs/src/changelog.md index c0fcf73b6..7efbec564 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -21,6 +21,10 @@ When releasing a new version, move the "Unreleased" changes to a new version sec ### Added +- `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)) + ### Changed ### Deprecated @@ -29,6 +33,8 @@ When releasing a new version, move the "Unreleased" changes to a new version sec ### Fixed +- `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)) + ### Performance ## [0.14.0](https://github.com/QuantumKitHub/MPSKit.jl/compare/v0.13.13...v0.14.0) - 2026-10-05 diff --git a/docs/src/man/algorithms.md b/docs/src/man/algorithms.md index 7b4a1e1fb..5d4029c87 100644 --- a/docs/src/man/algorithms.md +++ b/docs/src/man/algorithms.md @@ -123,10 +123,31 @@ Here, there are several algorithms available WI WII TaylorCluster +ClusterExpansion ``` See also [Time evolution accuracy](@ref). +#### Nonperturbative cluster expansion + +[`ClusterExpansion`](@ref) constructs evolution MPOs for finite and infinite +nearest-neighbor Hamiltonians by matching exact exponentials on clusters of up to +`N` sites [vanhecke2021cluster](@cite). Unlike [`TaylorCluster`](@ref), `N` specifies +the cluster size rather than a Taylor order. For nearest-neighbor interactions, +the one-step error is `O(dt^N)`; a finite chain is matched up to roundoff when +`N` reaches its length. + +```julia +U = make_time_mpo(H, dt, ClusterExpansion(; N = 3)) +ψ, envs, info = approximate(ψ, (U, ψ); trunc = truncrank(32)) +``` + +Matching small clusters to all orders can permit larger time steps and help grow +the bond dimension from weakly entangled states. The construction costs more than +the compact [`WII`](@ref) approximation, and increasing `N` requires exponentially +larger cluster exponentials and MPO bonds. Check convergence by varying `dt` and +`N`, as well as the MPS compression settings. + ## Excitations It might also be desirable to obtain information beyond the lowest energy state of a given system, and study the dispersion relation. diff --git a/src/MPSKit.jl b/src/MPSKit.jl index 4111185cc..ea209f449 100644 --- a/src/MPSKit.jl +++ b/src/MPSKit.jl @@ -36,7 +36,7 @@ export excitations export FiniteExcited, QuasiparticleAnsatz, ChepigaAnsatz, ChepigaAnsatz2 export time_evolve, timestep, timestep!, make_time_mpo export AlgorithmInfo, convergence_measure -export TDVP, TDVP2, BUG, WI, WII, TaylorCluster +export TDVP, TDVP2, BUG, WI, WII, TaylorCluster, ClusterExpansion export changebonds, changebonds! export VUMPSSvdCut, OptimalExpand, SvdCut, RandExpand, SketchedExpand export NoiseSchedule, FunctionalSchedule, ExponentialDecay, Warmup, DMRG3S @@ -174,6 +174,7 @@ include("algorithms/timestep/tdvp.jl") include("algorithms/timestep/bug.jl") include("algorithms/timestep/taylorcluster.jl") include("algorithms/timestep/wii.jl") +include("algorithms/timestep/clusterexpansion.jl") include("algorithms/timestep/integrators.jl") include("algorithms/timestep/time_evolve.jl") diff --git a/src/algorithms/timestep/clusterexpansion.jl b/src/algorithms/timestep/clusterexpansion.jl new file mode 100644 index 000000000..44d38b083 --- /dev/null +++ b/src/algorithms/timestep/clusterexpansion.jl @@ -0,0 +1,401 @@ +""" +$(TYPEDEF) + +Algorithm for constructing a finite or infinite time evolution MPO using the nonperturbative cluster expansion. +Exact exponentials are matched on clusters of up to `N` sites. +For nearest-neighbor Hamiltonians the error is `O(dt^N)`. + +# Fields + +$(TYPEDFIELDS) + +# Algorithm + +Exact evolution within clusters of up to `N` sites is included to all orders in +`dt`, together with products of disjoint clusters. The remaining error comes from +connected contributions extending beyond `N` sites: for nearest-neighbor +Hamiltonians the one-step error is `O(dt^N)`. Infinite Hamiltonians include clusters +crossing the unit-cell boundary; finite chains become exact up to roundoff when +`N` reaches the chain length. + +Compared with [`TaylorCluster`](@ref), which matches a specified Taylor order, +`ClusterExpansion` controls the spatial extent of the exact evolution. Its `N` +is a cluster size, not a Taylor order. Resumming smaller clusters can reduce the +error prefactor and permit larger time steps. Compared with the compact, local +construction of [`WII`](@ref), increasing `N` systematically incorporates larger +clusters, at the cost of exponentially growing cluster exponentials and MPO bond +dimensions. For full-rank uniform physical spaces of dimension `d`, the latter +scale as `1 + d^2 + ⋯ + d^(2 floor(N/2))`; mixed spaces or rank-deficient residuals +can require extra channels. + +This is useful for evolution with fewer, larger steps and for growing the MPS bond +dimension from weakly entangled states. Larger steps can also increase compression +costs, so check convergence and runtime by varying `dt` and `N`. MPS compression +introduces a separate error, and the approximate MPO need not be unitary. `tol` +only identifies shared SVD support: smaller singular components are retained in +complementary channels rather than discarded. + +# See also + +Used as the `algorithm` argument of [`make_time_mpo`](@ref). + +# References + +* [Vanhecke et al. Phys. Rev. A 103 (2021)](@cite vanhecke2021cluster) +""" +struct ClusterExpansion <: Algorithm + "maximum cluster size (positive integer, default: `2`)" + N::Int + "relative singular-value cutoff for shared SVD support (`0 < tol < 1`, default: `1.0e-12`)" + tol::Float64 + function ClusterExpansion(N::Integer; tol::Real = 1.0e-12) + N ≥ 1 || throw(ArgumentError("cluster size must be positive")) + isfinite(tol) && 0 < tol < 1 || throw(ArgumentError("tol must be between zero and one")) + return new(N, tol) + end +end +ClusterExpansion(; N::Integer = 2, kwargs...) = ClusterExpansion(N; kwargs...) + +function make_time_mpo(H::MPOHamiltonian, dt::Number, alg::ClusterExpansion) + # The only dispatch on a runtime cluster size. All subsequent stages and + # contraction shapes are determined by this Val parameter. + isempty(H) && throw(ArgumentError("Hamiltonian must contain at least one site")) + N = isfinite(H) ? min(alg.N, length(H)) : alg.N + return make_cluster_expansion_mpo(H, -im * dt, Val(N), alg.tol) +end + +function make_cluster_expansion_mpo(H::MPOHamiltonian, τ::Number, ::Val{N}, tol::Real) where {N} + storage_type = TensorKit.similarstoragetype(storagetype(H), promote_type(scalartype(H), typeof(τ))) + return make_cluster_expansion_mpo(H, τ, Val(N), tol, storage_type) +end +function make_cluster_expansion_mpo( + H::MPOHamiltonian, τ::Number, ::Val{N}, tol::Real, ::Type{A} + ) where {N, A} + tensors = map(parent(H)) do ham_tensor + return exponentiate_onsite_tensor(ham_tensor, τ, A) + end + # These cases are already exact products; avoid building channels from + # roundoff in differences of independently contracted exponentials. + if !(iszero(τ) || all(has_only_onsite_terms, parent(H))) + tensors = add_cluster_corrections(tensors, H, τ, Val(N), tol) + end + + return assemble_evolution_mpo(tensors) +end + +has_only_onsite_terms(ham_tensor) = size(ham_tensor, 1) ≤ 2 && size(ham_tensor, 4) ≤ 2 + +assemble_evolution_mpo(tensors::PeriodicVector) = InfiniteMPO(tensors) +function assemble_evolution_mpo(tensors::Vector) + tensors[1] = tensors[1][1, :, :, :] + tensors[end] = tensors[end][:, :, :, 1] + # Remove structurally unused levels only after all stages have finished. + return remove_orphans!(FiniteMPO(tensors); tol = 0) +end + +function exponentiate_onsite_tensor(ham_tensor, τ, ::Type{A}) where {A} + physical_space = physicalspace(ham_tensor) + levels = ⊞(oneunit(physical_space)) + blocktype = tensormaptype(spacetype(physical_space), 2, 2, A) + tensor = SparseBlockTensorMap{blocktype}(undef, levels ⊗ physical_space ← physical_space ⊗ levels) + onsite = removeunit(removeunit(ham_tensor[1, 1, 1, end], 4), 1) + tensor[1, 1, 1, 1] = add_util_leg(exp(τ * onsite)) + return tensor +end + +function add_cluster_corrections(tensors, H, τ, ::Val{N}, tol) where {N} + N == 1 && return tensors + previous = add_cluster_corrections(tensors, H, τ, Val(N - 1), tol) + starts = isfinite(H) ? (1:(length(H) - N + 1)) : (1:length(H)) + environment_length = (N - 1) ÷ 2 + level = N ÷ 2 + 1 + # Every translated window must see the completed previous stage. Compute + # all corrections before inserting any of them, including across the cell seam. + if isodd(N) + centers = [solve_center_correction(previous, start, evolution_cluster_residual(H, previous, start, τ, Val(N)), Val(N)) for start in starts] + for start in starts + previous[start + environment_length][level, 1, 1, level] = centers[start] + end + return previous + else + factors = [factor_two_site_correction(H, previous, start, τ, Val(N), tol) for start in starts] + expanded = expand_correction_bonds(previous, factors, Val(N)) + for start in starts + expanded[start + environment_length][level - 1, 1, 1, level] = factors[start][1] + expanded[start + environment_length + 1][level, 1, 1, level - 1] = factors[start][2] + end + return expanded + end +end + +function expand_correction_bonds(tensors::PeriodicVector, factors, ::Val{N}) where {N} + environment_length = (N - 1) ÷ 2 + bonds = PeriodicArray([right_virtualspace(factors[mod1(i - environment_length, length(tensors))][1]) for i in 1:length(tensors)]) + return PeriodicArray([append_virtual_spaces(tensors[i], bonds[i - 1], bonds[i]) for i in 1:length(tensors)]) +end + +function expand_correction_bonds(tensors::Vector, factors, ::Val{N}) where {N} + environment_length = (N - 1) ÷ 2 + empty_bond = zero(left_virtualspace(tensors[1])[1]) + # Only the middle bonds can support an N-site cluster. Empty levels keep + # the same level indices everywhere until the final boundary projection. + bonds = vcat( + fill(empty_bond, environment_length + 1), + [right_virtualspace(factor[1]) for factor in factors], + fill(empty_bond, environment_length + 1) + ) + return [append_virtual_spaces(tensors[i], bonds[i], bonds[i + 1]) for i in 1:length(tensors)] +end + +function append_virtual_spaces(tensor, left_bond, right_bond) + left_levels = left_virtualspace(tensor) ⊞ left_bond + right_levels = right_virtualspace(tensor) ⊞ right_bond + physical_space = physicalspace(tensor) + expanded = typeof(tensor)(undef, left_levels ⊗ physical_space ← physical_space ⊗ right_levels) + for (indices, block) in nonzero_pairs(tensor) + expanded[indices] = block + end + return expanded +end + +function evolution_cluster_residual(H, tensors, start::Int, τ, ::Val{N}) where {N} + boundary = size(H[start + N - 1], 4) + exact = exp(τ * contract_mpo_window(H, start, Val(N), boundary)) + return add_util_leg(exact - contract_mpo_window(tensors, start, Val(N), 1)) +end + +function factor_two_site_correction(H, tensors, start, τ, ::Val{N}, tol) where {N} + residual = evolution_cluster_residual(H, tensors, start, τ, Val(N)) + if N == 2 + left_factor, right_factor = factor_with_complementary_spaces(permute(residual, ((1, 2, 4), (3, 5, 6))), tol) + return permute(left_factor, ((1, 2), (3, 4))), permute(right_factor, ((1, 2), (3, 4))) + end + left_map, right_map = correction_environment_maps(tensors, start, Val(N)) + center = solve_center_correction(left_map, right_map, residual, Val(N)) + center_map = permute(center, ((1, 2, 4), (3, 5, 6))) + + # The minimum-norm solves restrict the center to the active environment + # subspaces. Complete only these directions, rather than the redundant + # virtual directions introduced by previous stages. Untruncated LQ/QR bases + # retain the entire physical environment support, including weak sectors. + _, left_basis = right_orth(left_map) + right_basis, _ = left_orth(right_map) + left_embedding = left_basis' ⊗ id(storagetype(center_map), codomain(center_map)[2] ⊗ codomain(center_map)[3]) + right_embedding = id(storagetype(center_map), domain(center_map)[1] ⊗ domain(center_map)[2]) ⊗ right_basis + left_factor, right_factor = factor_with_complementary_spaces(left_embedding' * center_map * right_embedding, tol) + return permute(left_embedding * left_factor, ((1, 2), (3, 4))), + permute(right_factor * right_embedding', ((1, 2), (3, 4))) +end + +""" + factor_with_complementary_spaces(tensor, tol) + +Factor `tensor` through its shared SVD support and separate left/right complementary +spaces. For exact rank `r_c`, the virtual multiplicity in sector `c` is +`m_c + n_c - r_c`, the minimum that gives both factors full environment support. +For exact null directions, the factors have the form + + A = [U√S γU⊥ 0], B = [√S Vᴴ; 0; γV⊥ᴴ]. + +At a finite rank cutoff, retain the remaining block in the complementary +channels, so the product still reconstructs `tensor` up to roundoff. +The scale `γ` balances these channels against the retained singular values. +""" +function factor_with_complementary_spaces(tensor::TensorMap, tol::Real) + tensor_norm = norm(tensor) + if iszero(tensor_norm) + # At zero rank, the two complementary spaces occupy disjoint channels. + left_space, right_space = fuse(codomain(tensor)), fuse(domain(tensor)) + storage_type = storagetype(tensor) + left_factor = catdomain( + isomorphism(storage_type, codomain(tensor) ← left_space), zeros(storage_type, codomain(tensor) ← right_space) + ) + right_factor = catcodomain( + zeros(storage_type, left_space ← domain(tensor)), isomorphism(storage_type, right_space ← domain(tensor)) + ) + return left_factor, right_factor + end + U, S, Vᴴ, _ = svd_trunc(tensor; trunc = trunctol(; rtol = tol, p = Inf)) + U_perp, V_perp = left_null(U), right_null(Vᴴ) + complement = U_perp' * tensor * V_perp' + complement_scale = sqrt(tensor_norm) + sqrt_values = sqrt(S) + left_factor = catdomain(catdomain(U * sqrt_values, complement_scale * U_perp), U_perp * complement / (2complement_scale)) + right_factor = catcodomain(catcodomain(sqrt_values * Vᴴ, complement * V_perp / (2complement_scale)), complement_scale * V_perp) + return left_factor, right_factor +end + +function contract_mpo_window(tensors, start::Int, ::Val{N}, boundary::Int) where {N} + N == 1 && return removeunit(removeunit(TensorMap(tensors[start][1, :, :, boundary]), 4), 1) + left = removeunit(TensorMap(tensors[start][1, :, :, :]), 1) + right = removeunit(TensorMap(tensors[start + N - 1][:, :, :, boundary]), 4) + sites = collect_bulk_tensors(tensors, start + 1, Val(N - 2)) + return contract_mpo_window(left, sites, right, Val(N)) +end + +function collect_bulk_tensors(tensors, start, ::Val{N}) where {N} + N == 0 && return () + return (TensorMap(tensors[start]), collect_bulk_tensors(tensors, start + 1, Val(N - 1))...) +end + +""" + contract_mpo_window(left, sites, right, ::Val{N}) + +Contract the open-boundary cluster, shown for `N = 4`. `L` and `R` are the +boundary tensors; `O₂` and `O₃` are the site-dependent bulk tensors. Negative labels are open +physical legs and positive labels are contracted virtual bonds. Domain legs +are above the tensors and codomain legs below. + +``` + -5 -6 -7 -8 + │ │ │ │ +┌─┴─┐ ┌─┴─┐ ┌─┴─┐ ┌─┴─┐ +│ L ├─1─┤O₂ ├─2─┤O₃ ├─3─┤ R │ +└─┬─┘ └─┬─┘ └─┬─┘ └─┬─┘ + │ │ │ │ + -1 -2 -3 -4 +``` +""" +@generated function contract_mpo_window( + left::AbstractTensorMap{<:Any, S, 1, 2}, sites::NTuple{K, AbstractTensorMap{<:Any, S, 2, 2}}, + right::AbstractTensorMap{<:Any, S, 2, 1}, ::Val{N} + ) where {S, N, K} + result = tensorexpr(:cluster, -(1:N), -((N + 1):(2N))) + first_site = tensorexpr(:left, -1, (-N - 1, 1)) + last_site = tensorexpr(:right, (N - 1, -N), -2N) + bulk = [tensorexpr(:(sites[$(i - 1)]), (i - 1, -i), (-N - i, i)) for i in 2:(N - 1)] + return macroexpand(@__MODULE__, :(return @plansor $result := *($first_site, $(bulk...), $last_site))) +end + +function left_correction_environment(tensors, start, ::Val{N}) where {N} + N == 1 && return tensors[start][1, 1, 1, 2] + left = left_correction_environment(tensors, start, Val(N - 1)) + return extend_left_environment(left, tensors[start + N - 1][N, 1, 1, N + 1]) +end + +function right_correction_environment(tensors, start, ::Val{N}) where {N} + N == 1 && return tensors[start][2, 1, 1, 1] + right = right_correction_environment(tensors, start + 1, Val(N - 1)) + return extend_right_environment(tensors[start][N + 1, 1, 1, N], right) +end + +""" + extend_left_environment(left, site) + +Append a site to the left environment, shown for `K = numout(left) = 2`. +For general `K`, each vertical leg of `left` is a bundle of `K - 1` physical +legs. The outer virtual legs remain open and bond `1` is contracted. Physical +domain legs are above the tensors and codomain legs below. + +``` + -4 -5 + │ │ + ┌──┴──┐ ┌──┴──┐ +-1───┤left ├─1─┤site ├───-6 + └──┬──┘ └──┬──┘ + │ │ + -2 -3 +``` +""" +@generated function extend_left_environment(left::AbstractTensorMap{<:Any, <:Any, K, K}, site::AbstractTensorMap) where {K} + result = tensorexpr(:next_left, -(1:(K + 1)), -((K + 2):(2K + 2))) + environment = tensorexpr(:left, -(1:K), (-((K + 2):(2K))..., 1)) + site_expression = tensorexpr(:site, (1, -K - 1), (-2K - 1, -2K - 2)) + return macroexpand(@__MODULE__, :(return @plansor $result := $environment * $site_expression)) +end + +""" + extend_right_environment(site, right) + +Prepend a site to the right environment, shown for `K = numout(right) = 2`. +For general `K`, each vertical leg of `right` is a bundle of `K - 1` physical +legs. The outer virtual legs remain open and bond `1` is contracted. Physical +domain legs are above the tensors and codomain legs below. + +``` + -4 -5 + │ │ + ┌──┴──┐ ┌──┴──┐ +-1───┤site ├─1─┤right├───-6 + └──┬──┘ └──┬──┘ + │ │ + -2 -3 +``` +""" +@generated function extend_right_environment(site::AbstractTensorMap, right::AbstractTensorMap{<:Any, <:Any, K, K}) where {K} + result = tensorexpr(:next_right, -(1:(K + 1)), -((K + 2):(2K + 2))) + site_expression = tensorexpr(:site, (-1, -2), (-K - 2, 1)) + environment = tensorexpr(:right, (1, -(3:(K + 1))...), -((K + 3):(2K + 2))) + return macroexpand(@__MODULE__, :(return @plansor $result := $site_expression * $environment)) +end + +function correction_environment_maps(tensors, start, ::Val{N}) where {N} + environment_length = (N - 1) ÷ 2 + left = left_correction_environment(tensors, start, Val(environment_length)) + right = right_correction_environment(tensors, start + N - environment_length, Val(environment_length)) + environment_indices = ntuple(identity, Val(2environment_length + 1)) + return permute(left, (environment_indices, (2environment_length + 2,))), permute(right, ((1,), environment_indices .+ 1)) +end + +function solve_center_correction(tensors, start::Int, residual::TensorMap, ::Val{N}) where {N} + left_map, right_map = correction_environment_maps(tensors, start, Val(N)) + return solve_center_correction(left_map, right_map, residual, Val(N)) +end + +""" + solve_center_correction(left_map, right_map, residual, ::Val{N}) + +Solve the environment equations from the left and right without forming inverse +maps. Completed environments have full row/column support in each symmetry +sector; rectangular solves choose the minimum-norm center. + +The grouped map layouts for `N = 3` are shown below. Domain bundles are above +and codomain bundles below. `pᵢ′` is the dual physical space, and `u_L`/`u_R` +are the trivial boundary spaces. + +First group the left physical legs of the residual into `B_L` and solve `L \\ B_L`: + +``` + V_L (p₂′, p₃′, p₂, p₃, u_R) + │ │ + ┌─────┴─────┐ ┌─────┴─────┐ + │ L │ │ B_L │ + └─────┬─────┘ └─────┬─────┘ + │ │ + (u_L, p₁, p₁′) (u_L, p₁, p₁′) +``` + +Regroup that solution into `B_R` and solve `B_R / R`: + +``` + (p₃′, p₃, u_R) (p₃′, p₃, u_R) + │ │ + ┌─────┴─────┐ ┌─────┴─────┐ + │ R │ │ B_R │ + └─────┬─────┘ └─────┬─────┘ + │ │ + V_R (V_L, p₂, p₂′) +``` + +Finally move `p₂′` back into the domain. The center has one physical leg per +side for odd `N`, or two for even `N`, and its two virtual legs. For larger +clusters, each outer physical space above is a bundle of `(N - 1) ÷ 2` legs. +""" +@generated function solve_center_correction( + left_map::TensorMap, right_map::TensorMap, residual::TensorMap, ::Val{N} + ) where {N} + environment_length = (N - 1) ÷ 2 + center_width = N - 2environment_length + # Group all left environment legs into the codomain for the left solve. + left_indices = (1, (2:(environment_length + 1))..., ((N + 2):(N + environment_length + 1))...) + remaining_indices = (((environment_length + 2):(N + 1))..., ((N + environment_length + 2):(2N + 2))...) + # After solving, group the right environment legs into the domain. + middle_indices = (1, (2:(center_width + 1))..., ((environment_length + center_width + 2):(environment_length + 2center_width + 1))...) + right_indices = (((center_width + 2):(environment_length + center_width + 1))..., ((environment_length + 2center_width + 2):(2environment_length + 2center_width + 2))...) + center_indices = (Tuple(1:(center_width + 1)), Tuple((center_width + 2):(2center_width + 2))) + return quote + left_solved = left_map \ permute(residual, $((left_indices, remaining_indices))) + solved = permute(left_solved, $((middle_indices, right_indices))) / right_map + return permute(solved, $center_indices) + end +end diff --git a/src/operators/mpohamiltonian.jl b/src/operators/mpohamiltonian.jl index 6a8c6e055..4b724a4eb 100644 --- a/src/operators/mpohamiltonian.jl +++ b/src/operators/mpohamiltonian.jl @@ -884,6 +884,8 @@ function isemptylevel(H::InfiniteMPOHamiltonian, i::Int) end function Base.convert(::Type{TensorMap}, H::FiniteMPOHamiltonian) + # A single site is both boundaries; strip its utility legs only once. + length(H) == 1 && return removeunit(removeunit(TensorMap(only(H)), 4), 1) L = removeunit(H[1], 1) R = removeunit(H[end], 4) M = Tuple(H[2:(end - 1)]) diff --git a/test/hamiltonian/finite.jl b/test/hamiltonian/finite.jl index e3c266249..14f7c9794 100644 --- a/test/hamiltonian/finite.jl +++ b/test/hamiltonian/finite.jl @@ -41,6 +41,10 @@ using TensorKit: ℙ end end + # A one-site Hamiltonian must not be contracted with itself. + Hsingle = FiniteMPOHamiltonian([V], 1 => O₁) + @test convert(TensorMap, Hsingle) ≈ O₁ + # check if constructor works by converting back to tensormap H1_tm = convert(TensorMap, H1) operators = vcat(fill(E, L - 1), O₁) diff --git a/test/timeevolution/clusterexpansion.jl b/test/timeevolution/clusterexpansion.jl new file mode 100644 index 000000000..43fc85756 --- /dev/null +++ b/test/timeevolution/clusterexpansion.jl @@ -0,0 +1,156 @@ +using Test +using MPSKit, TensorKit, LinearAlgebra, Random +using TensorKit: ℙ +using TensorKitTensors.SpinOperators: σˣ, σʸ, σᶻ + +# Independent dense windows also cover lengths that are not multiples of the unit cell. +function dense_mpo_window(O, start = 1, L = length(O)) + sites = [copy(O[start + i - 1]) for i in 1:L] + sites[1] = sites[1][1, :, :, :] + boundary = O isa MPOHamiltonian ? size(sites[end], 4) : 1 + sites[end] = sites[end][:, :, :, boundary] + window = O isa MPOHamiltonian ? FiniteMPOHamiltonian(sites) : MPSKit.DenseMPO(FiniteMPO(sites)) + return convert(TensorMap, window) +end + +function random_test_hamiltonian(rng, Ps; onsite_only = false, dimers = false) + infinite = Ps isa MPSKit.PeriodicVector + terms = Pair[] + for s in eachindex(Ps) + a = randn!(rng, zeros(ComplexF64, Ps[s] ← Ps[s])) + push!(terms, s => (a + a') / 2) + if (infinite || s < length(Ps)) && !onsite_only && (!dimers || isodd(s)) + b = randn!(rng, zeros(ComplexF64, Ps[s] ⊗ Ps[s + 1] ← Ps[s] ⊗ Ps[s + 1])) + push!(terms, (s, s + 1) => (b + b') / 2) + end + end + return infinite ? InfiniteMPOHamiltonian(Ps, terms...) : FiniteMPOHamiltonian(Ps, terms...) +end + +@testset "Nonperturbative cluster expansion" begin + X, Y, Z = σˣ(), σʸ(), σᶻ() + P = space(X, 1) + h = Z ⊗ Z + 0.3 * X ⊗ X + 0.2 * (X ⊗ Y + Y ⊗ X) + H = InfiniteMPOHamiltonian(PeriodicArray([P]), 1 => 0.7X, (1, 2) => h) + + @testset "Exact clusters and inference, N=$N" for N in 1:5 + U = @inferred MPSKit.make_cluster_expansion_mpo(H, -0.1im, Val(N), 1.0e-12) + @test U isa InfiniteMPO + @test dim(left_virtualspace(U, 1)) == sum(4^l for l in 0:(N ÷ 2)) + for L in 1:N + @test dense_mpo_window(U, 1, L) ≈ exp(-0.1im * dense_mpo_window(H, 1, L)) atol = 1.0e-10 + end + if N ≥ 3 + residual = @inferred MPSKit.evolution_cluster_residual(H, U, 1, -0.1im, Val(N)) + @inferred MPSKit.solve_center_correction(U, 1, residual, Val(N)) + end + if N in (3, 4) + larger = dense_mpo_window(U, 1, N + 1) + reflected = permute(larger, (Tuple((N + 1):-1:1), Tuple((2N + 2):-1:(N + 2)))) + @test larger ≈ reflected atol = 1.0e-10 + U2 = make_time_mpo(repeat(H, 2), 0.1, ClusterExpansion(N)) + @test dense_mpo_window(U2, 2, N + 1) ≈ larger atol = 1.0e-10 + end + end + + @testset "Leading error, N=$N" for N in 2:4 + for model in (H, open_boundary_conditions(H, 5)) + errors = [ + norm( + dense_mpo_window(make_time_mpo(model, dt, ClusterExpansion(N)), 1, N + 1) - + exp(-im * dt * dense_mpo_window(model, 1, N + 1)) + ) for dt in (0.04, 0.02) + ] + @test log2(errors[1] / errors[2]) ≥ N - 0.15 + end + end + + @testset "Finite boundaries and cluster-size cap, L=$L" for L in (1, 2, 3, 6) + finite = open_boundary_conditions(H, L) + U = make_time_mpo(finite, 0.1, ClusterExpansion(L + 3; tol = 0.5)) + @test U isa FiniteMPO + @test dim(left_virtualspace(U, 1)) == dim(right_virtualspace(U, L)) == 1 + @test dense_mpo_window(U) ≈ exp(-0.1im * convert(TensorMap, finite)) atol = 1.0e-10 + end + @testset "Finite/infinite consistency, N=$N" for N in (2, 3, 4, 5) + finite = make_time_mpo(open_boundary_conditions(H, 6), 0.1, ClusterExpansion(N)) + infinite = make_time_mpo(H, 0.1, ClusterExpansion(N)) + @test dense_mpo_window(finite) ≈ dense_mpo_window(infinite, 1, 6) atol = 1.0e-10 + end + + @test_throws ArgumentError ClusterExpansion(0) + @test_throws ArgumentError ClusterExpansion(2; tol = 0) + @test_throws ArgumentError ClusterExpansion(2; tol = Inf) + @test ClusterExpansion(; N = 3).N == 3 + empty = FiniteMPOHamiltonian(similar(parent(open_boundary_conditions(H, 2)), 0)) + @test_throws ArgumentError make_time_mpo(empty, 0.1, ClusterExpansion(2)) +end + +@testset "Mixed physical spaces and boundaries" begin + rng = MersenneTwister(25) + lattices = ( + [ℂ^2, ℂ^3], [ℂ^2, ℂ^3, ℂ^2], [ℝ^2, ℝ^3], [ℙ^2, ℙ^3], + [Rep[U₁](0 => 1, 1 => 1), dual(Rep[U₁](0 => 1, 1 => 1, 2 => 1))], + [dual(Rep[SU₂](1 // 2 => 1)), Rep[SU₂](1 => 1)], + [Vect[FermionParity](0 => 1, 1 => 1), Vect[FermionParity](0 => 2, 1 => 1)], + ) + for lattice in lattices, infinite in (true, false) + Ps = infinite ? PeriodicArray(lattice) : repeat(lattice, 3)[1:5] + H = random_test_hamiltonian(rng, Ps) + for N in (4, 5) + U = @inferred MPSKit.make_cluster_expansion_mpo(H, -0.1im, Val(N), 1.0e-12) + @test length(U) == length(H) + for s in eachindex(Ps) + @test physicalspace(U, s) == Ps[s] + (infinite || s < length(H)) && @test right_virtualspace(U, s) == left_virtualspace(U, s + 1) + width = infinite ? N : min(N, length(H) - s + 1) + @test dense_mpo_window(U, s, width) ≈ exp(-0.1im * dense_mpo_window(H, s, width)) atol = 1.0e-10 + end + end + end +end + +@testset "Time steps and rank-deficient environments" begin + rng = MersenneTwister(24) + Z = σᶻ() + P = space(Z, 1) + commuting = InfiniteMPOHamiltonian(PeriodicArray([P]), (1, 2) => Z ⊗ Z) + U = make_time_mpo(commuting, 0.1, ClusterExpansion(4)) + @test dense_mpo_window(U, 1, 5) ≈ exp(-0.1im * dense_mpo_window(commuting, 1, 5)) atol = 1.0e-10 + + for Ps in (PeriodicArray([ℂ^2, ℂ^3]), [ℂ^2, ℂ^3, ℂ^2, ℂ^3, ℂ^2]) + H = random_test_hamiltonian(rng, Ps) + onsite = random_test_hamiltonian(rng, Ps; onsite_only = true) + dimers = random_test_hamiltonian(rng, Ps; dimers = true) + for (model, dt) in ((onsite, 0.1), (dimers, 0.1), (H, 0.0), (H, -0.1im), (H, -0.1)) + U = make_time_mpo(model, dt, ClusterExpansion(5)) + width = isfinite(model) || model === H && !iszero(dt) ? 5 : 6 + starts = isfinite(model) ? (1:1) : eachindex(Ps) + for s in starts + @test dense_mpo_window(U, s, width) ≈ exp(-im * dt * dense_mpo_window(model, s, width)) atol = 1.0e-10 + end + model === onsite && @test dim(left_virtualspace(U, 1)) == 1 + end + # The SVD cutoff changes the channels, not the matched exponentials. + for N in (4, 5) + U = @inferred MPSKit.make_cluster_expansion_mpo(H, -0.1im, Val(N), 0.5) + starts = isfinite(H) ? (1:(length(H) - N + 1)) : eachindex(Ps) + for s in starts + @test dense_mpo_window(U, s, N) ≈ exp(-0.1im * dense_mpo_window(H, s, N)) atol = 1.0e-10 + end + end + end +end + +# Reconstruction alone does not ensure the factors retain all environment directions. +@testset "Completed SVD support" begin + L, R = ℂ^2, ℂ^3 + weak = TensorMap(ComplexF64[1 0 0; 0 1.0e-14 0], L ← R) + for C in (weak, zeros(ComplexF64, L ← R)) + A, B = @inferred MPSKit.factor_with_complementary_spaces(C, 1.0e-12) + @test norm(A * B - C) ≤ 1.0e-16 * norm(C) + @test A * pinv(A; rtol = 1.0e-12) ≈ id(storagetype(C), L) atol = 1.0e-10 + @test pinv(B; rtol = 1.0e-12) * B ≈ id(storagetype(C), R) atol = 1.0e-10 + @test dim(only(domain(A))) == (iszero(norm(C)) ? 5 : 4) + end +end