Support symmetric backends across the core algorithms

Project.toml

@@ -21,14 +21,15 @@ SimpleGraphConverter = "205b04f2-f585-4877-a239-566270b3f673"
 SplitApplyCombine = "03a91e81-4c3e-53e1-a0a4-9c0c8f19dd66"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+TensorAlgebra = "68bd88dc-f39d-4e12-b2ca-f046b68fcc6a"
 TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2"
 TypeParameterAccessors = "7e5a90cf-f82e-492e-a09b-e3e26432c138"
 
 [compat]
 Adapt = "4.3.0"
 Dictionaries = "0.4"
 Graphs = "1.8.0"
-ITensorBase = "0.10.2"
+ITensorBase = "0.10.9"
 KrylovKit = "0.10.2"
 LinearAlgebra = "1.11.0"
 MatrixAlgebraKit = "0.6.8"
@@ -40,6 +41,7 @@ SimpleGraphConverter = "0.1.0"
 SplitApplyCombine = "1.2.3"
 Statistics = "1.11.1"
 StatsBase = "0.34.4"
+TensorAlgebra = "0.17"
 TensorOperations = "5.2"
 TypeParameterAccessors = "0.3.10, 0.4"
 julia = "1.10"

examples/2dIsing_dynamics_Heisenbergpicture.jl

@@ -1,7 +1,7 @@
 using TensorNetworkQuantumSimulator
 using Graphs: center
 using TensorNetworkQuantumSimulator: setindex_preserve!, noprime
-using ITensors: ITensors, ITensor
+using TensorNetworkQuantumSimulator: ITensors, ITensor
 
 function main()
     nx, ny = 4, 4

examples/boundarymps.jl

@@ -1,4 +1,5 @@
 using TensorNetworkQuantumSimulator
+using TensorNetworkQuantumSimulator: ITensors
 
 using Random
 Random.seed!(1634)

examples/heavyhexIsing_dynamics.jl

@@ -1,5 +1,6 @@
 using TensorNetworkQuantumSimulator
 using Statistics
+using TensorNetworkQuantumSimulator: ITensor, ITensors
 
 function main()
     #Define the lattice

examples/hexagonal_heisenbergmodel_thermalstate.jl

@@ -1,6 +1,6 @@
 using TensorNetworkQuantumSimulator
 using TensorNetworkQuantumSimulator: scalar_factors_quotient, TensorNetworkQuantumSimulator, freenergy
-using ITensors: ITensors, ITensor
+using TensorNetworkQuantumSimulator: ITensors, ITensor
 
 function main()
     χ = 32

examples/loopcorrections.jl

@@ -1,4 +1,6 @@
 using TensorNetworkQuantumSimulator
+using TensorNetworkQuantumSimulator: ITensors
+
 using LinearAlgebra: norm
 
 using Random

src/Apply/apply_gates.jl

