Reference

TensorAlgebra.AbstractAlgorithm

Supertype for the algorithm objects operations dispatch on. An operation's own supertype subtypes this (for example ContractAlgorithm), which is what makes an instance pass through TensorAlgebra.select_algorithm unchanged.

ConjBroadcasted(parent)

Lazy conjugate node in the linear-combination broadcast fold, the LinearBroadcasted counterpart of ScaledBroadcasted/AddBroadcasted. Holds parent (any array-like operand or backend tensor, not necessarily an AbstractArray) and presents its axes conjugated (dualized); the op primitives absorb it by unwrapping to the parent and folding conj into their op (see the bipermutedimsopadd! method below). Produced internally by the conj lowering of a broadcast (linearbroadcasted(conj, a)), not a user-facing lazy-conjugate wrapper. Because it is not an AbstractArray, it works for non-array backends (e.g. a TensorMap) as well as dense arrays.

LinearBroadcasted

Abstract supertype for lazy linear broadcast expressions. Analogous to Base.Broadcast.Broadcasted but restricted to linear operations.

Materializes via the protocol: copy(lb) = copyto!(similar(lb), lb) copyto!(dest, lb) → scaleadd!(dest, lb, 1, 0)

PermutedDims(parent, perm)

Lazy permuted-dims wrapper storing parent and the permutation perm in fields (unlike Base.PermutedDimsArray, which encodes perm in a type parameter), so it constructs cheaply from a runtime permutation. Primarily for internal use to track permutations in linear broadcasting.

TensorOperationsContract(; backend = nothing, allocator = nothing)

Contract using TensorOperations, with backend selecting the contraction kernel and allocator the allocator for temporary tensors (e.g. TensorOperations.ManualAllocator()). A nothing field uses TensorOperations' default. Only usable with TensorOperations loaded.

add!(dest, src)

dest .+= src.

allocate_project(raw, axes_codomain, axes_domain) -> dest

Allocate the destination that projecting raw onto axes_codomain/axes_domain fills. This is a backend customization point (with projectto! and is_projected); the default is plain similar_map(raw, axes_codomain, axes_domain).

project projects into exactly the given axes, so raw must not have more axes than they account for. To append a derived flux-carrying auxiliary axis for a charge-shifting operator or a non-invariant state, use project_aux instead.

bipartition(t::Tuple, length1::Val) -> (t1, t2)
bipartition(t::Tuple, group1::Tuple, group2::Tuple) -> (p1, p2)

Split a flat tuple into two groups, returned as a pair of tuples.

The first form splits t in order, taking the first length1 entries as t1 and the remaining entries as t2. The second form gathers the entries of t at the two index groups group1 and group2, returning p1 = t[group1...] and p2 = t[group2...].

biperm(t, t1, t2) -> (p1, p2)

Locate the groups t1 and t2 within t, returning the positions of t1 as p1 and the positions of t2 as p2. The groups t1 and t2 must partition t, so the concatenation (p1..., p2...) is a permutation of eachindex(t) and the pair (p1, p2) is a bipartitioned permutation (a "biperm") splitting t into a codomain p1 and a domain p2.

bipermutedimsopadd!(dest, op, src, perm_codomain, perm_domain, α, β)

dest = β * dest + α * permutedims(op.(src), (perm_codomain..., perm_domain...)).

This is the primary overload point for downstream array types that want to implement op-aware bipartitioned permutation + accumulation (e.g., fuse conj into the copy, or use lazy wrappers like StridedView with op metadata).

The op is the conjugation flag expressed as a function — identity or conj, analogous to TensorOperations' boolean conjA/conjB. On graded axes conj dualizes; on dense axes it is a no-op. Transposition/permutation is carried by the perm arguments, not by op.

The default implementation flattens the bipartitioned permutation, applies op element-wise, permutes, then accumulates via broadcasting with Strided.jl optimization when possible.

contract(a1, labels1, a2, labels2; alg = nothing) -> a_dest, labels_dest

Contract the arrays over the labels they share, returning the result along with the labels of its dimensions. A label appearing on two operands is summed over, one appearing on a single operand survives, and labels_dest reports the surviving labels in the order the result carries them.

julia> using TensorAlgebra: contract

julia> a, b = randn(2, 3), randn(3, 4);

julia> ab, labels = contract(a, (:i, :j), b, (:j, :k));

julia> labels
2-element Vector{Symbol}:
 :i
 :k

julia> ab ≈ a * b
true

See also contractalign to name the output dimensions and get back just the array.

contractalign(labels_dest, a1, labels1, a2, labels2; alg = nothing) -> a_dest

Contract the input arrays over the shared labels, aligning the output array according to the specified destination labels labels_dest. labels_dest must match the uncontracted labels, i.e. issetequal(labels_dest, symdiff(labels1, labels2)) must be true.

julia> using TensorAlgebra: contractalign

julia> a, b = randn(2, 3), randn(3, 4);

julia> contractalign((:k, :i), a, (:i, :j), b, (:j, :k)) ≈ permutedims(a * b, (2, 1))
true

data(a)

The underlying storage of a: the value reached by following parent to its fixed point (an object that is its own parent). A wrapper returns the storage it ultimately wraps, and a plain array returns itself.

Examples

julia> import TensorAlgebra

julia> a = [1.0 2.0; 3.0 4.0];

julia> TensorAlgebra.data(transpose(a)) === a
true

julia> TensorAlgebra.data(a) === a
true

datatype(a) -> Type

The type of the underlying storage of a, i.e. typeof(data(a)), in contrast to scalartype/eltype, which give its element type alone.

Examples

julia> import TensorAlgebra

julia> TensorAlgebra.datatype([1.0 2.0; 3.0 4.0])
Matrix{Float64} (alias for Array{Float64, 2})

julia> TensorAlgebra.datatype(transpose([1.0 2.0; 3.0 4.0]))
Matrix{Float64} (alias for Array{Float64, 2})

TensorAlgebra.default_algorithm(f, args::Tuple)
TensorAlgebra.default_algorithm(f, Args::Type{<:Tuple})

The algorithm operation f runs with on args when the caller names none. The types form is the registration point for a storage type, and the values form defaults to it.

A storage type registers its choice per operation, so a backend that contracts its own way adds a method to default_algorithm(contract!, ::Type{<:Tuple{A_dest, A1, A2}}).

dual(a)

Returns the dual of the axis a. Falls back to returning a unchanged for AbstractUnitRange.

See also isdual.

eig_full(A, labels_A, labels_codomain, labels_domain; kwargs...) -> D, V
eig_full(A, perm_codomain, perm_domain; kwargs...) -> D, V
eig_full(A, ndims_codomain::Val; kwargs...) -> D, V

