From 0679931ddceebea6807a0212e106c403eb380250 Mon Sep 17 00:00:00 2001 From: Joey Date: Fri, 5 Apr 2024 10:48:20 -0400 Subject: [PATCH 01/23] Move ITensorsExt to external module --- src/ITensorNetworks.jl | 2 +- .../ITensorsExtensions.jl} | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) rename src/{ITensorsExt/itensorutils.jl => ITensorsExtensions/ITensorsExtensions.jl} (99%) diff --git a/src/ITensorNetworks.jl b/src/ITensorNetworks.jl index 5e58d1b7..6b11f791 100644 --- a/src/ITensorNetworks.jl +++ b/src/ITensorNetworks.jl @@ -38,7 +38,6 @@ include("contraction_tree_to_graph.jl") include("gauging.jl") include("utils.jl") include("tensornetworkoperators.jl") -include("ITensorsExt/itensorutils.jl") include("solvers/local_solvers/eigsolve.jl") include("solvers/local_solvers/exponentiate.jl") include("solvers/local_solvers/dmrg_x.jl") @@ -68,6 +67,7 @@ include("environment.jl") include("exports.jl") include("ModelHamiltonians/ModelHamiltonians.jl") include("ModelNetworks/ModelNetworks.jl") +include("ITensorsExtensions/ITensorsExtensions.jl") using PackageExtensionCompat: @require_extensions using Requires: @require diff --git a/src/ITensorsExt/itensorutils.jl b/src/ITensorsExtensions/ITensorsExtensions.jl similarity index 99% rename from src/ITensorsExt/itensorutils.jl rename to src/ITensorsExtensions/ITensorsExtensions.jl index e8ce2e26..d5add565 100644 --- a/src/ITensorsExt/itensorutils.jl +++ b/src/ITensorsExtensions/ITensorsExtensions.jl @@ -1,3 +1,4 @@ +module ITensorsExtensions using LinearAlgebra: pinv using ITensors.NDTensors: Block, From 9529c9c90300ae53ad9e056668d24ae03a20e513 Mon Sep 17 00:00:00 2001 From: Joey Date: Tue, 9 Apr 2024 11:22:22 -0400 Subject: [PATCH 02/23] Refactor apply interface, simplify. --- src/ITensorNetworks.jl | 3 +- src/ITensorsExtensions/ITensorsExtensions.jl | 48 ++++-- src/ModelHamiltonians/ModelHamiltonians.jl | 20 ++- src/apply.jl | 150 +++++++----------- src/gauging.jl | 11 +- src/itensornetwork.jl | 2 +- src/tensornetworkoperators.jl | 47 ------ src/treetensornetworks/opsum_to_ttn.jl | 1 + .../projttns/abstractprojttn.jl | 2 +- test/test_tno.jl | 64 -------- 10 files changed, 117 insertions(+), 231 deletions(-) delete mode 100644 src/tensornetworkoperators.jl delete mode 100644 test/test_tno.jl diff --git a/src/ITensorNetworks.jl b/src/ITensorNetworks.jl index 6b11f791..ceffe592 100644 --- a/src/ITensorNetworks.jl +++ b/src/ITensorNetworks.jl @@ -37,7 +37,7 @@ include("caches/beliefpropagationcache.jl") include("contraction_tree_to_graph.jl") include("gauging.jl") include("utils.jl") -include("tensornetworkoperators.jl") +include("ITensorsExtensions/ITensorsExtensions.jl") include("solvers/local_solvers/eigsolve.jl") include("solvers/local_solvers/exponentiate.jl") include("solvers/local_solvers/dmrg_x.jl") @@ -67,7 +67,6 @@ include("environment.jl") include("exports.jl") include("ModelHamiltonians/ModelHamiltonians.jl") include("ModelNetworks/ModelNetworks.jl") -include("ITensorsExtensions/ITensorsExtensions.jl") using PackageExtensionCompat: @require_extensions using Requires: @require diff --git a/src/ITensorsExtensions/ITensorsExtensions.jl b/src/ITensorsExtensions/ITensorsExtensions.jl index d5add565..8641b364 100644 --- a/src/ITensorsExtensions/ITensorsExtensions.jl +++ b/src/ITensorsExtensions/ITensorsExtensions.jl @@ -1,18 +1,37 @@ module ITensorsExtensions -using LinearAlgebra: pinv +using LinearAlgebra: LinearAlgebra, eigen, pinv +using ITensors: + ITensor, + Index, + commonind, + dag, + hasqns, + inds, + isdiag, + itensor, + map_diag, + noncommonind, + noprime, + replaceinds, + space, + sqrt_decomp using ITensors.NDTensors: + NDTensors, Block, Tensor, blockdim, blockoffsets, + denseblocks, diaglength, getdiagindex, nzblocks, setdiagindex!, + svd, tensor, DiagBlockSparseTensor, DenseTensor, BlockOffsets +using Observers: update!, insert_function! function NDTensors.blockoffsets(dense::DenseTensor) return BlockOffsets{ndims(dense)}([Block(ntuple(Returns(1), ndims(dense)))], [0]) @@ -26,19 +45,6 @@ NDTensors.blockdim(i::Index{Int}, b::Block) = blockdim(space(i), b) LinearAlgebra.isdiag(it::ITensor) = isdiag(tensor(it)) -function map_diag!(f::Function, it_destination::ITensor, it_source::ITensor) - return itensor(map_diag!(f, tensor(it_destination), tensor(it_source))) -end -map_diag(f::Function, it::ITensor) = map_diag!(f, copy(it), it) - -function map_diag!(f::Function, t_destination::Tensor, t_source::Tensor) - for i in 1:diaglength(t_destination) - setdiagindex!(t_destination, f(getdiagindex(t_source, i)), i) - end - return t_destination -end -map_diag(f::Function, t::Tensor) = map_diag!(f, copy(t), t) - # Convenience functions sqrt_diag(it::ITensor) = map_diag(sqrt, it) inv_diag(it::ITensor) = map_diag(inv, it) @@ -46,6 +52,18 @@ invsqrt_diag(it::ITensor) = map_diag(inv ∘ sqrt, it) pinv_diag(it::ITensor) = map_diag(pinv, it) pinvsqrt_diag(it::ITensor) = map_diag(pinv ∘ sqrt, it) +function map_itensor( + f::Function, A::ITensor, lind=first(inds(A)); regularization=nothing, kwargs... +) + USV = svd(A, lind; kwargs...) + U, S, V, spec, u, v = USV + S = map_diag(s -> f(s + regularization), S) + sqrtDL, δᵤᵥ, sqrtDR = sqrt_decomp(S, u, v) + sqrtDR = denseblocks(sqrtDR) * denseblocks(δᵤᵥ) + L, R = U * sqrtDL, V * sqrtDR + return L * R +end + # Analagous to `denseblocks`. # Extract the diagonal entries into a diagonal tensor. function diagblocks(D::Tensor) @@ -89,3 +107,5 @@ function group_commuting_itensors(its::Vector{ITensor}) return it_groups end + +end diff --git a/src/ModelHamiltonians/ModelHamiltonians.jl b/src/ModelHamiltonians/ModelHamiltonians.jl index f54faea3..eebc80ec 100644 --- a/src/ModelHamiltonians/ModelHamiltonians.jl +++ b/src/ModelHamiltonians/ModelHamiltonians.jl @@ -82,20 +82,26 @@ function heisenberg(g::AbstractGraph; J1=1, J2=0, h=0) (; J1, J2, h) = map(to_callable, (; J1, J2, h)) ℋ = OpSum() for e in edges(g) - ℋ += J1(e) / 2, "S+", src(e), "S-", dst(e) - ℋ += J1(e) / 2, "S-", src(e), "S+", dst(e) - ℋ += J1(e), "Sz", src(e), "Sz", dst(e) + if !iszero(J1(e)) + ℋ += J1(e) / 2, "S+", src(e), "S-", dst(e) + ℋ += J1(e) / 2, "S-", src(e), "S+", dst(e) + ℋ += J1(e), "Sz", src(e), "Sz", dst(e) + end end for v in vertices(g) for nn in next_nearest_neighbors(g, v) e = edgetype(g)(v, nn) - ℋ += J2(e) / 2, "S+", src(e), "S-", dst(e) - ℋ += J2(e) / 2, "S-", src(e), "S+", dst(e) - ℋ += J2(e), "Sz", src(e), "Sz", dst(e) + if !iszero(J2(e)) + ℋ += J2(e) / 2, "S+", src(e), "S-", dst(e) + ℋ += J2(e) / 2, "S-", src(e), "S+", dst(e) + ℋ += J2(e), "Sz", src(e), "Sz", dst(e) + end end end for v in vertices(g) - ℋ += h(v), "Sz", v + if !iszero(h(v)) + ℋ += h(v), "Sz", v + end end return ℋ end diff --git a/src/apply.jl b/src/apply.jl index 559438e2..b7b54b48 100644 --- a/src/apply.jl +++ b/src/apply.jl @@ -16,67 +16,16 @@ using ITensors: prime, replaceind, replaceinds, + factorize_svd, unioninds, uniqueinds using ITensors.ContractionSequenceOptimization: optimal_contraction_sequence using ITensors.ITensorMPS: siteinds using KrylovKit: linsolve using LinearAlgebra: eigen, norm, svd -using NamedGraphs: NamedEdge +using NamedGraphs: NamedEdge, has_edge using Observers: Observers -function sqrt_and_inv_sqrt( - A::ITensor; ishermitian=false, cutoff=nothing, regularization=nothing -) - if isdiag(A) - A = map_diag(x -> x + regularization, A) - sqrtA = sqrt_diag(A) - inv_sqrtA = inv_diag(sqrtA) - return sqrtA, inv_sqrtA - end - @assert ishermitian - D, U = eigen(A; ishermitian, cutoff) - D = map_diag(x -> x + regularization, D) - sqrtD = sqrt_diag(D) - # sqrtA = U * sqrtD * prime(dag(U)) - sqrtA = noprime(sqrtD * U) - inv_sqrtD = inv_diag(sqrtD) - # inv_sqrtA = U * inv_sqrtD * prime(dag(U)) - inv_sqrtA = noprime(inv_sqrtD * dag(U)) - return sqrtA, inv_sqrtA -end - -function symmetric_factorize( - A::ITensor, inds...; (observer!)=nothing, tags="", svd_kwargs... -) - if !isnothing(observer!) - Observers.insert_function!( - observer!, "singular_values" => (; singular_values) -> singular_values - ) - end - U, S, V = svd(A, inds...; lefttags=tags, righttags=tags, svd_kwargs...) - u = commonind(S, U) - v = commonind(S, V) - sqrtS = sqrt_diag(S) - Fu = U * sqrtS - Fv = V * sqrtS - if hasqns(A) - # Hack to make a generalized (non-diagonal) `δ` tensor. - # TODO: Make this easier, `ITensors.δ` doesn't work here. - δᵤᵥ = copy(S) - ITensors.data(δᵤᵥ) .= true - Fu *= dag(δᵤᵥ) - S = denseblocks(S) - S *= prime(dag(δᵤᵥ), u) - S = diagblocks(S) - else - Fu = replaceinds(Fu, v => u) - S = replaceinds(S, v => u') - end - Observers.update!(observer!; singular_values=S) - return Fu, Fv -end - function full_update_bp( o, ψ, @@ -85,7 +34,7 @@ function full_update_bp( nfullupdatesweeps=10, print_fidelity_loss=false, envisposdef=false, - (observer!)=nothing, + (singular_values!)=nothing, symmetrize=false, apply_kwargs..., ) @@ -116,8 +65,13 @@ function full_update_bp( apply_kwargs..., ) if symmetrize - Rᵥ₁, Rᵥ₂ = symmetric_factorize( - Rᵥ₁ * Rᵥ₂, inds(Rᵥ₁); tags=edge_tag(v⃗[1] => v⃗[2]), observer!, apply_kwargs... + Rᵥ₁, Rᵥ₂ = factorize_svd( + Rᵥ₁ * Rᵥ₂, + inds(Rᵥ₁); + ortho="none", + tags=edge_tag(v⃗[1] => v⃗[2]), + singular_values!, + apply_kwargs..., ) end ψᵥ₁ = Qᵥ₁ * Rᵥ₁ @@ -125,19 +79,25 @@ function full_update_bp( return ψᵥ₁, ψᵥ₂ end -function simple_update_bp_full(o, ψ, v⃗; envs, (observer!)=nothing, apply_kwargs...) +function simple_update_bp_full(o, ψ, v⃗; envs, (singular_values!)=nothing, apply_kwargs...) cutoff = 10 * eps(real(scalartype(ψ))) regularization = 10 * eps(real(scalartype(ψ))) envs_v1 = filter(env -> hascommoninds(env, ψ[v⃗[1]]), envs) envs_v2 = filter(env -> hascommoninds(env, ψ[v⃗[2]]), envs) - sqrt_and_inv_sqrt_envs_v1 = - sqrt_and_inv_sqrt.(envs_v1; ishermitian=true, cutoff, regularization) - sqrt_and_inv_sqrt_envs_v2 = - sqrt_and_inv_sqrt.(envs_v2; ishermitian=true, cutoff, regularization) - sqrt_envs_v1 = first.(sqrt_and_inv_sqrt_envs_v1) - inv_sqrt_envs_v1 = last.(sqrt_and_inv_sqrt_envs_v1) - sqrt_envs_v2 = first.(sqrt_and_inv_sqrt_envs_v2) - inv_sqrt_envs_v2 = last.(sqrt_and_inv_sqrt_envs_v2) + sqrt_envs_v1 = [ + ITensorsExtensions.map_itensor(sqrt, env; cutoff, regularization) for env in envs_v1 + ] + sqrt_envs_v2 = [ + ITensorsExtensions.map_itensor(sqrt, env; cutoff, regularization) for env in envs_v2 + ] + inv_sqrt_envs_v1 = [ + ITensorsExtensions.map_itensor(inv ∘ sqrt, env; cutoff, regularization) for + env in envs_v1 + ] + inv_sqrt_envs_v2 = [ + ITensorsExtensions.map_itensor(inv ∘ sqrt, env; cutoff, regularization) for + env in envs_v2 + ] ψᵥ₁ᵥ₂_tn = [ψ[v⃗[1]]; ψ[v⃗[2]]; sqrt_envs_v1; sqrt_envs_v2] ψᵥ₁ᵥ₂ = contract(ψᵥ₁ᵥ₂_tn; sequence=contraction_sequence(ψᵥ₁ᵥ₂_tn; alg="optimal")) oψ = apply(o, ψᵥ₁ᵥ₂) @@ -150,32 +110,38 @@ function simple_update_bp_full(o, ψ, v⃗; envs, (observer!)=nothing, apply_kwa v1_inds = [v1_inds; siteinds(ψ, v⃗[1])] v2_inds = [v2_inds; siteinds(ψ, v⃗[2])] e = v⃗[1] => v⃗[2] - ψᵥ₁, ψᵥ₂ = symmetric_factorize(oψ, v1_inds; tags=edge_tag(e), observer!, apply_kwargs...) + ψᵥ₁, ψᵥ₂ = factorize_svd( + oψ, v1_inds; ortho="none", tags=edge_tag(e), singular_values!, apply_kwargs... + ) for inv_sqrt_env_v1 in inv_sqrt_envs_v1 - # TODO: `dag` here? - ψᵥ₁ *= inv_sqrt_env_v1 + ψᵥ₁ *= dag(inv_sqrt_env_v1) end for inv_sqrt_env_v2 in inv_sqrt_envs_v2 - # TODO: `dag` here? - ψᵥ₂ *= inv_sqrt_env_v2 + ψᵥ₂ *= dag(inv_sqrt_env_v2) end return ψᵥ₁, ψᵥ₂ end # Reduced version -function simple_update_bp(o, ψ, v⃗; envs, (observer!)=nothing, apply_kwargs...) +function simple_update_bp(o, ψ, v⃗; envs, (singular_values!)=nothing, apply_kwargs...) cutoff = 10 * eps(real(scalartype(ψ))) regularization = 10 * eps(real(scalartype(ψ))) envs_v1 = filter(env -> hascommoninds(env, ψ[v⃗[1]]), envs) envs_v2 = filter(env -> hascommoninds(env, ψ[v⃗[2]]), envs) - sqrt_and_inv_sqrt_envs_v1 = - sqrt_and_inv_sqrt.(envs_v1; ishermitian=true, cutoff, regularization) - sqrt_and_inv_sqrt_envs_v2 = - sqrt_and_inv_sqrt.(envs_v2; ishermitian=true, cutoff, regularization) - sqrt_envs_v1 = first.(sqrt_and_inv_sqrt_envs_v1) - inv_sqrt_envs_v1 = last.(sqrt_and_inv_sqrt_envs_v1) - sqrt_envs_v2 = first.(sqrt_and_inv_sqrt_envs_v2) - inv_sqrt_envs_v2 = last.(sqrt_and_inv_sqrt_envs_v2) + sqrt_envs_v1 = [ + ITensorsExtensions.map_itensor(sqrt, env; cutoff, regularization) for env in envs_v1 + ] + sqrt_envs_v2 = [ + ITensorsExtensions.map_itensor(sqrt, env; cutoff, regularization) for env in envs_v2 + ] + inv_sqrt_envs_v1 = [ + ITensorsExtensions.map_itensor(inv ∘ sqrt, env; cutoff, regularization) for + env in envs_v1 + ] + inv_sqrt_envs_v2 = [ + ITensorsExtensions.map_itensor(inv ∘ sqrt, env; cutoff, regularization) for + env in envs_v2 + ] ψᵥ₁ = contract([ψ[v⃗[1]]; sqrt_envs_v1]) ψᵥ₂ = contract([ψ[v⃗[2]]; sqrt_envs_v2]) sᵥ₁ = siteinds(ψ, v⃗[1]) @@ -186,12 +152,16 @@ function simple_update_bp(o, ψ, v⃗; envs, (observer!)=nothing, apply_kwargs.. rᵥ₂ = commoninds(Qᵥ₂, Rᵥ₂) oR = apply(o, Rᵥ₁ * Rᵥ₂) e = v⃗[1] => v⃗[2] - Rᵥ₁, Rᵥ₂ = symmetric_factorize( - oR, unioninds(rᵥ₁, sᵥ₁); tags=edge_tag(e), observer!, apply_kwargs... + Rᵥ₁, Rᵥ₂ = factorize_svd( + oR, + unioninds(rᵥ₁, sᵥ₁); + ortho="none", + tags=edge_tag(e), + singular_values!, + apply_kwargs..., ) - # TODO: `dag` here? - Qᵥ₁ = contract([Qᵥ₁; inv_sqrt_envs_v1]) - Qᵥ₂ = contract([Qᵥ₂; inv_sqrt_envs_v2]) + Qᵥ₁ = contract([Qᵥ₁; dag.(inv_sqrt_envs_v1)]) + Qᵥ₂ = contract([Qᵥ₂; dag.(inv_sqrt_envs_v2)]) ψᵥ₁ = Qᵥ₁ * Rᵥ₁ ψᵥ₂ = Qᵥ₂ * Rᵥ₂ return ψᵥ₁, ψᵥ₂ @@ -206,7 +176,7 @@ function ITensors.apply( nfullupdatesweeps=10, print_fidelity_loss=false, envisposdef=false, - (observer!)=nothing, + (singular_values!)=nothing, variational_optimization_only=false, symmetrize=false, reduced=true, @@ -242,15 +212,15 @@ function ITensors.apply( nfullupdatesweeps, print_fidelity_loss, envisposdef, - observer!, + singular_values!, symmetrize, apply_kwargs..., ) else if reduced - ψᵥ₁, ψᵥ₂ = simple_update_bp(o, ψ, v⃗; envs, observer!, apply_kwargs...) + ψᵥ₁, ψᵥ₂ = simple_update_bp(o, ψ, v⃗; envs, singular_values!, apply_kwargs...) else - ψᵥ₁, ψᵥ₂ = simple_update_bp_full(o, ψ, v⃗; envs, observer!, apply_kwargs...) + ψᵥ₁, ψᵥ₂ = simple_update_bp_full(o, ψ, v⃗; envs, singular_values!, apply_kwargs...) end end if normalize @@ -366,13 +336,13 @@ function ITensors.apply(o, ψ::VidalITensorNetwork; normalize=false, apply_kwarg for vn in neighbors(ψ, src(e)) if (vn != dst(e)) - ψv1 = noprime(ψv1 * inv_diag(bond_tensor(ψ, vn => src(e)))) + ψv1 = noprime(ψv1 * ITensorsExtensions.inv_diag(bond_tensor(ψ, vn => src(e)))) end end for vn in neighbors(ψ, dst(e)) if (vn != src(e)) - ψv2 = noprime(ψv2 * inv_diag(bond_tensor(ψ, vn => dst(e)))) + ψv2 = noprime(ψv2 * ITensorsExtensions.inv_diag(bond_tensor(ψ, vn => dst(e)))) end end diff --git a/src/gauging.jl b/src/gauging.jl index 449a8cbd..e9cc305a 100644 --- a/src/gauging.jl +++ b/src/gauging.jl @@ -40,7 +40,7 @@ function ITensorNetwork( for e in edges(ψ) vsrc, vdst = src(e), dst(e) - root_S = sqrt_diag(bond_tensor(ψ_vidal, e)) + root_S = ITensorsExtensions.sqrt_diag(bond_tensor(ψ_vidal, e)) setindex_preserve_graph!(ψ, noprime(root_S * ψ[vsrc]), vsrc) setindex_preserve_graph!(ψ, noprime(root_S * ψ[vdst]), vdst) end @@ -88,11 +88,12 @@ function vidalitensornetwork_preserve_cache( Y_D, Y_U = eigen( only(message(cache, reverse(pe))); ishermitian=true, cutoff=message_cutoff ) - X_D, Y_D = map_diag(x -> x + regularization, X_D), - map_diag(x -> x + regularization, Y_D) + X_D, Y_D = ITensorsExtensions.map_diag(x -> x + regularization, X_D), + ITensorsExtensions.map_diag(x -> x + regularization, Y_D) - rootX_D, rootY_D = sqrt_diag(X_D), sqrt_diag(Y_D) - inv_rootX_D, inv_rootY_D = invsqrt_diag(X_D), invsqrt_diag(Y_D) + rootX_D, rootY_D = ITensorsExtensions.sqrt_diag(X_D), ITensorsExtensions.sqrt_diag(Y_D) + inv_rootX_D, inv_rootY_D = ITensorsExtensions.invsqrt_diag(X_D), + ITensorsExtensions.invsqrt_diag(Y_D) rootX = X_U * rootX_D * prime(dag(X_U)) rootY = Y_U * rootY_D * prime(dag(Y_U)) inv_rootX = X_U * inv_rootX_D * prime(dag(X_U)) diff --git a/src/itensornetwork.jl b/src/itensornetwork.jl index 5e84138e..f6789ca2 100644 --- a/src/itensornetwork.jl +++ b/src/itensornetwork.jl @@ -1,6 +1,6 @@ using DataGraphs: DataGraphs, DataGraph using Dictionaries: dictionary -using ITensors: ITensor +using ITensors: ITensor, factorize using NamedGraphs: NamedGraphs, NamedEdge, NamedGraph, vertextype struct Private end diff --git a/src/tensornetworkoperators.jl b/src/tensornetworkoperators.jl deleted file mode 100644 index 34c332dc..00000000 --- a/src/tensornetworkoperators.jl +++ /dev/null @@ -1,47 +0,0 @@ -using Graphs: has_edge -using ITensors: ITensors, commoninds, product -using LinearAlgebra: factorize - -""" -Take a vector of gates which act on different edges/ vertices of an Inds network and construct the tno which represents prod(gates). -""" -function gate_group_to_tno(s::IndsNetwork, gates::Vector{ITensor}) - - #Construct indsnetwork for TNO - s_O = union_all_inds(s, prime(s; links=[])) - - O = delta_network(s_O) - - for gate in gates - v⃗ = vertices(s)[findall(i -> (length(commoninds(s[i], inds(gate))) != 0), vertices(s))] - if length(v⃗) == 1 - O[v⃗[1]] = product(O[v⃗[1]], gate) - elseif length(v⃗) == 2 - e = v⃗[1] => v⃗[2] - if !has_edge(s, e) - error("Vertices where the gates are being applied must be neighbors for now.") - end - Osrc, Odst = factorize(gate, commoninds(O[v⃗[1]], gate)) - O[v⃗[1]] = product(O[v⃗[1]], Osrc) - O[v⃗[2]] = product(O[v⃗[2]], Odst) - else - error( - "Can only deal with gates acting on one or two sites for now. Physical indices of the gates must also match those in the IndsNetwork.", - ) - end - end - - return combine_linkinds(O) -end - -"""Take a series of gates acting on the physical indices specified by IndsNetwork convert into a series of tnos -whose product represents prod(gates). Useful for keeping the bond dimension of each tno low (as opposed to just building a single tno)""" -function get_tnos(s::IndsNetwork, gates::Vector{ITensor}) - tnos = ITensorNetwork[] - gate_groups = group_commuting_itensors(gates) - for group in gate_groups - push!(tnos, gate_group_to_tno(s, group)) - end - - return tnos -end diff --git a/src/treetensornetworks/opsum_to_ttn.jl b/src/treetensornetworks/opsum_to_ttn.jl index 4ffd1743..e4972977 100644 --- a/src/treetensornetworks/opsum_to_ttn.jl +++ b/src/treetensornetworks/opsum_to_ttn.jl @@ -14,6 +14,7 @@ using ITensors.LazyApply: coefficient using ITensors.LazyApply: Sum using ITensors.LazyApply: Prod using ITensors.NDTensors: Block +using ITensors.NDTensors: blockdim using ITensors.NDTensors: maxdim using ITensors.NDTensors: nblocks using ITensors.NDTensors: nnzblocks diff --git a/src/treetensornetworks/projttns/abstractprojttn.jl b/src/treetensornetworks/projttns/abstractprojttn.jl index 63ff4bf7..d86cc48a 100644 --- a/src/treetensornetworks/projttns/abstractprojttn.jl +++ b/src/treetensornetworks/projttns/abstractprojttn.jl @@ -1,6 +1,6 @@ using DataGraphs: DataGraphs, underlying_graph using Graphs: neighbors -using ITensors: ITensor, contract, order +using ITensors: ITensor, contract, order, product using ITensors.ITensorMPS: ITensorMPS, nsite using NamedGraphs: NamedGraphs, NamedEdge, incident_edges, vertextype diff --git a/test/test_tno.jl b/test/test_tno.jl deleted file mode 100644 index 4a3f5018..00000000 --- a/test/test_tno.jl +++ /dev/null @@ -1,64 +0,0 @@ -@eval module $(gensym()) -using Graphs: vertices -using ITensorNetworks: - apply, - contract_inner, - flatten_networks, - group_commuting_itensors, - gate_group_to_tno, - get_tnos, - random_tensornetwork, - siteinds -using ITensorNetworks.ModelHamiltonians: ModelHamiltonians -using ITensors: ITensor, noprime -using NamedGraphs: named_grid -using Test: @test, @testset - -@testset "TN operator Basics" begin - L = 3 - g = named_grid((L, L)) - s = siteinds("S=1/2", g) - - ℋ = ModelHamiltonians.ising(g; h=1.5) - gates = Vector{ITensor}(ℋ, s) - gate_groups = group_commuting_itensors(gates) - - @test typeof(gate_groups) == Vector{Vector{ITensor}} - - #Construct a number of tnos whose product is prod(gates) - tnos = get_tnos(s, gates) - @test length(tnos) == length(gate_groups) - - #Construct a single tno which represents prod(gates) - single_tno = gate_group_to_tno(s, gates) - - ψ = random_tensornetwork(s; link_space=2) - - ψ_gated = copy(ψ) - for gate in gates - ψ_gated = apply(gate, ψ_gated) - end - ψ_tnod = copy(ψ) - for tno in tnos - ψ_tnod = flatten_networks(ψ_tnod, tno) - for v in vertices(ψ_tnod) - ψ_tnod[v] = noprime(ψ_tnod[v]) - end - end - ψ_tno = copy(ψ) - ψ_tno = flatten_networks(ψ_tno, single_tno) - for v in vertices(ψ_tno) - ψ_tno[v] = noprime(ψ_tno[v]) - end - - z1 = contract_inner(ψ_gated, ψ_gated) - z2 = contract_inner(ψ_tnod, ψ_tnod) - z3 = contract_inner(ψ_tno, ψ_tno) - f12 = contract_inner(ψ_tnod, ψ_gated) / sqrt(z1 * z2) - f13 = contract_inner(ψ_tno, ψ_gated) / sqrt(z1 * z3) - f23 = contract_inner(ψ_tno, ψ_tnod) / sqrt(z2 * z3) - @test f12 * conj(f12) ≈ 1.0 - @test f13 * conj(f13) ≈ 1.0 - @test f23 * conj(f23) ≈ 1.0 -end -end From dc781c7684c13f033086c3b18b0d01ddd711f146 Mon Sep 17 00:00:00 2001 From: Joey Date: Fri, 12 Apr 2024 13:22:16 -0400 Subject: [PATCH 03/23] Refactor map_ITensor to be based on Eigen --- src/ITensorsExtensions/ITensorsExtensions.jl | 27 ++++++++++++++------ 1 file changed, 19 insertions(+), 8 deletions(-) diff --git a/src/ITensorsExtensions/ITensorsExtensions.jl b/src/ITensorsExtensions/ITensorsExtensions.jl index 8641b364..7c1f2b15 100644 --- a/src/ITensorsExtensions/ITensorsExtensions.jl +++ b/src/ITensorsExtensions/ITensorsExtensions.jl @@ -12,7 +12,9 @@ using ITensors: map_diag, noncommonind, noprime, + replaceind, replaceinds, + sim, space, sqrt_decomp using ITensors.NDTensors: @@ -53,15 +55,24 @@ pinv_diag(it::ITensor) = map_diag(pinv, it) pinvsqrt_diag(it::ITensor) = map_diag(pinv ∘ sqrt, it) function map_itensor( - f::Function, A::ITensor, lind=first(inds(A)); regularization=nothing, kwargs... + f::Function, A::ITensor, linds=Index[first(inds(A))]; regularization=nothing, kwargs... ) - USV = svd(A, lind; kwargs...) - U, S, V, spec, u, v = USV - S = map_diag(s -> f(s + regularization), S) - sqrtDL, δᵤᵥ, sqrtDR = sqrt_decomp(S, u, v) - sqrtDR = denseblocks(sqrtDR) * denseblocks(δᵤᵥ) - L, R = U * sqrtDL, V * sqrtDR - return L * R + rinds = setdiff(inds(A), linds) + + D, U = eigen(A, linds, rinds; kwargs...) + ul, ur = noncommonind(D, U), commonind(D, U) + ulnew = sim(noprime(ul)) + + Ur = dag(U) + Ul = replaceinds(U, (rinds..., ur), (linds..., ulnew)) + + D_mapped = map_diag(s -> f(s + regularization), D) + D_mapped = replaceind(D_mapped, ul, ulnew) + + sqrtD_mapped_L, δᵤᵥ, sqrtD_mapped_R = sqrt_decomp(D_mapped, ulnew, ur) + sqrtD_mapped_R = denseblocks(sqrtD_mapped_R) * denseblocks(δᵤᵥ) + + return Ul * sqrtD_mapped_L * sqrtD_mapped_R * Ur end # Analagous to `denseblocks`. From cb3616a1c5ef1aa5fa2f5b909a145c49d9894ceb Mon Sep 17 00:00:00 2001 From: Joey Date: Fri, 12 Apr 2024 17:29:12 -0400 Subject: [PATCH 04/23] Added testing --- src/ITensorsExtensions/ITensorsExtensions.jl | 33 ++++++------- src/apply.jl | 33 ++++++++----- src/tensornetworkoperators.jl | 0 src/treetensornetworks/opsum_to_ttn.jl | 2 +- test/test_itensorsextensions.jl | 49 ++++++++++++++++++++ 5 files changed, 88 insertions(+), 29 deletions(-) delete mode 100644 src/tensornetworkoperators.jl create mode 100644 test/test_itensorsextensions.jl diff --git a/src/ITensorsExtensions/ITensorsExtensions.jl b/src/ITensorsExtensions/ITensorsExtensions.jl index cd663d03..37d86757 100644 --- a/src/ITensorsExtensions/ITensorsExtensions.jl +++ b/src/ITensorsExtensions/ITensorsExtensions.jl @@ -1,5 +1,5 @@ module ITensorsExtensions -using LinearAlgebra: LinearAlgebra, eigen, pinv +using LinearAlgebra: LinearAlgebra, eigen, ishermitian, pinv using ITensors: ITensor, Index, @@ -55,24 +55,25 @@ pinv_diag(it::ITensor) = map_diag(pinv, it) pinvsqrt_diag(it::ITensor) = map_diag(pinv ∘ sqrt, it) function map_itensor( - f::Function, A::ITensor, linds=Index[first(inds(A))]; regularization=nothing, kwargs... + f::Function, + A::ITensor, + lind=first(inds(A)); + ishermitian=false, + regularization=nothing, + kwargs..., ) - rinds = setdiff(inds(A), linds) + USV = svd(A, lind; kwargs...) + U, S, V, spec, u, v = USV - D, U = eigen(A, linds, rinds; kwargs...) - ul, ur = noncommonind(D, U), commonind(D, U) - ulnew = sim(noprime(ul)) - - Ur = dag(U) - Ul = replaceinds(U, (rinds..., ur), (linds..., ulnew)) - - D_mapped = map_diag(s -> f(s + regularization), D) - D_mapped = replaceind(D_mapped, ul, ulnew) - - sqrtD_mapped_L, δᵤᵥ, sqrtD_mapped_R = sqrt_decomp(D_mapped, ulnew, ur) - sqrtD_mapped_R = denseblocks(sqrtD_mapped_R) * denseblocks(δᵤᵥ) + if !isnothing(regularization) + f = s -> f(s + regularization) + end - return Ul * sqrtD_mapped_L * sqrtD_mapped_R * Ur + S = map_diag(s -> f(s), S) + sqrtDL, δᵤᵥ, sqrtDR = sqrt_decomp(S, u, v) + sqrtDR = denseblocks(sqrtDR) * denseblocks(δᵤᵥ) + L, R = U * sqrtDL, V * sqrtDR + return L * R end # Analagous to `denseblocks`. diff --git a/src/apply.jl b/src/apply.jl index b7b54b48..ee32b2a6 100644 --- a/src/apply.jl +++ b/src/apply.jl @@ -10,6 +10,7 @@ using ITensors: contract, dag, denseblocks, + factorize, hasqns, isdiag, noprime, @@ -85,18 +86,22 @@ function simple_update_bp_full(o, ψ, v⃗; envs, (singular_values!)=nothing, ap envs_v1 = filter(env -> hascommoninds(env, ψ[v⃗[1]]), envs) envs_v2 = filter(env -> hascommoninds(env, ψ[v⃗[2]]), envs) sqrt_envs_v1 = [ - ITensorsExtensions.map_itensor(sqrt, env; cutoff, regularization) for env in envs_v1 + ITensorsExtensions.map_itensor(sqrt, env; ishermitian=true, cutoff, regularization) for + env in envs_v1 ] sqrt_envs_v2 = [ - ITensorsExtensions.map_itensor(sqrt, env; cutoff, regularization) for env in envs_v2 + ITensorsExtensions.map_itensor(sqrt, env; ishermitian=true, cutoff, regularization) for + env in envs_v2 ] inv_sqrt_envs_v1 = [ - ITensorsExtensions.map_itensor(inv ∘ sqrt, env; cutoff, regularization) for - env in envs_v1 + ITensorsExtensions.map_itensor( + inv ∘ sqrt, env; ishermitian=true, cutoff, regularization + ) for env in envs_v1 ] inv_sqrt_envs_v2 = [ - ITensorsExtensions.map_itensor(inv ∘ sqrt, env; cutoff, regularization) for - env in envs_v2 + ITensorsExtensions.map_itensor( + inv ∘ sqrt, env; ishermitian=true, cutoff, regularization + ) for env in envs_v2 ] ψᵥ₁ᵥ₂_tn = [ψ[v⃗[1]]; ψ[v⃗[2]]; sqrt_envs_v1; sqrt_envs_v2] ψᵥ₁ᵥ₂ = contract(ψᵥ₁ᵥ₂_tn; sequence=contraction_sequence(ψᵥ₁ᵥ₂_tn; alg="optimal")) @@ -129,18 +134,22 @@ function simple_update_bp(o, ψ, v⃗; envs, (singular_values!)=nothing, apply_k envs_v1 = filter(env -> hascommoninds(env, ψ[v⃗[1]]), envs) envs_v2 = filter(env -> hascommoninds(env, ψ[v⃗[2]]), envs) sqrt_envs_v1 = [ - ITensorsExtensions.map_itensor(sqrt, env; cutoff, regularization) for env in envs_v1 + ITensorsExtensions.map_itensor(sqrt, env; ishermitian=true, cutoff, regularization) for + env in envs_v1 ] sqrt_envs_v2 = [ - ITensorsExtensions.map_itensor(sqrt, env; cutoff, regularization) for env in envs_v2 + ITensorsExtensions.map_itensor(sqrt, env; ishermitian=true, cutoff, regularization) for + env in envs_v2 ] inv_sqrt_envs_v1 = [ - ITensorsExtensions.map_itensor(inv ∘ sqrt, env; cutoff, regularization) for - env in envs_v1 + ITensorsExtensions.map_itensor( + inv ∘ sqrt, env; ishermitian=true, cutoff, regularization + ) for env in envs_v1 ] inv_sqrt_envs_v2 = [ - ITensorsExtensions.map_itensor(inv ∘ sqrt, env; cutoff, regularization) for - env in envs_v2 + ITensorsExtensions.map_itensor( + inv ∘ sqrt, env; ishermitian=true, cutoff, regularization + ) for env in envs_v2 ] ψᵥ₁ = contract([ψ[v⃗[1]]; sqrt_envs_v1]) ψᵥ₂ = contract([ψ[v⃗[2]]; sqrt_envs_v2]) diff --git a/src/tensornetworkoperators.jl b/src/tensornetworkoperators.jl deleted file mode 100644 index e69de29b..00000000 diff --git a/src/treetensornetworks/opsum_to_ttn.jl b/src/treetensornetworks/opsum_to_ttn.jl index 9b555547..51d8da07 100644 --- a/src/treetensornetworks/opsum_to_ttn.jl +++ b/src/treetensornetworks/opsum_to_ttn.jl @@ -2,7 +2,7 @@ using Graphs: degree, is_tree using ITensors: flux, has_fermion_string, itensor, ops, removeqns, space, val using ITensors.ITensorMPS: ITensorMPS, cutoff, linkdims, truncate! using ITensors.LazyApply: Prod, Sum, coefficient -using ITensors.NDTensors: Block, maxdim, nblocks, nnzblocks +using ITensors.NDTensors: Block, blockdim, maxdim, nblocks, nnzblocks using ITensors.Ops: Op, OpSum using NamedGraphs: degrees, is_leaf, vertex_path using StaticArrays: MVector diff --git a/test/test_itensorsextensions.jl b/test/test_itensorsextensions.jl new file mode 100644 index 00000000..b966594c --- /dev/null +++ b/test/test_itensorsextensions.jl @@ -0,0 +1,49 @@ +using ITensors: + ITensors, + ITensor, + Index, + randomITensor, + dag, + replaceind, + noprime, + prime, + inds, + delta, + QN, + denseblocks, + replaceinds, + dir, + array +using ITensorNetworks.ITensorsExtensions: map_itensor +using ITensorNetworks: siteinds, random_tensornetwork +using NamedGraphs: named_grid +using Random +using Test: @test, @testset +using LinearAlgebra: ishermitian, isposdef + +Random.seed!(1234) +@testset "ITensorsExtensions" begin + for eltype in [Float64, ComplexF64] + for n in [2, 3, 5, 10] + i, j = Index(n, "i"), Index(n, "j") + A = randn(eltype, (n, n)) + A = A * A' + P = ITensor(A, i, j) + sqrtP = map_itensor(sqrt, P; ishermitian=true) + inv_P = dag(map_itensor(inv, P; ishermitian=true)) + inv_sqrtP = dag(map_itensor(inv ∘ sqrt, P; ishermitian=true)) + + sqrtPdag = replaceind(dag(sqrtP), i, i') + P2 = replaceind(sqrtP * sqrtPdag, i', j) + @test P2 ≈ P atol = 1e-12 + + invP = replaceind(inv_P, i, i') + I = invP * P + @test I ≈ delta(eltype, inds(I)) atol = 1e-12 + + inv_sqrtP = replaceind(inv_sqrtP, i, i') + I = inv_sqrtP * sqrtP + @test I ≈ delta(eltype, inds(I)) atol = 1e-12 + end + end +end From cc55807e9f2d70ebb73b10cb8597c197f0ec5c6c Mon Sep 17 00:00:00 2001 From: Joey Date: Fri, 12 Apr 2024 17:41:42 -0400 Subject: [PATCH 05/23] Restored ModelHamiltonians.jl --- src/ModelHamiltonians/ModelHamiltonians.jl | 20 +++++++------------- 1 file changed, 7 insertions(+), 13 deletions(-) diff --git a/src/ModelHamiltonians/ModelHamiltonians.jl b/src/ModelHamiltonians/ModelHamiltonians.jl index eebc80ec..f54faea3 100644 --- a/src/ModelHamiltonians/ModelHamiltonians.jl +++ b/src/ModelHamiltonians/ModelHamiltonians.jl @@ -82,26 +82,20 @@ function heisenberg(g::AbstractGraph; J1=1, J2=0, h=0) (; J1, J2, h) = map(to_callable, (; J1, J2, h)) ℋ = OpSum() for e in edges(g) - if !iszero(J1(e)) - ℋ += J1(e) / 2, "S+", src(e), "S-", dst(e) - ℋ += J1(e) / 2, "S-", src(e), "S+", dst(e) - ℋ += J1(e), "Sz", src(e), "Sz", dst(e) - end + ℋ += J1(e) / 2, "S+", src(e), "S-", dst(e) + ℋ += J1(e) / 2, "S-", src(e), "S+", dst(e) + ℋ += J1(e), "Sz", src(e), "Sz", dst(e) end for v in vertices(g) for nn in next_nearest_neighbors(g, v) e = edgetype(g)(v, nn) - if !iszero(J2(e)) - ℋ += J2(e) / 2, "S+", src(e), "S-", dst(e) - ℋ += J2(e) / 2, "S-", src(e), "S+", dst(e) - ℋ += J2(e), "Sz", src(e), "Sz", dst(e) - end + ℋ += J2(e) / 2, "S+", src(e), "S-", dst(e) + ℋ += J2(e) / 2, "S-", src(e), "S+", dst(e) + ℋ += J2(e), "Sz", src(e), "Sz", dst(e) end end for v in vertices(g) - if !iszero(h(v)) - ℋ += h(v), "Sz", v - end + ℋ += h(v), "Sz", v end return ℋ end From 6951df4d347bd8d2673e3edfacfa6ca93aba05a6 Mon Sep 17 00:00:00 2001 From: Joey Date: Fri, 12 Apr 2024 17:43:56 -0400 Subject: [PATCH 06/23] Remove utils.jl from include list --- src/ITensorNetworks.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/ITensorNetworks.jl b/src/ITensorNetworks.jl index d908924e..fec86bf8 100644 --- a/src/ITensorNetworks.jl +++ b/src/ITensorNetworks.jl @@ -35,7 +35,6 @@ include("formnetworks/quadraticformnetwork.jl") include("caches/beliefpropagationcache.jl") include("contraction_tree_to_graph.jl") include("gauging.jl") -include("utils.jl") include("ITensorsExtensions/ITensorsExtensions.jl") include("solvers/local_solvers/eigsolve.jl") include("solvers/local_solvers/exponentiate.jl") From fc687afc24b8a989904d21a7d2e84a5f30e424ad Mon Sep 17 00:00:00 2001 From: Joey Date: Fri, 12 Apr 2024 17:44:43 -0400 Subject: [PATCH 07/23] Add utils.jl back into include list --- src/ITensorNetworks.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/ITensorNetworks.jl b/src/ITensorNetworks.jl index fec86bf8..d908924e 100644 --- a/src/ITensorNetworks.jl +++ b/src/ITensorNetworks.jl @@ -35,6 +35,7 @@ include("formnetworks/quadraticformnetwork.jl") include("caches/beliefpropagationcache.jl") include("contraction_tree_to_graph.jl") include("gauging.jl") +include("utils.jl") include("ITensorsExtensions/ITensorsExtensions.jl") include("solvers/local_solvers/eigsolve.jl") include("solvers/local_solvers/exponentiate.jl") From 5b76d3451afaeb04373cf2464f0d853a32dbe18e Mon Sep 17 00:00:00 2001 From: Joey Date: Fri, 12 Apr 2024 17:51:05 -0400 Subject: [PATCH 08/23] Removed ishermitian flag --- src/ITensorsExtensions/ITensorsExtensions.jl | 1 - src/apply.jl | 16 ++++++++-------- 2 files changed, 8 insertions(+), 9 deletions(-) diff --git a/src/ITensorsExtensions/ITensorsExtensions.jl b/src/ITensorsExtensions/ITensorsExtensions.jl index 37d86757..8355f346 100644 --- a/src/ITensorsExtensions/ITensorsExtensions.jl +++ b/src/ITensorsExtensions/ITensorsExtensions.jl @@ -58,7 +58,6 @@ function map_itensor( f::Function, A::ITensor, lind=first(inds(A)); - ishermitian=false, regularization=nothing, kwargs..., ) diff --git a/src/apply.jl b/src/apply.jl index 68a9d83e..17416856 100644 --- a/src/apply.jl +++ b/src/apply.jl @@ -87,21 +87,21 @@ function simple_update_bp_full(o, ψ, v⃗; envs, (singular_values!)=nothing, ap envs_v1 = filter(env -> hascommoninds(env, ψ[v⃗[1]]), envs) envs_v2 = filter(env -> hascommoninds(env, ψ[v⃗[2]]), envs) sqrt_envs_v1 = [ - ITensorsExtensions.map_itensor(sqrt, env; ishermitian=true, cutoff, regularization) for + ITensorsExtensions.map_itensor(sqrt, env; cutoff, regularization) for env in envs_v1 ] sqrt_envs_v2 = [ - ITensorsExtensions.map_itensor(sqrt, env; ishermitian=true, cutoff, regularization) for + ITensorsExtensions.map_itensor(sqrt, env; cutoff, regularization) for env in envs_v2 ] inv_sqrt_envs_v1 = [ ITensorsExtensions.map_itensor( - inv ∘ sqrt, env; ishermitian=true, cutoff, regularization + inv ∘ sqrt, env; cutoff, regularization ) for env in envs_v1 ] inv_sqrt_envs_v2 = [ ITensorsExtensions.map_itensor( - inv ∘ sqrt, env; ishermitian=true, cutoff, regularization + inv ∘ sqrt, env; cutoff, regularization ) for env in envs_v2 ] ψᵥ₁ᵥ₂_tn = [ψ[v⃗[1]]; ψ[v⃗[2]]; sqrt_envs_v1; sqrt_envs_v2] @@ -135,21 +135,21 @@ function simple_update_bp(o, ψ, v⃗; envs, (singular_values!)=nothing, apply_k envs_v1 = filter(env -> hascommoninds(env, ψ[v⃗[1]]), envs) envs_v2 = filter(env -> hascommoninds(env, ψ[v⃗[2]]), envs) sqrt_envs_v1 = [ - ITensorsExtensions.map_itensor(sqrt, env; ishermitian=true, cutoff, regularization) for + ITensorsExtensions.map_itensor(sqrt, env; cutoff, regularization) for env in envs_v1 ] sqrt_envs_v2 = [ - ITensorsExtensions.map_itensor(sqrt, env; ishermitian=true, cutoff, regularization) for + ITensorsExtensions.map_itensor(sqrt, env; cutoff, regularization) for env in envs_v2 ] inv_sqrt_envs_v1 = [ ITensorsExtensions.map_itensor( - inv ∘ sqrt, env; ishermitian=true, cutoff, regularization + inv ∘ sqrt, env; cutoff, regularization ) for env in envs_v1 ] inv_sqrt_envs_v2 = [ ITensorsExtensions.map_itensor( - inv ∘ sqrt, env; ishermitian=true, cutoff, regularization + inv ∘ sqrt, env; cutoff, regularization ) for env in envs_v2 ] ψᵥ₁ = contract([ψ[v⃗[1]]; sqrt_envs_v1]) From 0dc60165e8fa48f35980318f6db1c85e0359b63a Mon Sep 17 00:00:00 2001 From: Joey Date: Fri, 12 Apr 2024 18:16:53 -0400 Subject: [PATCH 09/23] Fix Bug. Lower computational expense of test_apply.jl --- src/ITensorsExtensions/ITensorsExtensions.jl | 11 +++---- src/apply.jl | 32 ++++++++------------ test/test_apply.jl | 2 +- test/test_itensorsextensions.jl | 6 ++-- 4 files changed, 20 insertions(+), 31 deletions(-) diff --git a/src/ITensorsExtensions/ITensorsExtensions.jl b/src/ITensorsExtensions/ITensorsExtensions.jl index 8355f346..4680d466 100644 --- a/src/ITensorsExtensions/ITensorsExtensions.jl +++ b/src/ITensorsExtensions/ITensorsExtensions.jl @@ -55,20 +55,17 @@ pinv_diag(it::ITensor) = map_diag(pinv, it) pinvsqrt_diag(it::ITensor) = map_diag(pinv ∘ sqrt, it) function map_itensor( - f::Function, - A::ITensor, - lind=first(inds(A)); - regularization=nothing, - kwargs..., + f::Function, A::ITensor, lind=first(inds(A)); regularization=nothing, kwargs... ) USV = svd(A, lind; kwargs...) U, S, V, spec, u, v = USV if !isnothing(regularization) - f = s -> f(s + regularization) + S = map_diag(s -> f(s + regularization), S) + else + S = map_diag(s -> f(s), S) end - S = map_diag(s -> f(s), S) sqrtDL, δᵤᵥ, sqrtDR = sqrt_decomp(S, u, v) sqrtDR = denseblocks(sqrtDR) * denseblocks(δᵤᵥ) L, R = U * sqrtDL, V * sqrtDR diff --git a/src/apply.jl b/src/apply.jl index 17416856..0062b3f1 100644 --- a/src/apply.jl +++ b/src/apply.jl @@ -87,22 +87,18 @@ function simple_update_bp_full(o, ψ, v⃗; envs, (singular_values!)=nothing, ap envs_v1 = filter(env -> hascommoninds(env, ψ[v⃗[1]]), envs) envs_v2 = filter(env -> hascommoninds(env, ψ[v⃗[2]]), envs) sqrt_envs_v1 = [ - ITensorsExtensions.map_itensor(sqrt, env; cutoff, regularization) for - env in envs_v1 + ITensorsExtensions.map_itensor(sqrt, env; cutoff, regularization) for env in envs_v1 ] sqrt_envs_v2 = [ - ITensorsExtensions.map_itensor(sqrt, env; cutoff, regularization) for - env in envs_v2 + ITensorsExtensions.map_itensor(sqrt, env; cutoff, regularization) for env in envs_v2 ] inv_sqrt_envs_v1 = [ - ITensorsExtensions.map_itensor( - inv ∘ sqrt, env; cutoff, regularization - ) for env in envs_v1 + ITensorsExtensions.map_itensor(inv ∘ sqrt, env; cutoff, regularization) for + env in envs_v1 ] inv_sqrt_envs_v2 = [ - ITensorsExtensions.map_itensor( - inv ∘ sqrt, env; cutoff, regularization - ) for env in envs_v2 + ITensorsExtensions.map_itensor(inv ∘ sqrt, env; cutoff, regularization) for + env in envs_v2 ] ψᵥ₁ᵥ₂_tn = [ψ[v⃗[1]]; ψ[v⃗[2]]; sqrt_envs_v1; sqrt_envs_v2] ψᵥ₁ᵥ₂ = contract(ψᵥ₁ᵥ₂_tn; sequence=contraction_sequence(ψᵥ₁ᵥ₂_tn; alg="optimal")) @@ -135,22 +131,18 @@ function simple_update_bp(o, ψ, v⃗; envs, (singular_values!)=nothing, apply_k envs_v1 = filter(env -> hascommoninds(env, ψ[v⃗[1]]), envs) envs_v2 = filter(env -> hascommoninds(env, ψ[v⃗[2]]), envs) sqrt_envs_v1 = [ - ITensorsExtensions.map_itensor(sqrt, env; cutoff, regularization) for - env in envs_v1 + ITensorsExtensions.map_itensor(sqrt, env; cutoff, regularization) for env in envs_v1 ] sqrt_envs_v2 = [ - ITensorsExtensions.map_itensor(sqrt, env; cutoff, regularization) for - env in envs_v2 + ITensorsExtensions.map_itensor(sqrt, env; cutoff, regularization) for env in envs_v2 ] inv_sqrt_envs_v1 = [ - ITensorsExtensions.map_itensor( - inv ∘ sqrt, env; cutoff, regularization - ) for env in envs_v1 + ITensorsExtensions.map_itensor(inv ∘ sqrt, env; cutoff, regularization) for + env in envs_v1 ] inv_sqrt_envs_v2 = [ - ITensorsExtensions.map_itensor( - inv ∘ sqrt, env; cutoff, regularization - ) for env in envs_v2 + ITensorsExtensions.map_itensor(inv ∘ sqrt, env; cutoff, regularization) for + env in envs_v2 ] ψᵥ₁ = contract([ψ[v⃗[1]]; sqrt_envs_v1]) ψᵥ₂ = contract([ψ[v⃗[2]]; sqrt_envs_v2]) diff --git a/test/test_apply.jl b/test/test_apply.jl index d4c408e0..fab04ceb 100644 --- a/test/test_apply.jl +++ b/test/test_apply.jl @@ -19,7 +19,7 @@ using Test: @test, @testset @testset "apply" begin Random.seed!(5623) - g_dims = (2, 3) + g_dims = (2, 2) n = prod(g_dims) g = named_grid(g_dims) s = siteinds("S=1/2", g) diff --git a/test/test_itensorsextensions.jl b/test/test_itensorsextensions.jl index b966594c..5dc05ce4 100644 --- a/test/test_itensorsextensions.jl +++ b/test/test_itensorsextensions.jl @@ -29,9 +29,9 @@ Random.seed!(1234) A = randn(eltype, (n, n)) A = A * A' P = ITensor(A, i, j) - sqrtP = map_itensor(sqrt, P; ishermitian=true) - inv_P = dag(map_itensor(inv, P; ishermitian=true)) - inv_sqrtP = dag(map_itensor(inv ∘ sqrt, P; ishermitian=true)) + sqrtP = map_itensor(sqrt, P) + inv_P = dag(map_itensor(inv, P)) + inv_sqrtP = dag(map_itensor(inv ∘ sqrt, P)) sqrtPdag = replaceind(dag(sqrtP), i, i') P2 = replaceind(sqrtP * sqrtPdag, i', j) From 513921b3de05950b6b2017f50c94ad67a346e36d Mon Sep 17 00:00:00 2001 From: Joey Date: Mon, 15 Apr 2024 09:08:09 -0400 Subject: [PATCH 10/23] Switch to Eigen. No regularization. Better name --- src/ITensorsExtensions/ITensorsExtensions.jl | 30 +++++++++++--------- src/apply.jl | 22 +++++++------- src/caches/beliefpropagationcache.jl | 2 +- test/test_itensornetwork.jl | 2 +- test/test_itensorsextensions.jl | 8 +++--- 5 files changed, 35 insertions(+), 29 deletions(-) diff --git a/src/ITensorsExtensions/ITensorsExtensions.jl b/src/ITensorsExtensions/ITensorsExtensions.jl index 4680d466..939426bc 100644 --- a/src/ITensorsExtensions/ITensorsExtensions.jl +++ b/src/ITensorsExtensions/ITensorsExtensions.jl @@ -54,22 +54,26 @@ invsqrt_diag(it::ITensor) = map_diag(inv ∘ sqrt, it) pinv_diag(it::ITensor) = map_diag(pinv, it) pinvsqrt_diag(it::ITensor) = map_diag(pinv ∘ sqrt, it) -function map_itensor( - f::Function, A::ITensor, lind=first(inds(A)); regularization=nothing, kwargs... +function map_eigenvalues( + f::Function, A::ITensor, linds=Index[first(inds(A))]; ishermitian=false, kwargs... ) - USV = svd(A, lind; kwargs...) - U, S, V, spec, u, v = USV - - if !isnothing(regularization) - S = map_diag(s -> f(s + regularization), S) - else - S = map_diag(s -> f(s), S) + if isdiag(A) + return map_diag(s -> f(s), A) end - sqrtDL, δᵤᵥ, sqrtDR = sqrt_decomp(S, u, v) - sqrtDR = denseblocks(sqrtDR) * denseblocks(δᵤᵥ) - L, R = U * sqrtDL, V * sqrtDR - return L * R + @assert ishermitian + rinds = setdiff(inds(A), linds) + + D, U = eigen(A, linds, rinds; ishermitian, kwargs...) + ul, ur = noncommonind(D, U), commonind(D, U) + ulnew = sim(ul) + + Ul = replaceinds(U, (rinds..., ur), (linds..., ulnew)) + + D_mapped = map_diag(f, D) + D_mapped = replaceind(D_mapped, ul, ulnew) + + return Ul * D_mapped * dag(U) end # Analagous to `denseblocks`. diff --git a/src/apply.jl b/src/apply.jl index 0062b3f1..3f576a38 100644 --- a/src/apply.jl +++ b/src/apply.jl @@ -83,21 +83,22 @@ end function simple_update_bp_full(o, ψ, v⃗; envs, (singular_values!)=nothing, apply_kwargs...) cutoff = 10 * eps(real(scalartype(ψ))) - regularization = 10 * eps(real(scalartype(ψ))) envs_v1 = filter(env -> hascommoninds(env, ψ[v⃗[1]]), envs) envs_v2 = filter(env -> hascommoninds(env, ψ[v⃗[2]]), envs) sqrt_envs_v1 = [ - ITensorsExtensions.map_itensor(sqrt, env; cutoff, regularization) for env in envs_v1 + ITensorsExtensions.map_eigenvalues(sqrt, env; cutoff, ishermitian=true) for + env in envs_v1 ] sqrt_envs_v2 = [ - ITensorsExtensions.map_itensor(sqrt, env; cutoff, regularization) for env in envs_v2 + ITensorsExtensions.map_eigenvalues(sqrt, env; cutoff, ishermitian=true) for + env in envs_v2 ] inv_sqrt_envs_v1 = [ - ITensorsExtensions.map_itensor(inv ∘ sqrt, env; cutoff, regularization) for + ITensorsExtensions.map_eigenvalues(inv ∘ sqrt, env; cutoff, ishermitian=true) for env in envs_v1 ] inv_sqrt_envs_v2 = [ - ITensorsExtensions.map_itensor(inv ∘ sqrt, env; cutoff, regularization) for + ITensorsExtensions.map_eigenvalues(inv ∘ sqrt, env; cutoff, ishermitian=true) for env in envs_v2 ] ψᵥ₁ᵥ₂_tn = [ψ[v⃗[1]]; ψ[v⃗[2]]; sqrt_envs_v1; sqrt_envs_v2] @@ -127,21 +128,22 @@ end # Reduced version function simple_update_bp(o, ψ, v⃗; envs, (singular_values!)=nothing, apply_kwargs...) cutoff = 10 * eps(real(scalartype(ψ))) - regularization = 10 * eps(real(scalartype(ψ))) envs_v1 = filter(env -> hascommoninds(env, ψ[v⃗[1]]), envs) envs_v2 = filter(env -> hascommoninds(env, ψ[v⃗[2]]), envs) sqrt_envs_v1 = [ - ITensorsExtensions.map_itensor(sqrt, env; cutoff, regularization) for env in envs_v1 + ITensorsExtensions.map_eigenvalues(sqrt, env; cutoff, ishermitian=true) for + env in envs_v1 ] sqrt_envs_v2 = [ - ITensorsExtensions.map_itensor(sqrt, env; cutoff, regularization) for env in envs_v2 + ITensorsExtensions.map_eigenvalues(sqrt, env; cutoff, ishermitian=true) for + env in envs_v2 ] inv_sqrt_envs_v1 = [ - ITensorsExtensions.map_itensor(inv ∘ sqrt, env; cutoff, regularization) for + ITensorsExtensions.map_eigenvalues(inv ∘ sqrt, env; cutoff, ishermitian=true) for env in envs_v1 ] inv_sqrt_envs_v2 = [ - ITensorsExtensions.map_itensor(inv ∘ sqrt, env; cutoff, regularization) for + ITensorsExtensions.map_eigenvalues(inv ∘ sqrt, env; cutoff, ishermitian=true) for env in envs_v2 ] ψᵥ₁ = contract([ψ[v⃗[1]]; sqrt_envs_v1]) diff --git a/src/caches/beliefpropagationcache.jl b/src/caches/beliefpropagationcache.jl index 35769012..282b9ee8 100644 --- a/src/caches/beliefpropagationcache.jl +++ b/src/caches/beliefpropagationcache.jl @@ -8,7 +8,7 @@ using ITensors: dir using ITensors.ITensorMPS: ITensorMPS using NamedGraphs: boundary_partitionedges, partitionvertices, partitionedges -default_message(inds_e) = ITensor[denseblocks(delta(inds_e))] +default_message(inds_e) = ITensor[denseblocks(delta(i)) for i in inds_e] default_messages(ptn::PartitionedGraph) = Dictionary() default_message_norm(m::ITensor) = norm(m) function default_message_update(contract_list::Vector{ITensor}; kwargs...) diff --git a/test/test_itensornetwork.jl b/test/test_itensornetwork.jl index cc5c6217..7012a1b5 100644 --- a/test/test_itensornetwork.jl +++ b/test/test_itensornetwork.jl @@ -34,6 +34,7 @@ using ITensors: scalartype, sim, uniqueinds +using ITensors.NDTensors: NDTensors, dim using ITensorNetworks: ITensorNetworks, ⊗, @@ -53,7 +54,6 @@ using ITensorNetworks: ttn using LinearAlgebra: factorize using NamedGraphs: NamedEdge, incident_edges, named_comb_tree, named_grid -using NDTensors: NDTensors, dim using Random: Random, randn! using Test: @test, @test_broken, @testset diff --git a/test/test_itensorsextensions.jl b/test/test_itensorsextensions.jl index 5dc05ce4..8a7edf93 100644 --- a/test/test_itensorsextensions.jl +++ b/test/test_itensorsextensions.jl @@ -14,7 +14,7 @@ using ITensors: replaceinds, dir, array -using ITensorNetworks.ITensorsExtensions: map_itensor +using ITensorNetworks.ITensorsExtensions: map_eigenvalues using ITensorNetworks: siteinds, random_tensornetwork using NamedGraphs: named_grid using Random @@ -29,9 +29,9 @@ Random.seed!(1234) A = randn(eltype, (n, n)) A = A * A' P = ITensor(A, i, j) - sqrtP = map_itensor(sqrt, P) - inv_P = dag(map_itensor(inv, P)) - inv_sqrtP = dag(map_itensor(inv ∘ sqrt, P)) + sqrtP = map_eigenvalues(sqrt, P) + inv_P = dag(map_eigenvalues(inv, P)) + inv_sqrtP = dag(map_eigenvalues(inv ∘ sqrt, P)) sqrtPdag = replaceind(dag(sqrtP), i, i') P2 = replaceind(sqrtP * sqrtPdag, i', j) From 8e15f0b4a0b4f786270e8c26ddc3a0ecbb283df5 Mon Sep 17 00:00:00 2001 From: Joey Date: Mon, 15 Apr 2024 10:06:35 -0400 Subject: [PATCH 11/23] Separate logic out in map_eigvals and use default eigen inds setting --- src/ITensorsExtensions/ITensorsExtensions.jl | 26 ++++++------- src/apply.jl | 40 ++++++++++++-------- test/test_itensorsextensions.jl | 9 +++-- 3 files changed, 42 insertions(+), 33 deletions(-) diff --git a/src/ITensorsExtensions/ITensorsExtensions.jl b/src/ITensorsExtensions/ITensorsExtensions.jl index 939426bc..185c87e5 100644 --- a/src/ITensorsExtensions/ITensorsExtensions.jl +++ b/src/ITensorsExtensions/ITensorsExtensions.jl @@ -54,26 +54,26 @@ invsqrt_diag(it::ITensor) = map_diag(inv ∘ sqrt, it) pinv_diag(it::ITensor) = map_diag(pinv, it) pinvsqrt_diag(it::ITensor) = map_diag(pinv ∘ sqrt, it) -function map_eigenvalues( - f::Function, A::ITensor, linds=Index[first(inds(A))]; ishermitian=false, kwargs... -) - if isdiag(A) - return map_diag(s -> f(s), A) - end - +#TODO: Make this work for non-hermitian A +function eigendecomp(A::ITensor, linds, rinds; ishermitian=false, kwargs...) @assert ishermitian - rinds = setdiff(inds(A), linds) - D, U = eigen(A, linds, rinds; ishermitian, kwargs...) ul, ur = noncommonind(D, U), commonind(D, U) ulnew = sim(ul) - Ul = replaceinds(U, (rinds..., ur), (linds..., ulnew)) + Ul = replaceinds(U, vcat(rinds, ur), vcat(linds, ulnew)) + D = replaceind(D, ul, ulnew) + return Ul, D, dag(U) +end + +function map_eigvals(f::Function, A::ITensor, inds...; ishermitian=false, kwargs...) + if isdiag(A) + return map_diag(f, A) + end - D_mapped = map_diag(f, D) - D_mapped = replaceind(D_mapped, ul, ulnew) + Ul, D, Ur = eigendecomp(A, inds...; ishermitian, kwargs...) - return Ul * D_mapped * dag(U) + return Ul * map_diag(f, D) * Ur end # Analagous to `denseblocks`. diff --git a/src/apply.jl b/src/apply.jl index 3f576a38..1f70bee9 100644 --- a/src/apply.jl +++ b/src/apply.jl @@ -86,20 +86,24 @@ function simple_update_bp_full(o, ψ, v⃗; envs, (singular_values!)=nothing, ap envs_v1 = filter(env -> hascommoninds(env, ψ[v⃗[1]]), envs) envs_v2 = filter(env -> hascommoninds(env, ψ[v⃗[2]]), envs) sqrt_envs_v1 = [ - ITensorsExtensions.map_eigenvalues(sqrt, env; cutoff, ishermitian=true) for - env in envs_v1 + ITensorsExtensions.map_eigvals( + sqrt, env, first(inds(env)), last(inds(env)); cutoff, ishermitian=true + ) for env in envs_v1 ] sqrt_envs_v2 = [ - ITensorsExtensions.map_eigenvalues(sqrt, env; cutoff, ishermitian=true) for - env in envs_v2 + ITensorsExtensions.map_eigvals( + sqrt, env, first(inds(env)), last(inds(env)); cutoff, ishermitian=true + ) for env in envs_v2 ] inv_sqrt_envs_v1 = [ - ITensorsExtensions.map_eigenvalues(inv ∘ sqrt, env; cutoff, ishermitian=true) for - env in envs_v1 + ITensorsExtensions.map_eigvals( + inv ∘ sqrt, env, first(inds(env)), last(inds(env)); cutoff, ishermitian=true + ) for env in envs_v1 ] inv_sqrt_envs_v2 = [ - ITensorsExtensions.map_eigenvalues(inv ∘ sqrt, env; cutoff, ishermitian=true) for - env in envs_v2 + ITensorsExtensions.map_eigvals( + inv ∘ sqrt, env, first(inds(env)), last(inds(env)); cutoff, ishermitian=true + ) for env in envs_v2 ] ψᵥ₁ᵥ₂_tn = [ψ[v⃗[1]]; ψ[v⃗[2]]; sqrt_envs_v1; sqrt_envs_v2] ψᵥ₁ᵥ₂ = contract(ψᵥ₁ᵥ₂_tn; sequence=contraction_sequence(ψᵥ₁ᵥ₂_tn; alg="optimal")) @@ -131,20 +135,24 @@ function simple_update_bp(o, ψ, v⃗; envs, (singular_values!)=nothing, apply_k envs_v1 = filter(env -> hascommoninds(env, ψ[v⃗[1]]), envs) envs_v2 = filter(env -> hascommoninds(env, ψ[v⃗[2]]), envs) sqrt_envs_v1 = [ - ITensorsExtensions.map_eigenvalues(sqrt, env; cutoff, ishermitian=true) for - env in envs_v1 + ITensorsExtensions.map_eigvals( + sqrt, env, first(inds(env)), last(inds(env)); cutoff, ishermitian=true + ) for env in envs_v1 ] sqrt_envs_v2 = [ - ITensorsExtensions.map_eigenvalues(sqrt, env; cutoff, ishermitian=true) for - env in envs_v2 + ITensorsExtensions.map_eigvals( + sqrt, env, first(inds(env)), last(inds(env)); cutoff, ishermitian=true + ) for env in envs_v2 ] inv_sqrt_envs_v1 = [ - ITensorsExtensions.map_eigenvalues(inv ∘ sqrt, env; cutoff, ishermitian=true) for - env in envs_v1 + ITensorsExtensions.map_eigvals( + inv ∘ sqrt, env, first(inds(env)), last(inds(env)); cutoff, ishermitian=true + ) for env in envs_v1 ] inv_sqrt_envs_v2 = [ - ITensorsExtensions.map_eigenvalues(inv ∘ sqrt, env; cutoff, ishermitian=true) for - env in envs_v2 + ITensorsExtensions.map_eigvals( + inv ∘ sqrt, env, first(inds(env)), last(inds(env)); cutoff, ishermitian=true + ) for env in envs_v2 ] ψᵥ₁ = contract([ψ[v⃗[1]]; sqrt_envs_v1]) ψᵥ₂ = contract([ψ[v⃗[2]]; sqrt_envs_v2]) diff --git a/test/test_itensorsextensions.jl b/test/test_itensorsextensions.jl index 8a7edf93..bf5db060 100644 --- a/test/test_itensorsextensions.jl +++ b/test/test_itensorsextensions.jl @@ -14,7 +14,7 @@ using ITensors: replaceinds, dir, array -using ITensorNetworks.ITensorsExtensions: map_eigenvalues +using ITensorNetworks.ITensorsExtensions: map_eigvals using ITensorNetworks: siteinds, random_tensornetwork using NamedGraphs: named_grid using Random @@ -26,12 +26,13 @@ Random.seed!(1234) for eltype in [Float64, ComplexF64] for n in [2, 3, 5, 10] i, j = Index(n, "i"), Index(n, "j") + linds, rinds = Index[i], Index[j] A = randn(eltype, (n, n)) A = A * A' P = ITensor(A, i, j) - sqrtP = map_eigenvalues(sqrt, P) - inv_P = dag(map_eigenvalues(inv, P)) - inv_sqrtP = dag(map_eigenvalues(inv ∘ sqrt, P)) + sqrtP = map_eigvals(sqrt, P, linds, rinds; ishermitian=true) + inv_P = dag(map_eigvals(inv, P, linds, rinds; ishermitian=true)) + inv_sqrtP = dag(map_eigvals(inv ∘ sqrt, P, linds, rinds; ishermitian=true)) sqrtPdag = replaceind(dag(sqrtP), i, i') P2 = replaceind(sqrtP * sqrtPdag, i', j) From 5d0477dd01df2bf7913656cf43f979d6c63f43f3 Mon Sep 17 00:00:00 2001 From: Joey Date: Mon, 15 Apr 2024 11:00:07 -0400 Subject: [PATCH 12/23] Add QNS to test of map_eigvals --- src/ITensorsExtensions/ITensorsExtensions.jl | 4 +- test/test_itensorsextensions.jl | 72 ++++++++++++++------ 2 files changed, 51 insertions(+), 25 deletions(-) diff --git a/src/ITensorsExtensions/ITensorsExtensions.jl b/src/ITensorsExtensions/ITensorsExtensions.jl index 185c87e5..039ce20f 100644 --- a/src/ITensorsExtensions/ITensorsExtensions.jl +++ b/src/ITensorsExtensions/ITensorsExtensions.jl @@ -59,10 +59,8 @@ function eigendecomp(A::ITensor, linds, rinds; ishermitian=false, kwargs...) @assert ishermitian D, U = eigen(A, linds, rinds; ishermitian, kwargs...) ul, ur = noncommonind(D, U), commonind(D, U) - ulnew = sim(ul) + Ul = replaceinds(U, vcat(rinds, ur), vcat(linds, ul)) - Ul = replaceinds(U, vcat(rinds, ur), vcat(linds, ulnew)) - D = replaceind(D, ul, ulnew) return Ul, D, dag(U) end diff --git a/test/test_itensorsextensions.jl b/test/test_itensorsextensions.jl index bf5db060..e7be0a51 100644 --- a/test/test_itensorsextensions.jl +++ b/test/test_itensorsextensions.jl @@ -23,28 +23,56 @@ using LinearAlgebra: ishermitian, isposdef Random.seed!(1234) @testset "ITensorsExtensions" begin - for eltype in [Float64, ComplexF64] - for n in [2, 3, 5, 10] - i, j = Index(n, "i"), Index(n, "j") - linds, rinds = Index[i], Index[j] - A = randn(eltype, (n, n)) - A = A * A' - P = ITensor(A, i, j) - sqrtP = map_eigvals(sqrt, P, linds, rinds; ishermitian=true) - inv_P = dag(map_eigvals(inv, P, linds, rinds; ishermitian=true)) - inv_sqrtP = dag(map_eigvals(inv ∘ sqrt, P, linds, rinds; ishermitian=true)) - - sqrtPdag = replaceind(dag(sqrtP), i, i') - P2 = replaceind(sqrtP * sqrtPdag, i', j) - @test P2 ≈ P atol = 1e-12 - - invP = replaceind(inv_P, i, i') - I = invP * P - @test I ≈ delta(eltype, inds(I)) atol = 1e-12 - - inv_sqrtP = replaceind(inv_sqrtP, i, i') - I = inv_sqrtP * sqrtP - @test I ≈ delta(eltype, inds(I)) atol = 1e-12 + @testset "Test map eigvals without QNS" begin + for eltype in [Float64, ComplexF64] + for n in [2, 3, 5, 10] + i, j = Index(n, "i"), Index(n, "j") + linds, rinds = Index[i], Index[j] + A = randn(eltype, (n, n)) + A = A * A' + P = ITensor(A, i, j) + sqrtP = map_eigvals(sqrt, P, linds, rinds; ishermitian=true) + inv_P = dag(map_eigvals(inv, P, linds, rinds; ishermitian=true)) + inv_sqrtP = dag(map_eigvals(inv ∘ sqrt, P, linds, rinds; ishermitian=true)) + + sqrtPdag = replaceind(dag(sqrtP), i, i') + P2 = replaceind(sqrtP * sqrtPdag, i', j) + @test P2 ≈ P atol = 1e-12 + + invP = replaceind(inv_P, i, i') + I = invP * P + @test I ≈ delta(eltype, inds(I)) atol = 1e-12 + + inv_sqrtP = replaceind(inv_sqrtP, i, i') + I = inv_sqrtP * sqrtP + @test I ≈ delta(eltype, inds(I)) atol = 1e-12 + end + end + end + + @testset "Test map eigvals with QNS" begin + for eltype in [Float64, ComplexF64] + for n in [2, 3, 5, 10] + i, j = Index.(([QN() => n], [QN() => n])) + A = randomITensor(eltype, i, j) + P = A * prime(dag(A), i) + sqrtP = map_eigvals(sqrt, P, i, i'; ishermitian=true) + inv_P = dag(map_eigvals(inv, P, i, i'; ishermitian=true)) + inv_sqrtP = dag(map_eigvals(inv ∘ sqrt, P, i, i'; ishermitian=true)) + + new_ind = noprime(sim(i')) + sqrtPdag = replaceind(dag(sqrtP), i', new_ind) + P2 = replaceind(sqrtP * sqrtPdag, new_ind, i) + @test P2 ≈ P atol = 1e-12 + + inv_P = replaceind(inv_P, i', new_ind) + I = replaceind(inv_P * P, new_ind, i) + @test I ≈ op("I", i) atol = 1e-12 + + inv_sqrtP = replaceind(inv_sqrtP, i', new_ind) + I = replaceind(inv_sqrtP * sqrtP, new_ind, i) + @test I ≈ op("I", i) atol = 1e-12 + end end end end From 678284f9571dfe1a26e56929f98f8a43f73aaf0d Mon Sep 17 00:00:00 2001 From: Joey Date: Mon, 15 Apr 2024 13:00:13 -0400 Subject: [PATCH 13/23] Tidy namespace in test_itensorsextensions.jl --- test/test_itensorsextensions.jl | 17 ++++++----------- 1 file changed, 6 insertions(+), 11 deletions(-) diff --git a/test/test_itensorsextensions.jl b/test/test_itensorsextensions.jl index e7be0a51..0fa53a42 100644 --- a/test/test_itensorsextensions.jl +++ b/test/test_itensorsextensions.jl @@ -2,24 +2,19 @@ using ITensors: ITensors, ITensor, Index, - randomITensor, + QN, dag, - replaceind, + delta, + inds, noprime, prime, - inds, - delta, - QN, - denseblocks, + randomITensor, + replaceind, replaceinds, - dir, - array + sim using ITensorNetworks.ITensorsExtensions: map_eigvals -using ITensorNetworks: siteinds, random_tensornetwork -using NamedGraphs: named_grid using Random using Test: @test, @testset -using LinearAlgebra: ishermitian, isposdef Random.seed!(1234) @testset "ITensorsExtensions" begin From 27b002b97e99cf68454dbe7e3133aa57270f2820 Mon Sep 17 00:00:00 2001 From: Joseph Tindall <51231103+JoeyT1994@users.noreply.github.com> Date: Mon, 15 Apr 2024 13:26:19 -0400 Subject: [PATCH 14/23] Update test/test_itensorsextensions.jl Co-authored-by: Matt Fishman --- test/test_itensorsextensions.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_itensorsextensions.jl b/test/test_itensorsextensions.jl index 0fa53a42..95a6af0f 100644 --- a/test/test_itensorsextensions.jl +++ b/test/test_itensorsextensions.jl @@ -13,7 +13,7 @@ using ITensors: replaceinds, sim using ITensorNetworks.ITensorsExtensions: map_eigvals -using Random +using Random: Random using Test: @test, @testset Random.seed!(1234) From 8f3f16417a76f69ff290e33fd2561e5c88b6919d Mon Sep 17 00:00:00 2001 From: Joseph Tindall <51231103+JoeyT1994@users.noreply.github.com> Date: Mon, 15 Apr 2024 13:26:25 -0400 Subject: [PATCH 15/23] Update test/test_itensorsextensions.jl Co-authored-by: Matt Fishman --- test/test_itensorsextensions.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/test_itensorsextensions.jl b/test/test_itensorsextensions.jl index 95a6af0f..943b7a96 100644 --- a/test/test_itensorsextensions.jl +++ b/test/test_itensorsextensions.jl @@ -1,3 +1,4 @@ +@eval module $(gensym()) using ITensors: ITensors, ITensor, From 05d6a24a93b2d9c857efb25692d93e6bce6b6309 Mon Sep 17 00:00:00 2001 From: Joseph Tindall <51231103+JoeyT1994@users.noreply.github.com> Date: Mon, 15 Apr 2024 13:26:31 -0400 Subject: [PATCH 16/23] Update test/test_itensorsextensions.jl Co-authored-by: Matt Fishman --- test/test_itensorsextensions.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/test_itensorsextensions.jl b/test/test_itensorsextensions.jl index 943b7a96..1594e985 100644 --- a/test/test_itensorsextensions.jl +++ b/test/test_itensorsextensions.jl @@ -72,3 +72,4 @@ Random.seed!(1234) end end end +end From 4ff8d7740469e52b7adbd9793d555a1ccebefd5d Mon Sep 17 00:00:00 2001 From: Joseph Tindall <51231103+JoeyT1994@users.noreply.github.com> Date: Mon, 15 Apr 2024 13:26:36 -0400 Subject: [PATCH 17/23] Update test/test_itensorsextensions.jl Co-authored-by: Matt Fishman --- test/test_itensorsextensions.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_itensorsextensions.jl b/test/test_itensorsextensions.jl index 1594e985..1f582104 100644 --- a/test/test_itensorsextensions.jl +++ b/test/test_itensorsextensions.jl @@ -20,7 +20,7 @@ using Test: @test, @testset Random.seed!(1234) @testset "ITensorsExtensions" begin @testset "Test map eigvals without QNS" begin - for eltype in [Float64, ComplexF64] + for eltype in (Float64, ComplexF64) for n in [2, 3, 5, 10] i, j = Index(n, "i"), Index(n, "j") linds, rinds = Index[i], Index[j] From e03748b7392e281c55f78fb3aacae929484ace00 Mon Sep 17 00:00:00 2001 From: Joseph Tindall <51231103+JoeyT1994@users.noreply.github.com> Date: Mon, 15 Apr 2024 13:26:42 -0400 Subject: [PATCH 18/23] Update test/test_itensorsextensions.jl Co-authored-by: Matt Fishman --- test/test_itensorsextensions.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_itensorsextensions.jl b/test/test_itensorsextensions.jl index 1f582104..a06f1c23 100644 --- a/test/test_itensorsextensions.jl +++ b/test/test_itensorsextensions.jl @@ -47,7 +47,7 @@ Random.seed!(1234) end @testset "Test map eigvals with QNS" begin - for eltype in [Float64, ComplexF64] + for eltype in (Float64, ComplexF64) for n in [2, 3, 5, 10] i, j = Index.(([QN() => n], [QN() => n])) A = randomITensor(eltype, i, j) From bc3e500c58014d9c0dd5b8dc79a3ce97892e167d Mon Sep 17 00:00:00 2001 From: Joseph Tindall <51231103+JoeyT1994@users.noreply.github.com> Date: Mon, 15 Apr 2024 13:26:48 -0400 Subject: [PATCH 19/23] Update test/test_itensorsextensions.jl Co-authored-by: Matt Fishman --- test/test_itensorsextensions.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/test/test_itensorsextensions.jl b/test/test_itensorsextensions.jl index a06f1c23..5af4a931 100644 --- a/test/test_itensorsextensions.jl +++ b/test/test_itensorsextensions.jl @@ -19,7 +19,9 @@ using Test: @test, @testset Random.seed!(1234) @testset "ITensorsExtensions" begin - @testset "Test map eigvals without QNS" begin + @testset "Test map eigvals without QNS (eltype=$elt, dim=$n)" for elt in (Float32, Float64, Complex{Float32}, Complex{Float64}), n in [2, 3, 5, 10] + ... + end for eltype in (Float64, ComplexF64) for n in [2, 3, 5, 10] i, j = Index(n, "i"), Index(n, "j") From 1e99c3cb86ff1f60e91b3dfcbf6c7b4a87159ca5 Mon Sep 17 00:00:00 2001 From: Joey Date: Mon, 15 Apr 2024 13:37:06 -0400 Subject: [PATCH 20/23] Better test names in test_itensorsextensions --- test/test_itensorsextensions.jl | 89 ++++++++++++++++----------------- 1 file changed, 44 insertions(+), 45 deletions(-) diff --git a/test/test_itensorsextensions.jl b/test/test_itensorsextensions.jl index 5af4a931..b2438780 100644 --- a/test/test_itensorsextensions.jl +++ b/test/test_itensorsextensions.jl @@ -8,6 +8,7 @@ using ITensors: delta, inds, noprime, + op, prime, randomITensor, replaceind, @@ -19,59 +20,57 @@ using Test: @test, @testset Random.seed!(1234) @testset "ITensorsExtensions" begin - @testset "Test map eigvals without QNS (eltype=$elt, dim=$n)" for elt in (Float32, Float64, Complex{Float32}, Complex{Float64}), n in [2, 3, 5, 10] - ... - end - for eltype in (Float64, ComplexF64) - for n in [2, 3, 5, 10] - i, j = Index(n, "i"), Index(n, "j") - linds, rinds = Index[i], Index[j] - A = randn(eltype, (n, n)) - A = A * A' - P = ITensor(A, i, j) - sqrtP = map_eigvals(sqrt, P, linds, rinds; ishermitian=true) - inv_P = dag(map_eigvals(inv, P, linds, rinds; ishermitian=true)) - inv_sqrtP = dag(map_eigvals(inv ∘ sqrt, P, linds, rinds; ishermitian=true)) + @testset "Test map eigvals without QNS (eltype=$elt, dim=$n)" for elt in ( + Float32, Float64, Complex{Float32}, Complex{Float64} + ), + n in (2, 3, 5, 10) + + i, j = Index(n, "i"), Index(n, "j") + linds, rinds = Index[i], Index[j] + A = randn(elt, (n, n)) + A = A * A' + P = ITensor(A, i, j) + sqrtP = map_eigvals(sqrt, P, linds, rinds; ishermitian=true) + inv_P = dag(map_eigvals(inv, P, linds, rinds; ishermitian=true)) + inv_sqrtP = dag(map_eigvals(inv ∘ sqrt, P, linds, rinds; ishermitian=true)) - sqrtPdag = replaceind(dag(sqrtP), i, i') - P2 = replaceind(sqrtP * sqrtPdag, i', j) - @test P2 ≈ P atol = 1e-12 + sqrtPdag = replaceind(dag(sqrtP), i, i') + P2 = replaceind(sqrtP * sqrtPdag, i', j) + @test P2 ≈ P - invP = replaceind(inv_P, i, i') - I = invP * P - @test I ≈ delta(eltype, inds(I)) atol = 1e-12 + invP = replaceind(inv_P, i, i') + I = invP * P + @test I ≈ delta(elt, inds(I)) - inv_sqrtP = replaceind(inv_sqrtP, i, i') - I = inv_sqrtP * sqrtP - @test I ≈ delta(eltype, inds(I)) atol = 1e-12 - end - end + inv_sqrtP = replaceind(inv_sqrtP, i, i') + I = inv_sqrtP * sqrtP + @test I ≈ delta(elt, inds(I)) end - @testset "Test map eigvals with QNS" begin - for eltype in (Float64, ComplexF64) - for n in [2, 3, 5, 10] - i, j = Index.(([QN() => n], [QN() => n])) - A = randomITensor(eltype, i, j) - P = A * prime(dag(A), i) - sqrtP = map_eigvals(sqrt, P, i, i'; ishermitian=true) - inv_P = dag(map_eigvals(inv, P, i, i'; ishermitian=true)) - inv_sqrtP = dag(map_eigvals(inv ∘ sqrt, P, i, i'; ishermitian=true)) + @testset "Test map eigvals with QNS (eltype=$elt, dim=$n)" for elt in ( + Float32, Float64, Complex{Float32}, Complex{Float64} + ), + n in (2, 3, 5, 10) + + i, j = Index.(([QN() => n], [QN() => n])) + A = randomITensor(elt, i, j) + P = A * prime(dag(A), i) + sqrtP = map_eigvals(sqrt, P, i, i'; ishermitian=true) + inv_P = dag(map_eigvals(inv, P, i, i'; ishermitian=true)) + inv_sqrtP = dag(map_eigvals(inv ∘ sqrt, P, i, i'; ishermitian=true)) - new_ind = noprime(sim(i')) - sqrtPdag = replaceind(dag(sqrtP), i', new_ind) - P2 = replaceind(sqrtP * sqrtPdag, new_ind, i) - @test P2 ≈ P atol = 1e-12 + new_ind = noprime(sim(i')) + sqrtPdag = replaceind(dag(sqrtP), i', new_ind) + P2 = replaceind(sqrtP * sqrtPdag, new_ind, i) + @test P2 ≈ P - inv_P = replaceind(inv_P, i', new_ind) - I = replaceind(inv_P * P, new_ind, i) - @test I ≈ op("I", i) atol = 1e-12 + inv_P = replaceind(inv_P, i', new_ind) + I = replaceind(inv_P * P, new_ind, i) + @test I ≈ op("I", i) - inv_sqrtP = replaceind(inv_sqrtP, i', new_ind) - I = replaceind(inv_sqrtP * sqrtP, new_ind, i) - @test I ≈ op("I", i) atol = 1e-12 - end - end + inv_sqrtP = replaceind(inv_sqrtP, i', new_ind) + I = replaceind(inv_sqrtP * sqrtP, new_ind, i) + @test I ≈ op("I", i) end end end From 4fa0d164e5074e0720f89ae9a1234feaf12023c9 Mon Sep 17 00:00:00 2001 From: Joseph Tindall <51231103+JoeyT1994@users.noreply.github.com> Date: Mon, 22 Apr 2024 08:39:04 -0400 Subject: [PATCH 21/23] Update src/apply.jl Co-authored-by: Matt Fishman --- src/apply.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/apply.jl b/src/apply.jl index 1f70bee9..e77258e4 100644 --- a/src/apply.jl +++ b/src/apply.jl @@ -12,13 +12,13 @@ using ITensors: dag, denseblocks, factorize, + factorize_svd, hasqns, isdiag, noprime, prime, replaceind, replaceinds, - factorize_svd, unioninds, uniqueinds using ITensors.ContractionSequenceOptimization: optimal_contraction_sequence From 5b918fe0acd2944f86e3efb83accb06d38527656 Mon Sep 17 00:00:00 2001 From: Joey Date: Mon, 22 Apr 2024 08:47:01 -0400 Subject: [PATCH 22/23] Check envs have ndim == 2 in apply.jl --- src/apply.jl | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/src/apply.jl b/src/apply.jl index e77258e4..9405e550 100644 --- a/src/apply.jl +++ b/src/apply.jl @@ -85,24 +85,25 @@ function simple_update_bp_full(o, ψ, v⃗; envs, (singular_values!)=nothing, ap cutoff = 10 * eps(real(scalartype(ψ))) envs_v1 = filter(env -> hascommoninds(env, ψ[v⃗[1]]), envs) envs_v2 = filter(env -> hascommoninds(env, ψ[v⃗[2]]), envs) + @assert all(ndims(env) == 2 for env in vcat(envs_v1, envs_v2)) sqrt_envs_v1 = [ ITensorsExtensions.map_eigvals( - sqrt, env, first(inds(env)), last(inds(env)); cutoff, ishermitian=true + sqrt, env, inds(env)[1], inds(env)[2]; cutoff, ishermitian=true ) for env in envs_v1 ] sqrt_envs_v2 = [ ITensorsExtensions.map_eigvals( - sqrt, env, first(inds(env)), last(inds(env)); cutoff, ishermitian=true + sqrt, env, inds(env)[1], inds(env)[2]; cutoff, ishermitian=true ) for env in envs_v2 ] inv_sqrt_envs_v1 = [ ITensorsExtensions.map_eigvals( - inv ∘ sqrt, env, first(inds(env)), last(inds(env)); cutoff, ishermitian=true + inv ∘ sqrt, env, inds(env)[1], inds(env)[2]; cutoff, ishermitian=true ) for env in envs_v1 ] inv_sqrt_envs_v2 = [ ITensorsExtensions.map_eigvals( - inv ∘ sqrt, env, first(inds(env)), last(inds(env)); cutoff, ishermitian=true + inv ∘ sqrt, env, inds(env)[1], inds(env)[2]; cutoff, ishermitian=true ) for env in envs_v2 ] ψᵥ₁ᵥ₂_tn = [ψ[v⃗[1]]; ψ[v⃗[2]]; sqrt_envs_v1; sqrt_envs_v2] @@ -134,24 +135,25 @@ function simple_update_bp(o, ψ, v⃗; envs, (singular_values!)=nothing, apply_k cutoff = 10 * eps(real(scalartype(ψ))) envs_v1 = filter(env -> hascommoninds(env, ψ[v⃗[1]]), envs) envs_v2 = filter(env -> hascommoninds(env, ψ[v⃗[2]]), envs) + @assert all(ndims(env) == 2 for env in vcat(envs_v1, envs_v2)) sqrt_envs_v1 = [ ITensorsExtensions.map_eigvals( - sqrt, env, first(inds(env)), last(inds(env)); cutoff, ishermitian=true + sqrt, env, inds(env)[1], inds(env)[2]; cutoff, ishermitian=true ) for env in envs_v1 ] sqrt_envs_v2 = [ ITensorsExtensions.map_eigvals( - sqrt, env, first(inds(env)), last(inds(env)); cutoff, ishermitian=true + sqrt, env, inds(env)[1], inds(env)[2]; cutoff, ishermitian=true ) for env in envs_v2 ] inv_sqrt_envs_v1 = [ ITensorsExtensions.map_eigvals( - inv ∘ sqrt, env, first(inds(env)), last(inds(env)); cutoff, ishermitian=true + inv ∘ sqrt, env, inds(env)[1], inds(env)[2]; cutoff, ishermitian=true ) for env in envs_v1 ] inv_sqrt_envs_v2 = [ ITensorsExtensions.map_eigvals( - inv ∘ sqrt, env, first(inds(env)), last(inds(env)); cutoff, ishermitian=true + inv ∘ sqrt, env, inds(env)[1], inds(env)[2]; cutoff, ishermitian=true ) for env in envs_v2 ] ψᵥ₁ = contract([ψ[v⃗[1]]; sqrt_envs_v1]) From 80ab9fce0bcd28b1268f8e9dd804c1fbec743ec6 Mon Sep 17 00:00:00 2001 From: Joseph Tindall <51231103+JoeyT1994@users.noreply.github.com> Date: Mon, 22 Apr 2024 10:07:51 -0400 Subject: [PATCH 23/23] Update src/ITensorsExtensions/ITensorsExtensions.jl Co-authored-by: Matt Fishman --- src/ITensorsExtensions/ITensorsExtensions.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ITensorsExtensions/ITensorsExtensions.jl b/src/ITensorsExtensions/ITensorsExtensions.jl index 039ce20f..5b58e663 100644 --- a/src/ITensorsExtensions/ITensorsExtensions.jl +++ b/src/ITensorsExtensions/ITensorsExtensions.jl @@ -1,5 +1,5 @@ module ITensorsExtensions -using LinearAlgebra: LinearAlgebra, eigen, ishermitian, pinv +using LinearAlgebra: LinearAlgebra, eigen, pinv using ITensors: ITensor, Index,