Support symmetric backends across the core algorithms#1
Conversation
…e compat eigen Simple update gauges environments with `sqrth_invsqrth_safe` and `renyi_entropy` uses `sqrth_safe`, replacing the eigendecomposition helpers vendored in `utils.jl`. The compat `eigen` projects onto the Hermitian part with `MatrixAlgebraKit.project_hermitian` at the matricized level instead of a named-level add. The compat factorizations bipartition indices by name so graded bond duality does not drop legs. `Rxxyy` and `Rxxyyzz` are built as whole matrices from raising and lowering operators so they project onto symmetric backends, and the vendored gate set gains the `SWAP` matrix. Removes the singular value sign fix after gate application (MatrixAlgebraKit singular values are nonnegative, so it was a no-op). Fermionic sign handling is left for future work. Dense `expect` on networks with pruned trivial tensors is fixed by pruning at the network level. Adds `TensorAlgebra` as a direct dependency. Co-authored-by: Claude Fable 5 <noreply@anthropic.com>
The order-2 compat `delta` is now built with a checked `project` and carries exactly the indices it is given, matching the legacy convention where callers set the arrow directions. Dangling non-physical legs (for example the charge leg of a charged product state) pair directly between ket and bra, handled by a shared `bra_tensor` helper used by `norm_factors`, the forms, and the samplers. The compat `onehot` carries its charge on a derived auxiliary leg when needed, the compat `array` densifies any backend storage, and the boundary MPS message stitching mints its virtual legs as trivial sectors of the messages' backend.
…mple tests The loop-correction weights now match bond indices by name (a graded tensor carries the dual of its message's axis, so full Index comparisons never match), build the antiprojector identity from the actual message legs, and record dangling aux legs before the bond relabeling so the bra construction does not misread the relabeled bond as a charge leg. The compat layer gains the legacy `itensor` reshape semantics and a shared charged-state construction for state vectors, used by `onehot` and the product-state constructor. The example scripts run on the new backend and `test_examples.jl` is back in the test suite.
…ption The gauge transformation now works directly on the BP cache: per edge, take Hermitian square roots and pseudo-inverse roots of the two messages with sqrth_invsqrth_safe, absorb them into the endpoint tensors by renaming the state's bond indices onto the roots' domain names, and SVD the bond matrix formed by contracting the two roots' codomain legs. The factorization directions and names line up so no conjugation or arrow flips are needed on any backend (dense, graded, TensorKit). Small eigenvalues are clamped to zero through the pseudo-inverse instead of the legacy regularization shift. The safe spectral functions require Hermitian input, so this and the other square-root callsites (simple_update, entanglement) project with project_hermitian first. The compat sim now rebuilds the index name keeping tags and prime level, since uniquename mints bare names. Co-authored-by: Claude <noreply@anthropic.com>
…port is in review Co-authored-by: Claude Fable 5 <noreply@anthropic.com>
`sim` is now `uniquename` on the index, which keeps the tags and prime level. `onehot` and the state-vector constructor share a `project_aux` helper that carries a charge on a derived auxiliary leg only when the basis vector needs it. The order-2 identity map is spelled `id(codomain, domain)`, following TensorKit, and the direct identity-map callers (`tr`, `normalize_rdm`, `bilinearform`) use it, while `delta` keeps its own order-2 method for the count-dispatched copy-tensor cases. Factorization and message bond names mint bare `IndexName`s via the name type. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
The order-2 identity constructions (`tr`, `normalize_rdm`, `bilinearform`) now build their identity maps through ITensorBase's device-following `Base.one` instead of a CPU-materialized dense identity routed through `project`, so the result follows the input's device (fixing a GPU regression) and stays correct on graded backends. `tr` and `normalize_rdm` build the identity in a single `one(t, codomain, prime.(codomain))` over all prime pairs at once. The domain is `prime.(codomain)` rather than `inds(t; plev=1)`: `plev` filtering does not preserve the pairing order between the unprimed and primed groups, so pairing by prime keeps the trace independent of storage order. `bilinearform`'s ket tensor carries each site index but not its prime, so the tensor itself can't serve as the identity's shape prototype. A new `similar_map` helper allocates an identity-shaped tensor over the given codomain and domain following the prototype's backend and element type, which `one` then fills. Also drops the now-unused `id` (the direct identity-map callers were its only users) and the `map_diag`/`map_diag!` diagonal-map shims (no remaining callsites, and the spectral matrix functions cover those uses).
mtfishman
force-pushed
the
mf/itensorbase-backend
branch
from
1eeae5f to
6a04b9d
Compare
The example loaded the legacy `ITensors` package for `op` and `ITensor`. Bind them to the compatibility layer through `TensorNetworkQuantumSimulator` instead, matching the thermal-state example, so no example pulls in legacy `ITensors`.
The trace identity (`tr`, `normalize_rdm`) and the bilinear-form operator now build over the duals of the contracted legs, `dag.(codomain)` and `dag.(prime.(codomain))`, through `similar_map`. The device-following identity migration had dropped the `dag` that the previous construction carried. Dense backends ignore arrows, so that migration's tests passed, but a graded contraction requires arrow-opposite axes, so `reduced_density_matrix`, `tr`, and `inner` failed with an axis mismatch. The fix stays device-following and contracts on both dense and graded backends.
`vertices(::ITensor, ::TensorNetworkState)` resolved a gate's vertices with a full-`Index` `intersect`, which is arrow-sensitive on a graded axis: a gate's contraction leg is the dual of the site index, so the match came up empty and `apply_gates` hit a `BoundsError` in `simple_update`. The dense backend passed only because it ignores arrows. It now matches by `IndexName` through `nameisdisjoint`. This also reworks the compat index-set algebra into a generic small-collection layer keyed by a `by` transform (`smallintersect` and friends, a linear scan with no hashing), name-keyed wrappers built on it with `by = name`, and the `commoninds`/etc. named-tensor family delegating to those.
`op("I")` and the form BP `default_message` built their order-2 identity maps through `delta`/`project`. They now use the identity-map primitives instead: `op("I")` uses the from-scratch `id(eltype, codomain, domain)` (only the site index is in hand, no prototype tensor to follow), and `default_message(form)` uses the device-following `one(similar_map(...), ...)` pattern over the ket bonds (codomain) and the duals of the bra bonds (domain), which generalizes to any number of ket-to-bra bond pairs. With no order-2 identity caller left on `delta`, the `delta(eltype, ::Tuple{Index,Index})` method is deleted. The remaining `delta` callers are dense copy tensors that fall to the general `diagonaltensor` path.
lkdvos
left a comment
lkdvos
left a comment
There was a problem hiding this comment.
Left some very small comments, but overall looks good to me!
One overall comment - in the long run it might be useful to try and convert a lot of the prose that is currently in comments into docstrings and documentation, but obviously that is easy to handle afterwards.
| # (not `inds(t; plev=1)`) so each codomain index pairs with its own prime; `plev` filtering | ||
| # does not preserve that pairing order. Accessed qualified (`ITensors.tr`) so it doesn't shadow | ||
| # `LinearAlgebra.tr`, which TNQS still calls on plain matrices. | ||
| function tr(t::AbstractITensor) |
There was a problem hiding this comment.
note: this effectively makes this definition the supertrace instead of the trace for fermionic tensors, i.e. different from the sum of the diagonal elements.
There was a problem hiding this comment.
That's good to know, would be good to think about which one we really want here/what a good interface for distinguishing those two would be. Maybe as a followup when we work on fermions...
| # Project a raw vector as a state `ITensor` over `i`, adding an auxiliary leg only when the | ||
| # vector cannot live in the flux-zero space over `i` alone. The index axis selects the | ||
| # backend (dense, graded, `TensorMap`). Shared by the state constructors and `onehot`. | ||
| function project_aux(v::AbstractVector{<:Number}, i::Index) |
There was a problem hiding this comment.
I might be missing something but it is not entirely obvious to me how this works, is this silently assuming that there is only a single charge in the index? if not, I don't think I immediately understand how that would work (and I'm not really seeing it for TensorKit), although I would be interested to hear what the idea is here
There was a problem hiding this comment.
This is probably the most subtle part and something we should discuss more to decide if the design and implementation make sense. In TensorAlgebra as part of the project functionality there is a mode where it automatically tries to figure out the proper auxiliary space to cancel out the flux. For GradedArrays that's relatively easy (just search for a nonzero value and figure out which block/sector it is in), for TensorMap I did it in a pretty brute force way that probably needs to be improved (it loops over possible auxiliary sectors and tries to project into that sector). So basically what happens is if you pass in an array you want to project with one more dimension than the number of indices you want to project into, it will try to automatically infer the sector of the auxiliary space in that way. If the array shape matches the desired space it just does a straight projection into the specified space (i.e. assuming it is invariant).
There was a problem hiding this comment.
Let's discuss this tomorrow - there is similar functionality in TensorKit (axes(V, sector) gives you the range where a given sector contributes), but using that in anything non-abelian is tricky :)
There was a problem hiding this comment.
Ah, I see, that sounds helpful. Agreed, I was just realizing I wasn't sure if the current code works for non-abelian.
…and full_update qr, auxinds, and full_update computed index complements with Base setdiff, which builds a Set of Indexes. On a hash collision Julia falls back to isequal(::Index, ::Index), which iterates the index space, and a TensorKit GradedSpace is not iterable, so these threw on symmetric backends once a collision occurred (small index collections happened to work, larger ones did not). They now use the name-based namesetdiff, matching the rest of the compat name-keyed index-set algebra, which never hashes or iterates an Index.
Agreed, would be good to do as a followup. |
Also to clarify, the goal is to eventually (or soon, in some followup PRs) get rid of ITensorsITensorBaseCompat entirely, either by updating library code or moving functionality to TensorAlgebra/ITensorBase. So a lot of that code should be considered slightly hacky temporary workarounds, I'll keep working on that in followups. |
Replaces the multi-purpose compat `factorize_svd` with a direct balanced SVD inlined at the two apply call sites. `simple_update` and `full_update` now call `svd_trunc` and absorb the square root of the singular values into both factors, reading off the singular values and truncation error inline. That is all `factorize_svd` did for its only two callers, which both used the balanced (non-orthogonal) split. Replaces `_trunc_spec` with a plain `itensor_trunc(; maxdim, cutoff)` that builds a MatrixAlgebraKit truncation strategy (`maxdim` caps the rank, `cutoff` maps to a relative 2-norm error) and returns `notrunc()` when neither applies. `generic_apply`'s no-cap default becomes `nothing` rather than `typemax(Int)`, so no sentinel handling leaks into the utility, and the orthogonalizing `factorize` wrapper uses `itensor_trunc` too. Folds in a few smaller compat cleanups: `normalize_rdm` goes through `ITensors.tr`, the reshaping `itensor` constructor throws `DimensionMismatch` on a length mismatch, and the internal hermitian eigen helper is renamed and relies on `eigh_full`'s own hermitian tolerance instead of pre-projecting onto the hermitian part. Also trims over-long and outdated comments in the compat layer. The truncation error is still computed from norms, which loses precision when little weight is discarded. Surfacing MatrixAlgebraKit's own truncation error through `ITensorBase.svd_trunc` is left for a follow-up, noted at both call sites.
|
Thanks for the review @lkdvos, I think all comments are addressed (or deliberately left for future work). I'll go ahead and merge into JoeyT1994#78 so things are more consolidated. |
ITensorBase 0.10.9 and TensorAlgebra 0.17.0 are registered, so resolve against the released versions instead of the branch pins. The [compat] floors already require these versions.
mtfishman
changed the title
2 participants
Summary
Makes the core algorithms run on symmetric backends (graded arrays and TensorKit
TensorMaps), not just dense arrays, on top of the ITensorBase backend swap. Gauging, entanglement,inner,rdm,expect, sampling, loop corrections, and gate application all work when the site and bond spaces carry symmetry sectors.The bulk of the work is making the shared machinery symmetry-safe: factorizations and index-set operations match indices by name so graded bond duality never drops a leg, the gauge and square-root steps go through Hermitian-safe matrix functions, and state and gate construction carry their charges on the appropriate sectors.
symmetric_gaugeis rewritten to act directly on the BP cache.Fermionic sign handling is left for follow-up work.