Compute the eigenvalue decomposition of a generic N-dimensional array interpreted as a general (non-Hermitian) linear map from the domain to the codomain dimensions. The output eltype is always <:Complex. The partition is specified either via labels or directly through a bi-permutation.

See also MatrixAlgebraKit.eig_full!.

eig_trunc(A, labels_A, labels_codomain, labels_domain; trunc, kwargs...) -> D, V
eig_trunc(A, perm_codomain, perm_domain; trunc, kwargs...) -> D, V
eig_trunc(A, ndims_codomain::Val; trunc, kwargs...) -> D, V

Truncated general eigenvalue decomposition, like eig_full but keeping only the eigenvalues selected by the trunc strategy.

See also MatrixAlgebraKit.eig_trunc!.

eig_vals(A, labels_A, labels_codomain, labels_domain; kwargs...) -> D
eig_vals(A, perm_codomain, perm_domain; kwargs...) -> D
eig_vals(A, ndims_codomain::Val; kwargs...) -> D

Compute the eigenvalues of a generic N-dimensional array interpreted as a general (non-Hermitian) linear map from the domain to the codomain dimensions. The output is a vector of eigenvalues with <:Complex eltype.

See also MatrixAlgebraKit.eig_vals!.

eigh_full(A, labels_A, labels_codomain, labels_domain; kwargs...) -> D, V
eigh_full(A, perm_codomain, perm_domain; kwargs...) -> D, V
eigh_full(A, ndims_codomain::Val; kwargs...) -> D, V

Compute the eigenvalue decomposition of a generic N-dimensional array interpreted as a Hermitian linear map from the domain to the codomain dimensions. The partition is specified either via labels or directly through a bi-permutation.

See also MatrixAlgebraKit.eigh_full!.

eigh_trunc(A, labels_A, labels_codomain, labels_domain; trunc, kwargs...) -> D, V
eigh_trunc(A, perm_codomain, perm_domain; trunc, kwargs...) -> D, V
eigh_trunc(A, ndims_codomain::Val; trunc, kwargs...) -> D, V

Truncated Hermitian eigenvalue decomposition, like eigh_full but keeping only the eigenvalues selected by the trunc strategy.

See also MatrixAlgebraKit.eigh_trunc!.

eigh_vals(A, labels_A, labels_codomain, labels_domain; kwargs...) -> D
eigh_vals(A, perm_codomain, perm_domain; kwargs...) -> D
eigh_vals(A, ndims_codomain::Val; kwargs...) -> D

Compute the eigenvalues of a generic N-dimensional array interpreted as a Hermitian linear map from the domain to the codomain dimensions. The output is a vector of eigenvalues.

See also MatrixAlgebraKit.eigh_vals!.

zeros([T,] axes) -> A
ones([T,] axes) -> A
randn([rng,] [T,] axes) -> A
rand([rng,] [T,] axes) -> A
fill(v, axes) -> A

Axis-friendly counterparts of Base.zeros/Base.ones/Base.randn/Base.rand/Base.fill, taking the axes as a single tuple. Base.zeros/Base.ones/Base.fill already accept axes, but Base.randn/Base.rand accept only integer dims, so these fill that gap for dense Base.OneTo axes and otherwise forward to Base (so a graded-axis backend that extends Base.randn/rand on its axis type is picked up). These are the flat (non-map) companions of zeros_map.

zeros_map([T,] axes_codomain, axes_domain) -> M
ones_map([T,] axes_codomain, axes_domain) -> M
randn_map([rng,] [T,] axes_codomain, axes_domain) -> M
rand_map([rng,] [T,] axes_codomain, axes_domain) -> M
fill_map(v, axes_codomain, axes_domain) -> M

Construct an array shaped as a linear map from axes_domain to axes_codomain, filled with zeros (zeros_map), ones (ones_map), normally-distributed values (randn_map), uniformly-distributed values (rand_map), or the value v (fill_map), with element type T (defaulting to Float64; fill_map takes it from v). These are the value-filling companions of similar_map: the domain axes are given un-dualized (codomain facing) and stored dual, so the default flattens to the axis-friendly zeros/ones/randn/rand/fill over (axes_codomain..., conj.(axes_domain)...) (conj dualizes a graded axis and is a no-op on a dense one). Backends with map-shaped storage (e.g. a TensorMap) overload these to build the codomain/domain directly.

flattenlinear(bc::Broadcasted) -> LinearBroadcasted

Like tryflattenlinear, but throw an ArgumentError when the expression is not linear instead of returning nothing. The erroring counterpart to tryflattenlinear, following the parse/tryparse convention.

TensorAlgebra.has_bipartition(a) -> Bool

Whether a's type stores an intrinsic codomain/domain split, so that ndims_codomain reports a split a genuinely carries rather than the all-codomain fallback. A TensorMap returns true, an array false; a type that overloads ndims_codomain should also overload this.

Lets a consumer tell a genuinely all-codomain a apart from one with no notion of a split, which is what validating a caller-supplied split against a needs: there is nothing to validate it against unless a stores one.

infer_aux_space(raw, axes_codomain, axes_domain) -> aux

Derive the auxiliary axis the *_aux projection verbs append as the last domain axis, so the projected result is symmetry-allowed. raw carries the trailing slice axis whose space is derived. This is the backend customization point for that derivation: the generic (dense) method takes the space straight from raw, while a symmetric backend reads it from the sector structure (a graded backend derives per-slice sectors, the TensorMap backend scans the codomain ⊗ conj(domain) content).

invsqrth_safe(A, labels_A, labels_codomain, labels_domain; kwargs...) -> P
invsqrth_safe(A, perm_codomain, perm_domain; kwargs...) -> P
invsqrth_safe(A, ndims_codomain::Val; kwargs...) -> P

Pseudo-inverse square root of a generic N-dimensional array, interpreting it as a Hermitian positive semi-definite linear map from the domain to the codomain dimensions. The result carries the same codomain and domain axes as A. Eigenvalues below tolerance are clamped to zero (Moore-Penrose convention). The input must be Hermitian: project with project_hermitian first if it is Hermitian only up to numerical noise.

Keyword arguments

• alg: forwarded to MatrixAlgebraKit.eigh_full.

• atol::Real: absolute clamping threshold. Default 0.

