Add symmetric-backend matrix functions, projection, and factorizations#204
Merged
Conversation
…on glue Adds a Hermitian-safe matrix function family to `MatrixAlgebra` (`sqrth_safe`, `invsqrth_safe`, and `sqrth_invsqrth_safe`) that projects onto the Hermitian part and clamps eigenvalues below tolerance in the Moore-Penrose convention, with labeled, permutation, and `Val` tensor forms that work on dense, graded, and `TensorMap` matricizations. Adds a tensor-level `project_hermitian` matching the `MatrixAlgebraKit` function. Factorization outputs now unfuse through `unmatricize_codomain` and `unmatricize_domain`, which backends whose `matricize` does not fuse (TensorKit) override to validate and return the factor unchanged. `project` tolerates trailing length-1 axes on size-strict backends. Structured matrices such as `Diagonal` flow through `bipermutedimsopadd!` via `permuteddims` and the `isstrided` trait. Raises the `MatrixAlgebraKit` compat floor to `0.6`, the first version with `project_hermitian`. Co-authored-by: Claude Fable 5 <noreply@anthropic.com>
`trivialrange(r, n)` returns the n-dimensional identity range of a range family, defaulting to `Base.OneTo(n)`, with the TensorKit extension returning the direct sum of n unit spaces. Downstream packages overload the type-level method the way they already do for the length-1 form.
The AbstractMatrix method rebuilt the diagonal from a map whose clamp branch returned `zero(d)` and whose power branch returned `real(d)^p`, so the first below-tolerance entry of a complex diagonal widened the result to an abstract eltype. Convert both branches with `oftype`, matching what the in-place variant already does through its broadcast assignment.
powh_safe, sqrth_invsqrth_safe, and the gram functions call MatrixAlgebraKit.eigh_full directly instead of projecting onto the Hermitian part internally. Inputs that are Hermitian only up to numerical noise are the caller's responsibility: project with MatrixAlgebraKit.project_hermitian first, as the docstrings now state. Co-authored-by: Claude <noreply@anthropic.com>
pow_diag_safe now has one implementation for every backend: write the clamped powers back through MAK.diagview onto a copy of the input, dropping the separate AbstractMatrix method that rebuilt through MAK.diagonal. The result keeps the input's type and structure, and block-structured diagonal views (graded, TensorMap) come along through their blockwise copyto!. The TensorMap allocate_project aux-space scan tests each candidate slice with tryproject instead of a try/catch around the TensorMap constructor, so covariance is checked through the same verify path the rest of project uses rather than through TensorKit's constructor error types. Co-authored-by: Claude <noreply@anthropic.com>
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #204 +/- ##
==========================================
- Coverage 84.41% 80.20% -4.21%
==========================================
Files 24 23 -1
Lines 770 869 +99
==========================================
+ Hits 650 697 +47
- Misses 120 172 +52
Flags with carried forward coverage won't be shown. Click here to find out more. ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
This was referenced Jul 6, 2026
…overrides The generic factor-reconstruction helpers located a factor's new bond axis with `axes(X, 2)`, which addresses the bond only when `matricize` fuses the codomain into one axis. A non-fusing backend such as a `TensorMap` keeps the group's original legs plus the bond, so `axes(X, 2)` read a physical leg instead, and the TensorKit extension had to override both helpers. The bond is the factor's last axis (codomain factor) or first axis (domain factor) on every backend, so locating it with `axes(X, ndims(X))` makes the helpers backend-agnostic and lets the overrides go. Also trims the `allocate_project` aux-derivation comment. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…x-space inference Drops the `unmatricize_codomain`/`unmatricize_domain` helpers (both the generic definitions and the TensorKit-backend overrides). Each factorization now reconstructs its factors by calling `unmatricize` directly, reading the freshly created bond axis off the returned factor: the factor's last axis on a codomain factor, its first axis on a domain factor, on every backend. This matches how the SVD spectrum factor was already reconstructed inline, so the module uses one primitive throughout. Also splits the TensorKit `allocate_project` into a short dispatcher and an `infer_aux_space` helper. The dispatcher handles the passthrough and rank check, and `infer_aux_space` holds the sector scan that works out the space of a trailing auxiliary leg so the projected operator is symmetry-allowed. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
`pow_diag_safe` passes the whole matrix to `_pow_diag!` for each backend to overload, instead of passing `MAK.diagview(D)` and dispatching on what it returns. A `TensorMap`'s `diagview` is a `SectorVector` for a `DiagonalTensorMap` and a per-sector dict otherwise, which forced the shared code to carry a dict-recursion method. The TensorKit extension now clamps the diagonal per block, where each block's `diagview` is a plain vector view, and the generic path keeps the single-vector case for dense and graded diagonals. A shared helper holds the below-tolerance clamp and PSD-precondition rule for one entry so both backends use the same definition. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
`_pow_diag!` updates the diagonal with `map!(f, σ, σ)` instead of `copyto!(σ, map(f, σ))`, which allocated a temporary vector on every call. The generic path targets a dense, in-place-mappable diagonal view and the TensorKit extension maps over each block's plain vector view. A backend whose diagonal view is not in-place-mappable overloads `_pow_diag!`, so the generic path no longer needs the allocating map-and-copy form. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Renames the two locals in `infer_aux_space` to `aux_dim` (the surplus axis's position) and `aux_length` (its size), and folds a note into the comment that the slice-by-slice `tryproject` probe builds one `TensorMap` per candidate column. A faster derivation would read the covariant sectors from the fusion-tree block structure instead. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
The note in `infer_aux_space` claimed a faster derivation would read sectors "as the `TensorMap` constructor does", but `project_symmetric!` is given the tensor's space and only sorts dense data into its already-known blocks rather than discovering sectors. Rewords the note to say only that reading each column's sector from the fused content's block structure could avoid the per-slice projection. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
The projection verbs reshaped `raw`/`src` to the destination size with a bare `reshape`, which silently reinterprets same-length data of a genuinely different shape (a 2x3 input into a 3x2 destination, say). Adds an internal `check_project_shape` that throws unless the two sizes agree ignoring trailing length-1 axes, mirroring `Base.size(A, d)` for `d > ndims(A)`, and guards the `reshape` in the generic `projectto!` and `is_projected` and in the TensorKit `projectto!` with it. Trailing length-1 axes may still be added or dropped. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
`allocate_project` for `TensorMap` spaces only dispatched on a non-empty codomain, so an empty-codomain projection (a functional over the domain) skipped the auxiliary-leg derivation and fell through to the generic allocation. Splits it into two entries, codomain-led and domain-led, that read the elementary space type from whichever side is non-empty and share an `allocate_project_tensormap` builder, the same two-entry split `similar_map` uses, and threads that space type into `infer_aux_space`. Also renames the shape guard `check_project_shape` to `check_project_size` since it takes sizes. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Eliminate `TensorAlgebra.conjugate` and the `ConjArray` wrapper in favor of `ConjBroadcasted`, a lazy conjugation node that is a `LinearBroadcasted` rather than an `AbstractArray`. `ConjArray` could only wrap an `AbstractArray`, so a `TensorMap` had no lazy conjugation. `ConjBroadcasted` folds into `bipermutedimsopadd!` like the other linear-broadcast nodes and works on any backend, so `conj` now lowers to it and a nested `conj` cancels. The `conjed` name is dropped for now. It will return as a public interface for lazy conjugated contraction such as `X * conjed(X)`, backed by a proper conjugation wrapper.
Split `pow_diag_safe` into a public trio so a backend has a friendly overload point instead of the internal `_pow_diag!`. The in-place `pow_diag_safe!(Dp, D, p, tol)` is the method a backend overloads, with the keyword-tolerance `pow_diag_safe(D, p; atol, rtol)` and a positional-`tol` form on top. The generic kernel maps over `MAK.diagview`, and the TensorKit extension clamps per block.
mtfishman
added a commit
to ITensor/GradedArrays.jl
that referenced
this pull request
Jul 7, 2026
Clamp-power each reduced diagonal block through the new `pow_diag_safe!` hook (ITensor/TensorAlgebra.jl#204). This restores `gram_eigh_full_with_pinv` and the safe spectral functions on a graded matrix, whose diagonal view is stored per sector and cannot go through the generic `map!` path. Only the reduced degeneracy data is touched, which is correct even for a non-abelian sector type: a diagonal factor is `Diagonal(λ) ⊗ I` per sector, so the power passes to the reduced eigenvalues.
mtfishman
added a commit
to ITensor/ITensorBase.jl
that referenced
this pull request
Jul 7, 2026
Define `conj` on an `AbstractNamedTensor` as `conj.(a)`, so eager conjugation reuses the linear-combination broadcast machinery instead of a separate backend hook. `conj.` lowers to a `ConjBroadcasted` leaf (introduced in ITensor/TensorAlgebra.jl#204) that the op primitives absorb, materializing into an all-codomain destination with dualized axes, so a named tensor wrapping a `TensorMap` can be conjugated. The broadcast destination is now allocated from the flattened expression's axes rather than the prototype's, since an axis-dualizing operand makes them differ.
mtfishman
changed the title
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Extends the factorization, projection, and matrix-function layer to work on symmetric backends (graded arrays and TensorKit
TensorMaps), not just dense arrays. This is the foundation the next-generation symmetric ITensor backend builds on.Adds a Hermitian-safe matrix-function family (
sqrth_safe,invsqrth_safe,sqrth_invsqrth_safe,project_hermitian) at both the matrix and the labeled-tensor level, computing the square root and pseudo-inverse square root from a single eigendecomposition with a below-tolerance spectral clamp. The input must be Hermitian, as forMatrixAlgebraKit.eigh_full: callers project withproject_hermitianfirst when it is Hermitian only up to numerical noise. A backend supplies the underlying diagonal power through the publicpow_diag_safe!hook.Reworks the projection interface around three backend hooks (
allocate_project,projectto!,is_projected).projectis now the checked entry point that verifies nothing outside the symmetry-allowed structure was discarded,tryprojectis its nullable sibling, andunchecked_projectis the raw form. A dense array may omit trailing length-1 axes, and a symmetric backend derives the space of one trailing auxiliary axis (an operator's flux-carrying leg) from the data.Adds the TensorKit backend methods behind these interfaces, a length-
nform oftrivialrange, and the array-vocabulary verbs (size,scalar,sum) that a non-AbstractArraybackend overloads. Conjugation now lowers to a lazyConjBroadcastednode that the linear-combination primitives absorb, so a backend getsconjthrough broadcasting rather than a dedicated method.