diff --git a/docs/make.jl b/docs/make.jl
index e82688b..bc591bf 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -17,6 +17,7 @@ examples = [
   "./examples/fermi-hubbard-ks.jl",
   "./examples/haldane-shastry.jl",
   "./examples/electronic-structure.jl",
+  "./examples/symbolic-transverse-field-ising.jl",
   "./examples/fermi-hubbard-tc.jl",
 ]
 
@@ -53,6 +54,7 @@ makedocs(;
       "Real Space Fermi-Hubbard" => "examples/fermi-hubbard-rs.md",
       "Momentum Space Fermi-Hubbard" => "examples/fermi-hubbard-ks.md",
       "Electronic Structure" => "examples/electronic-structure.md",
+      "Symbolic Transverse Field Ising" => "examples/symbolic-transverse-field-ising.md",
       "Haldane-Shastry and Truncation" => "examples/haldane-shastry.md",
       "Challenge Problem" => "examples/fermi-hubbard-tc.md",
     ],
diff --git a/docs/src/documentation/MPO_new.md b/docs/src/documentation/MPO_new.md
index 49ddafc..509c3a6 100644
--- a/docs/src/documentation/MPO_new.md
+++ b/docs/src/documentation/MPO_new.md
@@ -9,11 +9,20 @@ The main exported function is `MPO_new`, which takes an `OpSum` and transforms i
 
 ```@docs
 MPO_new
+SymbolicMPO
+MPO_symbolic
 resume_MPO_construction!
 instantiate_MPO
+instantiate_MPO!
 sparsity
 block2_nnz
 ```
 
+## Symbolic MPO construction
+`MPO_symbolic` constructs the vertex-cover MPO structure once from an integer-labeled
+`OpIDSum`. The labels are positive one-based coefficient ids, and
+`instantiate_MPO` or `instantiate_MPO!` substitutes numeric values later without
+rebuilding the graph.
+
 ## A note on truncation
 `MPO_new` is designed to construct _numerically exact_ MPOs; the tolerance parameter should only be used to adjust the definition of "numerically exact", not to perform truncation. If truncation is desired, truncate the resulting MPO with `ITensorMPS.truncate!`. See [Haldane-Shastry and truncation](../examples/haldane-shastry.md) for an example.
diff --git a/docs/src/documentation/unexported.md b/docs/src/documentation/unexported.md
index 9f045f3..7c9c2a4 100644
--- a/docs/src/documentation/unexported.md
+++ b/docs/src/documentation/unexported.md
@@ -13,6 +13,9 @@ ITensorMPOConstruction.rewrite_in_operator_basis!
 ITensorMPOConstruction.op_sum_to_opID_sum
 ITensorMPOConstruction.check_os_for_errors
 ITensorMPOConstruction.prepare_opID_sum!
+ITensorMPOConstruction.SymbolicLocalMatrix
+ITensorMPOConstruction.SymbolicBlockSparseMatrix
+ITensorMPOConstruction.internalize_symbolic_ids!
 ITensorMPOConstruction.get_onsite_op
 ```
 
diff --git a/examples/electronic-structure.jl b/examples/electronic-structure.jl
index 876af39..6e7c1d9 100644
--- a/examples/electronic-structure.jl
+++ b/examples/electronic-structure.jl
@@ -55,12 +55,7 @@ function get_coefficients(N::Int)::Tuple{Array{Float64,2},Array{Float64,4}}
   return h, V
 end
 
-function electronic_structure_OpIDSum(
-  N::Int, h::Array{Float64,2}, V::Array{Float64,4}
-)::Tuple{Vector{<:Index},OpIDSum}
-  ↓ = false
-  ↑ = true
-
+function electronic_structure_sites_and_op_cache(N::Int)
   sites = siteinds("Electron", N; conserve_qns=true)
 
   operatorNames = [
@@ -85,6 +80,17 @@ function electronic_structure_OpIDSum(
   ]
 
   op_cache_vec = to_OpCacheVec(sites, operatorNames)
+  return sites, op_cache_vec
+end
+
+function electronic_structure_OpIDSum(
+  N::Int,
+  h::Array{Float64,2},
+  V::Array{Float64,4},
+  op_cache_vec::OpCacheVec,
+)::OpIDSum
+  ↓ = false
+  ↑ = true
 
   opC(k::Int, spin::Bool) = OpID{UInt8}(2 + spin, k)
   opCdag(k::Int, spin::Bool) = OpID{UInt8}(4 + spin, k)
@@ -129,40 +135,45 @@ function electronic_structure_OpIDSum(
     end
   end
 
-  return sites, os
+  return os
 end
 
-for alg in ("VC", "QR")
-  for N in [10, 20]
-    println("Constructing the electronic structure MPO for $N sites using $alg")
-
-    reset_timer!()
-    @time "Constructing OpIDSum" sites, os = electronic_structure_OpIDSum(
-      N, get_coefficients(N)...
-    )
-    @time "Constructing MPO" H = MPO_new(
-      os,
-      sites;
-      alg,
-      basis_op_cache_vec=os.op_cache_vec,
-      splitblocks=true,
-      check_for_errors=false,
-      checkflux=false,
-    ) # 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 = $percent_sparse%",
-    )
-    @assert maxlinkdim(H) == 2 * N^2 + 3 * N + 2
-
-    GC.gc(true)
-    println()
-  end
+function electronic_structure_OpIDSum(
+  N::Int, h::Array{Float64,2}, V::Array{Float64,4}
+)::Tuple{Vector{<:Index},OpIDSum}
+  sites, op_cache_vec = electronic_structure_sites_and_op_cache(N)
+  return sites, electronic_structure_OpIDSum(N, h, V, op_cache_vec)
 end
 
-nothing;
+# for alg in ("VC", "QR")
+#   for N in [10, 20]
+#     println("Constructing the electronic structure MPO for $N sites using $alg")
+
+#     reset_timer!()
+#     @time "Constructing OpIDSum" sites, os = electronic_structure_OpIDSum(
+#       N, get_coefficients(N)...
+#     )
+#     @time "Constructing MPO" H = MPO_new(
+#       os,
+#       sites;
+#       alg,
+#       basis_op_cache_vec=os.op_cache_vec,
+#       splitblocks=true,
+#       check_for_errors=false,
+#       checkflux=false,
+#     ) # 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 = $percent_sparse%",
+#     )
+#     @assert maxlinkdim(H) == 2 * N^2 + 3 * N + 2
+
+#     GC.gc(true)
+#     println()
+#   end
+# end
 
 # ## Results
 # Below are the runtime and sparsities for the QR decomposition algorithm and the vertex cover algorithm. Both algorithms produce MPOs of bond dimension ``2 N^2 + 3 N + 2``, which is optimal for a generic set of coefficients. These timings were taken with `julia -t8 --gcthreads=8,1` on a 2021 MacBook Pro with the M1 Max CPU and 32GB of memory.
@@ -192,3 +203,160 @@ nothing;
 # | ``\left[\text{Fe}_2 \text{S} (\text{C H}_3)(\text{SCH}_3)_4 \right]^{3-}`` | 36    | `VC`      | 2702           | 99.26%   |
 #
 # ### Thanks to [Huanchen Zhai](https://scholar.google.com/citations?user=HM_YBL0AAAAJ&hl=en) for providing the discussion and data motivating this section.
+
+# # Symbolic construction
+
+function electronic_structure_symbolic_OpIDSum(
+  N::Int,
+  op_cache_vec::OpCacheVec,
+)::Tuple{OpIDSum, Vector{Tuple{Int, Int}}, Vector{Tuple{Int, Int, Int, Int, Bool, Bool, Bool, Bool}}}
+  ↓ = false
+  ↑ = true
+
+  opC(k::Int, spin::Bool) = OpID(2 + spin, k)
+  opCdag(k::Int, spin::Bool) = OpID(4 + spin, k)
+
+  os = OpIDSum{4,SimpleWeight,Int}(2 * N^4 + 2 * N^2, op_cache_vec)
+
+  map_1e = Tuple{Int, Int}[]
+  for p in 1:N
+    for q in 1:N
+      push!(map_1e, (p, q))
+      id = length(map_1e)
+      for spin in (↓, ↑)
+        add!(os, SimpleWeight(id), opCdag(p, spin), opC(q, spin))
+      end
+    end
+  end
+
+  map_2e = Tuple{Int, Int, Int, Int, Bool, Bool, Bool, Bool}[]
+  for p in 1:N
+    for s_p in (↓, ↑)
+      for q in 1:N, s_q in (↓, ↑)
+        (q, s_q) <= (p, s_p) && continue
+
+        for r in 1:N, s_r in (↓, ↑)
+          for s in 1:N, s_s in (↓, ↑)
+            (s, s_s) <= (r, s_r) && continue
+
+            valid_entry = (s_r == s_q && s_s == s_p) ||
+              (s_s == s_q && s_r == s_p) ||
+              (s_r == s_p && s_s == s_q) ||
+              (s_s == s_p && s_r == s_q)
+
+            !valid_entry && continue
+
+            push!(map_2e, (p, q, r, s, s_p, s_q, s_r, s_s))
+
+            id = length(map_2e) + length(map_1e)
+            add!(
+              os,
+              SimpleWeight(id),
+              opCdag(p, s_p),
+              opCdag(q, s_q),
+              opC(r, s_r),
+              opC(s, s_s),
+            )
+          end
+        end
+      end
+    end
+  end
+
+  return os, map_1e, map_2e
+end
+
+function electronic_structure_symbolic_coefficients(
+  h::AbstractMatrix,
+  V::AbstractArray{<:Number,4},
+  map_1e::AbstractVector{<:Tuple{Int,Int}},
+  map_2e::AbstractVector{<:Tuple{Int,Int,Int,Int,Bool,Bool,Bool,Bool}},
+)::Vector{promote_type(eltype(h), eltype(V))}
+  coefficients = Vector{promote_type(eltype(h), eltype(V))}(
+    undef, length(map_1e) + length(map_2e)
+  )
+
+  for (id, (p, q)) in pairs(map_1e)
+    coefficients[id] = h[p, q]
+  end
+
+  offset = length(map_1e)
+  for (i, (p, q, r, s, s_p, s_q, s_r, s_s)) in pairs(map_2e)
+    coeff = zero(eltype(coefficients))
+    if s_r == s_q && s_s == s_p
+      coeff += V[p, s, q, r]
+    end
+    if s_s == s_q && s_r == s_p
+      coeff -= V[p, r, q, s]
+    end
+    if s_r == s_p && s_s == s_q
+      coeff -= V[q, s, p, r]
+    end
+    if s_s == s_p && s_r == s_q
+      coeff += V[q, r, p, s]
+    end
+
+    coefficients[offset + i] = coeff
+  end
+
+  return coefficients
+end
+
+function mpo_relative_difference(A::MPO, B::MPO)::Float64
+  AmB = add(A, -B; alg="directsum")
+  return real(inner(AmB, AmB)) / real(inner(B, B))
+end
+
+let
+  N = 6
+  println("Constructing a symbolic electronic structure MPO for $N sites")
+
+  h, V = get_coefficients(N)
+  sites, op_cache_vec = electronic_structure_sites_and_op_cache(N)
+  numeric_os = electronic_structure_OpIDSum(N, h, V, op_cache_vec)
+  symbolic_os, map_1e, map_2e = electronic_structure_symbolic_OpIDSum(N, op_cache_vec)
+  coefficients = electronic_structure_symbolic_coefficients(h, V, map_1e, map_2e)
+
+  @time "Constructing symbolic MPO" sym = MPO_symbolic(
+    symbolic_os,
+    sites;
+    basis_op_cache_vec=op_cache_vec,
+    check_for_errors=false,
+  )
+
+  H_symbolic = instantiate_MPO(
+    sym, coefficients; splitblocks=true, checkflux=false
+  )
+
+  @time "Instantiating symbolic MPO" H_symbolic = instantiate_MPO(
+    sym, coefficients; splitblocks=true, checkflux=false
+  )
+
+  instantiate_MPO!(
+    H_symbolic, sym, coefficients; checkflux=false
+  )
+
+  reset_timer!()
+
+  @time "Instantiating symbolic MPO 2" H_symbolic = instantiate_MPO!(
+    H_symbolic, sym, coefficients; checkflux=false
+  )
+
+  print_timer()
+
+  @time "Constructing reference MPO" H_numeric = MPO_new(
+    numeric_os,
+    sites;
+    alg="VC",
+    basis_op_cache_vec=op_cache_vec,
+    splitblocks=true,
+    check_for_errors=false,
+    checkflux=false,
+  )
+
+  relative_difference = mpo_relative_difference(H_symbolic, H_numeric)
+  println("Symbolic instantiation relative difference: $relative_difference")
+  @assert relative_difference < 1e-10
+end
+
+nothing;
diff --git a/examples/fermi-hubbard-tc.jl b/examples/fermi-hubbard-tc.jl
index 91d9f0e..ad51e83 100644
--- a/examples/fermi-hubbard-tc.jl
+++ b/examples/fermi-hubbard-tc.jl
@@ -431,7 +431,9 @@ using Serialization
 function call_back(
   n::Int,
   offsets::Vector{Vector{Int}},
-  block_sparse_matrices::Vector{Vector{ITensorMPOConstruction.BlockSparseMatrix{ValType}}},
+  block_sparse_matrices::Vector{
+    Vector{ITensorMPOConstruction.BlockSparseMatrix{Matrix{ValType}}}
+  },
   sites::Vector{<:Index},
   llinks::Vector{<:Index},
   g::ITensorMPOConstruction.MPOGraph,
diff --git a/examples/symbolic-transverse-field-ising.jl b/examples/symbolic-transverse-field-ising.jl
new file mode 100644
index 0000000..03e92ed
--- /dev/null
+++ b/examples/symbolic-transverse-field-ising.jl
@@ -0,0 +1,112 @@
+# # Symbolic all-to-all transverse-field Ising model
+#
+# The all-to-all transverse-field Ising model is
+# ```math
+# H = \sum_{i < j} J_{ij} Z_i Z_j + h \sum_i X_i .
+# ```
+# If the connectivity changes many times while the operator structure stays
+# fixed, the symbolic MPO path can construct the MPO once and instantiate
+# it repeatedly for different $J_{ij}$ and $h$.
+
+using ITensors, ITensorMPS, ITensorMPOConstruction
+using Random
+using Test
+
+function num_tfim_couplings(N::Int)::Int
+  return (N * (N - 1)) ÷ 2
+end
+
+function all_to_all_tfim_symbolic(N::Int)::SymbolicMPO
+  sites = siteinds("Qubit", N)
+  op_cache_vec = to_OpCacheVec(sites, [["I", "X", "Z"] for _ in 1:N])
+
+  X(n::Int) = OpID(2, n)
+  Z(n::Int) = OpID(3, n)
+
+  num_couplings = num_tfim_couplings(N)
+  os = OpIDSum{2,SimpleWeight,Int}(num_couplings + N, op_cache_vec)
+
+  label = 1
+  for i in 1:N
+    for j in (i + 1):N
+      add!(os, SimpleWeight(label), Z(i), Z(j))
+      label += 1
+    end
+  end
+
+  for i in 1:N
+    add!(os, SimpleWeight(label), X(i))
+  end
+
+  return MPO_symbolic(os, sites)
+end
+
+function all_to_all_tfim_numeric(
+  sites::Vector{<:Index}, coupling_weights::AbstractVector, h::Real
+)::MPO
+  N = length(sites)
+  num_couplings = num_tfim_couplings(N)
+
+  op_cache_vec = to_OpCacheVec(sites, [["I", "X", "Z"] for _ in 1:N])
+
+  X(n::Int) = OpID(2, n)
+  Z(n::Int) = OpID(3, n)
+
+  os = OpIDSum{2,Float64,Int}(num_couplings + N, op_cache_vec)
+
+  label = 1
+  for i in 1:N
+    for j in (i + 1):N
+      add!(os, coupling_weights[label], Z(i), Z(j))
+      label += 1
+    end
+  end
+
+  for i in 1:N
+    add!(os, h, X(i))
+  end
+
+  return MPO_new(os, sites; alg="VC", splitblocks=true) # TODO: remove splitbocks=true after warning
+end
+
+function tfim_symbolic_coefficients(coupling_weights::AbstractVector, h::Real)::Vector{Float64}
+  coefficients = Vector{Float64}(undef, length(coupling_weights) + 1)
+  coefficients[1:(end - 1)] .= coupling_weights
+  coefficients[end] = h
+  return coefficients
+end
+
+function test_symbolic_tfim_sample(
+  sym::SymbolicMPO, coupling_weights::AbstractVector, h::Real
+)::Float64
+  
+  symbolic_mpo = instantiate_MPO(
+    sym, tfim_symbolic_coefficients(coupling_weights, h)
+  )
+
+  numeric_mpo = all_to_all_tfim_numeric(sym.sites, coupling_weights, h)
+
+  return norm(add(symbolic_mpo, -numeric_mpo; alg="directsum"))
+end
+
+@testset "Symbolic TFIM" begin
+  N = 50
+  h = 1.0
+  num_samples = 5
+  
+  sym = all_to_all_tfim_symbolic(N)
+  num_couplings = num_tfim_couplings(N)
+
+  println(
+    "Constructed a symbolic all-to-all transverse-field Ising MPO template ",
+    "for $N sites",
+  )
+
+  for sample in 1:num_samples
+    error = test_symbolic_tfim_sample(sym, randn(num_couplings), h)
+    println("Sample $sample MPO error: $(round(error; sigdigits=3))")
+    @test iszero(error)
+  end
+end
+
+nothing
diff --git a/src/ITensorMPOConstruction.jl b/src/ITensorMPOConstruction.jl
index c81a4f9..d73ffb9 100644
--- a/src/ITensorMPOConstruction.jl
+++ b/src/ITensorMPOConstruction.jl
@@ -18,6 +18,8 @@ using ThreadsX
 #
 include("time_if.jl")
 include("OpIDSum.jl")