@@ -111,11 +111,13 @@ function apply_gate!(
     if length(v⃗) == 2
         v1, v2 = v⃗
         e = NamedEdge(v1 => v2)
-        ind2 = commonind(s_values, first(updated_tensors))
-        δuv = dag(copy(s_values))
-        δuv = replaceind(δuv, ind2, prime(ind2))
-        map_diag!(sign, δuv, δuv)
-        s_values = denseblocks(s_values) * denseblocks(δuv)
+        u = commonind(s_values, first(updated_tensors))
+        v = noncommonind(s_values, first(updated_tensors))
+        # The new messages are the singular values over the bond-and-prime pair.
+        # MatrixAlgebraKit singular values are nonnegative, so the legacy sign fix
+        # (`s * sign(s)` via `map_diag!`) was a no-op and is dropped; fermionic sign
+        # handling for this message construction is future work.
+        s_values = replaceind(s_values, v, prime(u))
         setmessage!(ψ_bpc, e, dag(s_values))
         setmessage!(ψ_bpc, reverse(e), s_values)
     end

src/Apply/full_update.jl

@@ -1,4 +1,5 @@
 using KrylovKit: linsolve
+using .ITensorsITensorBaseCompat: namesetdiff
 
 """
     full_update(o::ITensor, ψ::TensorNetworkState, v⃗; envs, kwargs...)
@@ -37,15 +38,20 @@ function full_update(
         apply_kwargs...,
     )
     if symmetrize
-        singular_values! = Ref{ITensor}()
-        Rᵥ₁, Rᵥ₂, spec = factorize_svd(
-            Rᵥ₁ * Rᵥ₂,
-            inds(Rᵥ₁);
-            ortho = "none",
-            singular_values!,
-            apply_kwargs...,
-        )
-        callback(; singular_values = singular_values![], truncation_error = spec.truncerr)
+        M = Rᵥ₁ * Rᵥ₂
+        codomain = inds(Rᵥ₁)
+        # Balanced SVD: split the singular values symmetrically (√S into each factor).
+        U, S, V = svd_trunc(M, codomain; trunc = itensor_trunc(; apply_kwargs...))
+        u = only(commoninds(U, S))
+        v = only(commoninds(S, V))
+        sqrtS = sqrth_safe(S, (u,), (v,); atol = 0, rtol = 0)
+        Rᵥ₁, Rᵥ₂ = U * replaceind(sqrtS, v, prime(u)), replaceind(sqrtS, u, prime(u)) * V
+        # Best-effort truncation error from norms; suffers catastrophic cancellation when little is
+        # discarded. TODO: expose MatrixAlgebraKit's `ϵ` from `ITensorBase.svd_trunc` and use it here.
+        total = abs2(norm(M))
+        truncation_error = iszero(total) ? zero(real(scalartype(M))) :
+            max(zero(real(scalartype(M))), 1 - abs2(norm(S)) / total)
+        callback(; singular_values = S, truncation_error)
     end
     ψᵥ₁ = Qᵥ₁ * Rᵥ₁
     ψᵥ₂ = Qᵥ₂ * Rᵥ₂
@@ -115,8 +121,8 @@ function optimise_p_q(
 
     fstart = print_fidelity_loss ? fidelity(envs, p_cur, q_cur, p, q, o) : 0
 
-    qs_ind = setdiff(inds(q_cur), collect(Iterators.flatten(inds.(vcat(envs, p_cur)))))
-    ps_ind = setdiff(inds(p_cur), collect(Iterators.flatten(inds.(vcat(envs, q_cur)))))
+    qs_ind = namesetdiff(inds(q_cur), collect(Iterators.flatten(inds.(vcat(envs, p_cur)))))
+    ps_ind = namesetdiff(inds(p_cur), collect(Iterators.flatten(inds.(vcat(envs, q_cur)))))
 
     function b(p::ITensor, q::ITensor, o::ITensor, envs::Vector{ITensor}, r::ITensor)
         ts = vcat(ITensor[p, q, o, dag(prime(r))], envs)

src/Apply/gate_definitions.jl

@@ -263,8 +263,13 @@ ITensors.op(o::OpName"xx_plus_yy", ::SiteType"Qubit"; θ::Number, β::Number) =
 
 Gate for rotation by XXYY at a given angle.
 """
-function ITensors.op(::OpName"Rxxyy", ::SiteType"S=1/2", s1::Index, s2::Index; θ = 1)
-    h = 0.5 * (op("X", s1) * op("X", s2) + op("Y", s1) * op("Y", s2))
+function ITensors.op(::OpName"Rxxyy", ::SiteType"S=1/2"; θ = 1)
+    # Built as one two-site matrix in the manifestly charge-conserving σ± form,
+    # ½(XX + YY) = σ⁺σ⁻ + σ⁻σ⁺, rather than from single-site `op("X", s)` factors:
+    # the gate conserves U(1) charge, but a standalone `X` does not, so the factored
+    # construction has no symmetric representation even though the sum does.
+    σp, σm = [0.0 1.0; 0.0 0.0], [0.0 0.0; 1.0 0.0]
+    h = kron(σp, σm) + kron(σm, σp)
     return exp(-im * θ * h)
 end
 ITensors.op(o::OpName"Rxxyy", ::SiteType"Qubit"; θ::Number) =
@@ -275,8 +280,11 @@ ITensors.op(o::OpName"Rxxyy", ::SiteType"Qubit"; θ::Number) =
 
 Gate for rotation by XXYYZZ at a given angle.
 """
-function ITensors.op(::OpName"Rxxyyzz", ::SiteType"S=1/2", s1::Index, s2::Index; θ = 1)
-    h = 0.5 * (op("X", s1) * op("X", s2) + op("Y", s1) * op("Y", s2) + op("Z", s1) * op("Z", s2))
+function ITensors.op(::OpName"Rxxyyzz", ::SiteType"S=1/2"; θ = 1)
+    # One two-site matrix in the σ± form, for the same reason as `Rxxyy` above:
+    # ½(XX + YY + ZZ) = σ⁺σ⁻ + σ⁻σ⁺ + ½ ZZ.
+    σp, σm, σz = [0.0 1.0; 0.0 0.0], [0.0 0.0; 1.0 0.0], [1.0 0.0; 0.0 -1.0]
+    h = kron(σp, σm) + kron(σm, σp) + 0.5 * kron(σz, σz)
     return exp(-im * θ * h)
 end
 ITensors.op(o::OpName"Rxxyyzz", ::SiteType"Qubit"; θ::Number) =

src/Apply/simple_update.jl

@@ -35,8 +35,14 @@ function simple_update(
         envs_v2 = filter(env -> hascommoninds(env, ψ⃗[2]), envs)
         @assert all(ndims(env) == 2 for env in vcat(envs_v1, envs_v2))
 
-        sqrt_inv_sqrt_envs_v1 = pseudo_sqrt_inv_sqrt.(envs_v1; cutoff = sqrt_cutoff)
-        sqrt_inv_sqrt_envs_v2 = pseudo_sqrt_inv_sqrt.(envs_v2; cutoff = sqrt_cutoff)
+        # The environments are hermitian only up to numerical noise, so project before
+        # the square roots (which require hermitian input).
+        sqrt_invsqrt = env -> sqrth_invsqrth_safe(
+            project_hermitian(env, (inds(env)[1],), (inds(env)[2],)),
+            (inds(env)[1],), (inds(env)[2],); atol = sqrt_cutoff, rtol = 0
+        )
+        sqrt_inv_sqrt_envs_v1 = map(sqrt_invsqrt, envs_v1)
+        sqrt_inv_sqrt_envs_v2 = map(sqrt_invsqrt, envs_v2)
         sqrt_envs_v1, inv_sqrt_envs_v1 = first.(sqrt_inv_sqrt_envs_v1), last.(sqrt_inv_sqrt_envs_v1)
         sqrt_envs_v2, inv_sqrt_envs_v2 = first.(sqrt_inv_sqrt_envs_v2), last.(sqrt_inv_sqrt_envs_v2)
 
@@ -49,16 +55,22 @@ function simple_update(
         rᵥ₁ = commoninds(Qᵥ₁, Rᵥ₁)
         rᵥ₂ = commoninds(Qᵥ₂, Rᵥ₂)
         oR = ITensors.apply(o, Rᵥ₁ * Rᵥ₂)
-        singular_values! = Ref{ITensor}()
-        Rᵥ₁, Rᵥ₂, spec = factorize_svd(
-            oR,
-            unioninds(rᵥ₁, sᵥ₁);
-            ortho = "none",
-            singular_values!,
-            apply_kwargs...,
-        )
-        err = spec.truncerr
-        s_values = singular_values![]
+        # Balanced SVD: split the singular values symmetrically (√S into each factor) so neither
+        # side is isometric. The bond stays on `prime(u)` (keeping `u`'s name), so once this
+        # function `noprime`s its result the bond becomes `u`, which the returned `s_values` (over
+        # `(u, v)`) still shares for `apply_gate!`'s bond-message construction.
+        U, S, V = svd_trunc(oR, unioninds(rᵥ₁, sᵥ₁); trunc = itensor_trunc(; apply_kwargs...))
+        u = only(commoninds(U, S))
+        v = only(commoninds(S, V))
+        sqrtS = sqrth_safe(S, (u,), (v,); atol = 0, rtol = 0)
+        Rᵥ₁, Rᵥ₂ = U * replaceind(sqrtS, v, prime(u)), replaceind(sqrtS, u, prime(u)) * V
+        s_values = S
+        # Best-effort truncation error from norms (SVD preserves the Frobenius norm); suffers
+        # catastrophic cancellation when little is discarded. TODO: expose MatrixAlgebraKit's `ϵ`
+        # from `ITensorBase.svd_trunc` and use it here instead.
+        total = abs2(norm(oR))
+        err = iszero(total) ? zero(real(scalartype(oR))) :
+            max(zero(real(scalartype(oR))), 1 - abs2(norm(S)) / total)
         Qᵥ₁ = contract([Qᵥ₁; dag.(inv_sqrt_envs_v1)])
         Qᵥ₂ = contract([Qᵥ₂; dag.(inv_sqrt_envs_v2)])
         updated_tensors = [Qᵥ₁ * Rᵥ₁, Qᵥ₂ * Rᵥ₂]

src/Forms/abstractform.jl

@@ -22,7 +22,9 @@ function virtualinds(form::AbstractForm, edge::NamedEdge)
 end
 
 function default_message(form::AbstractForm, edge::AbstractEdge)
-    return adapt_like(form, denseblocks(delta(virtualinds(form, edge))))
+    cod = virtualinds(ket(form), edge)
+    dom = dag.(bra_virtualinds(form, edge))
+    return one(similar_map(ket(form)[src(edge)], cod, dom), cod, dom)
 end
 
 function bp_factors(form::AbstractForm, verts::Vector)

src/Forms/bilinearform.jl

@@ -14,12 +14,15 @@ Base.copy(blf::BilinearForm) = BilinearForm(copy(blf.ket), copy(blf.operator), c
 
 #Constructor, bra is taken to be in the vector space of ket so the dual is taken
 function BilinearForm(ket::TensorNetworkState, bra::TensorNetworkState)
-    dtype = datatype(ket)
     @assert graph(ket) == graph(bra)
-    bra = map_tensors(t -> dag(prime(t)), bra)
     sinds = siteinds(ket)
     verts = collect(vertices(ket))
-    operator_tensors = [adapt(dtype)(reduce(*, ITensor[denseblocks(delta(sind, prime(dag(sind)))) for sind in sinds[v]])) for v in verts]
+    bra = TensorNetworkState(Dictionary(verts, [bra_tensor(bra, v) for v in verts]))
+    operator_tensors = [
+        let codomain = dag.(sinds[v]), domain = dag.(prime.(sinds[v]))
+            one(similar_map(ket[v], codomain, domain), codomain, domain)
+        end for v in verts
+    ]
     operator = TensorNetworkState(Dictionary(verts, operator_tensors))
     return BilinearForm(ket, operator, bra)
 end

src/Forms/quadraticform.jl

@@ -6,7 +6,7 @@ end
 ket(qf::QuadraticForm) = qf.ket
 operator(qf::QuadraticForm) = qf.operator
 bra(qf::QuadraticForm) = prime(dag(ket(qf)))
-bra_tensor(qf::QuadraticForm, v) = dag(prime(ket(qf)[v]))
+bra_tensor(qf::QuadraticForm, v) = bra_tensor(ket(qf), v)
 bra_virtualinds(qf::QuadraticForm, edge::NamedEdge) = dag.(prime.(virtualinds(ket(qf), edge)))
 
 Base.copy(qf::QuadraticForm) = QuadraticForm(copy(qf.ket), copy(qf.operator))

src/ITensorsITensorBaseCompat/ITensorsITensorBaseCompat.jl

@@ -26,15 +26,14 @@ include("ops.jl")
 
 export
     # Index access and set algebra
-    inds, commoninds, commonind, uniqueinds, noncommonind, noncommoninds, unioninds, hascommoninds,
+    inds, commoninds, commonind, uniqueinds, noncommonind, noncommoninds, unioninds,
+    hascommoninds,
     # Index operations
     sim, dag, prime, noprime, replaceind, replaceinds, dim, swapind,
     # ITensor construction
-    itensor, random_itensor, scalar, delta, onehot, combiner, combinedind,
+    itensor, random_itensor, scalar, delta, onehot,
     # Factorizations
-    qr, svd, eigen, factorize, factorize_svd,
-    # Diagonal manipulation
-    map_diag, map_diag!,
+    qr, svd, svd_trunc, eigen, factorize, itensor_trunc,
     # Storage / element-type accessors
     scalartype, datatype, array, data,
     # Dense / quantum-number no-ops