From 3ca0649f2f62d4537e339522cf1dc959bda024a3 Mon Sep 17 00:00:00 2001
From: Ben Corbett <corbett5@llnl.gov>
Date: Sun, 26 Apr 2026 00:15:16 -0600
Subject: [PATCH 1/4] Switched to vector of dict.

---
 src/large-graph-mpo-common.jl       |  21 ++---
 src/large-graph-mpo-qr.jl           |   9 ++-
 src/large-graph-mpo-vertex-cover.jl |   9 ++-
 src/sparse_tensor_construction.jl   | 121 +++++++++++++++-------------
 test/tensor-construction-tests.jl   |  36 +++++----
 5 files changed, 108 insertions(+), 88 deletions(-)

diff --git a/src/large-graph-mpo-common.jl b/src/large-graph-mpo-common.jl
index aec8010..09c8756 100644
--- a/src/large-graph-mpo-common.jl
+++ b/src/large-graph-mpo-common.jl
@@ -1,16 +1,15 @@
 """
     BlockSparseMatrix{C}
 
-Dictionary-backed block-sparse matrix representation used for intermediate MPO
+Vector-backed block-sparse matrix representation used for intermediate MPO
 tensor storage.
 
-Keys are `(left_link, right_link)` pairs and values are dense local operator
-matrices for that block. The right link ids are local to the component being
-processed; `at_site!` later returns offsets that place component-local
-right-link ids into the full outgoing MPO bond.
+The outer vector is indexed by component-local `right_link`. Each inner
+dictionary maps `left_link` to the dense local operator matrix for that block.
+`at_site!` later returns offsets that place component-local right-link ids into
+the full outgoing MPO bond.
 """