• rtol::Real: relative clamping threshold. Default eps(real(eltype(A)))^(3//4) when atol = 0, else 0.

See also sqrth_safe, sqrth_invsqrth_safe, and MatrixAlgebra.invsqrth_safe.

is_projected(dest, src, ndims_codomain::Val; kwargs...) -> Bool
is_projected(dest, src; kwargs...) -> Bool

Whether the projected dest still represents src within the isapprox tolerance, i.e. whether the projection that produced dest discarded only a negligible component of src. Keyword arguments are forwarded to isapprox. Compares src against unproject(dest, ndims_codomain), so a backend that changes basis in project is checked in the frame src was given in. The two-argument form uses the destination's own codomain rank.

Together with unchecked_project this is the backend customization point (project and tryproject derive from the two).

isdual(a) -> Bool

Returns true or false depending on if the axis a is dual. Falls back to false for AbstractUnitRange.

See also dual.

label_type(::Type{L}) -> Type

The label type to use when deriving a contraction.

Deriving a contraction makes several passes comparing labels, so a label type that is costly to compare can map here to a cheaper integer type: the labels are matched to integers by equality pattern, the bookkeeping runs on the integers, and the derived labels are mapped back. The default is the identity label_type(::Type{L}) = L, so the bookkeeping runs on the labels as-is. A label type with expensive equality opts in by defining e.g.

TensorAlgebra.label_type(::Type{MyLabel}) = Int

left_null(A, labels_A, labels_codomain, labels_domain; kwargs...) -> N
left_null(A, perm_codomain, perm_domain; kwargs...) -> N
left_null(A, ndims_codomain::Val; kwargs...) -> N

Compute the left nullspace of a generic N-dimensional array, by interpreting it as a linear map from the domain to the codomain dimensions. These can be specified either via their labels or directly through a bi-permutation. The output satisfies N' * A ≈ 0 and N' * N ≈ I.

Keyword arguments

• atol::Real=0: absolute tolerance for the nullspace computation.
• rtol::Real=0: relative tolerance for the nullspace computation.
• kind::Symbol: specify the kind of decomposition used to compute the nullspace. The options are :qr, :qrpos and :svd. The former two require 0 == atol == rtol. The default is :qrpos if atol == rtol == 0, and :svd otherwise.

left_orth(A, labels_A, labels_codomain, labels_domain; kwargs...) -> V, C
left_orth(A, perm_codomain, perm_domain; kwargs...) -> V, C
left_orth(A, ndims_codomain::Val; kwargs...) -> V, C

Compute the left orthogonal decomposition of a generic N-dimensional array, by interpreting it as a linear map from the domain to the codomain dimensions. These can be specified either via their labels or directly through a bi-permutation.

Keyword arguments

• Keyword arguments are passed on directly to MatrixAlgebraKit.

See also MatrixAlgebraKit.left_orth!.

left_polar(A, labels_A, labels_codomain, labels_domain; kwargs...) -> W, P
left_polar(A, perm_codomain, perm_domain; kwargs...) -> W, P
left_polar(A, ndims_codomain::Val; kwargs...) -> W, P

Compute the left polar decomposition of a generic N-dimensional array, by interpreting it as a linear map from the domain to the codomain dimensions. These can be specified either via their labels or directly through a bi-permutation.

Keyword arguments

• Keyword arguments are passed on directly to MatrixAlgebraKit.

See also MatrixAlgebraKit.left_polar!.

linearbroadcasted(f, args...)

Construct a LinearBroadcasted subtype from function f and arguments. Analogous to Base.Broadcast.broadcasted(f, args...).

Examples

linearbroadcasted(*, 2.0, a)   # ScaledBroadcasted(2.0, a)
linearbroadcasted(conj, a)     # ConjBroadcasted(a)
linearbroadcasted(+, a, b)     # AddBroadcasted(a, b)

lq_compact(A, labels_A, labels_codomain, labels_domain; kwargs...) -> L, Q
lq_compact(A, perm_codomain, perm_domain; kwargs...) -> L, Q
lq_compact(A, ndims_codomain::Val; kwargs...) -> L, Q

Compute the compact LQ decomposition of a generic N-dimensional array, by interpreting it as a linear map from the domain to the codomain dimensions, where L is square. The partition is specified either via labels or directly through a bi-permutation.

Keyword arguments

• positive::Bool=false: specify if the diagonal of L should be positive, leading to a unique decomposition.
• Other keywords are passed on directly to MatrixAlgebraKit.

See also MatrixAlgebraKit.lq_compact!.

lq_full(A, labels_A, labels_codomain, labels_domain; kwargs...) -> L, Q
lq_full(A, perm_codomain, perm_domain; kwargs...) -> L, Q
lq_full(A, ndims_codomain::Val; kwargs...) -> L, Q

Compute the full LQ decomposition of a generic N-dimensional array, by interpreting it as a linear map from the domain to the codomain dimensions, where Q is unitary. The partition is specified either via labels or directly through a bi-permutation.

Keyword arguments

• positive::Bool=false: specify if the diagonal of L should be positive, leading to a unique decomposition.
• Other keywords are passed on directly to MatrixAlgebraKit.

See also MatrixAlgebraKit.lq_full!.

matricize(a, perm_codomain, perm_domain)

matricizeop at identity. Has the same maybe-alias semantics.

matricizeop(op, a, perm_codomain, perm_domain)

Matricize a across the bipermutation with the element-wise operation op folded in, i.e. a matrix representing op.(permutedims(a, (perm_codomain..., perm_domain...))) with the codomain fused to rows and the domain to columns.

Has "maybe alias" semantics: the result may share a's memory or be fresh storage, depending on the array type. Treat it as read-only. Use matricizeopcopy for a matrix the caller owns, and matricizeopview (partial) for one guaranteed to alias.

TensorAlgebra.ndims_codomain(a) -> Int
TensorAlgebra.ndims_domain(a) -> Int

The codomain and domain ranks of a's intrinsic split, for a type that stores one. An array has no intrinsic split, so it defaults to all codomain and an empty domain, and a type that does store one overloads ndims_codomain (a TensorMap returns numout).

ndims_domain never needs overloading: it is whatever rank is left over, so the two agree by construction. A type that does overload ndims_codomain should also overload has_bipartition.

TensorAlgebra.ndims_codomain(a) -> Int
TensorAlgebra.ndims_domain(a) -> Int

The codomain and domain ranks of a's intrinsic split, for a type that stores one. An array has no intrinsic split, so it defaults to all codomain and an empty domain, and a type that does store one overloads ndims_codomain (a TensorMap returns numout).

ndims_domain never needs overloading: it is whatever rank is left over, so the two agree by construction. A type that does overload ndims_codomain should also overload has_bipartition.

TensorAlgebra.one(A, labels_A, labels_codomain, labels_domain) -> Id
TensorAlgebra.one(A, perm_codomain, perm_domain) -> Id
TensorAlgebra.one(A, ndims_codomain::Val) -> Id

Construct the identity operator tensor whose shape mirrors A, interpreted as a linear map from the domain to the codomain dimensions. The codomain and domain partition is specified either via labels or directly through a bi-permutation; fused codomain and domain sizes must match. A is treated as a shape prototype and is not mutated.

A tensor generalization in its own right, not an extension of Base.one, so it is neither exported nor imported. Qualify as TensorAlgebra.one(A, ...).

See also MatrixAlgebra.one!, the matrix-level fill this bottoms out on.

Examples

julia> using LinearAlgebra: I

julia> using TensorAlgebra: TensorAlgebra, matricize

julia> A = randn(2, 3, 2, 3);

julia> Id = TensorAlgebra.one(A, (:a, :b, :c, :d), (:a, :b), (:c, :d));

julia> matricize(Id, (1, 2), (3, 4)) ≈ I
true

zeros([T,] axes) -> A
ones([T,] axes) -> A
randn([rng,] [T,] axes) -> A
rand([rng,] [T,] axes) -> A
fill(v, axes) -> A

Axis-friendly counterparts of Base.zeros/Base.ones/Base.randn/Base.rand/Base.fill, taking the axes as a single tuple. Base.zeros/Base.ones/Base.fill already accept axes, but Base.randn/Base.rand accept only integer dims, so these fill that gap for dense Base.OneTo axes and otherwise forward to Base (so a graded-axis backend that extends Base.randn/rand on its axis type is picked up). These are the flat (non-map) companions of zeros_map.

zeros_map([T,] axes_codomain, axes_domain) -> M
ones_map([T,] axes_codomain, axes_domain) -> M
randn_map([rng,] [T,] axes_codomain, axes_domain) -> M
rand_map([rng,] [T,] axes_codomain, axes_domain) -> M
fill_map(v, axes_codomain, axes_domain) -> M

Construct an array shaped as a linear map from axes_domain to axes_codomain, filled with zeros (zeros_map), ones (ones_map), normally-distributed values (randn_map), uniformly-distributed values (rand_map), or the value v (fill_map), with element type T (defaulting to Float64; fill_map takes it from v). These are the value-filling companions of similar_map: the domain axes are given un-dualized (codomain facing) and stored dual, so the default flattens to the axis-friendly zeros/ones/randn/rand/fill over (axes_codomain..., conj.(axes_domain)...) (conj dualizes a graded axis and is a no-op on a dense one). Backends with map-shaped storage (e.g. a TensorMap) overload these to build the codomain/domain directly.

permuteddims(a, perm)

Lazy permutedims. For an AbstractArray this is a Base.PermutedDimsArray view; for any other operand it is a generic PermutedDims node. This is an extension hook: downstream types can overload it to return their own lazy permuted-dims type.

permutedims!(dest, a, perm)
permutedims!(dest, a, perm_codomain, perm_domain)

In-place counterpart of permutedims: write the permuted a into dest. Both forms forward to bipermutedimsopadd! with α, β = true, false; the flat form passes an empty domain permutation.

permutedims(a, perm)
permutedims(a, perm_codomain, perm_domain)

Out-of-place permutation of a, mirroring TensorKit.permute. The single-permutation form reorders every dimension into perm, giving an all-codomain result. The two-permutation form additionally splits the dimensions into a codomain/domain bipartition, with perm_codomain selecting the codomain dimensions and perm_domain the domain ones.

Allocates the destination with similar_map and materializes it through permutedims!, so any operand implementing the permutedimsopadd! / bipermutedimsopadd! interface (a dense array, a graded array, a TensorMap) is permuted with no dedicated method. A dense operand ignores the bipartition and stores the result flat.

permutedimsadd!(dest, src, perm, α, β)

dest = β * dest + α * permutedims(src, perm).

permutedimsop(op, src, perm_codomain, perm_domain)

Non-mutating version of bipermutedimsopadd!: returns op.(permutedims(src, (perm_codomain..., perm_domain...))).

permutedimsopadd!(dest, op, src, perm, α, β)

dest = β * dest + α * permutedims(op.(src), perm).

This is the single materialization primitive for LinearBroadcasted types. Downstream array types should implement bipermutedimsopadd! for the bipartitioned permutation version; this flat-permutation overload forwards to it with perm_domain = ().

project!(dest, src; kwargs...) -> dest

In-place checked projection: project src into the restricted space of dest via projectto! and verify with is_projected that only a negligible component was discarded, throwing an InexactError otherwise (keyword arguments are forwarded to the isapprox tolerance check). This is the checked sibling of the projectto! primitive, in the way copy! relates to copyto!; see project for the allocating form.

project(raw, axes_codomain, axes_domain; kwargs...) -> dest
project(raw, axes; kwargs...) -> dest

Project raw into a symmetry-restricted array shaped as a map from axes_domain to axes_codomain, verifying that only a negligible component of raw is discarded and throwing an InexactError otherwise (keyword arguments are forwarded to the isapprox tolerance check; the default tolerances are subject to change in future versions). See tryproject for a nullable version and unchecked_project for the unchecked projection this derives from.

raw must not have more axes than axes_codomain/axes_domain account for: project projects into exactly the given axes. To append a derived flux-carrying auxiliary axis (for a charge-shifting operator or a non-invariant state), use project_aux. The two-argument form takes a flat list of axes and is equivalent to an empty domain.

project_aux(raw, axes_codomain, axes_domain; kwargs...) -> dest
project_aux(raw, axes; kwargs...) -> dest

Project raw and append a derived auxiliary domain axis carrying its flux, giving a symmetry-allowed result whose squeezed data is the input (the flux-canceling MPO-virtual-leg idiom). Unlike project, which projects into exactly the given axes, project_aux derives the extra leg (see infer_aux_space). raw may have the physical rank (a single operator or state, its flux on a length-1 leg) or one trailing slice axis (an operator multiplet as laid out by stack). Like project, it verifies that only a negligible component is discarded; see unchecked_project_aux and tryproject_aux for the unchecked and nullable siblings.

project_hermitian(A, labels_A, labels_codomain, labels_domain; kwargs...) -> H
project_hermitian(A, perm_codomain, perm_domain; kwargs...) -> H
project_hermitian(A, ndims_codomain::Val; kwargs...) -> H

Hermitian part (M + M') / 2 of a generic N-dimensional array, interpreting it as a linear map M from the domain to the codomain dimensions. The result carries the same codomain and domain axes as A.

See also MatrixAlgebraKit.project_hermitian.

projectto!(dest, src) -> dest

Project src into the restricted space of dest without checking which components may have been projected out. The default reshapes src to size(dest) up to trailing length-1 axes (so a lower-rank src may omit them, e.g. an auxiliary flux-canceling leg a codomain/domain split introduces on a symmetric state) and copyto!s, throwing on a genuine shape mismatch rather than reinterpreting the data. A backend whose arrays are not copyto!-compatible with a dense array overloads this. This is the in-place fill primitive that unchecked_project allocates a destination for.

qr_compact(A, labels_A, labels_codomain, labels_domain; kwargs...) -> Q, R
qr_compact(A, perm_codomain, perm_domain; kwargs...) -> Q, R
qr_compact(A, ndims_codomain::Val; kwargs...) -> Q, R

Compute the compact QR decomposition of a generic N-dimensional array, by interpreting it as a linear map from the domain to the codomain dimensions, where R is square. The partition is specified either via labels or directly through a bi-permutation.

Keyword arguments

• positive::Bool=false: specify if the diagonal of R should be positive, leading to a unique decomposition.
• Other keywords are passed on directly to MatrixAlgebraKit.

See also MatrixAlgebraKit.qr_compact!.

qr_full(A, labels_A, labels_codomain, labels_domain; kwargs...) -> Q, R
qr_full(A, perm_codomain, perm_domain; kwargs...) -> Q, R
qr_full(A, ndims_codomain::Val; kwargs...) -> Q, R

Compute the full QR decomposition of a generic N-dimensional array, by interpreting it as a linear map from the domain to the codomain dimensions, where Q is unitary. The partition is specified either via labels or directly through a bi-permutation.

Keyword arguments

• positive::Bool=false: specify if the diagonal of R should be positive, leading to a unique decomposition.
• Other keywords are passed on directly to MatrixAlgebraKit.

See also MatrixAlgebraKit.qr_full!.

zeros([T,] axes) -> A
ones([T,] axes) -> A
randn([rng,] [T,] axes) -> A
rand([rng,] [T,] axes) -> A
fill(v, axes) -> A

Axis-friendly counterparts of Base.zeros/Base.ones/Base.randn/Base.rand/Base.fill, taking the axes as a single tuple. Base.zeros/Base.ones/Base.fill already accept axes, but Base.randn/Base.rand accept only integer dims, so these fill that gap for dense Base.OneTo axes and otherwise forward to Base (so a graded-axis backend that extends Base.randn/rand on its axis type is picked up). These are the flat (non-map) companions of zeros_map.

zeros_map([T,] axes_codomain, axes_domain) -> M
ones_map([T,] axes_codomain, axes_domain) -> M
randn_map([rng,] [T,] axes_codomain, axes_domain) -> M
rand_map([rng,] [T,] axes_codomain, axes_domain) -> M
fill_map(v, axes_codomain, axes_domain) -> M

Construct an array shaped as a linear map from axes_domain to axes_codomain, filled with zeros (zeros_map), ones (ones_map), normally-distributed values (randn_map), uniformly-distributed values (rand_map), or the value v (fill_map), with element type T (defaulting to Float64; fill_map takes it from v). These are the value-filling companions of similar_map: the domain axes are given un-dualized (codomain facing) and stored dual, so the default flattens to the axis-friendly zeros/ones/randn/rand/fill over (axes_codomain..., conj.(axes_domain)...) (conj dualizes a graded axis and is a no-op on a dense one). Backends with map-shaped storage (e.g. a TensorMap) overload these to build the codomain/domain directly.

zeros([T,] axes) -> A
ones([T,] axes) -> A
randn([rng,] [T,] axes) -> A
rand([rng,] [T,] axes) -> A
fill(v, axes) -> A

Axis-friendly counterparts of Base.zeros/Base.ones/Base.randn/Base.rand/Base.fill, taking the axes as a single tuple. Base.zeros/Base.ones/Base.fill already accept axes, but Base.randn/Base.rand accept only integer dims, so these fill that gap for dense Base.OneTo axes and otherwise forward to Base (so a graded-axis backend that extends Base.randn/rand on its axis type is picked up). These are the flat (non-map) companions of zeros_map.

zeros_map([T,] axes_codomain, axes_domain) -> M
ones_map([T,] axes_codomain, axes_domain) -> M
randn_map([rng,] [T,] axes_codomain, axes_domain) -> M
rand_map([rng,] [T,] axes_codomain, axes_domain) -> M
fill_map(v, axes_codomain, axes_domain) -> M

Construct an array shaped as a linear map from axes_domain to axes_codomain, filled with zeros (zeros_map), ones (ones_map), normally-distributed values (randn_map), uniformly-distributed values (rand_map), or the value v (fill_map), with element type T (defaulting to Float64; fill_map takes it from v). These are the value-filling companions of similar_map: the domain axes are given un-dualized (codomain facing) and stored dual, so the default flattens to the axis-friendly zeros/ones/randn/rand/fill over (axes_codomain..., conj.(axes_domain)...) (conj dualizes a graded axis and is a no-op on a dense one). Backends with map-shaped storage (e.g. a TensorMap) overload these to build the codomain/domain directly.

right_null(A, labels_A, labels_codomain, labels_domain; kwargs...) -> Nᴴ
right_null(A, perm_codomain, perm_domain; kwargs...) -> Nᴴ
right_null(A, ndims_codomain::Val::Val; kwargs...) -> Nᴴ

Compute the right nullspace of a generic N-dimensional array, by interpreting it as a linear map from the domain to the codomain dimensions. These can be specified either via their labels or directly through a bi-permutation. The output satisfies A * Nᴴ' ≈ 0 and Nᴴ * Nᴴ' ≈ I.

Keyword arguments

• atol::Real=0: absolute tolerance for the nullspace computation.
• rtol::Real=0: relative tolerance for the nullspace computation.
• kind::Symbol: specify the kind of decomposition used to compute the nullspace. The options are :lq, :lqpos and :svd. The former two require 0 == atol == rtol. The default is :lqpos if atol == rtol == 0, and :svd otherwise.

right_orth(A, labels_A, labels_codomain, labels_domain; kwargs...) -> C, V
right_orth(A, perm_codomain, perm_domain; kwargs...) -> C, V
right_orth(A, ndims_codomain::Val; kwargs...) -> C, V

Compute the right orthogonal decomposition of a generic N-dimensional array, by interpreting it as a linear map from the domain to the codomain dimensions. These can be specified either via their labels or directly through a bi-permutation.

Keyword arguments

• Keyword arguments are passed on directly to MatrixAlgebraKit.

See also MatrixAlgebraKit.right_orth!.

right_polar(A, labels_A, labels_codomain, labels_domain; kwargs...) -> P, W
right_polar(A, perm_codomain, perm_domain; kwargs...) -> P, W
right_polar(A, ndims_codomain::Val; kwargs...) -> P, W

Compute the right polar decomposition of a generic N-dimensional array, by interpreting it as a linear map from the domain to the codomain dimensions. These can be specified either via their labels or directly through a bi-permutation.

Keyword arguments

• Keyword arguments are passed on directly to MatrixAlgebraKit.

See also MatrixAlgebraKit.right_polar!.

scalar(a)

The single scalar held by a rank-0 (zero-dimensional) a, i.e. a[].

Examples

julia> import TensorAlgebra

julia> TensorAlgebra.scalar(fill(3.0))
3.0

scale!(a::AbstractArray, β::Number) -> a

In-place scaling: multiply every entry of a by β.

scaleadd!(dest, src, α, β)

dest = β * dest + α * src.

TensorAlgebra.select_algorithm(f, alg, args::Tuple)

Resolve the algorithm operation f should run with on args. An alg of nothing defers to TensorAlgebra.default_algorithm, an AbstractAlgorithm passes through unchanged, and anything else is an error.

alg is positional so each operation can dispatch on the algorithm type. The user-facing entry points take it as a keyword and hand it here.

The selection-relevant arguments are packed into a tuple rather than spliced, so that the value and type domains stay disjoint. (Float64, Int) is a pair of arguments that happen to be types, while Tuple{Float64, Int} names the types of a pair of arguments.

similar_map(prototype, [T,] axes_codomain, axes_domain) -> M

Allocate an array shaped as a linear map from axes_domain to axes_codomain with element type T (defaulting to eltype(prototype)), using prototype to determine the array backend. The domain axes are given un-dualized (codomain facing) and stored dual, so the default is similar(prototype, T, (axes_codomain..., conj.(axes_domain)...)). conj dualizes a graded axis and is a no-op on a dense axis. Backends with map-shaped storage (e.g. a TensorMap) overload this to build the codomain/domain directly.

Examples

julia> using TensorAlgebra: similar_map

julia> cod, dom = (Base.OneTo(2), Base.OneTo(3)), (Base.OneTo(4), Base.OneTo(5));

julia> M = similar_map(randn(3), Float32, cod, dom);

julia> eltype(M), size(M)
(Float32, (2, 3, 4, 5))

sqrth_invsqrth_safe(A, labels_A, labels_codomain, labels_domain; kwargs...) -> P, Pinv
sqrth_invsqrth_safe(A, perm_codomain, perm_domain; kwargs...) -> P, Pinv
sqrth_invsqrth_safe(A, ndims_codomain::Val; kwargs...) -> P, Pinv

Square root and pseudo-inverse square root of a generic N-dimensional array (see sqrth_safe and invsqrth_safe), from a single eigendecomposition. Both results carry the same codomain and domain axes as A.

Keyword arguments

• alg: forwarded to MatrixAlgebraKit.eigh_full.

• atol::Real: absolute clamping threshold. Default 0.

• rtol::Real: relative clamping threshold. Default eps(real(eltype(A)))^(3//4) when atol = 0, else 0.

See also MatrixAlgebra.sqrth_invsqrth_safe.

sqrth_safe(A, labels_A, labels_codomain, labels_domain; kwargs...) -> P
sqrth_safe(A, perm_codomain, perm_domain; kwargs...) -> P
sqrth_safe(A, ndims_codomain::Val; kwargs...) -> P

Square root of a generic N-dimensional array, interpreting it as a Hermitian positive semi-definite linear map from the domain to the codomain dimensions. The result carries the same codomain and domain axes as A. Eigenvalues below tolerance are clamped to zero. The input must be Hermitian: project with project_hermitian first if it is Hermitian only up to numerical noise.

Keyword arguments

• alg: forwarded to MatrixAlgebraKit.eigh_full.

• atol::Real: absolute clamping threshold. Default 0.

• rtol::Real: relative clamping threshold. Default eps(real(eltype(A)))^(3//4) when atol = 0, else 0.

See also invsqrth_safe, sqrth_invsqrth_safe, and MatrixAlgebra.sqrth_safe.

svd_compact(A, labels_A, labels_codomain, labels_domain; kwargs...) -> U, S, Vᴴ
svd_compact(A, perm_codomain, perm_domain; kwargs...) -> U, S, Vᴴ
svd_compact(A, ndims_codomain::Val; kwargs...) -> U, S, Vᴴ

Compute the compact (thin) SVD of a generic N-dimensional array, by interpreting it as a linear map from the domain to the codomain dimensions, where U and Vᴴ are isometric. The partition is specified either via labels or directly through a bi-permutation.

See also MatrixAlgebraKit.svd_compact!.

svd_full(A, labels_A, labels_codomain, labels_domain; kwargs...) -> U, S, Vᴴ
svd_full(A, perm_codomain, perm_domain; kwargs...) -> U, S, Vᴴ
svd_full(A, ndims_codomain::Val; kwargs...) -> U, S, Vᴴ

Compute the full (thick) SVD of a generic N-dimensional array, by interpreting it as a linear map from the domain to the codomain dimensions, where U and Vᴴ are unitary. The partition is specified either via labels or directly through a bi-permutation.

See also MatrixAlgebraKit.svd_full!.

svd_trunc(A, labels_A, labels_codomain, labels_domain; trunc, kwargs...) -> U, S, Vᴴ, ϵ
svd_trunc(A, perm_codomain, perm_domain; trunc, kwargs...) -> U, S, Vᴴ, ϵ
svd_trunc(A, ndims_codomain::Val; trunc, kwargs...) -> U, S, Vᴴ, ϵ

Compute the truncated SVD of a generic N-dimensional array, by interpreting it as a linear map from the domain to the codomain dimensions. The partition is specified either via labels or directly through a bi-permutation. In addition to the factors, returns the truncation error ϵ, the 2-norm of the discarded singular values.

Keyword arguments

• trunc: truncation strategy, passed on to MatrixAlgebraKit.svd_trunc.
• Other keywords are passed on directly to MatrixAlgebraKit.

Examples

julia> using TensorAlgebra: svd_trunc, contractalign

julia> A = randn(4, 4);

julia> U, S, Vᴴ, ϵ = svd_trunc(A, (:i, :j), (:i,), (:j,));

julia> SV = contractalign((:u, :j), S, (:u, :v), Vᴴ, (:v, :j));

julia> contractalign((:i, :j), U, (:i, :u), SV, (:u, :j)) ≈ A
true

julia> isapprox(ϵ, 0; atol = 1e-10)
true

See also MatrixAlgebraKit.svd_trunc!.

svd_vals(A, labels_A, labels_codomain, labels_domain) -> S
svd_vals(A, perm_codomain, perm_domain) -> S
svd_vals(A, ndims_codomain::Val) -> S

Compute the singular values of a generic N-dimensional array, by interpreting it as a linear map from the domain to the codomain dimensions. The partition is specified either via labels or directly through a bi-permutation. The output is a vector of singular values.

See also MatrixAlgebraKit.svd_vals!.

TensorAlgebra.to_range(space)

Convert a description of a space into a default range type, returning the range unchanged if it already is one. This lets range and axis constructors accept a space uniformly instead of each reimplementing the conversion.

• an Integer length -> Base.OneTo
• an existing AbstractUnitRange -> itself (idempotent passthrough)

Downstream packages extend this for richer spaces; for example, GradedArrays adds a method that turns a vector of sector-to-multiplicity pairs into a graded range.

TensorAlgebra.tr(A, labels_A, labels_codomain, labels_domain)
TensorAlgebra.tr(A, perm_codomain, perm_domain)
TensorAlgebra.tr(A, ndims_codomain::Val)

Trace of a generic N-dimensional array A interpreted as a linear map from its domain to its codomain dimensions. The map is matricized into its square matrix, then the matrix trace is taken, so the backend's own matrix tr (dense, graded, or TensorMap) does the work. The partition is specified via labels, a bi-permutation, or directly as the codomain rank, matching the factorization entry points.

This is TensorAlgebra's own function, distinct from LinearAlgebra.tr; the two-argument and higher forms take a codomain/domain partition rather than a bare matrix.

Examples

julia> import TensorAlgebra

julia> A = randn(2, 2, 2, 2);

julia> TensorAlgebra.tr(A, (:i, :j, :k, :l), (:i, :k), (:j, :l)) ≈
       sum(A[i, i, k, k] for i in 1:2, k in 1:2)
true

TensorAlgebra.trivialrange(R::Type{<:AbstractUnitRange}[, n::Integer])
TensorAlgebra.trivialrange(r::AbstractUnitRange[, n::Integer])

Return the identity range for fusing ranges of type R: a one-dimensional range t for which fusing t with any other range of the same family leaves that range unchanged. Defaults to Base.OneTo(1). Downstream packages overload the type-level methods to return their own identity (for example, a charge-0 one-dimensional sector for a graded range).

With a length n, return the n-dimensional analogue: n copies of the identity range stacked into one range (for a graded range, a charge-0 sector of dimension n).

tryflattenlinear(bc::Broadcasted) -> LinearBroadcasted or nothing

Recursively convert a Broadcasted tree to a LinearBroadcasted tree. Returns nothing if any node is not linear (as determined by islinearbroadcast).

Analogous to Broadcast.flatten for Broadcasted trees, but converts to LinearBroadcasted subtypes via linearbroadcasted.

Downstream styles call this from Base.copy(::Broadcasted{MyStyle}) to opt into linear broadcasting at materialization time.

tryproject(raw, axes_codomain, axes_domain; kwargs...) -> Union{dest, Nothing}
tryproject(raw, axes; kwargs...) -> Union{dest, Nothing}

Like project, but return nothing instead of throwing when more than a negligible component of raw would be discarded. Useful for branching on whether raw is symmetry-allowed in the given axes, e.g. projecting a state as invariant and falling back to deriving an auxiliary flux-carrying leg:

@something tryproject(v, (cod,)) project_aux(v, (cod,))

Keyword arguments are forwarded to the isapprox tolerance check.

tryproject_aux(raw, axes_codomain, axes_domain; kwargs...) -> Union{dest, Nothing}
tryproject_aux(raw, axes; kwargs...) -> Union{dest, Nothing}

The nullable sibling of project_aux: derive and append the auxiliary axis, returning nothing instead of throwing when more than a negligible component of raw would be discarded.

unchecked_project(raw, axes_codomain, axes_domain) -> dest
unchecked_project(raw, axes) -> dest

Project raw into a symmetry-restricted array shaped as a map from axes_domain to axes_codomain, without checking which components are discarded: entries of raw outside the symmetry-allowed structure are dropped without inspection. Most callers want project, which verifies that nothing was discarded, or tryproject, its nullable sibling. All three derive from the backend customization points: this one is projectto!(allocate_project(raw, axes_codomain, axes_domain), raw). The two-argument form takes a flat list of axes and is equivalent to an empty domain.

unchecked_project_aux(raw, axes_codomain, axes_domain) -> dest
unchecked_project_aux(raw, axes) -> dest

The unchecked sibling of project_aux: derive and append the auxiliary axis, then project without verifying which components are discarded.

TensorAlgebra.ungrade(r)

Return the plain range underlying an axis, keeping only its extent and stripping any added structure. On a plain AbstractUnitRange this is the identity (there is nothing to strip, and any offset is preserved). Downstream packages extend it for richer axes: GradedArrays maps a graded range to the Base.OneTo of its total dimension, and a native TensorKit space maps to the Base.OneTo of its dimension.

Examples

julia> import TensorAlgebra

julia> TensorAlgebra.ungrade(2:5)
2:5

unproject(a, ndims_codomain::Val) -> raw

Inverse of project: recover the dense array that project maps to a, given the codomain/domain split ndims_codomain as a Val. The default is convert(Array, a); a backend that changes basis in project overloads this to undo that change, so that

unproject(project(raw, axes_codomain, axes_domain), Val(length(axes_codomain))) ≈ raw

zero!(a::AbstractArray) -> a

In-place zero: set every entry of a to zero.

zeros([T,] axes) -> A
ones([T,] axes) -> A
randn([rng,] [T,] axes) -> A
rand([rng,] [T,] axes) -> A
fill(v, axes) -> A

Axis-friendly counterparts of Base.zeros/Base.ones/Base.randn/Base.rand/Base.fill, taking the axes as a single tuple. Base.zeros/Base.ones/Base.fill already accept axes, but Base.randn/Base.rand accept only integer dims, so these fill that gap for dense Base.OneTo axes and otherwise forward to Base (so a graded-axis backend that extends Base.randn/rand on its axis type is picked up). These are the flat (non-map) companions of zeros_map.

zeros_map([T,] axes_codomain, axes_domain) -> M
ones_map([T,] axes_codomain, axes_domain) -> M
randn_map([rng,] [T,] axes_codomain, axes_domain) -> M
rand_map([rng,] [T,] axes_codomain, axes_domain) -> M
fill_map(v, axes_codomain, axes_domain) -> M

Construct an array shaped as a linear map from axes_domain to axes_codomain, filled with zeros (zeros_map), ones (ones_map), normally-distributed values (randn_map), uniformly-distributed values (rand_map), or the value v (fill_map), with element type T (defaulting to Float64; fill_map takes it from v). These are the value-filling companions of similar_map: the domain axes are given un-dualized (codomain facing) and stored dual, so the default flattens to the axis-friendly zeros/ones/randn/rand/fill over (axes_codomain..., conj.(axes_domain)...) (conj dualizes a graded axis and is a no-op on a dense one). Backends with map-shaped storage (e.g. a TensorMap) overload these to build the codomain/domain directly.

invsqrt_diag_safe(D; atol=0, rtol=eps(real(eltype(D)))^(3//4)) -> D^(-1//2)

Inverse square root of a diagonal-structured matrix D, treating diagonal entries below tolerance as zero (Moore-Penrose convention). Equivalent to pow_diag_safe(D, -1//2; atol, rtol).

Keyword arguments

• atol::Real: absolute clamping threshold. Default 0.
• rtol::Real: relative clamping threshold. Default eps(real(eltype(D)))^(3//4) when atol = 0, else 0.

invsqrth_safe(M; alg=nothing, atol=0, rtol=eps(real(eltype(M)))^(3//4)) -> M^(-1//2)

Inverse square root of a Hermitian positive semi-definite matrix. Equivalent to powh_safe(M, -1//2; alg, atol, rtol).

Keyword arguments

• alg: forwarded to MatrixAlgebraKit.eigh_full.

• atol::Real: absolute clamping threshold. Default 0.

• rtol::Real: relative clamping threshold. Default eps(real(eltype(M)))^(3//4) when atol = 0, else 0.

MatrixAlgebra.one!(m) -> m

Fill m with the identity in place. The matrix-level identity fill the tensor-level TensorAlgebra.one/one! bottom out on, and the customization point a backend overloads when its matricization is a type MatrixAlgebraKit.one! does not handle (TensorKit owns its own one! generic rather than extending MatrixAlgebraKit's, so a TensorMap fills through that).

pow_diag_safe!(Dp, D, p, tol) -> Dp

In-place kernel behind pow_diag_safe: write the clamped powers of D's diagonal onto Dp, where entries d with abs(d) < tol become zero and the rest become d^p. Dp must be diagonal-structured with the same structure as D (in the allocating pow_diag_safe path it is copy(D)).

This is the backend overload point. The generic method maps over MAK.diagview, valid when the diagonal view is a single in-place-mappable vector (dense). A backend whose diagonal is stored per sector (a TensorMap, or a graded matrix) overloads this method to clamp each reduced block.

pow_diag_safe(D, p; atol=0, rtol=eps(real(eltype(D)))^(3//4)) -> D^p
pow_diag_safe(D, p, tol) -> D^p

Raise a diagonal-structured matrix D to the power p. Diagonal entries d with abs(d) < tol are clamped to zero before exponentiation, where tol = max(atol, rtol * norm(D, Inf)) (the largest-magnitude entry, which is the largest-magnitude diagonal entry for a diagonal-structured matrix). The three-argument form takes tol directly. Negative d above tol cause d^p to error for fractional p (e.g. p = 1//2) and pass through for integer p, so the operation itself enforces the PSD precondition per-power. Errors if isdiag(D) is false.

The clamped powers are written onto a copy of D via pow_diag_safe!, so the result keeps the input's type and structure. This drives sqrt_diag_safe, invsqrt_diag_safe, and the powh_safe family, and works for any diagonal-structured backend (dense, graded, or a TensorMap).

Keyword arguments

• atol::Real: absolute clamping threshold. Default 0.
• rtol::Real: relative clamping threshold. Default eps(real(eltype(D)))^(3//4) when atol = 0, else 0.

powh_safe(M, p; alg=nothing, atol=0, rtol=eps(real(eltype(M)))^(3//4)) -> M^p

Raise a Hermitian positive semi-definite matrix to the power p. For diagonal-structured M (isdiag(M) == true), dispatches to pow_diag_safe and skips the eigendecomposition. Otherwise computes via M = V * D * V' as V * pow_diag_safe(D, p; atol, rtol) * V'.

The input must be Hermitian (as for MatrixAlgebraKit.eigh_full): project with MatrixAlgebraKit.project_hermitian first if it is Hermitian only up to numerical noise.

Keyword arguments

• alg: forwarded to MatrixAlgebraKit.eigh_full.

• atol::Real: absolute clamping threshold. Default 0.

• rtol::Real: relative clamping threshold. Default eps(real(eltype(M)))^(3//4) when atol = 0, else 0.

sqrt_diag_safe(D; atol=0, rtol=eps(real(eltype(D)))^(3//4)) -> D^(1//2)

Square root of a diagonal-structured matrix D, equivalent to pow_diag_safe(D, 1//2; atol, rtol).

Keyword arguments

• atol::Real: absolute clamping threshold. Default 0.
• rtol::Real: relative clamping threshold. Default eps(real(eltype(D)))^(3//4) when atol = 0, else 0.

sqrth_invsqrth_safe(M; alg=nothing, atol=0, rtol=eps(real(eltype(M)))^(3//4)) -> M^(1//2), M^(-1//2)

Square root and pseudo-inverse square root of a Hermitian positive semi-definite matrix, from a single eigendecomposition. Equivalent to (sqrth_safe(M; ...), invsqrth_safe(M; ...)) but with the eigendecomposition computed once. Eigenvalues below tolerance are clamped to zero in both factors (Moore-Penrose convention for the inverse).

The input must be Hermitian (as for MatrixAlgebraKit.eigh_full): project with MatrixAlgebraKit.project_hermitian first if it is Hermitian only up to numerical noise.

Keyword arguments

• alg: forwarded to MatrixAlgebraKit.eigh_full.

• atol::Real: absolute clamping threshold. Default 0.

• rtol::Real: relative clamping threshold. Default eps(real(eltype(M)))^(3//4) when atol = 0, else 0.

sqrth_safe(M; alg=nothing, atol=0, rtol=eps(real(eltype(M)))^(3//4)) -> M^(1//2)

Square root of a Hermitian positive semi-definite matrix. Equivalent to powh_safe(M, 1//2; alg, atol, rtol).

Keyword arguments

• alg: forwarded to MatrixAlgebraKit.eigh_full.

• atol::Real: absolute clamping threshold. Default 0.

• rtol::Real: relative clamping threshold. Default eps(real(eltype(M)))^(3//4) when atol = 0, else 0.

truncdegen(trunc::TruncationStrategy; atol::Real=0, rtol::Real=0)

Modify a truncation strategy so that if the truncation falls within a degenerate subspace, the entire subspace gets truncated as well. A value val is considered degenerate if norm(val - truncval) ≤ max(atol, rtol * norm(truncval)) where truncval is the largest value truncated by the original truncation strategy trunc.

For now, this truncation strategy assumes the spectrum being truncated has already been reverse sorted and the strategy being wrapped outputs a contiguous subset of values including the largest one. It also only truncates for now, so may not respect if a minimum dimension was requested in the strategy being wrapped. These restrictions may be lifted in the future or provided through a different truncation strategy.