+include("SymbolicWeight.jl")
+include("SymbolicLocalMatrix.jl")
 include("ops.jl")
 include("BipartiteGraph.jl")
 include("connected-components.jl")
@@ -35,7 +37,10 @@ include("MPOConstruction.jl")
 # OpIDSum.jl
 export OpInfo, OpCacheVec, to_OpCacheVec, OpID, OpIDSum
 
+export SimpleWeight
+
 # MPOConstruction.jl
-export resume_MPO_construction!, MPO_new, instantiate_MPO, sparsity, block2_nnz
+export resume_MPO_construction!,
+  MPO_new, SymbolicMPO, MPO_symbolic, instantiate_MPO, instantiate_MPO!, sparsity, block2_nnz
 
 end
diff --git a/src/MPOConstruction.jl b/src/MPOConstruction.jl
index 5a8edd8..b4523c0 100644
--- a/src/MPOConstruction.jl
+++ b/src/MPOConstruction.jl
@@ -45,7 +45,7 @@ $_resume_MPO_kwargs
 function resume_MPO_construction!(
   n_init::Int,
   offsets::Vector{Vector{Int}},
-  block_sparse_matrices::Vector{Vector{BlockSparseMatrix{ValType}}},
+  block_sparse_matrices::Vector{Vector{BlockSparseMatrix{MatrixType}}},
   sites::Vector{<:Index},
   llinks::Vector{<:Index},
   g::MPOGraph,
@@ -57,14 +57,14 @@ function resume_MPO_construction!(
   call_back::Function=(
     cur_site::Int,
     offsets::Vector{Vector{Int}},
-    block_sparse_matrices::Vector{Vector{BlockSparseMatrix{ValType}}},
+    block_sparse_matrices::Vector{Vector{BlockSparseMatrix{MatrixType}}},
     sites::Vector{<:Index},
     llinks::Vector{<:Index},
     cur_graph::MPOGraph,
     op_cache_vec::OpCacheVec,
   ) -> nothing,
   output_level::Int=0,
-)::Nothing where {ValType<:Number}
+)::Nothing where {MatrixType}
   @assert !ITensors.using_auto_fermion() # TODO: This should be fixed.
   @assert tol >= 0
 