-BlockSparseMatrix{C} = Dict{Tuple{Int,Int},Matrix{C}}
-# TODO: consider changing to Vector{Dict{Int, Matrix{C}}} from right_link => (left_link => matrix) and use Dictionaries.jl
+BlockSparseMatrix{C} = Vector{Dict{Int,Matrix{C}}}
 
 """
     MPOGraph{N,C,Ti}
@@ -332,6 +331,7 @@ function merge_qn_sectors(
   return new_order, new_qi
 end
 
+# TODO: Move this into large-graph-mpo-qr.jl
 """
     process_single_left_vertex_cc!(
       matrix_of_cc, rank_of_cc, next_edges_of_cc, g, ccs, cc, n, sites, op_cache_vec
@@ -362,7 +362,10 @@ function process_single_left_vertex_cc!(
   lv = left_vertex(g, lv_id)
   local_op = op_cache_vec[n][lv.op_id].matrix
 
-  matrix_element = get!(matrix_of_cc[cc], (lv.link, rank)) do
+  matrix = [Dict{Int,Matrix{ValType}}() for _ in 1:rank]
+  matrix_of_cc[cc] = matrix
+
+  matrix_element = get!(matrix[rank], lv.link) do
     return zeros(ValType, dim(sites[n]), dim(sites[n]))
   end
 
@@ -371,7 +374,7 @@ function process_single_left_vertex_cc!(
   if n == length(sites)
     scaling = only(g.edge_weights_from_left[lv_id])
 
-    for block in values(matrix_of_cc[cc])
+    for column in matrix, block in values(column)
       block .*= scaling
     end
 
diff --git a/src/large-graph-mpo-qr.jl b/src/large-graph-mpo-qr.jl
index d632b17..38bffe0 100644
--- a/src/large-graph-mpo-qr.jl
+++ b/src/large-graph-mpo-qr.jl
@@ -115,7 +115,6 @@ blocks.
       continue
     end
 
-    matrix = matrix_of_cc[cc]
     W, left_map, right_map = get_cc_matrix(g, ccs, cc; clear_edges=true)
 
     ## Compute the decomposition and then free W
@@ -124,6 +123,9 @@ blocks.
 
     rank_of_cc[cc] = rank
 
+    matrix = [Dict{Int,Matrix{ValType}}() for _ in 1:rank]
+    matrix_of_cc[cc] = matrix
+
     ## Form the local transformation tensor.
     for_non_zeros_batch(Q, rank) do weights::AbstractVector{C}, m::Int
       for (i, weight) in enumerate(weights)
@@ -133,7 +135,7 @@ blocks.
         lv = left_vertex(g, left_map[prow[i]])
         local_op = op_cache_vec[n][lv.op_id].matrix
 
-        matrix_element = get!(matrix, (lv.link, m)) do
+        matrix_element = get!(matrix[m], lv.link) do
           return zeros(ValType, dim(sites[n]), dim(sites[n]))
         end
 
@@ -148,7 +150,7 @@ blocks.
     ## If we are at the last site, then R will be a 1x1 matrix containing an overall scaling.
     if n == length(sites)
       scaling = only(R)
-      for block in values(matrix)
+      for column in matrix, block in values(column)
         block .*= scaling
       end
 
@@ -208,4 +210,3 @@ blocks.
 
   return nothing
 end
-
diff --git a/src/large-graph-mpo-vertex-cover.jl b/src/large-graph-mpo-vertex-cover.jl
index 2909a5d..69fff9d 100644
--- a/src/large-graph-mpo-vertex-cover.jl
+++ b/src/large-graph-mpo-vertex-cover.jl
@@ -34,7 +34,6 @@ otherwise this also builds `next_edges_of_cc` for the graph at site `n + 1`.
 
   # TODO: Consider nested multithreading
   Threads.@threads for cc in 1:num_connected_components(ccs)
-    matrix = matrix_of_cc[cc]
     lvs_of_component::Vector{Int} = ccs.lvs_of_component[cc]
     position_of_rvs_in_component = ccs.position_of_rvs_in_component
     rv_size_of_component = ccs.rv_size_of_component[cc]
@@ -45,6 +44,9 @@ otherwise this also builds `next_edges_of_cc` for the graph at site `n + 1`.
     rank = length(left_cover) + length(right_cover)
     rank_of_cc[cc] = rank
 
+    matrix = [Dict{Int,Matrix{ValType}}() for _ in 1:rank]
+    matrix_of_cc[cc] = matrix
+
     ## Construct the tensor from the left cover.
     @inbounds for m in eachindex(left_cover)
       lv_id = lvs_of_component[left_cover[m]]
@@ -53,7 +55,7 @@ otherwise this also builds `next_edges_of_cc` for the graph at site `n + 1`.
 
       matrix_element = zeros(ValType, site_dim, site_dim)
       add_to_local_matrix!(matrix_element, one(ValType), local_op, lv.needs_JW_string)
-      matrix[lv.link, m] = matrix_element
+      matrix[m][lv.link] = matrix_element
     end
 
     ## Construct the tensor from the right cover.
@@ -85,7 +87,7 @@ otherwise this also builds `next_edges_of_cc` for the graph at site `n + 1`.
         for (rv_id, weight) in weighted_edge_iterator(g, lv_id)
           m = right_cover_m[position_of_rvs_in_component[rv_id]]
 
-          matrix_element = get!(matrix, (lv.link, m)) do
+          matrix_element = get!(matrix[m], lv.link) do
             zeros(ValType, site_dim, site_dim)
           end
 
@@ -158,4 +160,3 @@ otherwise this also builds `next_edges_of_cc` for the graph at site `n + 1`.
 
   return nothing
 end
-
diff --git a/src/sparse_tensor_construction.jl b/src/sparse_tensor_construction.jl
index 492e078..344cff4 100644
--- a/src/sparse_tensor_construction.jl
+++ b/src/sparse_tensor_construction.jl
@@ -88,21 +88,23 @@ function _fill_splitblocks!(
   blocks::Vector{Block{4}},
   block_values::Vector{C},
   offsets::Vector{Int},
-  block_sparse_matrices::Vector{Dict{Tuple{Int,Int},Matrix{C}}},
+  block_sparse_matrices::Vector{BlockSparseMatrix{C}},
 ) where {C}
   entry = 0
 
   for (offset, matrix) in zip(offsets, block_sparse_matrices)
-    for ((left_link, right_link), block) in matrix
+    for (right_link, column) in enumerate(matrix)
       shifted_right_link = right_link + offset
 
-      @inbounds for i in axes(block, 1), j in axes(block, 2)
-        value = block[i, j]
-        iszero(value) && continue
+      for (left_link, block) in column
+        @inbounds for site_prime_link in axes(block, 1), site_link in axes(block, 2)
+          value = block[site_prime_link, site_link]
+          iszero(value) && continue
 
-        entry += 1
-        blocks[entry] = Block((left_link, shifted_right_link, i, j))
-        block_values[entry] = value
+          entry += 1
+          blocks[entry] = Block((left_link, shifted_right_link, site_prime_link, site_link))
+          block_values[entry] = value
+        end
       end
     end
   end
@@ -123,7 +125,7 @@ are shifted to global right-link ids by `right_link + offsets[i]`.
 """
 function _to_ITensor_splitblocks(
   offsets::Vector{Int},
-  block_sparse_matrices::Vector{Dict{Tuple{Int,Int},Matrix{C}}},
+  block_sparse_matrices::Vector{BlockSparseMatrix{C}},
   llink::ITensors.QNIndex,
   rlink::ITensors.QNIndex,
   site::ITensors.QNIndex;
@@ -136,7 +138,7 @@ function _to_ITensor_splitblocks(
 
   num_nonzero_entries = sum(
     count(value -> !iszero(value), block) for matrix in block_sparse_matrices for
-    block in values(matrix)
+    column in matrix for block in values(column)
   )
 
   blocks = Vector{Block{4}}(undef, num_nonzero_entries)
@@ -157,9 +159,9 @@ end
 Assemble a block-sparse ITensor MPO tensor from component-local block matrices.
 
 `block_sparse_matrices[i]` stores the local MPO matrix blocks for component `i`
-as `(left_link, right_link) => local_operator_matrix`. The corresponding
-`offsets[i]` shifts those component-local right-link ids into the global
-outgoing bond space. The returned tensor has indices
+as `right_link => (left_link => local_operator_matrix)`. The corresponding
+`offsets[i]` shifts those component-local right-link ids into the global outgoing
+bond space. The returned tensor has indices
 `(dag(llink), rlink, prime(site), dag(site))`.
 
 When `splitblocks=true`, construction delegates to `_to_ITensor_splitblocks`.
@@ -169,7 +171,7 @@ writes nonzero entries into the appropriate block views.
 """
 function to_ITensor(
   offsets::Vector{Int},
-  block_sparse_matrices::Vector{Dict{Tuple{Int,Int},Matrix{C}}},
+  block_sparse_matrices::Vector{BlockSparseMatrix{C}},
   llink::ITensors.QNIndex,
   rlink::ITensors.QNIndex,
   site::ITensors.QNIndex;
@@ -191,30 +193,33 @@ function to_ITensor(
   num_site_blocks = nblocks(inds[4])
 
   block_ids = Set{Int}()
-  sizehint!(block_ids, sum(length, block_sparse_matrices))
+  sizehint!(block_ids, sum(matrix -> sum(length, matrix), block_sparse_matrices))
 
   for (offset, matrix) in zip(offsets, block_sparse_matrices)
-    for ((left_link, right_link), block) in matrix
+    for (right_link, column) in enumerate(matrix)
       shifted_right_link = right_link + offset
-      left_block = left_blocks[left_link]
       right_block = right_blocks[shifted_right_link]
 
-      @inbounds for i in axes(block, 1), j in axes(block, 2)
-        value = block[i, j]
-        iszero(value) && continue
-
-        push!(
-          block_ids,
-          _blockid(
-            left_block,
-            right_block,
-            site_prime_blocks[i],
-            site_blocks[j],
-            num_right_blocks,
-            num_site_prime_blocks,
-            num_site_blocks,
-          ),
-        )
+      for (left_link, block) in column
+        left_block = left_blocks[left_link]
+
+        @inbounds for site_prime_link in axes(block, 1), site_link in axes(block, 2)
+          value = block[site_prime_link, site_link]
+          iszero(value) && continue
+
+          push!(
+            block_ids,
+            _blockid(
+              left_block,
+              right_block,
+              site_prime_blocks[site_prime_link],
+              site_blocks[site_link],
+              num_right_blocks,
+              num_site_prime_blocks,
+              num_site_blocks,
+            ),
+          )
+        end
       end
     end
   end
@@ -245,29 +250,35 @@ function to_ITensor(
     offset = offsets[i]
     matrix = block_sparse_matrices[i]
 
-    for ((left_link, right_link), block) in matrix
+    for (right_link, column) in enumerate(matrix)
       shifted_right_link = right_link + offset
-      left_block = left_blocks[left_link]
       right_block = right_blocks[shifted_right_link]
-      left_offset = left_offsets[left_link]
       right_offset = right_offsets[shifted_right_link]
 
-      @inbounds for i in axes(block, 1), j in axes(block, 2)
-        value = block[i, j]
-        iszero(value) && continue
-
-        block_id = _blockid(
-          left_block,
-          right_block,
-          site_prime_blocks[i],
-          site_blocks[j],
-          num_right_blocks,
-          num_site_prime_blocks,
-          num_site_blocks,
-        )
-        block_views[block_id][
-          left_offset, right_offset, site_prime_offsets[i], site_offsets[j]
-        ] = value
+      for (left_link, block) in column
+        left_block = left_blocks[left_link]
+        left_offset = left_offsets[left_link]
+
+        @inbounds for site_prime_link in axes(block, 1), site_link in axes(block, 2)
+          value = block[site_prime_link, site_link]
+          iszero(value) && continue
+
+          block_id = _blockid(
+            left_block,
+            right_block,
+            site_prime_blocks[site_prime_link],
+            site_blocks[site_link],
+            num_right_blocks,
+            num_site_prime_blocks,
+            num_site_blocks,
+          )
+          block_views[block_id][
+            left_offset,
+            right_offset,
+            site_prime_offsets[site_prime_link],
+            site_offsets[site_link],
+          ] = value
+        end
       end
     end
   end
@@ -299,8 +310,10 @@ function to_ITensor(
   Threads.@threads for i in eachindex(offsets)
     offset = offsets[i]
     matrix = block_sparse_matrices[i]
-    for ((left_link, right_link), block) in matrix
-      tensor[left_link, right_link + offset, :, :] = block
+    for (right_link, column) in enumerate(matrix)
+      for (left_link, block) in column
+        tensor[left_link, right_link + offset, :, :] = block
+      end
     end
   end
 
diff --git a/test/tensor-construction-tests.jl b/test/tensor-construction-tests.jl
index b401610..fd6e42b 100644
--- a/test/tensor-construction-tests.jl
+++ b/test/tensor-construction-tests.jl
@@ -4,14 +4,14 @@ using Test
 
 function reference_to_ITensor(
   offsets::Vector{Int},
-  block_sparse_matrices,
+  block_sparse_matrices::Vector{ITensorMPOConstruction.BlockSparseMatrix{C}},
   llink::ITensors.QNIndex,
   rlink::ITensors.QNIndex,
   site::ITensors.QNIndex;
   splitblocks=false,
   tol=0.0,
   checkflux=true,
-)::ITensor
+)::ITensor where {C}
   if splitblocks
     llink = ITensors.splitblocks(llink)
     rlink = ITensors.splitblocks(rlink)
@@ -19,17 +19,19 @@ function reference_to_ITensor(
   end
 
   T = ITensors.BlockSparseTensor(
-    eltype(first(first(values(first(block_sparse_matrices))))),
+    C,
     Block{4}[],
     (dag(llink), rlink, prime(site), dag(site)),
   )
 
   for (offset, matrix) in zip(offsets, block_sparse_matrices)
-    for ((left_link, right_link), block) in matrix
-      for i in axes(block, 1)
-        for j in axes(block, 2)
-          if abs(block[i, j]) > tol
-            T[left_link, right_link + offset, i, j] = block[i, j]
+    for (right_link, column) in enumerate(matrix)
+      for (left_link, block) in column
+        for i in axes(block, 1)
+          for j in axes(block, 2)
+            if abs(block[i, j]) > tol
+              T[left_link, right_link + offset, i, j] = block[i, j]
+            end
           end
         end
       end
@@ -47,14 +49,14 @@ function test_sparse_to_ITensor_matches_reference(splitblocks::Bool=false)
 
   offsets = [0, 1]
   block_sparse_matrices = [
-    Dict(
-      (1, 1) => [1.0 -2.0 0.0; 3.0 4.0 0.0; 0.0 0.0 5.0],
-      (3, 2) => [0.0 0.0 0.0; -6.0 0.0 0.0; 0.0 0.0 -1.5],
-    ),
-    Dict(
-      (2, 1) => [0.0 0.0 7.0; 0.0 0.0 -8.0; 0.0 0.0 0.0],
-      (3, 2) => [0.25 0.5 0.0; -0.75 1.0 0.0; 0.0 0.0 -2.0],
-    ),
+    [
+      Dict(1 => [1.0 -2.0 0.0; 3.0 4.0 0.0; 0.0 0.0 5.0]),
+      Dict(3 => [0.0 0.0 0.0; -6.0 0.0 0.0; 0.0 0.0 -1.5]),
+    ],
+    [
+      Dict(2 => [0.0 0.0 7.0; 0.0 0.0 -8.0; 0.0 0.0 0.0]),
+      Dict(3 => [0.25 0.5 0.0; -0.75 1.0 0.0; 0.0 0.0 -2.0]),
+    ],
   ]
 
   expected = reference_to_ITensor(
@@ -74,7 +76,7 @@ function test_sparse_to_ITensor_all_zero_blocks(splitblocks::Bool=false)
   site = Index([QN("N", 0) => 2, QN("N", 1) => 1]; tags="Site", dir=ITensors.Out)
 
   offsets = [0]
-  block_sparse_matrices = [Dict((1, 1) => zeros(ComplexF64, dim(site), dim(site)))]
+  block_sparse_matrices = [[Dict(1 => zeros(ComplexF64, dim(site), dim(site)))]]
 
   expected = reference_to_ITensor(
     offsets, block_sparse_matrices, llink, rlink, site; splitblocks

From 40e76e5ba96d44d981caa886393dea1b4454f082 Mon Sep 17 00:00:00 2001
From: Ben Corbett <corbett5@llnl.gov>
Date: Sun, 26 Apr 2026 14:15:19 -0600
Subject: [PATCH 2/4] Use Dictionaries.jl

---
 Project.toml                        |  2 ++
 examples/electronic-structure.jl    | 24 +++++++++++++-----------
 examples/fermi-hubbard-ks.jl        | 16 ++++++++--------
 examples/fermi-hubbard-tc.jl        | 15 ++++++++++++---
 src/ITensorMPOConstruction.jl       |  1 +
 src/large-graph-mpo-common.jl       |  6 +++---
 src/large-graph-mpo-qr.jl           |  2 +-
 src/large-graph-mpo-vertex-cover.jl |  4 ++--
 src/sparse_tensor_construction.jl   |  8 ++++----
 test/tensor-construction-tests.jl   | 13 +++++++------
 10 files changed, 53 insertions(+), 38 deletions(-)

diff --git a/Project.toml b/Project.toml
index 4ea64e2..690c757 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,6 +4,7 @@ authors = ["Ben Corbett and contributors"]
 version = "0.2.1"
 
 [deps]
+Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4"
 FunctionWrappers = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e"
 ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2"
 ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"
@@ -15,6 +16,7 @@ ThreadsX = "ac1d9e8a-700a-412c-b207-f0111f4b6c0d"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 
 [compat]
+Dictionaries = "0.4.6"
 FunctionWrappers = "1.1.3"
 ITensorMPS = "0.1, 0.2, 0.3"
 ITensors = "0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9"
diff --git a/examples/electronic-structure.jl b/examples/electronic-structure.jl
index f7621de..876af39 100644
--- a/examples/electronic-structure.jl
+++ b/examples/electronic-structure.jl
@@ -151,11 +151,13 @@ for alg in ("VC", "QR")
     ) # TODO: remove splitblocks=true after warning
     N > 5 && print_timer()
 
+    percent_sparse = round(100 * sparsity(H); digits=2)
     println(
-      "The maximum bond dimension is $(maxlinkdim(H)), sparsity = $(ITensorMPOConstruction.sparsity(H))",
+      "The maximum bond dimension is $(maxlinkdim(H)), sparsity = $percent_sparse%",
     )
     @assert maxlinkdim(H) == 2 * N^2 + 3 * N + 2
 
+    GC.gc(true)
     println()
   end
 end
@@ -167,16 +169,16 @@ nothing;
 #
 # | ``N`` | Minimum Vertex Cover    | QR Decomposition      |
 # |-------|-------------------------|-----------------------|
-# | 10    | 0.02s / 98.23%          | 0.13s / 95.69%        |
-# | 20    | 0.90s / 98.84%          | 1.20s / 97.01%        |
-# | 30    | 1.04s / 99.14%          | 8.37s / 97.54%        |
-# | 40    | 3.30s / 99.32%          | 16.3s / 97.72%        |
-# | 50    | 9.13s / 99.44%          | 56.1s / 97.83%        |
-# | 60    | 21.1s / 99.52%          | 169s / 97.90%         |
-# | 70    | 50.7s / 99.58%          | 516s / 97.95%         |
-# | 80    | 106s / 99.63%           | N/A                   |
-# | 90    | 221s / 99.67%           | N/A                   |
-# | 100   | 744s / 99.70%           | N/A                   |
+# | 10    | 0.01s / 98.23%          | 0.04s / 95.69%        |
+# | 20    | 0.81s / 98.84%          | 0.54s / 97.01%        |
+# | 30    | 0.94s / 99.14%          | 3.49s / 97.54%        |
+# | 40    | 2.88s / 99.32%          | 14.4s / 97.73%        |
+# | 50    | 7.42s / 99.44%          | 50.8s / 97.83%        |
+# | 60    | 17.2s / 99.52%          | 159s / 97.90%         |
+# | 70    | 40.2s / 99.58%          | 452s / 97.95%         |
+# | 80    | 95.0s / 99.63%          | N/A                   |
+# | 90    | 207s / 99.67%           | N/A                   |
+# | 100   | 739s / 99.70%           | N/A                   |
 #
 # Not only does the vertex cover algorithm produce an MPO much faster than the QR algorithm, but the resulting MPO has almost five times fewer entries (note: if `splitblocks=false` then the sparsities of the MPOs are equal). This story remains mostly unchanged when we construct the Hamiltonian for real systems.
 
diff --git a/examples/fermi-hubbard-ks.jl b/examples/fermi-hubbard-ks.jl
index 6cc15fa..a5a74e6 100644
--- a/examples/fermi-hubbard-ks.jl
+++ b/examples/fermi-hubbard-ks.jl
@@ -83,14 +83,14 @@ end
 # | ``N`` | ITensorMPS    | ITensorMPOConstruction: QR  |
 # |-------|---------------|-----------------------------|
 # | 10    | 0.32s / 92.7% | 0.01s / 99.32%              |
-# | 20    | 30.6s / 92.6% | 0.06s / 99.70%              |
-# | 30    | 792s / 92.6%  | 0.53s / 99.81%              |
-# | 40    | N/A           | 0.55s / 99.86%              |
-# | 50    | N/A           | 0.38s / 99.89%              |
-# | 100   | N/A           | 2.81s / 99.94%              |
-# | 200   | N/A           | 30.3s / 99.97%              |
+# | 20    | 30.6s / 92.6% | 0.04s / 99.70%              |
+# | 30    | 792s / 92.6%  | 0.45s / 99.81%              |
+# | 40    | N/A           | 0.48s / 99.86%              |
+# | 50    | N/A           | 0.29s / 99.89%              |
+# | 100   | N/A           | 2.46s / 99.94%              |
+# | 200   | N/A           | 39.1s / 99.97%              |
 # | 300   | N/A           | 136s / 99.982%              |
-# | 400   | N/A           | 415s / 99.986%              |
-# | 500   | N/A           | 1274s / 99.989%             |
+# | 400   | N/A           | 438s / 99.986%              |
+# | 500   | N/A           | 1313s / 99.989%             |
 #
 # Additionally, we constructed the MPO using the vertex cover algorithm, which turns out to be particularly ill-suited to this problem. Unlike the two other algorithms, the vertex cover algorithm produces MPOs of bond dimension ``1.5 N^2 + 4 N + 2``, far from the optimal of ``10 N - 4``. Granted, with `splitblocks=true` and `N = 100`, the MPO has a sparsity of 99.994%, but with a bond dimension of 15402 it still contains 25 times more nonzero entries than the MPO produced with the QR decomposition.
diff --git a/examples/fermi-hubbard-tc.jl b/examples/fermi-hubbard-tc.jl
index 4191211..91d9f0e 100644
--- a/examples/fermi-hubbard-tc.jl
+++ b/examples/fermi-hubbard-tc.jl
@@ -389,7 +389,7 @@ end
 using TimerOutputs
 for grid_size in ((2, 2), (6, 6))
   let t = 1, U = 4, J = -0.5, mapping = bipartite_mapping(grid_size)
-    @time "Constructing OpIDSum" sites, os = transcorrelated_fermi_hubbard(t, U, J, mapping)
+    @time "Constructing OpIDSum" sites, os = transcorrelated_fermi_hubbard(t, U, J, mapping; conserve_momentum=true)
     reset_timer!()
     @time "Constructing MPO" H = MPO_new(
       os,
@@ -400,10 +400,19 @@ for grid_size in ((2, 2), (6, 6))
       check_for_errors=false,
     )
 
+    prod(grid_size) > 4 && print_timer()
+
+    percent_sparse = round(100 * sparsity(H); digits=3)
     println(
-      "Constructed the MPO of bond dimension $(maxlinkdim(H)) and sparsity $(sparsity(H))"
+      "The maximum bond dimension is $(maxlinkdim(H)), sparsity = $percent_sparse%",
     )
-    grid_size != (2, 2) && print_timer()
+
+    if iszero(J) && length(grid_size) == 1 && mapping == standard_mapping(grid_size)
+      @assert maxlinkdim(H) == 10 * N - 4
+    end
+
+    GC.gc(true)
+    println()
   end
 end
 
diff --git a/src/ITensorMPOConstruction.jl b/src/ITensorMPOConstruction.jl
index 9f4d46f..c81a4f9 100644
--- a/src/ITensorMPOConstruction.jl
+++ b/src/ITensorMPOConstruction.jl
@@ -5,6 +5,7 @@ module ITensorMPOConstruction
 #
 using ITensorMPS
 using ITensors
+using Dictionaries
 using LinearAlgebra
 using SparseArrays
 using TimerOutputs
diff --git a/src/large-graph-mpo-common.jl b/src/large-graph-mpo-common.jl
index 09c8756..91d3b8e 100644
--- a/src/large-graph-mpo-common.jl
+++ b/src/large-graph-mpo-common.jl
@@ -5,11 +5,11 @@ Vector-backed block-sparse matrix representation used for intermediate MPO
 tensor storage.
 
 The outer vector is indexed by component-local `right_link`. Each inner
-dictionary maps `left_link` to the dense local operator matrix for that block.
+`Dictionary` maps `left_link` to the dense local operator matrix for that block.
 `at_site!` later returns offsets that place component-local right-link ids into
 the full outgoing MPO bond.
 """
-BlockSparseMatrix{C} = Vector{Dict{Int,Matrix{C}}}
+BlockSparseMatrix{C} = Vector{Dictionary{Int,Matrix{C}}}
 
 """
     MPOGraph{N,C,Ti}
@@ -362,7 +362,7 @@ function process_single_left_vertex_cc!(
   lv = left_vertex(g, lv_id)
   local_op = op_cache_vec[n][lv.op_id].matrix
 
-  matrix = [Dict{Int,Matrix{ValType}}() for _ in 1:rank]
+  matrix = [Dictionary{Int,Matrix{ValType}}() for _ in 1:rank]
   matrix_of_cc[cc] = matrix
 
   matrix_element = get!(matrix[rank], lv.link) do
diff --git a/src/large-graph-mpo-qr.jl b/src/large-graph-mpo-qr.jl
index 38bffe0..6442245 100644
--- a/src/large-graph-mpo-qr.jl
+++ b/src/large-graph-mpo-qr.jl
@@ -123,7 +123,7 @@ blocks.
 
     rank_of_cc[cc] = rank
 
-    matrix = [Dict{Int,Matrix{ValType}}() for _ in 1:rank]
+    matrix = [Dictionary{Int,Matrix{ValType}}() for _ in 1:rank]
     matrix_of_cc[cc] = matrix
 
     ## Form the local transformation tensor.
diff --git a/src/large-graph-mpo-vertex-cover.jl b/src/large-graph-mpo-vertex-cover.jl
index 69fff9d..3a548af 100644
--- a/src/large-graph-mpo-vertex-cover.jl
+++ b/src/large-graph-mpo-vertex-cover.jl
@@ -44,7 +44,7 @@ otherwise this also builds `next_edges_of_cc` for the graph at site `n + 1`.
     rank = length(left_cover) + length(right_cover)
     rank_of_cc[cc] = rank
 
-    matrix = [Dict{Int,Matrix{ValType}}() for _ in 1:rank]
+    matrix = [Dictionary{Int,Matrix{ValType}}() for _ in 1:rank]
     matrix_of_cc[cc] = matrix
 
     ## Construct the tensor from the left cover.
@@ -55,7 +55,7 @@ otherwise this also builds `next_edges_of_cc` for the graph at site `n + 1`.
 
       matrix_element = zeros(ValType, site_dim, site_dim)
       add_to_local_matrix!(matrix_element, one(ValType), local_op, lv.needs_JW_string)
-      matrix[m][lv.link] = matrix_element
+      set!(matrix[m], lv.link, matrix_element)
     end
 
     ## Construct the tensor from the right cover.
diff --git a/src/sparse_tensor_construction.jl b/src/sparse_tensor_construction.jl
index 344cff4..566ffae 100644
--- a/src/sparse_tensor_construction.jl
+++ b/src/sparse_tensor_construction.jl
@@ -96,7 +96,7 @@ function _fill_splitblocks!(
     for (right_link, column) in enumerate(matrix)
       shifted_right_link = right_link + offset
 
-      for (left_link, block) in column
+      for (left_link, block) in pairs(column)
         @inbounds for site_prime_link in axes(block, 1), site_link in axes(block, 2)
           value = block[site_prime_link, site_link]
           iszero(value) && continue
@@ -200,7 +200,7 @@ function to_ITensor(
       shifted_right_link = right_link + offset
       right_block = right_blocks[shifted_right_link]
 
-      for (left_link, block) in column
+      for (left_link, block) in pairs(column)
         left_block = left_blocks[left_link]
 
         @inbounds for site_prime_link in axes(block, 1), site_link in axes(block, 2)
@@ -255,7 +255,7 @@ function to_ITensor(
       right_block = right_blocks[shifted_right_link]
       right_offset = right_offsets[shifted_right_link]
 
-      for (left_link, block) in column
+      for (left_link, block) in pairs(column)
         left_block = left_blocks[left_link]
         left_offset = left_offsets[left_link]
 
@@ -311,7 +311,7 @@ function to_ITensor(
     offset = offsets[i]
     matrix = block_sparse_matrices[i]
     for (right_link, column) in enumerate(matrix)
-      for (left_link, block) in column
+      for (left_link, block) in pairs(column)
         tensor[left_link, right_link + offset, :, :] = block
       end
     end
diff --git a/test/tensor-construction-tests.jl b/test/tensor-construction-tests.jl
index fd6e42b..f9d6d5a 100644
--- a/test/tensor-construction-tests.jl
+++ b/test/tensor-construction-tests.jl
@@ -1,4 +1,5 @@
 using ITensorMPOConstruction
+using Dictionaries
 using ITensors
 using Test
 
@@ -26,7 +27,7 @@ function reference_to_ITensor(
 
   for (offset, matrix) in zip(offsets, block_sparse_matrices)
     for (right_link, column) in enumerate(matrix)
-      for (left_link, block) in column
+      for (left_link, block) in pairs(column)
         for i in axes(block, 1)
           for j in axes(block, 2)
             if abs(block[i, j]) > tol
@@ -50,12 +51,12 @@ function test_sparse_to_ITensor_matches_reference(splitblocks::Bool=false)
   offsets = [0, 1]
   block_sparse_matrices = [
     [
-      Dict(1 => [1.0 -2.0 0.0; 3.0 4.0 0.0; 0.0 0.0 5.0]),
-      Dict(3 => [0.0 0.0 0.0; -6.0 0.0 0.0; 0.0 0.0 -1.5]),
+      Dictionary([1], [[1.0 -2.0 0.0; 3.0 4.0 0.0; 0.0 0.0 5.0]]),
+      Dictionary([3], [[0.0 0.0 0.0; -6.0 0.0 0.0; 0.0 0.0 -1.5]]),
     ],
     [
-      Dict(2 => [0.0 0.0 7.0; 0.0 0.0 -8.0; 0.0 0.0 0.0]),
-      Dict(3 => [0.25 0.5 0.0; -0.75 1.0 0.0; 0.0 0.0 -2.0]),
+      Dictionary([2], [[0.0 0.0 7.0; 0.0 0.0 -8.0; 0.0 0.0 0.0]]),
+      Dictionary([3], [[0.25 0.5 0.0; -0.75 1.0 0.0; 0.0 0.0 -2.0]]),
     ],
   ]
 
@@ -76,7 +77,7 @@ function test_sparse_to_ITensor_all_zero_blocks(splitblocks::Bool=false)
   site = Index([QN("N", 0) => 2, QN("N", 1) => 1]; tags="Site", dir=ITensors.Out)
 
   offsets = [0]
-  block_sparse_matrices = [[Dict(1 => zeros(ComplexF64, dim(site), dim(site)))]]
+  block_sparse_matrices = [[Dictionary([1], [zeros(ComplexF64, dim(site), dim(site))])]]
 
   expected = reference_to_ITensor(
     offsets, block_sparse_matrices, llink, rlink, site; splitblocks

From 2b276a936ca342fffdbe41b6ee8d5e9e651b95e7 Mon Sep 17 00:00:00 2001
From: Ben Corbett <corbett5@llnl.gov>
Date: Sun, 26 Apr 2026 14:20:16 -0600
Subject: [PATCH 3/4] Fixed test dependencies.

---
 test/Project.toml | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/test/Project.toml b/test/Project.toml
index ec6a311..09c7811 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,4 +1,5 @@
 [deps]
+Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2"
 ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"
@@ -7,6 +8,7 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
+Dictionaries = "0.4.6"
 Graphs = "1.14.0"
 ITensorMPS = "0.1, 0.2, 0.3"
 ITensors = "0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9"

From 8e6b270eaf80ed497c3181c819a9e2a6e58f28ac Mon Sep 17 00:00:00 2001
From: Ben Corbett <corbett5@llnl.gov>
Date: Sun, 26 Apr 2026 14:32:43 -0600
Subject: [PATCH 4/4] Move qr specifics out of common.

---
 src/large-graph-mpo-common.jl | 68 -----------------------------------
 src/large-graph-mpo-qr.jl     | 68 +++++++++++++++++++++++++++++++++++
 2 files changed, 68 insertions(+), 68 deletions(-)

diff --git a/src/large-graph-mpo-common.jl b/src/large-graph-mpo-common.jl
index 91d3b8e..839c79c 100644
--- a/src/large-graph-mpo-common.jl
+++ b/src/large-graph-mpo-common.jl
@@ -331,74 +331,6 @@ function merge_qn_sectors(
   return new_order, new_qi
 end
 
-# TODO: Move this into large-graph-mpo-qr.jl
-"""
-    process_single_left_vertex_cc!(
-      matrix_of_cc, rank_of_cc, next_edges_of_cc, g, ccs, cc, n, sites, op_cache_vec
-    ) -> Nothing
-
-Handle the common connected-component case with exactly one left vertex.
-
-This fills the local MPO tensor contribution for the component, records rank
-`1`, and either applies the terminal scaling on the last site or builds the
-outgoing edges for the next site. For nonterminal sites, the consumed left
-adjacency list is cleared after the outgoing edges are built.
-"""
-function process_single_left_vertex_cc!(
-  matrix_of_cc::Vector{BlockSparseMatrix{ValType}},
-  rank_of_cc::Vector{Int},
-  next_edges_of_cc::Vector{Matrix{Tuple{Vector{Int},Vector{C}}}},
-  g::MPOGraph{N,C,Ti},
-  ccs::BipartiteGraphConnectedComponents,
-  cc::Int,
-  n::Int,
-  sites::Vector{<:Index},
-  op_cache_vec::OpCacheVec,
-)::Nothing where {ValType<:Number,N,C,Ti}
-  lv_id = only(ccs.lvs_of_component[cc])
-  rank = 1
-  rank_of_cc[cc] = rank
-
-  lv = left_vertex(g, lv_id)
-  local_op = op_cache_vec[n][lv.op_id].matrix
-
-  matrix = [Dictionary{Int,Matrix{ValType}}() for _ in 1:rank]
-  matrix_of_cc[cc] = matrix
-
-  matrix_element = get!(matrix[rank], lv.link) do
-    return zeros(ValType, dim(sites[n]), dim(sites[n]))
-  end
-
-  add_to_local_matrix!(matrix_element, one(ValType), local_op, lv.needs_JW_string)
-
-  if n == length(sites)
-    scaling = only(g.edge_weights_from_left[lv_id])
-
-    for column in matrix, block in values(column)
-      block .*= scaling
-    end
-
-    return nothing
-  end
-
-  next_edges = Matrix{Tuple{Vector{Int},Vector{C}}}(
-    undef, rank, length(op_cache_vec[n + 1])
-  )
-  for i in eachindex(next_edges)
-    next_edges[i] = (Int[], C[])
-  end
-
-  build_next_edges_specialization!(
-    next_edges, g, n, g.right_vertex_ids_from_left[lv_id], g.edge_weights_from_left[lv_id]
-  )
-
-  clear_edges_from_left!(g, lv_id)
-
-  next_edges_of_cc[cc] = next_edges
-
-  return nothing
-end
-
 """
     at_site!(ValType, g, n, sites, tol, absolute_tol, op_cache_vec, alg;
              combine_qn_sectors, output_level=0)
diff --git a/src/large-graph-mpo-qr.jl b/src/large-graph-mpo-qr.jl
index 6442245..db2e1ac 100644
--- a/src/large-graph-mpo-qr.jl
+++ b/src/large-graph-mpo-qr.jl
@@ -80,6 +80,74 @@ function for_non_zeros_batch(
   end
 end
 
+"""
+    process_single_left_vertex_cc!(
+      matrix_of_cc, rank_of_cc, next_edges_of_cc, g, ccs, cc, n, sites, op_cache_vec
+    ) -> Nothing
+
+Handle the common connected-component case with exactly one left vertex.
+
+This fills the local MPO tensor contribution for the component, records rank
+`1`, and either applies the terminal scaling on the last site or builds the
+outgoing edges for the next site. For nonterminal sites, the consumed left
+adjacency list is cleared after the outgoing edges are built.
+"""
+function process_single_left_vertex_cc!(
+  matrix_of_cc::Vector{BlockSparseMatrix{ValType}},
+  rank_of_cc::Vector{Int},
+  next_edges_of_cc::Vector{Matrix{Tuple{Vector{Int},Vector{C}}}},
+  g::MPOGraph{N,C,Ti},
+  ccs::BipartiteGraphConnectedComponents,
+  cc::Int,
+  n::Int,
+  sites::Vector{<:Index},
+  op_cache_vec::OpCacheVec,
+)::Nothing where {ValType<:Number,N,C,Ti}
+  lv_id = only(ccs.lvs_of_component[cc])
+  rank = 1
+  rank_of_cc[cc] = rank
+
+  lv = left_vertex(g, lv_id)
+  local_op = op_cache_vec[n][lv.op_id].matrix
+
+  matrix = [Dictionary{Int,Matrix{ValType}}() for _ in 1:rank]
+  matrix_of_cc[cc] = matrix
+
+  matrix_element = get!(matrix[rank], lv.link) do
+    return zeros(ValType, dim(sites[n]), dim(sites[n]))
+  end
+
+  add_to_local_matrix!(matrix_element, one(ValType), local_op, lv.needs_JW_string)
+
+  if n == length(sites)
+    scaling = only(g.edge_weights_from_left[lv_id])
+
+    for column in matrix, block in values(column)
+      block .*= scaling
+    end
+
+    return nothing
+  end
+
+  next_edges = Matrix{Tuple{Vector{Int},Vector{C}}}(
+    undef, rank, length(op_cache_vec[n + 1])
+  )
+  for i in eachindex(next_edges)
+    next_edges[i] = (Int[], C[])
+  end
+
+  build_next_edges_specialization!(
+    next_edges, g, n, g.right_vertex_ids_from_left[lv_id], g.edge_weights_from_left[lv_id]
+  )
+
+  clear_edges_from_left!(g, lv_id)
+
+  next_edges_of_cc[cc] = next_edges
+
+  return nothing
+end
+
+
 """
     process_qr(
       matrix_of_cc, rank_of_cc, next_edges_of_cc, g, ccs, n, sites, tol, absolute_tol, op_cache_vec