From 28af89ddc64ccc39f6c4b4825458038ba4dff42e Mon Sep 17 00:00:00 2001 From: Ben Corbett Date: Tue, 19 Nov 2024 18:09:56 -0700 Subject: [PATCH 1/3] Moved changes over from ITensor. --- src/mpo.jl | 45 +++++++++++++++++++++++++++++++------------ test/base/test_mpo.jl | 24 ++++++++++++----------- 2 files changed, 46 insertions(+), 23 deletions(-) diff --git a/src/mpo.jl b/src/mpo.jl index 9acd561b..2e085724 100644 --- a/src/mpo.jl +++ b/src/mpo.jl @@ -630,6 +630,24 @@ function ITensors.contract(A::MPO, ψ::MPS; alg=nothing, method=alg, kwargs...) return contract(Algorithm(alg), A, ψ; kwargs...) end +function zipup_docstring(isMPOMPO::Bool)::String + rhsTypeString = isMPOMPO ? "MPO" : "MPS" + rhsString = isMPOMPO ? "B" : "ψ" + return """ - "zipup": The MPO and $rhsTypeString tensors are contracted then truncated at each site without enforcing + the appropriate orthogonal gauge. Once this sweep is complete a call to `truncate!` occurs. + Because the initial truncation is not locally optimal it is recommended to use a loose + `cutoff` and `maxdim` and then pass the desired truncation parameters to the locally optimal + `truncate!` sweep via the additional keyword argument `truncate_kwargs`. + A set of parameters suggested in [^Paeckel2019] is + `contract(A, $rhsString; method="zipup", cutoff=cutoff / 10, maxdim=2 * maxdim, truncate_kwargs=(; cutoff, maxdim))`.""" +end + +function Paeckel2019_citation_docstring()::String + return """ + [^Paeckel2019]: Time-evolution methods for matrix-product states. Sebastian Paeckel et al. [arXiv:1901.05824](https://arxiv.org/abs/1901.05824) + """ +end + contract_mpo_mps_doc = """ contract(ψ::MPS, A::MPO; kwargs...) -> MPS *(::MPS, ::MPO; kwargs...) -> MPS @@ -663,13 +681,15 @@ Choose the method with the `method` keyword, for example - `mindim::Int=1`: the minimal bond dimension of the resulting MPS. - `normalize::Bool=false`: whether or not to normalize the resulting MPS. - `method::String="densitymatrix"`: the algorithm to use for the contraction. - Currently the options are "densitymatrix", where the network formed by the - MPO and MPS is squared and contracted down to a density matrix which is - diagonalized iteratively at each site, and "naive", where the MPO and MPS - tensor are contracted exactly at each site and then a truncation of the - resulting MPS is performed. + - "densitymatrix": The network formed by the MPO and MPS is squared and contracted down to + a density matrix which is diagonalized iteratively at each site. + - "naive": The MPO and MPS tensor are contracted exactly at each site and then a truncation + of the resulting MPS is performed. + $(zipup_docstring(false)) See also [`apply`](@ref). + +$(Paeckel2019_citation_docstring()) """ @doc """ @@ -823,10 +843,11 @@ end function ITensors.contract( ::Algorithm"zipup", A::MPO, - B::MPO; + B::AbstractMPS; cutoff=1e-14, maxdim=maxlinkdim(A) * maxlinkdim(B), mindim=1, + truncate_kwargs=(; cutoff, maxdim, mindim), kwargs..., ) if hassameinds(siteinds, A, B) @@ -838,13 +859,13 @@ function ITensors.contract( N != length(B) && throw(DimensionMismatch("lengths of MPOs A ($N) and B ($(length(B))) do not match")) # Special case for a single site - N == 1 && return MPO([A[1] * B[1]]) + N == 1 && return typeof(B)([A[1] * B[1]]) A = orthogonalize(A, 1) B = orthogonalize(B, 1) A = sim(linkinds, A) sA = siteinds(uniqueinds, A, B) sB = siteinds(uniqueinds, B, A) - C = MPO(N) + C = typeof(B)(N) lCᵢ = Index[] R = ITensor(true) for i in 1:(N - 2) @@ -875,7 +896,7 @@ function ITensors.contract( mindim, kwargs..., ) - truncate!(C; kwargs...) + truncate!(C; truncate_kwargs...) return C end @@ -946,16 +967,16 @@ C = apply(A, B; alg="naive", truncate=false) in general you should set a `cutoff` value. - `maxdim::Int=maxlinkdim(A) * maxlinkdim(B))`: the maximal bond dimension of the results MPS. - `mindim::Int=1`: the minimal bond dimension of the resulting MPS. -- `alg="zipup"`: Either `"zipup"` or `"naive"`. `"zipup"` contracts pairs of - site tensors and truncates with SVDs in a sweep across the sites, while `"naive"` - first contracts pairs of tensor exactly and then truncates at the end if `truncate=true`. - `truncate=true`: Enable or disable truncation. If `truncate=false`, ignore other truncation parameters like `cutoff` and `maxdim`. This is most relevant for the `"naive"` version, if you just want to contract the tensors pairwise exactly. This can be useful if you are contracting MPOs that have diverging norms, such as MPOs originating from sums of local operators. + $(zipup_docstring(true)) See also [`apply`](@ref) for details about the arguments available. + +$(Paeckel2019_citation_docstring()) """ @doc """ diff --git a/test/base/test_mpo.jl b/test/base/test_mpo.jl index 047e7386..1ecca7d8 100644 --- a/test/base/test_mpo.jl +++ b/test/base/test_mpo.jl @@ -242,14 +242,14 @@ end @test_throws DimensionMismatch error_contract(phi, K, badpsi) end - @testset "contract" begin + @testset "contract(::MPO, ::MPS)" for method in ["densitymatrix", "naive", "zipup"] phi = random_mps(sites) K = random_mpo(sites) @test maxlinkdim(K) == 1 psi = random_mps(sites) - psi_out = contract(K, psi; maxdim=1) + psi_out = contract(K, psi; method, maxdim=1) @test inner(phi', psi_out) ≈ inner(phi', K, psi) - psi_out = contract(psi, K; maxdim=1) + psi_out = contract(psi, K; method, maxdim=1) @test inner(phi', psi_out) ≈ inner(phi', K, psi) psi_out = psi * K @test inner(phi', psi_out) ≈ inner(phi', K, psi) @@ -257,7 +257,7 @@ end badsites = [Index(2, "Site") for n in 1:(N + 1)] badpsi = random_mps(badsites) - @test_throws DimensionMismatch contract(K, badpsi) + @test_throws DimensionMismatch contract(K, badpsi; method) # make bigger random MPO... for link_dim in 2:5 @@ -289,7 +289,9 @@ end orthogonalize!(psi, 1; maxdim=link_dim) orthogonalize!(K, 1; maxdim=link_dim) orthogonalize!(phi, 1; normalize=true, maxdim=link_dim) - psi_out = contract(deepcopy(K), deepcopy(psi); maxdim=10 * link_dim, cutoff=0.0) + psi_out = contract( + deepcopy(K), deepcopy(psi); method, maxdim=10 * link_dim, cutoff=0.0 + ) @test inner(phi', psi_out) ≈ inner(phi', K, psi) end end @@ -356,14 +358,14 @@ end @test maxlinkdim(H) ≤ maxlinkdim(H₁) + maxlinkdim(H₂) end - @testset "contract(::MPO, ::MPO)" begin + @testset "contract(::MPO, ::MPO)" for alg in ["naive", "zipup"] psi = random_mps(sites) K = random_mpo(sites) L = random_mpo(sites) @test maxlinkdim(K) == 1 @test maxlinkdim(L) == 1 - KL = contract(prime(K), L; maxdim=1) - psi_kl_out = contract(prime(K), contract(L, psi; maxdim=1); maxdim=1) + KL = contract(prime(K), L; alg, maxdim=1) + psi_kl_out = contract(prime(K), contract(L, psi; alg, maxdim=1); alg, maxdim=1) @test inner(psi'', KL, psi) ≈ inner(psi'', psi_kl_out) atol = 5e-3 # where both K and L have differently labelled sites @@ -375,15 +377,15 @@ end replaceind!(K[ii], sites[ii]', othersitesk[ii]) replaceind!(L[ii], sites[ii]', othersitesl[ii]) end - KL = contract(K, L; maxdim=1) + KL = contract(K, L; alg, maxdim=1) psik = random_mps(othersitesk) psil = random_mps(othersitesl) - psi_kl_out = contract(K, contract(L, psil; maxdim=1); maxdim=1) + psi_kl_out = contract(K, contract(L, psil; alg, maxdim=1); alg, maxdim=1) @test inner(psik, KL, psil) ≈ inner(psik, psi_kl_out) atol = 5e-3 badsites = [Index(2, "Site") for n in 1:(N + 1)] badL = random_mpo(badsites) - @test_throws DimensionMismatch contract(K, badL) + @test_throws DimensionMismatch contract(K, badL; alg) end @testset "*(::MPO, ::MPO)" begin From 7167279432f10eb2e125b8959197d09c539a7b2b Mon Sep 17 00:00:00 2001 From: Ben Corbett Date: Mon, 8 Dec 2025 14:58:06 -0800 Subject: [PATCH 2/3] Formatting changes from runic. --- src/mpo.jl | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/src/mpo.jl b/src/mpo.jl index 3b04da07..e125efd6 100644 --- a/src/mpo.jl +++ b/src/mpo.jl @@ -631,21 +631,21 @@ function ITensors.contract(A::MPO, ψ::MPS; alg = nothing, method = alg, kwargs. end function zipup_docstring(isMPOMPO::Bool)::String - rhsTypeString = isMPOMPO ? "MPO" : "MPS" - rhsString = isMPOMPO ? "B" : "ψ" - return """ - "zipup": The MPO and $rhsTypeString tensors are contracted then truncated at each site without enforcing - the appropriate orthogonal gauge. Once this sweep is complete a call to `truncate!` occurs. - Because the initial truncation is not locally optimal it is recommended to use a loose - `cutoff` and `maxdim` and then pass the desired truncation parameters to the locally optimal - `truncate!` sweep via the additional keyword argument `truncate_kwargs`. - A set of parameters suggested in [^Paeckel2019] is - `contract(A, $rhsString; method="zipup", cutoff=cutoff / 10, maxdim=2 * maxdim, truncate_kwargs=(; cutoff, maxdim))`.""" + rhsTypeString = isMPOMPO ? "MPO" : "MPS" + rhsString = isMPOMPO ? "B" : "ψ" + return """ - "zipup": The MPO and $rhsTypeString tensors are contracted then truncated at each site without enforcing + the appropriate orthogonal gauge. Once this sweep is complete a call to `truncate!` occurs. + Because the initial truncation is not locally optimal it is recommended to use a loose + `cutoff` and `maxdim` and then pass the desired truncation parameters to the locally optimal + `truncate!` sweep via the additional keyword argument `truncate_kwargs`. + A set of parameters suggested in [^Paeckel2019] is + `contract(A, $rhsString; method="zipup", cutoff=cutoff / 10, maxdim=2 * maxdim, truncate_kwargs=(; cutoff, maxdim))`.""" end function Paeckel2019_citation_docstring()::String - return """ - [^Paeckel2019]: Time-evolution methods for matrix-product states. Sebastian Paeckel et al. [arXiv:1901.05824](https://arxiv.org/abs/1901.05824) - """ + return """ + [^Paeckel2019]: Time-evolution methods for matrix-product states. Sebastian Paeckel et al. [arXiv:1901.05824](https://arxiv.org/abs/1901.05824) + """ end contract_mpo_mps_doc = """ @@ -847,7 +847,7 @@ function ITensors.contract( cutoff = 1.0e-14, maxdim = maxlinkdim(A) * maxlinkdim(B), mindim = 1, - truncate_kwargs=(; cutoff, maxdim, mindim), + truncate_kwargs = (; cutoff, maxdim, mindim), kwargs..., ) if hassameinds(siteinds, A, B) From 7ad7173eeb68b28a2e2edb28d5499206bec877a4 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Mon, 23 Mar 2026 16:08:34 -0400 Subject: [PATCH 3/3] Edits to PR #97 --- Project.toml | 2 +- src/mpo.jl | 70 +++++++++++++++++++++++++------------------ test/base/test_mpo.jl | 17 ++++++----- 3 files changed, 52 insertions(+), 37 deletions(-) diff --git a/Project.toml b/Project.toml index cfbf7cff..b96d1682 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "ITensorMPS" uuid = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2" -version = "0.3.44" +version = "0.3.45" authors = ["Matthew Fishman ", "Miles Stoudenmire "] [workspace] diff --git a/src/mpo.jl b/src/mpo.jl index 0e1833ab..34ec5d2d 100644 --- a/src/mpo.jl +++ b/src/mpo.jl @@ -637,19 +637,18 @@ function ITensors.contract(A::MPO, ψ::MPS; alg = nothing, method = alg, kwargs. return contract(Algorithm(alg), A, ψ; kwargs...) end -function zipup_docstring(isMPOMPO::Bool)::String - rhsTypeString = isMPOMPO ? "MPO" : "MPS" - rhsString = isMPOMPO ? "B" : "ψ" - return """ - "zipup": The MPO and $rhsTypeString tensors are contracted then truncated at each site without enforcing - the appropriate orthogonal gauge. Once this sweep is complete a call to `truncate!` occurs. - Because the initial truncation is not locally optimal it is recommended to use a loose - `cutoff` and `maxdim` and then pass the desired truncation parameters to the locally optimal - `truncate!` sweep via the additional keyword argument `truncate_kwargs`. - A set of parameters suggested in [^Paeckel2019] is - `contract(A, $rhsString; method="zipup", cutoff=cutoff / 10, maxdim=2 * maxdim, truncate_kwargs=(; cutoff, maxdim))`.""" +function _zipup_docstring(rhs_type::String) + rhs_name = rhs_type == "MPO" ? "B" : "ψ" + return """ - `"zipup"`: The MPO and $rhs_type tensors are contracted then truncated at each site + without enforcing the appropriate orthogonal gauge. Once this sweep is complete a call + to `truncate!` occurs. Because the initial truncation is not locally optimal it is + recommended to use a loose `cutoff` and `maxdim` and then pass the desired truncation + parameters to the locally optimal `truncate!` sweep via the additional keyword argument + `truncate_kwargs`. A set of parameters suggested in [^Paeckel2019] is + `contract(A, $rhs_name; alg="zipup", cutoff=cutoff / 10, maxdim=2 * maxdim, truncate_kwargs=(; cutoff, maxdim))`.""" end -function Paeckel2019_citation_docstring()::String +function _paeckel2019_citation_docstring() return """ [^Paeckel2019]: Time-evolution methods for matrix-product states. Sebastian Paeckel et al. [arXiv:1901.05824](https://arxiv.org/abs/1901.05824) """ @@ -677,8 +676,7 @@ convenience function `apply`: apply(A, x; kwargs...) = replaceprime(contract(A, x; kwargs...), 2 => 1)`. ``` -Choose the method with the `method` keyword, for example -`"densitymatrix"` and `"naive"`. +Choose the algorithm with the `alg` keyword argument (or the deprecated `method` alias). # Keywords - `cutoff::Float64=1e-13`: the cutoff value for truncating the density matrix @@ -687,16 +685,18 @@ Choose the method with the `method` keyword, for example - `maxdim::Int=maxlinkdim(A) * maxlinkdim(ψ))`: the maximal bond dimension of the results MPS. - `mindim::Int=1`: the minimal bond dimension of the resulting MPS. - `normalize::Bool=false`: whether or not to normalize the resulting MPS. -- `method::String="densitymatrix"`: the algorithm to use for the contraction. - - "densitymatrix": The network formed by the MPO and MPS is squared and contracted down to - a density matrix which is diagonalized iteratively at each site. - - "naive": The MPO and MPS tensor are contracted exactly at each site and then a truncation - of the resulting MPS is performed. - $(zipup_docstring(false)) +- `alg::String="densitymatrix"`: the algorithm to use for the contraction. Options: + - `"densitymatrix"`: The network formed by the MPO and MPS is squared and contracted down + to a density matrix which is diagonalized iteratively at each site. + - `"naive"`: The MPO and MPS tensors are contracted exactly at each site and then a + truncation of the resulting MPS is performed. +$(_zipup_docstring("MPS")) +- `truncate_kwargs::NamedTuple`: (for `alg="zipup"` only) truncation parameters passed to + the final `truncate!` sweep. Defaults to `(; cutoff, maxdim, mindim)`. See also [`apply`](@ref). -$(Paeckel2019_citation_docstring()) +$(_paeckel2019_citation_docstring()) """ @doc """ @@ -865,11 +865,16 @@ function ITensors.contract( maxdim = maxlinkdim(A) * maxlinkdim(B), mindim = 1, truncate_kwargs = (; cutoff, maxdim, mindim), - kwargs..., + kwargs... ) if hassameinds(siteinds, A, B) error( - "In `contract(A::MPO, B::AbstractMPS)`, A and B have the same site indices. The indices of the MPOs in the contraction are taken literally, and therefore they should only share one site index per site so the contraction results in an MPO or MPS. You may want to use `replaceprime(contract(A', B), 2 => 1)` or `apply(A, B)` which automatically adjusts the prime levels assuming the inputs have pairs of primed and unprimed indices.", + "In `contract(A::MPO, B::AbstractMPS)`, A and B have the same site indices. " * + "The indices in the contraction are taken literally, and therefore they should " * + "only share one site index per site so the contraction results in an MPO or MPS. " * + "You may want to use `replaceprime(contract(A', B), 2 => 1)` or `apply(A, B)` " * + "which automatically adjusts the prime levels assuming the inputs have pairs of " * + "primed and unprimed indices." ) end N = length(A) @@ -984,16 +989,23 @@ C = apply(A, B; alg="naive", truncate=false) in general you should set a `cutoff` value. - `maxdim::Int=maxlinkdim(A) * maxlinkdim(B))`: the maximal bond dimension of the results MPS. - `mindim::Int=1`: the minimal bond dimension of the resulting MPS. -- `truncate=true`: Enable or disable truncation. If `truncate=false`, ignore - other truncation parameters like `cutoff` and `maxdim`. This is most relevant - for the `"naive"` version, if you just want to contract the tensors pairwise - exactly. This can be useful if you are contracting MPOs that have diverging - norms, such as MPOs originating from sums of local operators. - $(zipup_docstring(true)) +- `alg::String="zipup"`: the algorithm to use for the contraction. Options: + - `"zipup"`: Contracts pairs of site tensors and truncates with SVDs in a sweep across the + sites, then performs a final `truncate!` sweep. + - `"naive"`: Contracts pairs of tensors exactly and then truncates at the end if + `truncate=true`. +$(_zipup_docstring("MPO")) +- `truncate_kwargs::NamedTuple`: (for `alg="zipup"` only) truncation parameters passed to + the final `truncate!` sweep. Defaults to `(; cutoff, maxdim, mindim)`. +- `truncate=true`: (for `alg="naive"` only) Enable or disable truncation. If + `truncate=false`, ignore other truncation parameters like `cutoff` and `maxdim`. This is + most relevant if you just want to contract the tensors pairwise exactly. This can be useful + if you are contracting MPOs that have diverging norms, such as MPOs originating from sums + of local operators. See also [`apply`](@ref) for details about the arguments available. -$(Paeckel2019_citation_docstring()) +$(_paeckel2019_citation_docstring()) """ @doc """ diff --git a/test/base/test_mpo.jl b/test/base/test_mpo.jl index c852bfda..6e629e59 100644 --- a/test/base/test_mpo.jl +++ b/test/base/test_mpo.jl @@ -254,22 +254,23 @@ end @test_throws DimensionMismatch error_contract(phi, K, badpsi) end - @testset "contract(::MPO, ::MPS)" for method in ["densitymatrix", "naive", "zipup"] + @testset "contract(::MPO, ::MPS; alg=$alg)" for alg in + ["densitymatrix", "naive", "zipup"] phi = random_mps(sites) K = random_mpo(sites) @test maxlinkdim(K) == 1 psi = random_mps(sites) - psi_out = contract(K, psi; method, maxdim = 1) + psi_out = contract(K, psi; alg, maxdim = 1) @test inner(phi', psi_out) ≈ inner(phi', K, psi) - psi_out = contract(psi, K; method, maxdim = 1) + psi_out = contract(psi, K; alg, maxdim = 1) @test inner(phi', psi_out) ≈ inner(phi', K, psi) psi_out = psi * K @test inner(phi', psi_out) ≈ inner(phi', K, psi) - @test_throws MethodError contract(K, psi; method = "fakemethod") + @test_throws MethodError contract(K, psi; alg = "fakemethod") badsites = [Index(2, "Site") for n in 1:(N + 1)] badpsi = random_mps(badsites) - @test_throws DimensionMismatch contract(K, badpsi; method) + @test_throws DimensionMismatch contract(K, badpsi; alg) # make bigger random MPO... for link_dim in 2:5 @@ -301,7 +302,9 @@ end orthogonalize!(psi, 1; maxdim = link_dim) orthogonalize!(K, 1; maxdim = link_dim) orthogonalize!(phi, 1; normalize = true, maxdim = link_dim) - psi_out = contract(deepcopy(K), deepcopy(psi); method, maxdim = 10 * link_dim, cutoff = 0.0) + psi_out = contract( + deepcopy(K), deepcopy(psi); alg, maxdim = 10 * link_dim, cutoff = 0.0 + ) @test inner(phi', psi_out) ≈ inner(phi', K, psi) end end @@ -368,7 +371,7 @@ end @test maxlinkdim(H) ≤ maxlinkdim(H₁) + maxlinkdim(H₂) end - @testset "contract(::MPO, ::MPO)" for alg in ["naive", "zipup"] + @testset "contract(::MPO, ::MPO; alg=$alg)" for alg in ["naive", "zipup"] psi = random_mps(sites) K = random_mpo(sites) L = random_mpo(sites)