@@ -73,7 +73,7 @@ function resume_MPO_construction!(
       "At site $n/$(length(sites)) the graph takes up $(Base.format_bytes(Base.summarysize(g)))",
     )
     @time_if output_level 1 "at_site!" g, offsets[n], block_sparse_matrices[n], llinks[n + 1] = at_site!(
-      ValType,
+      MatrixType,
       g,
       n,
       sites,
@@ -91,6 +91,26 @@ function resume_MPO_construction!(
   return nothing
 end
 
+function _remove_symbolic_mpo_boundary_links!(H::MPO, llinks::Vector{<:Index})::Nothing
+  @assert ndims(H[1]) == 4
+  @assert ndims(H[end]) == 4
+
+  L = ITensor(llinks[1])
+  L[end] = 1.0
+  H[1] *= L
+
+  R = ITensor(dag(llinks[end]))
+  R[1] = 1.0
+  H[end] *= R
+
+  @assert ndims(H[1]) == 3
+  @assert ndims(H[end]) == 3
+
+  return nothing
+end
+
+# TODO: Move splitblocks warning to all instantiate_MPO calls
+
 """
     instantiate_MPO(offsets, block_sparse_matrices, sites, llinks; splitblocks, checkflux) -> MPO
 
@@ -103,7 +123,7 @@ Afterward, the dummy left and right boundary links are contracted away. If
 """
 function instantiate_MPO(
   offsets::Vector{Vector{Int}},
-  block_sparse_matrices::Vector{Vector{BlockSparseMatrix{ValType}}},
+  block_sparse_matrices::Vector{Vector{BlockSparseMatrix{Matrix{ValType}}}},
   sites::Vector{<:Index},
   llinks::Vector{<:Index};
   splitblocks::Bool,
@@ -111,24 +131,429 @@ function instantiate_MPO(
 )::MPO where {ValType<:Number}
   H = MPO(sites)
 
-  @timeit "to_ITensor" Threads.@threads for n in 1:length(sites)
+  @timeit "to_ITensor" for n in 1:length(sites)
     H[n] = to_ITensor(
       offsets[n], block_sparse_matrices[n], llinks[n], llinks[n + 1], sites[n]; splitblocks
     )
   end
 
-  # Remove the dummy link indices from the left and right.
-  L = ITensor(llinks[1])
-  L[end] = 1.0
+  _remove_symbolic_mpo_boundary_links!(H, llinks)
 
-  R = ITensor(dag(llinks[end]))
-  R[1] = 1.0
+  if checkflux
+    @timeit "checkflux" Threads.@threads for n in 1:length(sites)
+      ITensors.checkflux(H[n])
+    end
+  end
+
+  return H
+end
+
+"""
+    SymbolicMPO
+
+Intermediate symbolic MPO representation produced by [`MPO_symbolic`](@ref).
+
+The stored `block_sparse_matrices` have the same structural layout as numeric
+MPO construction, but each block entry is a symbolic local matrix term list.
+`max_user_label` records the largest positive coefficient label needed for
+later numeric substitution.
+"""
+struct SymbolicMPO{
+  W<:AbstractWeight,Ti<:Integer,Sites<:AbstractVector{<:Index},Links<:AbstractVector{<:Index}
+}
+  offsets::Vector{Vector{Int}}
+  block_sparse_matrices::Vector{Vector{BlockSparseMatrix{SymbolicLocalMatrix{W, Ti}}}}
+  sites::Sites
+  llinks::Links
+  op_cache_vec::OpCacheVec
+end
+
+function _symbolic_mpo_eltype(sym::SymbolicMPO, coefficients::AbstractVector)::Type
+  val_type = promote_type(Float64, eltype(coefficients))
+  for op_cache in sym.op_cache_vec
+    for op_info in op_cache
+      val_type = promote_type(val_type, eltype(op_info.matrix))
+    end
+  end
+  return val_type
+end
+
+function _evaluate_symbolic_block_sparse_matrix(
+  symbolic_matrix::BlockSparseMatrix{SymbolicLocalMatrix{W, Ti}},
+  coefficients::AbstractVector,
+  op_cache::Vector{OpInfo},
+  ::Type{ValType},
+)::BlockSparseMatrix{Matrix{ValType}} where {W,Ti,ValType}
+  matrix = BlockSparseMatrix(Matrix{ValType}, length(symbolic_matrix))
+
+  for right_link in eachindex(symbolic_matrix)
+    for (left_link, terms) in pairs(symbolic_matrix[right_link])
+      block = _evaluate_symbolic_local_matrix(ValType, terms, coefficients, op_cache)
+      set!(matrix[right_link], left_link, block)
+    end
+  end
+
+  return matrix
+end
+
+function _evaluate_symbolic_block_sparse_matrices(
+  symbolic_matrices::Vector{BlockSparseMatrix{SymbolicLocalMatrix{W, Ti}}},
+  coefficients::AbstractVector,
+  op_cache::Vector{OpInfo},
+  ::Type{ValType},
+)::Vector{BlockSparseMatrix{Matrix{ValType}}} where {W, Ti,ValType}
+  matrices = Vector{BlockSparseMatrix{Matrix{ValType}}}(undef, length(symbolic_matrices))
+  for cc in eachindex(symbolic_matrices)
+    matrices[cc] = _evaluate_symbolic_block_sparse_matrix(
+      symbolic_matrices[cc], coefficients, op_cache, ValType
+    )
+  end
+  return matrices
+end
+
+"""
+    instantiate_MPO(sym::SymbolicMPO, coefficients; splitblocks=true, checkflux=false) -> MPO
 
+Evaluate a symbolic MPO with a numeric coefficient vector and return a fresh
+numeric `MPO`.
+
+The coefficient vector is indexed by the original positive labels used to build
+`sym`. Each site's symbolic block storage is first evaluated into ordinary
+`BlockSparseMatrix` storage, then converted to ITensors through the same
+assembly path used by numeric MPO construction.
+"""
+function instantiate_MPO(
+  sym::SymbolicMPO, coefficients::AbstractVector; splitblocks::Bool=true, checkflux::Bool=false
+)::MPO
+  C = _symbolic_mpo_eltype(sym, coefficients)
+
+  H = MPO(sym.sites)
+
+  @timeit "to_ITensor" for n in eachindex(sym.sites)
+    @timeit "_evaluate_symbolic_block_sparse_matrices" block_sparse_matrices = _evaluate_symbolic_block_sparse_matrices(
+      sym.block_sparse_matrices[n], coefficients, sym.op_cache_vec[n], C
+    )
+
+    @timeit "to_ITensor" H[n] = to_ITensor(
+      sym.offsets[n],
+      block_sparse_matrices,
+      sym.llinks[n],
+      sym.llinks[n + 1],
+      sym.sites[n],
+      ;
+      splitblocks,
+    )
+  end
+
+  _remove_symbolic_mpo_boundary_links!(H, sym.llinks)
+
+  if checkflux
+    @timeit "checkflux" Threads.@threads for n in eachindex(sym.sites)
+      ITensors.checkflux(H[n])
+    end
+  end
+
+  return H
+end
+
+function _add_symbolic_mpo_boundary_links!(H::MPO, sym::SymbolicMPO)::Nothing
+  L = ITensor(dag(sym.llinks[1]))
+  L[1] = 1.0
   H[1] *= L
-  H[length(sites)] *= R
+
+  R = ITensor(sym.llinks[end])
+  R[1] = 1.0
+  H[end] *= R
+
+  return nothing
+end
+
+function _fill_symbolic_mpo_template_tensor!(
+  tensor::NDTensors.DenseTensor{C,4},
+  sym::SymbolicMPO,
+  coefficients::AbstractVector,
+  n::Int,
+)::Nothing where {C}
+  op_cache = sym.op_cache_vec[n]
+  tensor_data = array(tensor)
+
+  fill!(tensor_data, zero(C))
+  @timeit "evaluate and fill entries" for (offset, matrix) in zip(
+    sym.offsets[n], sym.block_sparse_matrices[n]
+  )
+    for (right_link, column) in enumerate(matrix)
+      shifted_right_link = right_link + offset
+      for (left_link, terms) in pairs(column)
+        local_matrix = _evaluate_symbolic_local_matrix(C, terms, coefficients, op_cache)
+        @inbounds for site_prime_link in axes(local_matrix, 1),
+          site_link in axes(local_matrix, 2)
+
+          value = local_matrix[site_prime_link, site_link]
+          iszero(value) && continue
+          tensor_data[left_link, shifted_right_link, site_prime_link, site_link] = value
+        end
+      end
+    end
+  end
+
+  return nothing
+end
+
+@inline function _block_key(block::Block{4})::NTuple{4,Int}
+  return (Int(block[1]), Int(block[2]), Int(block[3]), Int(block[4]))
+end
+
+@inline function _block_storage_index(
+  block_offset::Int, entry_offsets::NTuple{4,Int}, block_dims::NTuple{4,Int}
+)::Int
+  storage_index = block_offset + entry_offsets[1]
+  stride = block_dims[1]
+  for dim in 2:4
+    storage_index += (entry_offsets[dim] - 1) * stride
+    stride *= block_dims[dim]
+  end
+  return storage_index
+end
+
+function _try_fill_symbolic_mpo_splitblocks_sequential!(
+  tensor::NDTensors.BlockSparseTensor{C,4},
+  sym::SymbolicMPO,
+  coefficients::AbstractVector,
+  n::Int,
+  left_blocks::Vector{Int},
+  right_blocks::Vector{Int},
+  site_prime_blocks::Vector{Int},
+  site_blocks::Vector{Int},
+  tensor_data,
+)::Bool where {C}
+  length(tensor_data) == NDTensors.nnzblocks(tensor) || return false
+
+  op_cache = sym.op_cache_vec[n]
+  block_iterator = ITensors.eachnzblock(tensor)
+  block_state = iterate(block_iterator)
+  current_block_key = isnothing(block_state) ? nothing : _block_key(block_state[1])
+  storage_index = 1
+
+  fill!(tensor_data, zero(C))
+
+  local_matrix = zeros(C, size(op_cache[1].matrix))
+  @timeit "sequential splitblock refill" for (offset, matrix) in zip(
+    sym.offsets[n], sym.block_sparse_matrices[n]
+  )
+    for (right_link, column) in enumerate(matrix)
+      shifted_right_link = right_link + offset
+      right_block = right_blocks[shifted_right_link]
+      for (left_link, terms) in pairs(column)
+        left_block = left_blocks[left_link]
+        _evaluate_symbolic_local_matrix!(local_matrix, terms, coefficients, op_cache)
+
+        @inbounds for site_prime_link in axes(local_matrix, 1),
+          site_link in axes(local_matrix, 2)
+
+          expected_block_key = (
+            left_block,
+            right_block,
+            site_prime_blocks[site_prime_link],
+            site_blocks[site_link],
+          )
+          value = local_matrix[site_prime_link, site_link]
+
+          if !isnothing(current_block_key) && expected_block_key == current_block_key
+            tensor_data[storage_index] = value
+            storage_index += 1
+            block_state = iterate(block_iterator, block_state[2])
+            current_block_key = isnothing(block_state) ? nothing : _block_key(block_state[1])
+          elseif !iszero(value)
+            return false
+          end
+        end
+      end
+    end
+  end
+
+  return isnothing(current_block_key)
+end
+
+function _fill_symbolic_mpo_template_tensor!(
+  tensor::NDTensors.BlockSparseTensor{C,4},
+  sym::SymbolicMPO,
+  coefficients::AbstractVector,
+  n::Int,
+)::Nothing where {C}
+  tensor_inds = inds(tensor)
+
+  left_blocks, left_offsets = _qn_position_map(tensor_inds[1])
+  right_blocks, right_offsets = _qn_position_map(tensor_inds[2])
+  site_prime_blocks, site_prime_offsets = _qn_position_map(tensor_inds[3])
+  site_blocks, site_offsets = _qn_position_map(tensor_inds[4])
+
+  tensor_data = NDTensors.data(tensor)
+  if _try_fill_symbolic_mpo_splitblocks_sequential!(
+    tensor,
+    sym,
+    coefficients,
+    n,
+    left_blocks,
+    right_blocks,
+    site_prime_blocks,
+    site_blocks,
+    tensor_data,
+  )
+    return nothing
+  end
+
+  block_offsets = sizehint!(Dict{NTuple{4,Int},Int}(), NDTensors.nnzblocks(tensor))
+
+  
+  fill!(tensor_data, zero(C))
+
+  @timeit "index template blocks" begin
+    for block in ITensors.eachnzblock(tensor)
+      block_offsets[_block_key(block)] = NDTensors.offset(NDTensors.blockoffsets(tensor), block)
+    end
+  end
+
+  @timeit "evaluate and fill entries" _fill_symbolic_mpo_tensor_entries!(
+    tensor,
+    sym,
+    coefficients,
+    n,
+    left_blocks,
+    right_blocks,
+    site_prime_blocks,
+    site_blocks,
+    left_offsets,
+    right_offsets,
+    site_prime_offsets,
+    site_offsets,
+    tensor_inds,
+    tensor_data,
+    block_offsets,
+  )
+
+  return nothing
+end
+
+function _fill_symbolic_mpo_tensor_entries!(
+  tensor::NDTensors.BlockSparseTensor{C,4},
+  sym::SymbolicMPO,
+  coefficients::AbstractVector,
+  n::Int,
+  left_blocks::Vector{Int},
+  right_blocks::Vector{Int},
+  site_prime_blocks::Vector{Int},
+  site_blocks::Vector{Int},
+  left_offsets::Vector{Int},
+  right_offsets::Vector{Int},
+  site_prime_offsets::Vector{Int},
+  site_offsets::Vector{Int},
+  tensor_inds,
+  tensor_data,
+  block_offsets::Dict{NTuple{4,Int},Int},
+)::Nothing where {C}
+  op_cache = sym.op_cache_vec[n]
+
+  for (offset, matrix) in zip(sym.offsets[n], sym.block_sparse_matrices[n])
+    for (right_link, column) in enumerate(matrix)
+      shifted_right_link = right_link + offset
+      right_block = right_blocks[shifted_right_link]
+      right_offset = right_offsets[shifted_right_link]
+      for (left_link, terms) in pairs(column)
+        left_block = left_blocks[left_link]
+        left_offset = left_offsets[left_link]
+        local_matrix = _evaluate_symbolic_local_matrix(C, terms, coefficients, op_cache)
+
+        @inbounds for site_prime_link in axes(local_matrix, 1),
+          site_link in axes(local_matrix, 2)
+
+          value = local_matrix[site_prime_link, site_link]
+          iszero(value) && continue
+
+          site_prime_block = site_prime_blocks[site_prime_link]
+          site_block = site_blocks[site_link]
+          block_id = (left_block, right_block, site_prime_block, site_block)
+
+          block_offset = get(block_offsets, block_id, nothing)
+          if isnothing(block_offset)
+            tensor[left_link, shifted_right_link, site_prime_link, site_link] = value
+          else
+            site_prime_offset = site_prime_offsets[site_prime_link]
+            site_offset = site_offsets[site_link]
+            entry_offsets = (left_offset, right_offset, site_prime_offset, site_offset)
+            block_dims = ntuple(i -> blockdim(tensor_inds[i], block_id[i]), Val(4))
+            tensor_data[_block_storage_index(block_offset, entry_offsets, block_dims)] = value
+          end
+        end
+      end
+    end
+  end
+
+  return nothing
+end
+
+function _check_symbolic_mpo_template(H::MPO, sym::SymbolicMPO)::Nothing
+  if length(H) != length(sym.sites)
+    throw(
+      ArgumentError(
+        "Template MPO length $(length(H)) is incompatible with SymbolicMPO length $(length(sym.sites)).",
+      ),
+    )
+  end
+
+  for n in eachindex(sym.sites)
+    if !hasinds(H[n], prime(sym.sites[n]), dag(sym.sites[n]))
+      throw(ArgumentError("Template MPO site index at site $n is incompatible with SymbolicMPO."))
+    end
+  end
+
+  for n in 2:(length(sym.llinks) - 1)
+    template_link = linkind(H, n - 1)
+    if isnothing(template_link) || template_link != sym.llinks[n]
+      throw(ArgumentError("Template MPO link index at bond $(n - 1) is incompatible with SymbolicMPO."))
+    end
+  end
+
+  return nothing
+end
+
+"""
+    instantiate_MPO!(H::MPO, sym::SymbolicMPO, coefficients; checkflux=false) -> MPO
+
+Evaluate `sym` with `coefficients` into the existing MPO `H`.
+
+The target MPO fixes the tensor block layout, so this overload does not accept
+`splitblocks`. Site indices must match `sym`, and internal link dimensions must
+match the symbolic link dimensions. The dummy boundary links are temporarily
+restored so each tensor can be zeroed and refilled through the same four-index
+layout used by numeric tensor assembly, then contracted away again. The MPO
+object `H` is updated in place and returned.
+"""
+function instantiate_MPO!(
+  H::MPO, sym::SymbolicMPO, coefficients::AbstractVector; checkflux::Bool=false
+)::MPO
+  _check_symbolic_mpo_template(H, sym)
+
+  _add_symbolic_mpo_boundary_links!(H, sym)
+
+  @timeit "refill template tensors" for n in eachindex(sym.sites)
+    t = H[n]
+    llink = dag(sym.llinks[n])
+    rlink = sym.llinks[n + 1]
+    site = sym.sites[n]
+    t = ITensors.permute(
+      t, llink, rlink, prime(site), dag(site); allow_alias=true
+    )
+
+    _fill_symbolic_mpo_template_tensor!(t.tensor, sym, coefficients, n)
+    H[n] = t
+  end
+
+  _remove_symbolic_mpo_boundary_links!(
+    H, sym.llinks
+  )
 
   if checkflux
-    @timeit "checkflux" Threads.@threads for n in 1:length(sites)
+    @timeit "checkflux" Threads.@threads for n in eachindex(sym.sites)
       ITensors.checkflux(H[n])
     end
   end
@@ -136,6 +561,74 @@ function instantiate_MPO(
   return H
 end
 
+"""
+    instantiate_MPO(H_template::MPO, sym::SymbolicMPO, coefficients; checkflux=false) -> MPO
+
+Copy a compatible MPO template and evaluate `sym` into the copied layout.
+
+The template fixes the tensor block layout, so this overload does not accept
+`splitblocks`.
+"""
+function instantiate_MPO(
+  H_template::MPO, sym::SymbolicMPO, coefficients::AbstractVector; checkflux::Bool=false
+)::MPO
+  return instantiate_MPO!(deepcopy(H_template), sym, coefficients; checkflux)
+end
+
+"""
+    MPO_symbolic(os::OpIDSum, sites; basis_op_cache_vec=nothing,
+      check_for_errors=true, combine_qn_sectors=false, output_level=0) -> SymbolicMPO
+
+Run the vertex-cover MPO construction once for an integer-labeled `OpIDSum`,
+storing symbolic local matrix terms that can be numerically instantiated later.
+
+The input `OpIDSum` is preprocessed in place, matching `MPO_new`: positive user
+coefficient labels are remapped to internal symbolic ids, an optional basis
+rewrite is applied, and duplicate symbolic graph entries are preserved.
+"""
+function MPO_symbolic(
+  os::OpIDSum{N,W,Ti},
+  sites::Vector{<:Index};
+  basis_op_cache_vec=nothing,
+  check_for_errors::Bool=true,
+  combine_qn_sectors::Bool=false,
+  output_level::Int=0,
+  kwargs...,
+)::SymbolicMPO where {N,W<:AbstractWeight,Ti}
+  prepare_opID_sum!(os, to_OpCacheVec(sites, basis_op_cache_vec))
+  check_for_errors && check_os_for_errors(os)
+
+  label = "Constructing symbolic MPOGraph from $(length(os)) terms"
+  @time_if output_level 0 label g = MPOGraph(os)
+
+  llinks = Vector{Index}(undef, length(sites) + 1)
+  if hasqns(sites)
+    llinks[1] = Index(QN() => 1; tags="Link,l=0", dir=ITensors.Out)
+  else
+    llinks[1] = Index(1; tags="Link,l=0")
+  end
+
+  offsets = Vector{Vector{Int}}(undef, length(sites))
+  block_sparse_matrices = Vector{Vector{BlockSparseMatrix{SymbolicLocalMatrix{W, Ti}}}}(
+    undef, length(sites)
+  )
+
+  @time_if output_level 0 "Constructing symbolic MPO terms" resume_MPO_construction!(
+    1,
+    offsets,
+    block_sparse_matrices,
+    sites,
+    llinks,
+    g,
+    os.op_cache_vec;
+    alg="VC",
+    combine_qn_sectors,
+    output_level,
+  )
+
+  return SymbolicMPO(offsets, block_sparse_matrices, sites, llinks, os.op_cache_vec)
+end
+
 @doc """
     MPO_new(ValType, os::OpIDSum, sites; basis_op_cache_vec=nothing,
       check_for_errors=true, checkflux=true, splitblocks=true, output_level=0,
@@ -193,7 +686,9 @@ function MPO_new(
   end
 
   offsets = Vector{Vector{Int}}(undef, length(sites))
-  block_sparse_matrices = Vector{Vector{BlockSparseMatrix{ValType}}}(undef, length(sites))
+  block_sparse_matrices = Vector{Vector{BlockSparseMatrix{Matrix{ValType}}}}(
+    undef, length(sites)
+  )
 
   @time_if output_level 0 "Constructing MPO terms" resume_MPO_construction!(
     1,
diff --git a/src/OpIDSum.jl b/src/OpIDSum.jl
index e02eef5..8a8667b 100644
--- a/src/OpIDSum.jl
+++ b/src/OpIDSum.jl
@@ -39,8 +39,6 @@ function OpInfo(local_op::Op, site::Index)
   )
 end
 
-# TODO: Create a symbolic sum.
-
 """
     OpCacheVec
 
@@ -192,7 +190,7 @@ mutable struct OpIDSum{N,C,Ti}
   op_cache_vec::OpCacheVec
   abs_tol::Float64
   modify!::FunctionWrappers.FunctionWrapper{
-    C,
+    Float64,
     Tuple{
       SubArray{
         OpID{Ti},
@@ -229,7 +227,7 @@ function OpIDSum{N,C,Ti}(
   data = Vector{NTuple{N,OpID{Ti}}}(undef, max_terms)
   terms = reinterpret(reshape, OpID{Ti}, data)
   wrapped_modify! = FunctionWrappers.FunctionWrapper{
-    C,
+    Float64,
     Tuple{
       SubArray{
         OpID{Ti},
@@ -265,7 +263,7 @@ Terms with `abs(scalar) <= abs_tol` are ignored when added.
 function OpIDSum{N,C,Ti}(
   max_terms::Int, op_cache_vec::OpCacheVec; abs_tol::Real=0
 )::OpIDSum{N,C,Ti} where {N,C,Ti}
-  return OpIDSum{N,C,Ti}(max_terms, op_cache_vec, ops -> 1; abs_tol)
+  return OpIDSum{N,C,Ti}(max_terms, op_cache_vec, ops -> 1.0; abs_tol)
 end
 
 """
@@ -308,7 +306,7 @@ term-modification callback associated with `os` is applied. Terms below
 merge duplicate terms or check capacity; callers must construct `os` with enough
 room for all accepted terms.
 """
-function ITensorMPS.add!(os::OpIDSum, scalar::Number, ops)::Nothing
+function ITensorMPS.add!(os::OpIDSum, scalar, ops)::Nothing
   abs(scalar) <= os.abs_tol && return nothing
 
   num_terms = Threads.atomic_add!(os.num_terms, 1) + 1
@@ -340,7 +338,7 @@ end
 
 Vararg convenience method for `add!`.
 """
-function ITensorMPS.add!(os::OpIDSum, scalar::Number, ops::OpID...)::Nothing
+function ITensorMPS.add!(os::OpIDSum, scalar, ops::OpID...)::Nothing
   return add!(os, scalar, ops)
 end
 
@@ -386,7 +384,9 @@ function for_equal_sites(f::Function, ops::AbstractVector{<:OpID})::Nothing
 end
 
 """
-    rewrite_in_operator_basis!(os::OpIDSum, basis_op_cache_vec::OpCacheVec)
+    rewrite_in_operator_basis!(
+      os::OpIDSum, basis_op_cache_vec::OpCacheVec; symbolic_coefficients=false
+    ) -> Nothing
 
 Rewrite each same-site operator product in `os` into a single operator drawn
 from `basis_op_cache_vec`.
@@ -399,10 +399,15 @@ matrix, the whole term is set to zero. An error is thrown if a nonzero product
 cannot be represented by a single basis operator. After rewriting, each
 nonzero term is resorted by site and `os.op_cache_vec` is replaced by
 `basis_op_cache_vec`.
+
+When `symbolic_coefficients=true`, the stored coefficients are interpreted as
+internal signed symbolic ids. In that mode, rewrite factors must be exactly
+`+1` or `-1`; unsupported factors throw an `ArgumentError` instead of
+converting the symbolic id into an ordinary numeric coefficient.
 """
 @timeit function rewrite_in_operator_basis!(
   os::OpIDSum{N,C,Ti}, basis_op_cache_vec::OpCacheVec
-) where {N,C,Ti}
+)::Nothing where {N,C,Ti}
   op_cache_vec = os.op_cache_vec
 
   function scale_by_first_nz!(matrix::Matrix)
@@ -500,7 +505,7 @@ nonzero term is resorted by site and `os.op_cache_vec` is replaced by
 
   os.op_cache_vec = basis_op_cache_vec
 
-  return os
+  return nothing
 end
 
 """
@@ -622,15 +627,21 @@ term has even fermion parity.
 end
 
 """
-    prepare_opID_sum!(os::OpIDSum, basis_op_cache_vec::Union{Nothing,OpCacheVec}) -> Nothing
+    prepare_opID_sum!(
+      os::OpIDSum, basis_op_cache_vec::Union{Nothing,OpCacheVec};
+      symbolic_coefficients=false
+    ) -> Nothing
 
 Apply optional preprocessing to `os` before MPO construction.
 
 Currently this rewrites terms into `basis_op_cache_vec` when a basis cache is
-provided and otherwise leaves `os` unchanged.
+provided and otherwise leaves `os` unchanged. Set `symbolic_coefficients=true`
+when integer coefficients are internal signed symbolic ids rather than numeric
+scalars.
 """
 function prepare_opID_sum!(
-  os::OpIDSum, basis_op_cache_vec::Union{Nothing,OpCacheVec}
+  os::OpIDSum,
+  basis_op_cache_vec::Union{Nothing,OpCacheVec}
 )::Nothing
   if !isnothing(basis_op_cache_vec)
     rewrite_in_operator_basis!(os, basis_op_cache_vec)
diff --git a/src/SymbolicLocalMatrix.jl b/src/SymbolicLocalMatrix.jl
new file mode 100644
index 0000000..d1d5b15
--- /dev/null
+++ b/src/SymbolicLocalMatrix.jl
@@ -0,0 +1,54 @@
+"""
+    SymbolicLocalMatrix{Ti}
+
+Symbolic weighted sum of cached local operators for one MPO block entry.
+
+Each term is `(signed_weight_id, signed_local_op_id)`. The signed weight id is
+an internal symbolic coefficient id, where `+1` and `-1` represent constants,
+and `+k` or `-k` for `k > 1` map to user coefficient `coefficients[k - 1]` with
+the stored sign. The absolute value of `signed_local_op_id` is the local
+operator id; negative local operator ids indicate that a Jordan-Wigner string
+should be applied while evaluating the cached operator matrix.
+"""
+const SymbolicLocalMatrix{Weight<:AbstractWeight, Ti<:Integer} = Vector{Tuple{Weight, Ti}}
+
+function _append_symbolic_local_matrix_term!(
+  terms::SymbolicLocalMatrix{W, Ti},
+  weight::W,
+  signed_op_id::Ti,
+)::Nothing where {W, Ti}
+  push!(terms, (weight, signed_op_id))
+  return nothing
+end
+
+function _evaluate_symbolic_local_matrix!(
+  result::Matrix{C}, terms::SymbolicLocalMatrix, coefficients::AbstractVector, op_cache::Vector{OpInfo}
+)::Nothing where {C}
+  result .= 0
+
+  needs_jw = 0
+  for (signed_weight_id, signed_local_op_id) in terms
+    local_op_id = abs(Int(signed_local_op_id))
+    local_op = op_cache[local_op_id].matrix
+    weight = substitute_weight(signed_weight_id, coefficients)
+
+    needs_jw += signed_local_op_id < 0
+    add_to_local_matrix!(
+      result, weight, local_op, false
+    )
+  end
+
+  @assert needs_jw ∈ (0, length(terms))
+  needs_jw > 0 && apply_jw_string!(result)
+
+  return nothing
+end
+
+# TODO: Make sure this isn't used anywhere
+function _evaluate_symbolic_local_matrix(
+  ::Type{C}, terms::SymbolicLocalMatrix, coefficients::AbstractVector, op_cache::Vector{OpInfo}
+)::Matrix{C} where {C}
+  result = Matrix{C}(undef, size(op_cache[1].matrix))
+  _evaluate_symbolic_local_matrix!(result, terms, coefficients, op_cache)
+  return result
+end
diff --git a/src/SymbolicWeight.jl b/src/SymbolicWeight.jl
new file mode 100644
index 0000000..ae1ac62
--- /dev/null
+++ b/src/SymbolicWeight.jl
@@ -0,0 +1,28 @@
+# TODO Rename file to AbstractWeight
+abstract type AbstractWeight end
+
+# Abs is used to drop terms in OpIDSum which have small coefficients.
+Base.abs(::AbstractWeight) = Inf64
+
+struct SimpleWeight <: AbstractWeight
+  id::Int
+end
+
+Base.one(::Type{SimpleWeight}) = SimpleWeight(typemax(Int))
+
+function Base.:*(s::SimpleWeight, x)::SimpleWeight
+  x == +1 && return s
+  x == -1 && return SimpleWeight(-s.id)
+  error("SimpleWeight can only be multiplied by ±1, not $x")
+end
+
+function substitute_weight(weight::SimpleWeight, coefficients::Vector{ValType})::ValType where {ValType}
+  abs_weight_id = abs(weight.id)
+
+  value = one(ValType)
+  if abs_weight_id != typemax(Int)
+    value = coefficients[abs_weight_id]
+  end
+
+  return weight.id > 0 ? value : -value
+end
\ No newline at end of file
diff --git a/src/large-graph-mpo-common.jl b/src/large-graph-mpo-common.jl
index 839c79c..591a60a 100644
--- a/src/large-graph-mpo-common.jl
+++ b/src/large-graph-mpo-common.jl
@@ -1,15 +1,21 @@
 """
-    BlockSparseMatrix{C}
+    BlockSparseMatrix{MatrixType}
 
 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 one local block value. Numeric construction
+uses dense local operator matrices as `BlockSparseMatrix{Matrix{C}}`, while the
+symbolic path stores term lists as `BlockSparseMatrix{SymbolicLocalMatrix}`.
 `at_site!` later returns offsets that place component-local right-link ids into
 the full outgoing MPO bond.
 """
-BlockSparseMatrix{C} = Vector{Dictionary{Int,Matrix{C}}}
+BlockSparseMatrix{MatrixType} = Vector{Dictionary{Int,MatrixType}}
+
+function BlockSparseMatrix(::Type{MatrixType}, ncols::Int)::BlockSparseMatrix{MatrixType} where {MatrixType}
+  return [Dictionary{Int,MatrixType}() for _ in 1:ncols]
+end
 
 """
     MPOGraph{N,C,Ti}
@@ -200,7 +206,7 @@ function add_to_next_graph!(
 end
 
 """
-    MPOGraph(os::OpIDSum{N,C,Ti}) -> MPOGraph{N,C,Ti}
+    MPOGraph(os::OpIDSum{N,C,Ti}; symbolic_coefficients=false) -> MPOGraph{N,C,Ti}
 
 Convert an `OpIDSum` into the initial bipartite graph used by the
 MPO construction algorithm.
@@ -208,14 +214,18 @@ MPO construction algorithm.
 Operators within each term are put in reverse order (by decreasing site), then
 the terms are sorted along with the scalars. This sorting puts terms that share
 a terminating sequence of operators (which is now at the front of the term) nearby.
-Duplicate terms are then combined, and resulting terms with a weight below
-`os.abs_tol` are dropped. The returned graph is split before the first site:
-left vertices represent the identity incoming link and the operator emitted at
-site 1, while right vertices carry the remaining operator tuples.
+For numeric construction, duplicate terms are then combined and resulting terms
+with a weight below `os.abs_tol` are dropped. When `symbolic_coefficients=true`,
+scalars are interpreted as signed internal symbolic ids, so duplicate operator
+tuples are preserved as separate edge entries instead of being added together.
+The returned graph is split before the first site: left vertices represent the
+identity incoming link and the operator emitted at site 1, while right vertices
+carry the remaining operator tuples.
 """
-@timeit function MPOGraph(os::OpIDSum{N,C,Ti})::MPOGraph{N,C,Ti} where {N,C,Ti}
+@timeit function MPOGraph(
+  os::OpIDSum{N,C,Ti}
+)::MPOGraph{N,C,Ti} where {N,C,Ti}
   @assert size(os.terms, 1) == N
-
   ## Reverse the terms in the sum, ignoring trailing identity operators.
   Threads.@threads for i in 1:length(os)
     for j in N:-1:1
@@ -234,17 +244,26 @@ site 1, while right vertices carry the remaining operator tuples.
     alg=(Threads.nthreads() > 1) ? ThreadsX.QuickSort : Base.QuickSort,
   )
 
-  ## Combine duplicate terms and remove terms which are below the tolerance.
+  ## Combine duplicate numeric terms and remove terms which are below the tolerance.
+  ## Symbolic ids are labels, not numeric weights, so preserve duplicate entries.
   nnz = 0
-  for i in eachindex(os)
-    if i < length(os) && os._data[i] == os._data[i + 1]
-      os.scalars[i + 1] += os.scalars[i]
-      os.scalars[i] = 0
-    elseif abs(os.scalars[i]) > os.abs_tol
+  if C <: AbstractWeight
+    for i in eachindex(os)
       nnz += 1
       os.scalars[nnz] = os.scalars[i]
       os._data[nnz] = os._data[i]
     end
+  else
+    for i in eachindex(os)
+      if i < length(os) && os._data[i] == os._data[i + 1]
+        os.scalars[i + 1] += os.scalars[i]
+        os.scalars[i] = 0
+      elseif abs(os.scalars[i]) > os.abs_tol
+        nnz += 1
+        os.scalars[nnz] = os.scalars[i]
+        os._data[nnz] = os._data[i]
+      end
+    end
   end
 
   os.num_terms[] = nnz
@@ -267,6 +286,24 @@ site 1, while right vertices carry the remaining operator tuples.
   return g
 end
 
+# TODO: Replace add_to_local_matrix! with this
+function apply_jw_string!(a::Matrix)::Nothing
+  if size(a, 1) == 2
+    @inbounds for i in 1:2
+      a[i, 2] *= -1
+    end
+  elseif size(a, 1) == 4
+    @inbounds for i in 1:4
+      a[i, 2] *= -1
+      a[i, 3] *= -1
+    end
+  else
+    error("Unknown fermionic site.")
+  end
+
+  return nothing
+end
+
 """
     add_to_local_matrix!(a, weight, local_op, needs_JW_string) -> Nothing
 
@@ -331,10 +368,20 @@ function merge_qn_sectors(
   return new_order, new_qi
 end
 
+function _check_qr_block_storage(
+  ::Type{BlockSparseMatrix{SymbolicLocalMatrix{Ti}}}
+)::Nothing where {Ti}
+  throw(ArgumentError("QR construction is only supported for numeric block storage."))
+end
+
+function _check_qr_block_storage(::Type)::Nothing
+  return nothing
+end
+
 """
-    at_site!(ValType, g, n, sites, tol, absolute_tol, op_cache_vec, alg;
+    at_site!(MatrixType, g, n, sites, tol, absolute_tol, op_cache_vec, alg;
              combine_qn_sectors, output_level=0)
-        -> Tuple{MPOGraph,Vector{Int},Vector{BlockSparseMatrix{ValType}},Index}
+        -> Tuple{MPOGraph,Vector{Int},Vector{MatrixType},Index}
 
 Process one site of the MPO construction algorithm.
 
@@ -352,7 +399,7 @@ merges adjacent outgoing link sectors with the same QN after component ranks are
 known.
 """
 @timeit function at_site!(
-  ::Type{ValType},
+  ::Type{MatrixType},
   g::MPOGraph{N,C,Ti},
   n::Int,
   sites::Vector{<:Index},
@@ -362,9 +409,7 @@ known.
   alg::String;
   combine_qn_sectors::Bool,
   output_level::Int=0,
-)::Tuple{
-  MPOGraph{N,C,Ti},Vector{Int},Vector{BlockSparseMatrix{ValType}},Index
-} where {ValType<:Number,N,C,Ti}
+)::Tuple{MPOGraph{N,C,Ti},Vector{Int},Vector{BlockSparseMatrix{MatrixType}},Index} where {MatrixType,N,C,Ti}
   has_qns = hasqns(sites)
 
   workspace = combine_duplicate_adjacent_right_vertices!(g, terms_eq_from(n + 1))
@@ -377,7 +422,7 @@ known.
   rank_of_cc = zeros(Int, nccs)
 
   ## The MPO tensor for each component.
-  matrix_of_cc = [BlockSparseMatrix{ValType}() for _ in 1:nccs]
+  matrix_of_cc = [BlockSparseMatrix{MatrixType}() for _ in 1:nccs]
 
   ## The QN of each component
   qi_of_cc = Vector{Pair{QN,Int}}(undef, nccs)
@@ -398,6 +443,7 @@ known.
   )
 
   if alg == "QR"
+    _check_qr_block_storage(MatrixType)
     process_qr(
       matrix_of_cc,
       rank_of_cc,
diff --git a/src/large-graph-mpo-qr.jl b/src/large-graph-mpo-qr.jl
index db2e1ac..edef29a 100644
--- a/src/large-graph-mpo-qr.jl
+++ b/src/large-graph-mpo-qr.jl
@@ -93,7 +93,7 @@ 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}},
+  matrix_of_cc::Vector{BlockSparseMatrix{Matrix{ValType}}},
   rank_of_cc::Vector{Int},
   next_edges_of_cc::Vector{Matrix{Tuple{Vector{Int},Vector{C}}}},
   g::MPOGraph{N,C,Ti},
@@ -162,7 +162,7 @@ On the final site, the scalar `R` factor is applied directly to the local
 blocks.
 """
 @timeit function process_qr(
-  matrix_of_cc::Vector{BlockSparseMatrix{ValType}},
+  matrix_of_cc::Vector{BlockSparseMatrix{Matrix{ValType}}},
   rank_of_cc::Vector{Int},
   next_edges_of_cc::Vector{Matrix{Tuple{Vector{Int},Vector{C}}}},
   g::MPOGraph{N,C,Ti},
diff --git a/src/large-graph-mpo-vertex-cover.jl b/src/large-graph-mpo-vertex-cover.jl
index 3a548af..ffc7654 100644
--- a/src/large-graph-mpo-vertex-cover.jl
+++ b/src/large-graph-mpo-vertex-cover.jl
@@ -1,3 +1,42 @@
+
+function _signed_symbolic_local_op_id(lv::LeftVertex, ::Type{Ti})::Ti where {Ti<:Integer}
+  signed_local_op_id = lv.needs_JW_string ? -Int(lv.op_id) : Int(lv.op_id)
+  return Ti(signed_local_op_id)
+end
+
+function _add_vertex_cover_term!(
+  matrix::BlockSparseMatrix{Matrix{ValType}},
+  m::Int,
+  lv::LeftVertex,
+  weight::Number,
+  op_cache::Vector{OpInfo},
+  site_dim::Int,
+)::Nothing where {ValType}
+  local_op = op_cache[lv.op_id].matrix
+  matrix_element = get!(matrix[m], lv.link) do
+    zeros(ValType, site_dim, site_dim)
+  end
+  add_to_local_matrix!(matrix_element, weight, local_op, lv.needs_JW_string)
+  return nothing
+end
+
+function _add_vertex_cover_term!(
+  matrix::BlockSparseMatrix{SymbolicLocalMatrix{W, Ti}},
+  m::Int,
+  lv::LeftVertex,
+  weight::W,
+  op_cache::Vector{OpInfo},
+  site_dim::Int,
+)::Nothing where {W, Ti}
+  terms = get!(matrix[m], lv.link) do
+    SymbolicLocalMatrix{W, Ti}()
+  end
+  _append_symbolic_local_matrix_term!(
+    terms, weight, _signed_symbolic_local_op_id(lv, Ti)
+  )
+  return nothing
+end
+
 """
     process_vertex_cover!(
       matrix_of_cc, rank_of_cc, next_edges_of_cc, g, ccs, n, sites, op_cache_vec
@@ -12,7 +51,7 @@ stored edge weights. On the final site, only the local tensor blocks are needed;
 otherwise this also builds `next_edges_of_cc` for the graph at site `n + 1`.
 """
 @timeit function process_vertex_cover!(
-  matrix_of_cc::Vector{BlockSparseMatrix{ValType}},
+  matrix_of_cc::Vector{BlockSparseMatrix{MatrixType}},
   rank_of_cc::Vector{Int},
   next_edges_of_cc::Vector{Matrix{Tuple{Vector{Int},Vector{C}}}},
   g::MPOGraph{N,C,Ti},
@@ -20,7 +59,7 @@ otherwise this also builds `next_edges_of_cc` for the graph at site `n + 1`.
   n::Int,
   sites::Vector{<:Index},
   op_cache_vec::OpCacheVec,
-)::Nothing where {ValType<:Number,N,C,Ti}
+)::Nothing where {MatrixType,N,C,Ti}
   site_dim = dim(sites[n])
   op_cache = op_cache_vec[n]
 
@@ -32,7 +71,6 @@ otherwise this also builds `next_edges_of_cc` for the graph at site `n + 1`.
     end
   end
 
-  # TODO: Consider nested multithreading
   Threads.@threads for cc in 1:num_connected_components(ccs)
     lvs_of_component::Vector{Int} = ccs.lvs_of_component[cc]
     position_of_rvs_in_component = ccs.position_of_rvs_in_component
@@ -44,22 +82,18 @@ 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 = [Dictionary{Int,Matrix{ValType}}() for _ in 1:rank]
+    matrix = BlockSparseMatrix(MatrixType, 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]]
       lv = left_vertex(g, lv_id)
-      local_op = op_cache[lv.op_id].matrix
 
-      matrix_element = zeros(ValType, site_dim, site_dim)
-      add_to_local_matrix!(matrix_element, one(ValType), local_op, lv.needs_JW_string)
-      set!(matrix[m], lv.link, matrix_element)
+      _add_vertex_cover_term!(matrix, m, lv, one(C), op_cache, site_dim)
     end
 
     ## Construct the tensor from the right cover.
-    # TODO: Merge this loop with the one above, using `in_left_cover`
     let
       in_left_cover = falses(length(lvs_of_component))
       @inbounds for local_id in left_cover
@@ -82,16 +116,11 @@ otherwise this also builds `next_edges_of_cc` for the graph at site `n + 1`.
 
       @inbounds for lv_id in uncovered_left_ids
         lv = left_vertex(g, lv_id)
-        local_op = op_cache[lv.op_id].matrix
 
         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[m], lv.link) do
-            zeros(ValType, site_dim, site_dim)
-          end
-
-          add_to_local_matrix!(matrix_element, weight, local_op, lv.needs_JW_string)
+          _add_vertex_cover_term!(matrix, m, lv, weight, op_cache, site_dim)
         end
       end
     end
diff --git a/src/sparse_tensor_construction.jl b/src/sparse_tensor_construction.jl
index 566ffae..b9241ea 100644
--- a/src/sparse_tensor_construction.jl
+++ b/src/sparse_tensor_construction.jl
@@ -88,7 +88,7 @@ function _fill_splitblocks!(
   blocks::Vector{Block{4}},
   block_values::Vector{C},
   offsets::Vector{Int},
-  block_sparse_matrices::Vector{BlockSparseMatrix{C}},
+  block_sparse_matrices::Vector{BlockSparseMatrix{Matrix{C}}},
 ) where {C}
   entry = 0
 
@@ -125,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{BlockSparseMatrix{C}},
+  block_sparse_matrices::Vector{BlockSparseMatrix{Matrix{C}}},
   llink::ITensors.QNIndex,
   rlink::ITensors.QNIndex,
   site::ITensors.QNIndex;
@@ -171,7 +171,7 @@ writes nonzero entries into the appropriate block views.
 """
 function to_ITensor(
   offsets::Vector{Int},
-  block_sparse_matrices::Vector{BlockSparseMatrix{C}},
+  block_sparse_matrices::Vector{BlockSparseMatrix{Matrix{C}}},
   llink::ITensors.QNIndex,
   rlink::ITensors.QNIndex,
   site::ITensors.QNIndex;
@@ -200,11 +200,11 @@ function to_ITensor(
       shifted_right_link = right_link + offset
       right_block = right_blocks[shifted_right_link]
 
-      for (left_link, block) in pairs(column)
+      for (left_link, block_matrix) in pairs(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]
+        @inbounds for site_prime_link in axes(block_matrix, 1), site_link in axes(block_matrix, 2)
+          value = block_matrix[site_prime_link, site_link]
           iszero(value) && continue
 
           push!(
@@ -255,12 +255,12 @@ function to_ITensor(
       right_block = right_blocks[shifted_right_link]
       right_offset = right_offsets[shifted_right_link]
 
-      for (left_link, block) in pairs(column)
+      for (left_link, block_matrix) in pairs(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]
+        @inbounds for site_prime_link in axes(block_matrix, 1), site_link in axes(block_matrix, 2)
+          value = block_matrix[site_prime_link, site_link]
           iszero(value) && continue
 
           block_id = _blockid(
@@ -300,7 +300,7 @@ effect for dense indices.
 """
 function to_ITensor(
   offsets::Vector{Int},
-  block_sparse_matrices::Vector{BlockSparseMatrix{C}},
+  block_sparse_matrices::Vector{BlockSparseMatrix{Matrix{C}}},
   llink::Index,
   rlink::Index,
   site::Index;
@@ -311,8 +311,8 @@ function to_ITensor(
     offset = offsets[i]
     matrix = block_sparse_matrices[i]
     for (right_link, column) in enumerate(matrix)
-      for (left_link, block) in pairs(column)
-        tensor[left_link, right_link + offset, :, :] = block
+      for (left_link, block_matrix) in pairs(column)
+        tensor[left_link, right_link + offset, :, :] = block_matrix
       end
     end
   end
diff --git a/test/runtests.jl b/test/runtests.jl
index 3541f69..61c9eb7 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,5 +1,7 @@
 include("tensor-construction-tests.jl")
 include("ops-tests.jl")
+include("symbolic-local-matrix-tests.jl")
+include("symbolic-mpo-tests.jl")
 include("large-graph-tests.jl")
 include("test-MPOConstruction.jl")
 include("mponew_h1_accuracy.jl")
diff --git a/test/symbolic-local-matrix-tests.jl b/test/symbolic-local-matrix-tests.jl
new file mode 100644
index 0000000..e1b6cd2
--- /dev/null
+++ b/test/symbolic-local-matrix-tests.jl
@@ -0,0 +1,293 @@
+using ITensorMPOConstruction
+using ITensorMPOConstruction:
+  MPOGraph,
+  SimpleWeight,
+  SymbolicLocalMatrix,
+  _append_symbolic_local_matrix_term!,
+  _evaluate_symbolic_local_matrix,
+  combine_duplicate_adjacent_right_vertices!,
+  left_size,
+  left_vertex,
+  prepare_opID_sum!,
+  substitute_weight
+using ITensorMPOConstruction: right_size, terms_eq_from
+using ITensors
+using ITensorMPS
+using Test
+
+function test_symbolic_local_matrix_operations()::Nothing
+  coefficients = [10.0, 20.0]
+
+  @test substitute_weight(one(SimpleWeight), coefficients) == 1.0
+  @test substitute_weight(SimpleWeight(1), coefficients) == coefficients[1]
+  @test substitute_weight(SimpleWeight(-1), coefficients) == -coefficients[1]
+  @test substitute_weight(SimpleWeight(2), coefficients) == coefficients[2]
+
+  @test SimpleWeight(3) * 1 == SimpleWeight(3)
+  @test SimpleWeight(3) * -1 == SimpleWeight(-3)
+  @test SimpleWeight(-3) * -1 == SimpleWeight(3)
+  @test_throws ErrorException SimpleWeight(3) * 2
+  @test_throws ErrorException SimpleWeight(3) * im
+
+  terms = SymbolicLocalMatrix{SimpleWeight,Int}()
+  _append_symbolic_local_matrix_term!(terms, SimpleWeight(1), -2)
+  @test terms == [(SimpleWeight(1), -2)]
+
+  qubit_sites = siteinds("Qubit", 1)
+  qubit_op_cache_vec = to_OpCacheVec(qubit_sites, [["I", "X"]])
+  duplicate_terms = SymbolicLocalMatrix{SimpleWeight,Int}([
+    (SimpleWeight(2), 2), (SimpleWeight(2), 2)
+  ])
+  @test _evaluate_symbolic_local_matrix(
+    Float64, duplicate_terms, coefficients, qubit_op_cache_vec[1]
+  ) ≈ 2 * coefficients[2] * qubit_op_cache_vec[1][2].matrix
+
+  fermion_sites = siteinds("Fermion", 1)
+  fermion_op_cache_vec = to_OpCacheVec(fermion_sites, [["I", "C", "Cdag", "N"]])
+  plain_matrix = _evaluate_symbolic_local_matrix(
+    Float64,
+    SymbolicLocalMatrix{SimpleWeight,Int}([(SimpleWeight(2), 2)]),
+    coefficients,
+    fermion_op_cache_vec[1],
+  )
+  jw_matrix = _evaluate_symbolic_local_matrix(
+    Float64,
+    SymbolicLocalMatrix{SimpleWeight,Int}([(SimpleWeight(2), -2)]),
+    coefficients,
+    fermion_op_cache_vec[1],
+  )
+
+  expected_plain_matrix = coefficients[2] * fermion_op_cache_vec[1][2].matrix
+  expected_jw_matrix = coefficients[2] * copy(fermion_op_cache_vec[1][2].matrix)
+  expected_jw_matrix[:, 2] .*= -1
+
+  @test plain_matrix ≈ expected_plain_matrix
+  @test jw_matrix ≈ expected_jw_matrix
+
+  return nothing
+end
+
+function symbolic_basis_rewrite_cache_vecs(factor)
+  identity_matrix = [1.0 0.0; 0.0 1.0]
+  x_matrix = [0.0 1.0; 1.0 0.0]
+
+  op_cache_vec = [[
+    OpInfo("I", identity_matrix, false, QN()),
+    OpInfo("scaled I", factor * identity_matrix, false, QN()),
+    OpInfo("X", x_matrix, false, QN()),
+  ]]
+  basis_op_cache_vec = [[
+    OpInfo("I", identity_matrix, false, QN()), OpInfo("X", x_matrix, false, QN())
+  ]]
+
+  return op_cache_vec, basis_op_cache_vec
+end
+
+function symbolic_rewrite_opID_sum(factor, user_label::Int=7)
+  op_cache_vec, basis_op_cache_vec = symbolic_basis_rewrite_cache_vecs(factor)
+
+  scaled_identity(n) = OpID(2, n)
+  X(n) = OpID(3, n)
+
+  os = OpIDSum{2,SimpleWeight,Int}(1, op_cache_vec)
+  add!(os, SimpleWeight(user_label), scaled_identity(1), X(1))
+
+  return os, basis_op_cache_vec
+end
+
+function test_symbolic_basis_rewrite_sign_factors()::Nothing
+  user_label = 7
+
+  os, basis_op_cache_vec = symbolic_rewrite_opID_sum(-1.0, user_label)
+  prepare_opID_sum!(os, basis_op_cache_vec)
+
+  scalar, ops = os[1]
+  nonzero_ops = [op for op in ops if op != zero(op)]
+  @test scalar == SimpleWeight(-user_label)
+  @test nonzero_ops == [OpID(2, 1)]
+  @test count(op -> op == zero(op), ops) == 1
+  @test os.op_cache_vec === basis_op_cache_vec
+
+  os, basis_op_cache_vec = symbolic_rewrite_opID_sum(1.0, user_label)
+  prepare_opID_sum!(os, basis_op_cache_vec)
+
+  scalar, ops = os[1]
+  nonzero_ops = [op for op in ops if op != zero(op)]
+  @test scalar == SimpleWeight(user_label)
+  @test nonzero_ops == [OpID(2, 1)]
+  @test count(op -> op == zero(op), ops) == 1
+  @test os.op_cache_vec === basis_op_cache_vec
+
+  return nothing
+end
+
+function test_symbolic_basis_rewrite_rejects_unsupported_factors()::Nothing
+  for factor in (im, 2im, 0.5, 2.0)
+    os, basis_op_cache_vec = symbolic_rewrite_opID_sum(factor)
+    err = nothing
+    try
+      prepare_opID_sum!(os, basis_op_cache_vec)
+    catch e
+      err = e
+    end
+
+    @test !isnothing(err)
+    error_message = sprint(showerror, err)
+    @test occursin("SimpleWeight can only be multiplied", error_message)
+  end
+
+  return nothing
+end
+
+function two_site_qubit_symbolic_fixture()
+  sites = siteinds("Qubit", 2)
+  op_cache_vec = to_OpCacheVec(sites, [["I", "X"] for _ in eachindex(sites)])
+  return sites, op_cache_vec
+end
+
+function symbolic_terms_from_first_left_vertex(g)::SymbolicLocalMatrix{SimpleWeight,Int}
+  @test left_size(g) == 1
+  lv = left_vertex(g, 1)
+  signed_local_op_id = lv.needs_JW_string ? -lv.op_id : lv.op_id
+  terms = SymbolicLocalMatrix{SimpleWeight,Int}()
+  for signed_weight_id in g.edge_weights_from_left[1]
+    _append_symbolic_local_matrix_term!(terms, signed_weight_id, Int(signed_local_op_id))
+  end
+  return terms
+end
+
+function compact_duplicate_symbolic_right_vertices!(g)::Nothing
+  combine_duplicate_adjacent_right_vertices!(g, terms_eq_from(2))
+  @test right_size(g) == 1
+  @test g.right_vertex_ids_from_left == [[1 for _ in g.edge_weights_from_left[1]]]
+  return nothing
+end
+
+function test_numeric_mpo_graph_still_compacts_duplicate_terms()::Nothing
+  sites = siteinds("Qubit", 2)
+  op_cache_vec = to_OpCacheVec(sites, [["I", "X", "Z"] for _ in eachindex(sites)])
+
+  X(n) = OpID(2, n)
+  Z(n) = OpID(3, n)
+
+  os = OpIDSum{2,Float64,Int}(4, op_cache_vec; abs_tol=0.25)
+  add!(os, 1.0, X(1), X(2))
+  add!(os, 2.0, X(1), X(2))
+  add!(os, 0.5, Z(1), Z(2))
+  add!(os, -0.4, Z(1), Z(2))
+
+  g = MPOGraph(os)
+
+  @test left_size(g) == 1
+  @test right_size(g) == 1
+  @test left_vertex(g, 1).op_id == 2
+  @test g.right_vertex_ids_from_left == [[1]]
+  @test g.edge_weights_from_left == [[3.0]]
+
+  return nothing
+end
+
+function test_symbolic_mpo_graph_preserves_duplicate_signed_ids()::Nothing
+  _, op_cache_vec = two_site_qubit_symbolic_fixture()
+
+  X(n) = OpID(2, n)
+
+  os = OpIDSum{2,SimpleWeight,Int}(2, op_cache_vec)
+  add!(os, SimpleWeight(1), X(1), X(2))
+  add!(os, SimpleWeight(2), X(1), X(2))
+
+  g = MPOGraph(os)
+  @test length(g.edge_weights_from_left[1]) == 2
+  @test sort(getfield.(g.edge_weights_from_left[1], :id)) == [1, 2]
+
+  compact_duplicate_symbolic_right_vertices!(g)
+  terms = symbolic_terms_from_first_left_vertex(g)
+  @test terms == [(SimpleWeight(1), 2), (SimpleWeight(2), 2)]
+
+  coefficients = [11.0, 17.0]
+  evaluated = _evaluate_symbolic_local_matrix(Float64, terms, coefficients, op_cache_vec[1])
+  @test evaluated ≈ sum(coefficients) * op_cache_vec[1][2].matrix
+
+  return nothing
+end
+
+function test_symbolic_mpo_graph_preserves_duplicate_label_multiplicity()::Nothing
+  _, op_cache_vec = two_site_qubit_symbolic_fixture()
+
+  X(n) = OpID(2, n)
+
+  os = OpIDSum{2,SimpleWeight,Int}(2, op_cache_vec)
+  add!(os, SimpleWeight(1), X(1), X(2))
+  add!(os, SimpleWeight(1), X(1), X(2))
+
+  g = MPOGraph(os)
+  compact_duplicate_symbolic_right_vertices!(g)
+  terms = symbolic_terms_from_first_left_vertex(g)
+  @test terms == [(SimpleWeight(1), 2), (SimpleWeight(1), 2)]
+
+  coefficient = 13.0
+  evaluated = _evaluate_symbolic_local_matrix(Float64, terms, [coefficient], op_cache_vec[1])
+  @test evaluated ≈ 2 * coefficient * op_cache_vec[1][2].matrix
+
+  return nothing
+end
+
+function signed_symbolic_rewrite_cache_vecs()
+  identity_matrix = [1.0 0.0; 0.0 1.0]
+  x_matrix = [0.0 1.0; 1.0 0.0]
+
+  op_cache_vec = [
+    [
+      OpInfo("I", identity_matrix, false, QN()),
+      OpInfo("X", x_matrix, false, QN()),
+      OpInfo("-X", -x_matrix, false, QN()),
+    ],
+    [OpInfo("I", identity_matrix, false, QN()), OpInfo("X", x_matrix, false, QN())],
+  ]
+  basis_op_cache_vec = [
+    [OpInfo("I", identity_matrix, false, QN()), OpInfo("X", x_matrix, false, QN())],
+    [OpInfo("I", identity_matrix, false, QN()), OpInfo("X", x_matrix, false, QN())],
+  ]
+
+  return op_cache_vec, basis_op_cache_vec
+end
+
+function test_symbolic_mpo_graph_cancels_opposite_signed_ids()::Nothing
+  op_cache_vec, basis_op_cache_vec = signed_symbolic_rewrite_cache_vecs()
+
+  X(n) = OpID(2, n)
+  negX(n) = OpID(3, n)
+
+  os = OpIDSum{2,SimpleWeight,Int}(2, op_cache_vec)
+  add!(os, SimpleWeight(1), X(1), X(2))
+  add!(os, SimpleWeight(1), negX(1), X(2))
+  prepare_opID_sum!(os, basis_op_cache_vec)
+
+  g = MPOGraph(os)
+  compact_duplicate_symbolic_right_vertices!(g)
+  terms = symbolic_terms_from_first_left_vertex(g)
+  @test terms == [(SimpleWeight(1), 2), (SimpleWeight(-1), 2)]
+
+  evaluated = _evaluate_symbolic_local_matrix(Float64, terms, [7.0], basis_op_cache_vec[1])
+  @test iszero(evaluated)
+
+  return nothing
+end
+
+@testset "SymbolicLocalMatrix" begin
+  @testset "symbolic local matrix terms" begin
+    test_symbolic_local_matrix_operations()
+  end
+
+  @testset "basis rewrite" begin
+    test_symbolic_basis_rewrite_sign_factors()
+    test_symbolic_basis_rewrite_rejects_unsupported_factors()
+  end
+
+  @testset "symbolic graph duplicate terms" begin
+    test_numeric_mpo_graph_still_compacts_duplicate_terms()
+    test_symbolic_mpo_graph_preserves_duplicate_signed_ids()
+    test_symbolic_mpo_graph_preserves_duplicate_label_multiplicity()
+    test_symbolic_mpo_graph_cancels_opposite_signed_ids()
+  end
+end
diff --git a/test/symbolic-mpo-tests.jl b/test/symbolic-mpo-tests.jl
new file mode 100644
index 0000000..534bd55
--- /dev/null
+++ b/test/symbolic-mpo-tests.jl
@@ -0,0 +1,298 @@
+using ITensorMPOConstruction
+using ITensorMPOConstruction: MPO_symbolic, SimpleWeight, SymbolicMPO
+using ITensors
+using ITensorMPS
+using Test
+
+function two_site_qubit_symbolic_mpo_fixture()
+  sites = siteinds("Qubit", 2)
+  op_cache_vec = to_OpCacheVec(sites, [["I", "X"] for _ in eachindex(sites)])
+  return sites, op_cache_vec
+end
+
+function transverse_field_ising_symbolic_fixture()
+  sites = siteinds("Qubit", 2)
+  op_cache_vec = to_OpCacheVec(sites, [["I", "X", "Z"] for _ in eachindex(sites)])
+
+  X(n) = OpID(2, n)
+  Z(n) = OpID(3, n)
+
+  os = OpIDSum{2,SimpleWeight,Int}(3, op_cache_vec)
+  add!(os, SimpleWeight(1), X(1), X(2))
+  add!(os, SimpleWeight(2), Z(1))
+  add!(os, SimpleWeight(3), Z(2))
+
+  return os, sites
+end
+
+function transverse_field_ising_numeric_opidsum(sites, coefficients)
+  op_cache_vec = to_OpCacheVec(sites, [["I", "X", "Z"] for _ in eachindex(sites)])
+
+  X(n) = OpID(2, n)
+  Z(n) = OpID(3, n)
+
+  os = OpIDSum{2,eltype(coefficients),Int}(3, op_cache_vec)
+  add!(os, coefficients[1], X(1), X(2))
+  add!(os, coefficients[2], Z(1))
+  add!(os, coefficients[3], Z(2))
+
+  return os
+end
+
+function qn_number_symbolic_fixture()
+  sites = siteinds("Fermion", 2; conserve_qns=true)
+  op_cache_vec = to_OpCacheVec(sites, [["I", "N"] for _ in eachindex(sites)])
+
+  Nop(n) = OpID(2, n)
+
+  os = OpIDSum{2,SimpleWeight,Int}(3, op_cache_vec)
+  add!(os, SimpleWeight(1), Nop(1))
+  add!(os, SimpleWeight(2), Nop(2))
+  add!(os, SimpleWeight(3), Nop(1), Nop(2))
+
+  return os, sites
+end
+
+function qn_number_numeric_opidsum(sites, coefficients)
+  op_cache_vec = to_OpCacheVec(sites, [["I", "N"] for _ in eachindex(sites)])
+
+  Nop(n) = OpID(2, n)
+
+  os = OpIDSum{2,eltype(coefficients),Int}(3, op_cache_vec)
+  add!(os, coefficients[1], Nop(1))
+  add!(os, coefficients[2], Nop(2))
+  add!(os, coefficients[3], Nop(1), Nop(2))
+
+  return os
+end
+
+function mpo_relative_error(A::MPO, B::MPO)::Float64
+  AmB = add(A, -B; alg="directsum")
+  numerator = real(inner(AmB, AmB))
+  denominator = real(inner(A, A))
+  return numerator / denominator
+end
+
+function test_mpos_approx_equal(A::MPO, B::MPO; tol::Real=1e-10)::Nothing
+  relative_error = mpo_relative_error(A, B)
+  @test relative_error < tol
+  return nothing
+end
+
+function symbolic_mpo_terms(sym::SymbolicMPO)::Vector
+  terms = []
+  for site_matrices in sym.block_sparse_matrices
+    for matrix_of_component in site_matrices
+      for block in matrix_of_component
+        for block_terms in values(block)
+          append!(terms, block_terms)
+        end
+      end
+    end
+  end
+  return terms
+end
+
+function test_symbolic_mpo_construction_metadata()::Nothing
+  os, sites = transverse_field_ising_symbolic_fixture()
+
+  sym = MPO_symbolic(os, sites)
+
+  @test sym isa SymbolicMPO
+  @test length(sym.offsets) == length(sites)
+  @test length(sym.block_sparse_matrices) == length(sites)
+  @test length(sym.llinks) == length(sites) + 1
+  @test sym.op_cache_vec === os.op_cache_vec
+  @test !isempty(symbolic_mpo_terms(sym))
+
+  return nothing
+end
+
+function test_symbolic_mpo_fresh_instantiation_matches_numeric()::Nothing
+  os, sites = transverse_field_ising_symbolic_fixture()
+  sym = MPO_symbolic(os, sites)
+
+  for coefficients in ([1.25, -0.5, 0.75], [-2.0, 0.25, 1.5])
+    symbolic_mpo = instantiate_MPO(sym, coefficients; splitblocks=true)
+    numeric_mpo = MPO_new(
+      transverse_field_ising_numeric_opidsum(sites, coefficients),
+      sites;
+      alg="VC",
+      splitblocks=true,
+      checkflux=false,
+    )
+
+    test_mpos_approx_equal(symbolic_mpo, numeric_mpo)
+  end
+
+  return nothing
+end
+
+function test_symbolic_mpo_duplicate_labels_instantiate_like_numeric()::Nothing
+  sites, op_cache_vec = two_site_qubit_symbolic_mpo_fixture()
+
+  X(n) = OpID(2, n)
+
+  symbolic_os = OpIDSum{2,SimpleWeight,Int}(2, op_cache_vec)
+  add!(symbolic_os, SimpleWeight(1), X(1), X(2))
+  add!(symbolic_os, SimpleWeight(2), X(1), X(2))
+
+  sym = MPO_symbolic(symbolic_os, sites)
+  coefficients = [2.0, 3.0]
+  symbolic_mpo = instantiate_MPO(sym, coefficients; splitblocks=true)
+
+  numeric_os = OpIDSum{2,Float64,Int}(2, op_cache_vec)
+  add!(numeric_os, coefficients[1], X(1), X(2))
+  add!(numeric_os, coefficients[2], X(1), X(2))
+  numeric_mpo = MPO_new(numeric_os, sites; alg="VC", splitblocks=true, checkflux=false)
+
+  test_mpos_approx_equal(symbolic_mpo, numeric_mpo)
+
+  return nothing
+end
+
+function test_symbolic_mpo_qn_fresh_instantiation()::Nothing
+  os, sites = qn_number_symbolic_fixture()
+  sym = MPO_symbolic(os, sites)
+
+  coefficients = [0.8, -1.1, 0.3]
+  for splitblocks in (false, true)
+    symbolic_mpo = instantiate_MPO(sym, coefficients; splitblocks, checkflux=true)
+    numeric_mpo = MPO_new(
+      qn_number_numeric_opidsum(sites, coefficients),
+      sites;
+      alg="VC",
+      splitblocks,
+      checkflux=true,
+    )
+
+    test_mpos_approx_equal(symbolic_mpo, numeric_mpo)
+  end
+
+  nonzero_mpo = instantiate_MPO(sym, [1.0, 2.0, 3.0]; splitblocks=true)
+  zero_mpo = instantiate_MPO(sym, [0.0, 0.0, 0.0]; splitblocks=true)
+
+  @test linkdims(zero_mpo) == linkdims(nonzero_mpo)
+  @test_throws BoundsError instantiate_MPO(sym, [1.0, 2.0])
+
+  return nothing
+end
+
+function test_symbolic_mpo_missing_coefficients_rejected()::Nothing
+  os, sites = transverse_field_ising_symbolic_fixture()
+  sym = MPO_symbolic(os, sites)
+  coefficients = [1.0, 2.0, 3.0]
+  short_coefficients = [1.0, 2.0]
+
+  @test_throws BoundsError instantiate_MPO(sym, short_coefficients; splitblocks=true)
+
+  return nothing
+end
+
+function test_symbolic_mpo_inplace_instantiation()::Nothing
+  os, sites = qn_number_symbolic_fixture()
+  sym = MPO_symbolic(os, sites)
+
+  coefficients1 = [0.8, -1.1, 0.3]
+  coefficients2 = [-0.2, 1.4, 0.6]
+
+  H = instantiate_MPO(sym, coefficients1; splitblocks=true)
+  original_objectid = objectid(H)
+  original_linkdims = linkdims(H)
+
+  result = instantiate_MPO!(H, sym, coefficients2; checkflux=true)
+  expected = instantiate_MPO(sym, coefficients2; splitblocks=true, checkflux=true)
+
+  @test result === H
+  @test objectid(H) == original_objectid
+  @test linkdims(H) == original_linkdims
+  test_mpos_approx_equal(H, expected)
+
+  coefficients3 = [0.0, 1.4, 0.6]
+  result = instantiate_MPO!(H, sym, coefficients3; checkflux=true)
+  expected = instantiate_MPO(sym, coefficients3; splitblocks=true, checkflux=true)
+
+  @test result === H
+  @test objectid(H) == original_objectid
+  @test linkdims(H) == original_linkdims
+  test_mpos_approx_equal(H, expected)
+
+  @test_throws MethodError instantiate_MPO!(H, sym, coefficients2; splitblocks=false)
+
+  return nothing
+end
+
+function test_symbolic_mpo_template_instantiation()::Nothing
+  os, sites = qn_number_symbolic_fixture()
+  sym = MPO_symbolic(os, sites)
+
+  coefficients1 = [0.8, -1.1, 0.3]
+  coefficients2 = [-0.2, 1.4, 0.6]
+
+  H_template = instantiate_MPO(sym, coefficients1; splitblocks=false)
+  H_template_expected = instantiate_MPO(sym, coefficients1; splitblocks=false)
+  H = instantiate_MPO(H_template, sym, coefficients2; checkflux=true)
+  expected = instantiate_MPO(sym, coefficients2; splitblocks=false, checkflux=true)
+
+  @test H !== H_template
+  @test linkdims(H) == linkdims(H_template)
+  test_mpos_approx_equal(H, expected)
+  test_mpos_approx_equal(H_template, H_template_expected)
+
+  @test_throws MethodError instantiate_MPO(H_template, sym, coefficients2; splitblocks=true)
+
+  return nothing
+end
+
+function test_symbolic_mpo_incompatible_template_rejected()::Nothing
+  os, sites = qn_number_symbolic_fixture()
+  sym = MPO_symbolic(os, sites)
+  coefficients = [0.8, -1.1, 0.3]
+
+  bad_os, bad_sites = qn_number_symbolic_fixture()
+  bad_sym = MPO_symbolic(bad_os, bad_sites)
+  H_bad = instantiate_MPO(bad_sym, coefficients; splitblocks=false)
+  @test_throws ArgumentError instantiate_MPO(H_bad, sym, coefficients)
+
+  return nothing
+end
+
+function test_symbolic_mpo_uses_negative_local_op_id_for_jw_terms()::Nothing
+  sites = siteinds("Fermion", 2)
+  op_cache_vec = to_OpCacheVec(sites, [["I", "C", "Cdag", "N"] for _ in eachindex(sites)])
+
+  C(n) = OpID(2, n)
+  Cdag(n) = OpID(3, n)
+
+  os = OpIDSum{2,SimpleWeight,Int}(1, op_cache_vec)
+  add!(os, SimpleWeight(1), Cdag(1), C(2))
+
+  sym = MPO_symbolic(os, sites)
+  terms = symbolic_mpo_terms(sym)
+  @test any(term -> term[2] < 0, terms)
+
+  return nothing
+end
+
+@testset "SymbolicMPO" begin
+  @testset "fresh instantiation" begin
+    test_symbolic_mpo_construction_metadata()
+    test_symbolic_mpo_fresh_instantiation_matches_numeric()
+    test_symbolic_mpo_duplicate_labels_instantiate_like_numeric()
+    test_symbolic_mpo_qn_fresh_instantiation()
+    test_symbolic_mpo_uses_negative_local_op_id_for_jw_terms()
+  end
+
+  @testset "public API rejections" begin
+    test_symbolic_mpo_missing_coefficients_rejected()
+    test_symbolic_mpo_incompatible_template_rejected()
+  end
+
+  @testset "in-place instantiation" begin
+    test_symbolic_mpo_inplace_instantiation()
+  end
+
+  @testset "template-assisted instantiation" begin
+    test_symbolic_mpo_template_instantiation()
+  end
+end
diff --git a/test/tensor-construction-tests.jl b/test/tensor-construction-tests.jl
index f9d6d5a..16847bd 100644
--- a/test/tensor-construction-tests.jl
+++ b/test/tensor-construction-tests.jl
@@ -5,7 +5,7 @@ using Test
 
 function reference_to_ITensor(
   offsets::Vector{Int},
-  block_sparse_matrices::Vector{ITensorMPOConstruction.BlockSparseMatrix{C}},
+  block_sparse_matrices::Vector{ITensorMPOConstruction.BlockSparseMatrix{Matrix{C}}},
   llink::ITensors.QNIndex,
   rlink::ITensors.QNIndex,
   site::ITensors.QNIndex;