Upgrade tree-DMRG to the ITensorBase 0.8 and TensorAlgebra 0.11 stack

The DMRG code predated the INN v0.7 upgrade. Move it onto the released stack:
named arrays and operators now come from ITensorBase rather than NamedDimsArrays
(state, operator, apply, lazy, dimnames, replacedimnames, uniquename), the QR and
SVD factorizations go through MatrixAlgebraKit (qr_compact, svd_trunc), and the
tensor network types are ITensorNetwork and AbstractITensorNetwork. The
VectorInterface methods that KrylovKit needs were dropped here and moved into
ITensorBase (ITensor/ITensorBase.jl#186), so the Aqua
piracy check is back to the full default.

Author: mtfishman

Project.toml

@@ -24,6 +24,10 @@ SimpleTraits = "699a6c99-e7fa-54fc-8d76-47d257e15c1d"
 SplitApplyCombine = "03a91e81-4c3e-53e1-a0a4-9c0c8f19dd66"
 TensorAlgebra = "68bd88dc-f39d-4e12-b2ca-f046b68fcc6a"
 
+[sources.ITensorBase]
+rev = "mf/vectorinterface-itensor"
+url = "https://github.com/ITensor/ITensorBase.jl"
+
 [compat]
 Adapt = "4.3"
 AlgorithmsInterface = "0.1"

src/ITensorNetworksNext.jl

@@ -21,7 +21,6 @@ include("beliefpropagation/normnetwork.jl")
 
 include("apply/apply_operators.jl")
 
-include("nameddimsarrays_extensions.jl")
 include("tensornetworkoperator.jl")
 include("quadraticformnetwork.jl")
 include("orthogonalize.jl")

src/dmrg.jl

@@ -2,10 +2,11 @@ import .AlgorithmsInterfaceExtensions as AIE
 import AlgorithmsInterface as AI
 using Base: @kwdef
 using Graphs: dst, edges, edgetype, src
+using ITensorBase:
+    ITensorBase as ITB, dimnames, operator, replacedimnames, state, uniquename
 using KrylovKit: eigsolve
-using NamedDimsArrays: NamedDimsArrays as NDA, dimnames, replacedimnames
+using MatrixAlgebraKit: svd_trunc
 using NamedGraphs.GraphsExtensions: edge_path
-using TensorAlgebra: TensorAlgebra as TA
 
 # ============================ Top-level entry point =====================================
 
@@ -35,7 +36,7 @@ function dmrg(
         trunc = nothing,
         regions = default_dmrg_regions(ket),
         link_index_map = Dict(
-            ln => randname(ln) for e in edges(ket) for ln in linknames(ket, e)
+            ln => uniquename(ln) for e in edges(ket) for ln in linknames(ket, e)
         ),
         alg = nothing,
         kwargs...
@@ -361,44 +362,42 @@ function region_update_nsite!(::Val{2}, alg::TwoSiteEigsolve, qf, env, region)
     return qf, energy
 end
 
-function svd_split(T, rows; trunc)
-    return isnothing(trunc) ? TA.svd(T, rows) : TA.svd(T, rows; trunc)
-end
+svd_split(T, rows; trunc) = svd_trunc(T, rows; trunc)
 
 # Lowest eigenpair of the effective Hamiltonian operator `H_eff` acting on the named tensor
-# `T` (via `NamedDimsArrays.apply`), using `KrylovKit.eigsolve`.
+# `T` (via `ITensorBase.apply`), using `KrylovKit.eigsolve`.
 function eigsolve_named(H_eff, T, which)
-    vals, vecs = eigsolve(x -> NDA.apply(H_eff, x), T, 1, which; ishermitian = true)
+    vals, vecs = eigsolve(x -> ITB.apply(H_eff, x), T, 1, which; ishermitian = true)
     return real(vals[1]), vecs[1]
 end
 
 # The effective Hamiltonian for a `region` (a vector of vertices) is the projected operator
 # obtained by contracting the operator tensors on the region with the incoming environment
 # messages on the region's boundary. Contracting those gives a single tensor whose codomain
 # (output) names are the region's bra names and whose domain (input) names are its ket names
-# — i.e. exactly a `NamedDimsArrays` operator. So `effective_hamiltonian` returns that
-# operator, and its action on a region ket tensor `T` is `NamedDimsArrays.apply(H, T)` (which
-# contracts the ket names and renames the resulting bra names back to ket names).
+# — i.e. exactly an ITensor operator. So `effective_hamiltonian` returns that operator, and
+# its action on a region ket tensor `T` is `ITensorBase.apply(H, T)` (which contracts the ket
+# names and renames the resulting bra names back to ket names).
 
 """
-    effective_hamiltonian(qf::QuadraticFormNetwork, env, region) -> NamedDimsOperator
+    effective_hamiltonian(qf::QuadraticFormNetwork, env, region) -> ITensorOperator
 
-Effective (projected) Hamiltonian for `region` (a vector of vertices) as a `NamedDimsArrays`
+Effective (projected) Hamiltonian for `region` (a vector of vertices) as an ITensor
 operator: its domain (input) names are the region's ket names and its codomain (output)
 names are the matching bra names. Apply it to a region ket tensor `T` with
-`NamedDimsArrays.apply`. The environment `env` is a `MessageCache` of
+`ITensorBase.apply`. The environment `env` is a `MessageCache` of
 [`quadratic_form_environments`](@ref).
 """
 function effective_hamiltonian(qf::QuadraticFormNetwork, env, region)
     operators = [operator_tensor(qf, v) for v in region]
-    boundary = [NDA.state(m) for m in incoming_edge_data(env, region)]
+    boundary = [state(m) for m in incoming_edge_data(env, region)]
     h = contract_network([operators; boundary])
     sitemap = site_index_map(qf.operator)
     ketnames = [
         n for n in dimnames(h) if haskey(sitemap, n) || haskey(qf.link_index_map, n)
     ]
     branames = [bra_name_map(qf)[n] for n in ketnames]
-    return NDA.operator(h, branames, ketnames)
+    return operator(h, branames, ketnames)
 end
 
 # ============================ Energy-based convergence ==================================

src/nameddimsarrays_extensions.jl

@@ -1,7 +1,7 @@
 using Graphs: dst, src
-using NamedDimsArrays: dimnames, replacedimnames
+using ITensorBase: dimnames, replacedimnames
+using MatrixAlgebraKit: qr_compact
 using NamedGraphs.GraphsExtensions: edge_path, post_order_dfs_edges
-using TensorAlgebra: TensorAlgebra as TA
 
 # Isometric (QR) gauge on a tree. The orthogonality center is *not* stored on the network
 # — it is tracked by the caller. `orthogonalize(state, center)` canonicalizes the whole
@@ -19,7 +19,7 @@ function gauge_move!(state, v, w)
     ln = only(linknames(state, v => w))
     tv = state[v]
     rows = collect(setdiff(dimnames(tv), [ln]))
-    Q, R = TA.qr(tv, rows)
+    Q, R = qr_compact(tv, rows)
     r = only(setdiff(dimnames(Q), rows))
     new_w = R * state[w]
     setindex_preserve_graph!(state, replacedimnames(Q, r => ln), v)

src/orthogonalize.jl

@@ -1,17 +1,16 @@
-using .LazyNamedDimsArrays: lazy
 using DataGraphs: DataGraphs, underlying_graph
 using Graphs: dst, edgetype, neighbors, src, vertices
-using NamedDimsArrays: NamedDimsArrays as NDA, dimnames, replacedimnames
+using ITensorBase: dimnames, lazy, operator, replacedimnames, state
 using NamedGraphs.GraphsExtensions: vertextype
 
-# A lazy `⟨ψ|H|ψ⟩` network: a ket `TensorNetwork`, a `TensorNetworkOperator` (which carries
+# A lazy `⟨ψ|H|ψ⟩` network: a ket `ITensorNetwork`, a `TensorNetworkOperator` (which carries
 # the ket → bra *site* name map), and a forward ket → bra map for the *link* names. The bra
 # layer is derived from the ket (`conj` + index renaming), never stored, so updating a ket
-# tensor is reflected in the bra. As an `AbstractTensorNetwork`, the data on vertex `v` is
+# tensor is reflected in the bra. As an `AbstractITensorNetwork`, the data on vertex `v` is
 # the lazy product `lazy(ket) * lazy(operator) * lazy(bra)`, so the existing
 # `contract_network` / `MessageCache` machinery treats it like any other tensor network.
 struct QuadraticFormNetwork{V, VD, Ket, Operator, LinkMap} <:
-    AbstractTensorNetwork{V, VD}
+    AbstractITensorNetwork{V, VD}
     ket::Ket
     operator::Operator
     link_index_map::LinkMap
@@ -79,22 +78,22 @@ end
 
 function incoming_subtree_messages(messages, graph, v, w)
     return [
-        NDA.state(messages[edgetype(graph)(u, v)]) for
+        state(messages[edgetype(graph)(u, v)]) for
             u in neighbors(graph, v) if u != w
     ]
 end
 
 function environment_operator(message, link_index_map)
     ketnames = [n for n in dimnames(message) if haskey(link_index_map, n)]
     branames = [link_index_map[n] for n in ketnames]
-    return NDA.operator(message, branames, ketnames)
+    return operator(message, branames, ketnames)
 end
 
 """
     quadratic_form_environments(qf::QuadraticFormNetwork; root) -> MessageCache
 
 Exact projected-Hamiltonian environments of `⟨ψ|H|ψ⟩` on a tree, as a `MessageCache` of
-`NamedDimsArrays` operators keyed by directed edges. The message on `v → w` is the
+ITensor operators keyed by directed edges. The message on `v → w` is the
 contraction of the `⟨ψ|H|ψ⟩` subtree on `v`'s side of `(v, w)`, wrapped as an operator
 recording the bra ↔ ket link correspondence (see [`environment_operator`](@ref)).
 """

src/quadraticformnetwork.jl

@@ -1,32 +1,32 @@
 using DataGraphs: DataGraphs, underlying_graph
-using NamedDimsArrays: NamedDimsArrays as NDA
+using ITensorBase: ITensorBase as ITB
 using NamedGraphs.GraphsExtensions: vertextype
 
-# A tensor-network operator: an operator `TensorNetwork` together with the map between its
-# bra-side and ket-side physical names. The tensor-network analogue of a `NamedDimsArrays`
-# operator, with `state` the underlying operator network (mirroring `NamedDimsArrays.state`
-# on a `NamedDimsOperator`) and the tensor-network interface forwarded to it.
-struct TensorNetworkOperator{V, VD, State, SiteMap} <: AbstractTensorNetwork{V, VD}
+# A tensor-network operator: an operator `ITensorNetwork` together with the map between its
+# bra-side and ket-side physical names. The tensor-network analogue of an ITensor operator,
+# with `state` the underlying operator network (mirroring `ITensorBase.state` on an
+# `ITensorOperator`) and the tensor-network interface forwarded to it.
+struct TensorNetworkOperator{V, VD, State, SiteMap} <: AbstractITensorNetwork{V, VD}
     state::State
     site_index_map::SiteMap
 end
 
-NDA.state(o::TensorNetworkOperator) = o.state
+ITB.state(o::TensorNetworkOperator) = o.state
 site_index_map(o::TensorNetworkOperator) = o.site_index_map
 
-DataGraphs.underlying_graph(o::TensorNetworkOperator) = underlying_graph(NDA.state(o))
+DataGraphs.underlying_graph(o::TensorNetworkOperator) = underlying_graph(ITB.state(o))
 function DataGraphs.is_vertex_assigned(o::TensorNetworkOperator, v)
-    return DataGraphs.is_vertex_assigned(NDA.state(o), v)
+    return DataGraphs.is_vertex_assigned(ITB.state(o), v)
 end
 function DataGraphs.is_edge_assigned(o::TensorNetworkOperator, e)
-    return DataGraphs.is_edge_assigned(NDA.state(o), e)
+    return DataGraphs.is_edge_assigned(ITB.state(o), e)
 end
 function DataGraphs.get_vertex_data(o::TensorNetworkOperator, v)
-    return DataGraphs.get_vertex_data(NDA.state(o), v)
+    return DataGraphs.get_vertex_data(ITB.state(o), v)
 end
 
 function Base.copy(o::TensorNetworkOperator)
-    return TensorNetworkOperator(copy(NDA.state(o)), site_index_map(o))
+    return TensorNetworkOperator(copy(ITB.state(o)), site_index_map(o))
 end
 
 function TensorNetworkOperator(state, site_index_map)

src/tensornetworkoperator.jl

@@ -3,9 +3,5 @@ using ITensorNetworksNext: ITensorNetworksNext
 using Test: @testset
 
 @testset "Code quality (Aqua.jl)" begin
-    # Piracy is checked separately as `broken`: `dmrg.jl` temporarily pirates a few
-    # `VectorInterface` methods on `AbstractNamedDimsArray` (needed by `KrylovKit.eigsolve`).
-    # These are slated to move into `NamedDimsArrays`; drop the `broken` marker once they do.
-    Aqua.test_all(ITensorNetworksNext; persistent_tasks = false, piracies = false)
-    Aqua.test_piracies(ITensorNetworksNext; broken = true)
+    Aqua.test_all(ITensorNetworksNext; persistent_tasks = false)
 end

test/test_aqua.jl

@@ -1,9 +1,8 @@
 using Graphs: add_edge!, edges, vertices
-using ITensorBase: Index
-using ITensorNetworksNext: TensorNetwork, TensorNetworkOperator, dmrg, insertlink!
+using ITensorBase: Index, name, nameddims, setname, uniquename
+using ITensorNetworksNext: ITensorNetwork, TensorNetworkOperator, dmrg, insertlink!
 using LinearAlgebra: eigen
 using MatrixAlgebraKit: truncrank
-using NamedDimsArrays: name, nameddims, randname, setname
 using NamedGraphs.NamedGraphGenerators: named_path_graph
 using NamedGraphs: NamedGraph
 using StableRNGs: StableRNG
@@ -34,7 +33,7 @@ function tfim_path_operator(N, sites, sitemap; J, h)
     verts = collect(1:N)
     bond_edges = [(verts[i], verts[i + 1]) for i in 1:(N - 1)]
     bonds = Dict(
-        e => setname(Index(Base.OneTo(3)), randname(name(sites[verts[1]])))
+        e => setname(Index(Base.OneTo(3)), uniquename(name(sites[verts[1]])))
             for e in bond_edges
     )
     W = ising_mpo_tensor(; J, h)
@@ -58,7 +57,7 @@ function tfim_path_operator(N, sites, sitemap; J, h)
             W, (name(bonds[bond_edges[left]]), name(bonds[bond_edges[right]]), k, b)
         )
     end
-    tn = TensorNetwork(g) do v
+    tn = ITensorNetwork(g) do v
         return tensor(v)
     end
     return TensorNetworkOperator(tn, sitemap)
@@ -83,7 +82,7 @@ function tfim_star_operator(sites, sitemap; J, h)
         add_edge!(g, 1, v)
     end
     bonds = Dict(
-        v => setname(Index(Base.OneTo(3)), randname(name(sites[1]))) for v in leaves
+        v => setname(Index(Base.OneTo(3)), uniquename(name(sites[1]))) for v in leaves
     )
     # Center tensor: one bond per leaf, with `-hX` once and `-JZ`/`I` on each leaf channel.
     center = zeros(3, 3, 3, 2, 2)
@@ -105,7 +104,7 @@ function tfim_star_operator(sites, sitemap; J, h)
     end
     tensors =
         Dict(1 => nameddims(center, center_names), (v => leaf_tensor(v) for v in leaves)...)
-    tn = TensorNetwork(g) do v
+    tn = ITensorNetwork(g) do v
         return tensors[v]
     end
     return TensorNetworkOperator(tn, sitemap)
@@ -124,7 +123,7 @@ end
 
 function random_ket(rng, g)
     sites = Dict(v => Index(Base.OneTo(2)) for v in vertices(g))
-    ket = TensorNetwork(NamedGraph(collect(vertices(g)))) do v
+    ket = ITensorNetwork(NamedGraph(collect(vertices(g)))) do v
         return randn(rng, Float64, (sites[v],))
     end
     for e in edges(g)
@@ -138,7 +137,7 @@ end
         rng = StableRNG(8)
         g = named_path_graph(N)
         ket0, sites = random_ket(rng, g)
-        sitemap = Dict(name(sites[v]) => randname(name(sites[v])) for v in vertices(g))
+        sitemap = Dict(name(sites[v]) => uniquename(name(sites[v])) for v in vertices(g))
         operator = tfim_path_operator(N, sites, sitemap; J = 1.0, h = 0.7)
         exact = minimum(eigen(tfim_path_dense(N; J = 1.0, h = 0.7)).values)
 
@@ -151,7 +150,7 @@ end
         rng = StableRNG(8)
         g = named_path_graph(N)
         ket0, sites = random_ket(rng, g)
-        sitemap = Dict(name(sites[v]) => randname(name(sites[v])) for v in vertices(g))
+        sitemap = Dict(name(sites[v]) => uniquename(name(sites[v])) for v in vertices(g))
         operator = tfim_path_operator(N, sites, sitemap; J = 1.0, h = 0.7)
         exact = minimum(eigen(tfim_path_dense(N; J = 1.0, h = 0.7)).values)
 
@@ -165,7 +164,7 @@ end
         rng = StableRNG(8)
         g = named_path_graph(N)
         ket0, sites = random_ket(rng, g)
-        sitemap = Dict(name(sites[v]) => randname(name(sites[v])) for v in vertices(g))
+        sitemap = Dict(name(sites[v]) => uniquename(name(sites[v])) for v in vertices(g))
         operator = tfim_path_operator(N, sites, sitemap; J = 1.0, h = 0.7)
         exact = minimum(eigen(tfim_path_dense(N; J = 1.0, h = 0.7)).values)
 
@@ -180,7 +179,7 @@ end
         rng = StableRNG(8)
         g = named_path_graph(N)
         ket0, sites = random_ket(rng, g)
-        sitemap = Dict(name(sites[v]) => randname(name(sites[v])) for v in vertices(g))
+        sitemap = Dict(name(sites[v]) => uniquename(name(sites[v])) for v in vertices(g))
         operator = tfim_path_operator(N, sites, sitemap; J = 1.0, h = 0.7)
         exact = minimum(eigen(tfim_path_dense(N; J = 1.0, h = 0.7)).values)
 
@@ -197,7 +196,7 @@ end
             add_edge!(g, 1, v)
         end
         ket0, sites = random_ket(rng, g)
-        sitemap = Dict(name(sites[v]) => randname(name(sites[v])) for v in vertices(g))
+        sitemap = Dict(name(sites[v]) => uniquename(name(sites[v])) for v in vertices(g))
         operator = tfim_star_operator(sites, sitemap; J = 1.0, h = 0.6)
         exact = minimum(eigen(tfim_star_dense(; J = 1.0, h = 0.6)).values)
 
@@ -209,7 +208,7 @@ end
         rng = StableRNG(8)
         g = named_path_graph(4)
         ket0, sites = random_ket(rng, g)
-        sitemap = Dict(name(sites[v]) => randname(name(sites[v])) for v in vertices(g))
+        sitemap = Dict(name(sites[v]) => uniquename(name(sites[v])) for v in vertices(g))
         operator = tfim_path_operator(4, sites, sitemap; J = 1.0, h = 0.7)
         @test_throws ArgumentError dmrg(operator, ket0)
     end