Corbett/symbolic

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",
     ],

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.

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
 ```
 

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;

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,