Internals

DeclaredScales(w, t)

The two node channels a block may ask for beyond the values: sensitivity_scales seeded with the input scales w, and the live_nodes flags. Each is built by the first block that asks for it and shared by every block after.

Deferring the sweeps is what makes them free on a program that never needs them. Deciding in advance would take its own O(L) pass over the recorded blocks, which on a tape whose blocks all keep a usable value pair costs about as much as the cover sweep it precedes — so every program would pay for the possibility of a fallback none of its blocks reach.

DeferredScale(store, index)

A declared scale named but not yet computed: slot index of store, resolved by declared_scale when a block first reaches the fallback.

ElementalSite(val, operands)

The CoverSite of an elemental in a forward propagation. Without a tape to number the elemental, only its value and operand values identify it, so those are what the site reports.

operands is held as given and read only when an error is raised, so naming a site costs nothing on the path that succeeds.

ForwardChannels{T, S, M, SM}

The two vector channels an FVar carries beyond its value and its declared scale: the intermediate gradient g = ∇v (storage S, entries T), and the squared cover G (storage SM, entries M = boundtype(T)) maintaining sqrt.(G) * sqrt.(G)' >= magnitude.(∇²v) entrywise.

G is squared rather than rooted because the recursion that maintains it is linear in the parents' G, which is what keeps a per-node re-merge from compounding with depth.

The channels live in the FVar's own fields; channels assembles this view of them on demand, and it costs nothing once inlined.

ForwardContext{T, S, M, SM}

State shared by every variable in a forward propagation: the degeneracy threshold local covers split at, whether each local cover is verified, and the scratch a node writes through while assembling its cover.

S is the storage of the intermediate-gradient channel, one entry of T per input, and SM that of the cover channel, one entry of M = boundtype(T) per input. forward_cover derives both from the point it is given, so a Vector input propagates in Vectors and a StaticVector input in StaticVectors. zerogradient and zerochannel are zeros of the two storages, which fresh channels are copied from; the input variables share zerochannel as their cover channel. scratch is the atom accumulator, written in place when SM is mutable and left alone when it is not; and rho holds the slot weights of the elemental being covered.

It is a reference type so that every variable of the propagation holds a pointer to the one context: storing an immutable into a variable's field would copy it, at one allocation per node. Identity is then also what === compares, which is what distinguishes two propagations.

ForwardScale(sigma, live)

The declared channel of one slot of a forward propagation: the scale sigma that stands in for a parent value carrying none, and whether the slot's intermediate gradient is identically zero.

Liveness travels separately from sigma because sigma cannot report it. It is a sum of products of first partials, so it underflows on a long chain of small factors while the intermediate gradient it bounds is still there, and a slot dropped on that evidence would leave a real Hessian entry uncovered. The flag is carried by the same recursion live_nodes runs over a tape — a slot is live when some parent is live and reaches it through a nonzero partial — so both backends drop exactly the same entries.

MaxMerge <: MergeRule

b_j = max_k Ũ_kj/a_k, the exact minimal b for the given a. Non-smooth, so its marginal alternation can settle into a limit cycle: the weights keep changing by a fixed amount every round, which no larger budget reduces. balance_cover reports that as a stall.

MergeRule

How balance_cover's b-step combines a column of the scaled surrogate into one input weight: given the entries Ũ_kj/a_k over the outputs k, it returns b_j. Every rule that dominates each of its terms yields a valid cover, so the choice trades tightness, never validity.

SumMerge, QuadMerge and MaxMerge are the p = 1, 2, ∞ members of one family, so a parameterized member (a p-norm merge) is reachable through the same b-step without an API break. They are types rather than names so that the b-step's doubly-nested loop specializes on the choice, and so that the non-convergence hint can dispatch per rule.

SignAveragedMerge is the exception: it runs the whole alternation, not only the final b-step, in quadrature on R, and the pair it returns is a sign-averaged cover, which does not cover J.

arule, the other strategy axis of balance_cover, is a Symbol instead, because it is closed: see that function's docstring.

QuadMerge <: MergeRule

b_j = √(∑_k (Ũ_kj/a_k)²), combining the same terms in quadrature.

SignAveragedMerge <: MergeRule

The sign-averaged pair: balance_cover's marginal alternation run in quadrature on R (quadrature_surrogate),

a_k ← sqrt((1/ρ_k) ∑_j R_kj² / b_j²),    b_j ← sqrt((1/γ_j) ∑_k R_kj² / a_k²),

followed by the final b-step b_j = sqrt(∑_k R_kj² / a_k²). Squared, this is the SumMerge alternation applied to the squared channel Ũ⁽²⁾ = R.^2 with factors (a², b²), the same targets and round budget, and the squared start. The iteration works with a, b and R themselves and accumulates every sum of squares with hypot, so it never forms the square of a Jacobian-scale quantity, and the pair is representable wherever the plain pair is. On a tape each half-step is one sweep combining local partials in quadrature, so a round costs O(L) and R is never formed; the final b-step equals gramvariance(t, outidxs, w) with w_k = 1/a_k² on the live outputs and zero on the others, where those squares are representable.

The pair is a sign-averaged cover, not a cover of J. The final b-step gives every live column ∑_k R_kj² / (a_k b_j)² = 1, so a_k b_j >= R_kj, which is what check = true verifies. But R <= Ũ, with equality only where an entry has a single path, so a_k b_j can fall below Ũ_kj, and below abs(J_kj) where several paths add with the same sign. If each parent slot's local partial carries an independent, equiprobable sign, R_kj² = E[J_kj²], and at convergence the pair standardizes the signed Jacobian in expectation: the normalized iterate (a, b ./ sqrt.(γ)) has expected row sums of squares ∑_j E[J_kj²] / (a_k² b_j² / γ_j) = ρ_k and column sums γ_j, and the returned b has expected column sums of squares 1.

SumMerge <: MergeRule

b_j = ∑_k Ũ_kj/a_k, the loosest merge and, with SignAveragedMerge, the only one the tape form can run as O(L) sweeps.

UndecidedAbs(x::CoverVar)

abs(x) for a non-constant x whose sign is not decided (see decided_sign). It records nothing. A constant power UndecidedAbs(x)^p with p >= 2 records |x|^p as one node on x; every recorder throws a DomainError when handed one.

_flow_targets(t::Tape, outidxs, li::LiveInfo) -> (r, c)

Row and column targets from the two flow weightings of the live slot graph, mixed so that neither side dominates: r = ρ_fwd/nrow + ρ_rev/ncol and c = γ_fwd/nrow + γ_rev/ncol, with nrow the live output positions and ncol the live inputs. Marginals are linear in the weight matrix, so the mixture is itself the marginal pair of a nonnegative matrix with Ũ's pattern, which is what makes the targets realizable.

Two of the four marginals are constant by construction — forward flow reaches every live node with mass 1, and reverse flow leaves every live input with mass 1 — so only two sweeps run. O(L), forming no matrix; every weight is at most 1, so nothing can overflow. The two totals are 2 each, to roundoff.

_structural_pattern(t::Tape, outidxs, li::LiveInfo) -> P::Matrix{UInt64}

The structural sparsity pattern of Ũ packed 64 inputs to a word: bit (j - 1) % 64 of P[k, cld(j, 64)] is set exactly when output position k reaches input j along a path of nonzero local derivatives. cld(n, 64) integer passes over the tape and m * cld(n, 64) words, both small beside the m * n surrogate this is compared against.

_support_counts(t::Tape, outidxs, li::LiveInfo) -> (r, c)

Exact reachability counts: r[k] live inputs reaching output position k, c[j] live output positions reached by input j — the row and column nonzero counts of Ũ, obtained without forming it. Word-parallel bit propagation over the live slot graph: cld(n, 64) forward passes and cld(m, 64) reverse passes, each O(L) integer work.

abs_forward(t::Tape, theta::AbstractVector) -> s

One forward sweep in absolute-value arithmetic, the transpose of abs_reverse: s[k] = theta[k] + sum_i abs(∂φ_k/∂v_{p_i}) * s[p_i] over the parents p_i of node k, in tape order. Seeding theta[j] = 1 at a single input j gives, at every node, the absolute sensitivity of that node to input j — the j-th column of the surrogate Ũ read at the output nodes — and the sweep is linear in the seeds, so any weighted combination of inputs costs the same one O(L) traversal.

abs_reverse(t::Tape, Theta::AbstractMatrix) -> Abar

Vector-mode form of abs_reverse: each row of Theta is an independent seed vector over the nodes (columns), and all rows are propagated in one traversal. Abar[:, a] for an input a holds that input's value for every channel. The seeds are nonnegative magnitudes, as group_seeds returns them in boundtype of the tape's element type, and the sweep computes in their type.

abs_reverse(t::Tape, theta) -> abar

One reverse sweep in absolute-value arithmetic: abar[j] = theta[j] + sum_{k : j in P_k} magnitude(∂φ_k/∂v_j) * abar[k].

With theta from local_covers, abar[1:t.ninputs] is the sum-merge cover — the whole merge in a single gradient-cost sweep, because sum_k s_k = sum_j theta[j] * utilde_j is linear in the seeds.

The sweep computes in the element type of theta, so seeds in boundtype give a sweep that rounds as the tape asks.

balance_cover(t::Tape, outidxs; kwargs...) -> (; a, b, iters, delta, roots, nties)
balance_cover(r, x0::AbstractVector; kwargs...)
balance_cover(Utilde::AbstractMatrix; kwargs...)

Output and input weights covering the absolute-sensitivity surrogate: a[k] * b[j] >= Ũ_kj, and so >= abs(J_kj), for every merge except SignAveragedMerge.

By default the weights come from a scale-covariant start, rounds rounds of a marginal alternation toward row and column targets the sparsity pattern makes realizable, and a final MergeRule b-step for the a those rounds reached. The b-step is what makes the pair a cover, so the result covers whatever the rounds did or did not settle. The accompanying manuscript derives the weighting rules and their convergence; the keywords below say what each one does without it.

Keywords

• start = :covariant: products on the support then transform exactly under a rescaling of the inputs and outputs — and, on the tape path, the intermediates — at every round count. start = :equal starts at a = 1 instead, which is permutation equivariant at every round count but covariant only in the limit. Both reach the same fixed point.
• a0: starting weights of your own, positive on every live output and zero on the others. Combining it with an explicit start is an error.
• targets: :flow (tape default; O(L); maximizes a flow-weighted mean log-tightness, which can be noticeably looser than :support when outputs differ widely in the number of paths reaching them), :support (matrix default; maximizes the geometric-mean tightness of the sum-merge cover; on a tape it costs cld(n, 64) + cld(m, 64) bit passes), an explicit (r, c), or :uniform for unit row and column sums. :uniform needs the pattern to have total support — every nonzero entry on some positive diagonal — and without it the weights drift, separating without bound while every iterate still covers, so a longer run returns a looser cover. Applies to arule = :marginal only, and :flow needs a tape.
• rounds = 4: run exactly this many rounds, with no convergence test and no error. 0 merges the start as it is, without computing targets, and reports delta = Inf. nothing runs to delta <= rtol instead and throws if maxiter rounds do not reach it; naming maxiter or rtol and no rounds asks for that too. Such a run overrelaxes the alternation toward targets, in log space with a factor adapted to the observed contraction, which reaches the same fixed point in far fewer rounds on slowly converging patterns; a fixed budget always runs the plain rounds.
• maxiter, rtol: bound and tolerance of a rounds = nothing run. maxiter = nothing picks 1000 rounds for arule = :marginal and 10000 for :geomean. Passing maxiter together with a rounds count is an error.
• merge: the final b-step, SumMerge (default), QuadMerge or MaxMerge. SignAveragedMerge instead runs the whole alternation in quadrature on R and returns a sign-averaged pair, which does not cover J; the matrix form then takes R (quadrature_surrogate) in place of Ũ, and check verifies a_k b_j >= R_kj.
• arule: :marginal (default) alternates marginal_weights with the b-step; :geomean fits log Ũ_kj ≈ log a_k + log b_j instead and prescribes no targets.
• root, cgiters: passed to tape_potentials for the tape-potential start. Naming them on a call that does not reach that start is an error.
• check = true: verify the returned pair against the surrogate.

Returns

a and b by output position and input, iters rounds run, delta the quantity the convergence test applies to (the gauge-invariant product change max_kj |log((a_k b_j)_new/(a_k b_j)_old)|, or the relative row-marginal residual max_k |∑_j Ũ_kj/(a_k b_j r_k) - 1| for non-uniform targets), and roots, nties from tape_potentials when the tape-potential start was used, nothing otherwise. A zero row of Ũ gets a[k] = 0 and a zero column b[j] = 0. With SignAveragedMerge(), delta is measured on the pair itself and R: the product change is that of a_k b_j, and the residual is max_k |sqrt(∑_j R_kj²/b_j²) / (a_k sqrt(r_k)) - 1|.

Cost

The tape form with merge = SumMerge() or SignAveragedMerge() and arule = :marginal runs as O(L) sweeps per round and never forms Ũ. Every other combination, and check = true, forms it (or R) at O(m*L).

Examples

julia> r(x) = [x[1] * x[2], exp(x[1]) - x[3], x[2] + 2x[3]];

julia> bal = DerivativeCovers.balance_cover(r, [0.5, 2.0, -1.0]);

julia> jc = cover_jacobian(WithDiagnostics(), r, [0.5, 2.0, -1.0]; a = bal.a);

julia> bal.iters, all(jc.b_max .<= bal.b)
(4, true)

balanced_split(c11, c12, c22, z1, z2) -> (rho1, rho2)

Rank-1 cover of an off-diagonal-active 2×2 block [c11 c12; c12 c22], which requires c12^2 > c11*c22 (see offdiag_active). Returns rho1 >= sqrt(c11) and rho2 >= sqrt(c22) with rho1*rho2 == c12: the cover is exactly tight on the off-diagonal, so no cover of the block is smaller in both entries.

Those constraints leave a one-parameter family. The scale references z1, z2 — one per parent, supplied by split_block: the parent value magnitudes wherever they carry a scale, the weighted absolute sensitivities where they do not — pick a member, giving rho_i the units of sqrt(f)/[v_{p_i}] and making the cover scale-covariant. The result is continuous in all five arguments except where a vanishing scale reference meets a degenerate block; see the extended help.

Extended help

The three constraints confine rho1 to the segment [sqrt(c11), c12/sqrt(c22)] — nonempty precisely because c12^2 > c11*c22 — and fix rho2 = c12/rho1. Every returned cover lies on this segment, which is what makes it entrywise unimprovable: a competitor smaller in both entries would have rho1*rho2 < c12.

The position along the segment is the point that equalizes the value-weighted scales, rho1*z1 == rho2*z2, clamped to the segment when that point falls outside. Covariance follows: under a rescaling v_i -> d_i*v_i of the intermediate variables, c_ij -> c_ij/(d_i*d_j) and z_i -> d_i*z_i, so both endpoints and the equalizing point scale as 1/d_1; clamp is positively homogeneous, so rho_1 -> rho_1/d_1 and rho_2 -> rho_2/d_2.

Clamping is also what keeps the result continuous as a parent value vanishes: the equalizing point runs off to an endpoint, and the endpoint absorbs it. Nothing computes those limits separately, because floating-point arithmetic already carries them. sqrt(c12*z2/z1) is Inf at z1 == 0 and 0 at z2 == 0, which clamp resolves to the upper and the lower endpoint; c12/sqrt(c22) is Inf at c22 == 0, placing the upper endpoint there. Only z1 == z2 == 0 leaves a NaN that no clamp resolves, and that is what reaches the fallback.

That reasoning needs the endpoint being approached to be a positive finite scale, and it is not when the corresponding diagonal entry vanishes: c22 == 0 puts the upper endpoint at infinity, c11 == 0 puts the lower endpoint at zero. Nor is there any limit to take when both values vanish at once, since the equalizing point then depends on the ratio z2/z1 alone. In those cases the geometric mean of the endpoints is used, which is covariant and permutation-equivariant.

So the result is continuous except at three families of inputs, where it is the ratio z2/z1 that carries the information and the ratio is unavailable:

• z1 == 0 with c22 == 0,
• z2 == 0 with c11 == 0,
• z1 == z2 == 0, for any block.

A multiplication node with a vanishing operand has c11 == c22 == 0 and meets all three at once. Among the other elementals of overloads.jl, a quotient u/w meets the first at u == 0, since c22 = 2u/w^3 vanishes with the numerator, and x^y meets the second at a zero exponent, where c11 == 0 and z2 == 0 while c22 > 0.

The continuity above is continuity in the values with the block held fixed. At those two elementals the block moves with the value: the diagonal entry that bounds the feasible segment is proportional to the vanishing value, so the endpoint the clamp settles on moves too, and the seed of the slot whose value vanishes grows like |value|^(-1/2). The split stays covariant and valid, but the cover grows without bound as the value approaches zero while the derivatives it covers stay finite, and at the zero itself the block falls into the families above and the cover jumps to a finite one.

These zero-reference cases concern the value channel, and they are exactly the ones value_bounds_split detects: split_block substitutes the declared scales there, and a declared scale vanishes only at a slot whose intermediate gradient is zero, which cover_slots! covers diagonally instead of splitting. So a split reaches this function with a pair that is not degenerate, and a pair carrying a ratio reaches the geometric mean only when c12*z2/z1 leaves the range of the element type at a block whose feasible segment has a degenerate endpoint. There the equalizing point is recomputed from the separated exponents, which represents it whenever it is representable at all; see split_fallback. A split that is genuinely out of range throws, since every pair returned here is a cover. Entries that are already non-finite are returned as IEEE arithmetic gives them, which is what check=false asks for.

channels(x::FVar) -> ForwardChannels

The intermediate-gradient and cover channels behind x. A constant has neither: every elemental folds a constant operand into one that propagates, so reaching here with a constant means a caller bypassed that folding.

check_input_scales(w, ninputs::Integer) -> w

Validate declared input scales: nothing for unit scales, or one strictly positive finite entry per input.

Positivity is not a formality. The declared scales stand in for the parent values at blocks whose values carry none, so a zero entry would blank a live slot's scale, leaving its split with no ratio to read and, at a block whose feasible segment is degenerate too, no cover to return.

cyclic_partition(t::Tape, ngroups::Integer) -> groups

Assign node k of t to group 1 + (k-1) % ngroups, a node-granular partition for cover_hybrid.

Looseness comes from summing contributions whose covers overlap, and overlapping contributions are typically siblings — the many product terms of a quadratic form, the fan-out of a reduction — which a tape lists consecutively. Dealing them round-robin therefore spreads them across groups, where quadrature combines them. Binning by depth instead (see level_partition) puts every sibling in the same group, so on such programs it recovers little of the quadrature merge.

Siblings absorbed into a single array primitive are siblings this partition cannot separate, since every atom of a node lands in the node's group. Dealing atoms instead — cover_hybrid's integer form, see group_seeds — spreads them whether the program packs them into one primitive or spells them out.

declared_scale(x) -> sigma

The declared scale behind one slot's reference: the absolute sensitivity that stands in for a parent value carrying no scale of its own.

dependency_sweep(t::Tape, theta, theta2) -> (abar, vbar)

The staged form's propagation of the hard and squared channels, one reverse sweep each:

abar[j] = theta[j]  + sum_{k : j in P_k} magnitude(∂φ_k/∂v_j)   * abar[k]
vbar[j] = theta2[j] + sum_{k : j in P_k} magnitude(∂φ_k/∂v_j)^2 * vbar[k]

The hard channel is abs_reverse: seeded with theta from local_covers, its value at the inputs is the sum merge. The squared channel squares both the seeds and the coefficients, so it collects the same path contributions in quadrature; seeded with the matching theta2 from squared_seeds, sqrt.(vbar) at the inputs is the sign-averaged cover, an estimate of the quadrature merge from below. Both channels hold magnitudes: vbar here is the squared channel, not an adjoint.

cover_hessian forms both channels within its single fused sweep instead; this form seeds them independently.

dispatchcurvature(f, C::Curvature, ::Type{T})

Call f(C), branching on the built-in Curvature forms over the element type T first.

Tape stores its curvatures in a container with abstract element type, so an interface method invoked on one is a dynamic dispatch — once per node per sweep, which is a large fraction of the sweep's cost. Testing the built-in forms by isa first lets the compiler resolve those calls statically. A form defined elsewhere falls through to the same call, dispatched the ordinary way: this is an optimization for the common cases, not a restriction on which forms are allowed.

T is the element type of the tape the curvature came from, which is the one the elementals of overloads.jl record in; a DeclaredCurvature holds magnitudes, so its form is resolved over magtype(T). A curvature over some other element or matrix type is correct here, just not statically resolved.

C is @nospecialized because specializing on it would defeat the point: a caller that knows the concrete form would compile an instance in which every branch below but one is dead, having already paid the dispatch this ladder exists to avoid. One instance per call site and element type, branching at run time, is the intended shape.

group_seeds(t::Tape, rhos, groups::AbstractVector{<:Integer}) -> Theta
group_seeds(t::Tape, rhos, ngroups::Integer) -> Theta

Seed weights for the vector-mode reverse sweep, Theta[g, j] = sum over the contributions in group g of the weight each places on node j. Channels index the rows so that the weights of one contribution are contiguous.

The two forms differ in what a group holds. groups[k] assigns tape node k, keeping all of its atoms together, and maximum(groups) sets the channel count. An integer ngroups assigns atoms (see atom_slots): they are numbered in tape order and dealt round-robin, so a node exposing several atoms spreads them over several groups. The channel count is ngroups, capped at natoms(t) since further channels would stay empty.

Since a node's atoms partition its slots, both forms deposit the same total weight and differ only in how they distribute it over channels. Dealing atoms makes the distribution independent of how many atoms one elemental carries, so an array primitive that absorbs several atoms into one node still has them spread across channels. Which atoms share a channel remains a matter of the order they appear in.

With rhos from local_covers, column j collects exactly the seeds that abs_reverse propagates into the group partial sums. Theta is in boundtype of the tape's element type.

hessian_bound(t::Tape, vbar, utilde) -> B

The elementwise bound B = sum_k magnitude(vbar[k]) * sum_{i,j} magnitude(C^k[i,j]) * utilde_i * utilde_j', which satisfies B >= abs.(H) entrywise, in boundtype of the tape's element type. Diagnostic only: B is n × n, and the two inequalities it discards — cancellation along paths (inside each utilde) and cancellation across elementals — are exactly the looseness this bound trades for its cost.

Each node's term is accumulated by bound_contribution!, which every Curvature supplies in the form its structure allows. A DeclaredCurvature node contributes its declared bound s[i]*s[j] >= abs(C[i, j]), so B >= abs.(H) remains valid conditional on the declaration.

intermediate_gradients(t::Tape) -> u

The signed intermediate gradients: u[k][a] = ∂v_k/∂x_a, obtained by propagating u[k] = sum_{i in P_k} (∂φ_k/∂v_i) * u[i] from u[x_a] = e_a.

The signed counterpart of sensitivities, which runs the same recursion in absolute-value arithmetic and so returns a bound utilde >= abs.(u) rather than the derivatives themselves. The gap between the two is the cancellation that occurs along paths from an input to a node.

Costs O(n * L) time and memory. Only first partials enter, so a declared sub-map node (see submapnode) propagates exactly like any other.

jacobian_surrogate(t::Tape, outidxs::AbstractVector{<:Integer}) -> Utilde

The m × n surrogate of absolute sensitivities of the outputs outidxs with respect to the inputs of t: Utilde[k, j] ≥ |J_kj|, tight when no signed path cancellation occurs between output k and input j. Formed by one m-channel absolute-value reverse sweep seeded at the output nodes (see abs_reverse), so it costs O(m*L) — a diagnostic, like the Hessian bound B. A repeated entry in outidxs yields an identical row. The sweep runs in boundtype of the tape's element type and the result is returned in its magtype.

level_partition(t::Tape, ngroups::Integer) -> groups

Assign each node of t to one of ngroups groups by binning node_levels into contiguous depth bands. Nodes at the same depth always land in the same group, and when the tape is shallower than ngroups each depth gets a group of its own (leaving the surplus groups empty). Inputs carry no second-order contribution, so their group is irrelevant; they are assigned group 1.

This separates nodes along a chain but never separates siblings, which limits how much of the quadrature merge it can recover; cyclic_partition is the default for that reason.

live_nodes(t::Tape) -> live

Which nodes of t have an intermediate gradient that is not identically zero. live[k] is false exactly when every path from the inputs to node k passes through a recorded first partial that is exactly zero (whose magnitude is zero), and then ∇v_k is the zero vector in any arithmetic: every term of every path product carries an exact zero factor.

The test is on the recorded partials rather than on a propagated magnitude, so it neither underflows a live node to dead nor depends on the declared input scales.

Costs one O(L) forward pass, paid once per cover and only by a program carrying a block whose off-diagonal is active.

local_cover_slots(t::Tape{T}, k::Integer, av; check::Bool=true) -> Union{Nothing,Vector{M}}
local_cover_slots(ratio, t::Tape, k::Integer, av, scales; check::Bool=true)

The rank-1 cover of a single node, slot-aligned with its parents: nonnegative rho with rho[i]*rho[j] >= av*magnitude(C^k[i,j]) wherever both slots have a nonzero intermediate gradient (see dead_slot), and av the absolute adjoint magnitude(vbar[k]). rho and av are in M = magtype(T). Returns nothing when the node contributes no second-order term (av == 0, or a NoCurvature node). This is the per-node kernel shared by local_covers, which accumulates the results into (rhos, theta), and the CoverOnly sweep, which deposits them directly into the absolute adjoint channel. It reads the node's Curvature off the tape and forwards to cover_slots, which holds the construction for each form. check verifies the covering invariant on recorded nodes.

The first form takes the default split, resolving the threshold against the tape's element type. The second takes a resolved ratio together with scales = split_scales(w, t), both computed once by the caller and shared across nodes.

This is the tape's adapter onto cover_slots: it names the node with a TapeSite and gathers scales[t.parents[k]] into slot order, both lazily.

A node whose covering entries are not finite is reported by check_finite_partials rather than by the invariant test.

local_covers(t::Tape, vbar; check::Bool=true) -> (rhos, theta)

Per-node rank-1 covers. For node k, rhos[k][i] >= 0 satisfies rhos[k][i] * rhos[k][j] >= c^k[i,j] where c^k[i,j] = magnitude(vbar[k]) * magnitude(C^k[i,j]), except at entries a dead_slot empties; nodes with no second-order contribution get an empty rhos[k]. theta[j] = sum_k rhos[k][i] collects the weight each node receives as a parent, which seeds the reverse sweep of abs_reverse; it is in boundtype of the tape's element type, so the sweep it seeds rounds as the tape asks.

Each node's cover comes from cover_slots, dispatching on its Curvature. A DeclaredCurvature node (a declared sub-map output, see submapnode) brings its own rank-1 cover, and there is no recorded C to check it against, so check cannot audit those nodes — validity is conditional on the declaration.

scales declares the input scales the split falls back on and ratio the threshold at which it does so; see split_block and split_ratio. check verifies the covering invariant at every recorded node.

marginal_weights(t::Tape, outidxs::AbstractVector{<:Integer}, b::AbstractVector) -> a
marginal_weights(Utilde::AbstractMatrix, b::AbstractVector) -> a

The output weights a_k = ∑_j Ũ_kj / b_j — the row marginal of the surrogate weighted by the current input cover b. Each term of the sum already certifies its own entry (a_k b_j ≥ Ũ_kj), so (a, b) is a valid cover of the surrogate for any positive b.

This is the covariant output-weighting rule: under an output rescaling r → diag(c)*r it picks up exactly diag(c), and under an input rescaling x → diag(d)*x (with b → diag(d)⁻¹*b alongside) it is exactly invariant — each term Ũ_kj/b_j carries the units of r_k alone. The bare row marginal ∑_j Ũ_kj lacks the second property, and the value rule a_k = |r_k| degenerates wherever residuals vanish.

The tape method costs one O(L) abs_forward sweep seeded 1/b_j at the inputs, never forming Ũ. A zero b_j is accepted only for an input no output depends on: it contributes nothing to any marginal. A zero b_j on an input that does reach an output is an error — no output weighting can restore validity there.

maxpath_forward(t::Tape, d) -> TopTwo

Forward (max, ×) sweep with per-input seeds 1 / d[l]: for every node k, the two largest values of U[k, l] / d[l] over inputs l, with their inputs, where U[k, l] is the largest product of absolute first partials along any path from input l to node k. An input with d[l] == 0 seeds nothing.

O(L) time and memory.

maxpath_rowcorrelation(t::Tape, outidx::Integer, d; check=true, scales=nothing, ratio=nothing) -> m

The (max, ×) row statistic of the Hessian at outidx relative to the scales d:

m[i] = max over atoms of smax[i] * max_{l != i} smax[l] / d[l] / d[i],

where an atom's smax[l] = max_{slots} rho[slot] * U[parent, l] takes each slot along its largest path (U as in maxpath_forward) and the largest of its slots, and rho are the local covers of local_cover_slots under the same scales and ratio the cover sweeps use.

This is the max-merge counterpart of the largest normalized off-diagonal entry in row i of the quadrature merge's surrogate sum_atoms s_atom s_atom', and is at most that entry. It is exact for a sum of squared linear forms written with each input once per form. It is a statistic for starting a metric fit, not a bound, and not a cover. With d covariant (for instance a cover), m is invariant under rescaling of the inputs and of the output. A coordinate with d[i] == 0 neither seeds the forward sweep nor receives a value: m[i] = 0.

Two sweeps: the forward (max, ×) sweep, and one reverse sweep that carries the adjoint and the (max, ×) channel together, O(L) time and memory.

merge_atoms(t::Tape, rhos, u) -> (h_sum, h_quad)

Both merges of the cover atoms, from per-node cover slots rhos and an intermediate-gradient channel u. An atom's cover vector is s_atom = sum_{i in atom} rhos[k][i] * u[P_k[i]], and the merges are h_sum = sum_atoms s_atom and h_quad[i] = sqrt(sum_atoms s_atom[i]^2).

u[j] must be nonnegative and dominate abs.(∇v_j) entrywise for the results to cover; both sensitivities and the magnitudes of intermediate_gradients qualify. The merges are in boundtype of the tape's element type. Which slots share an atom comes from atom_slots, and the loop depends on its partitioning the slots: a slot in no atom would lose its contribution, one named twice would be counted twice.

Costs O(n * L) time and O(n) memory beyond u.

nls_partition(N, sqidx; ngroups=nothing) -> groups

Node partition for the residuals-as-groups hybrid merge of a least-squares loss tape with N nodes and residual-square nodes sqidx. Each residual square is placed in a group of its own (ngroups === nothing, one group per residual, giving the exact Gauss–Newton diagonal in quadrature) or spread cyclically across ngroups groups (a cheaper upper bound). Every other node — the residual-curvature subgraph — is placed in one further group, so its second-order contribution is summed there and combined in quadrature with the Gauss–Newton groups.

That last group is always occupied, so maximum(groups) - 1 recovers the number of residual groups even when ngroups exceeds the residual count.

node_levels(t::Tape) -> levels

Depth of each node: 0 for the inputs, and 1 + maximum(levels[parents]) otherwise. Parents always precede their children on the tape, so one forward pass suffices.

offdiag_active(c11, c12, c22) -> Bool

Whether the diagonal cover (sqrt(c11), sqrt(c22)) of the nonnegative block [c11 c12; c12 c22] falls short of its off-diagonal, i.e. whether c12^2 > c11*c22, in which case balanced_split is needed.

The test is c12 > sqrt(c11)*sqrt(c22), never the squared form: c12^2 underflows to zero for c12 below the square root of the smallest subnormal (about 1e-162 in Float64), which would misread a block with c11 == c22 == 0 and a tiny but nonzero cross term as diagonal and cover it with (0, 0). A long running product of small factors produces exactly such blocks. Square roots keep every representable entry inside the representable range, so this form neither underflows nor overflows.

quadrature_surrogate(t::Tape, outidxs::AbstractVector{<:Integer}) -> R

The m × n quadrature counterpart of jacobian_surrogate: R[k, j] is the 2-norm, over the paths from input j to output k, of the product of the local partial magnitudes along each path, R_kj = sqrt(∑_paths ∏ abs(d)^2). A 2-norm is at most the 1-norm, so R <= Utilde entrywise, with equality where an entry has a single path. If each parent slot's local partial carries an independent, equiprobable sign, R[k, j]^2 = E[J_kj^2]. This is the matrix balance_cover takes with merge = SignAveragedMerge(). Formed by one m-channel reverse sweep, O(m*L), that accumulates with hypot. No path product is squared, so an entry is accurate to rounding wherever the path products themselves are representable, as for Utilde.

sensitivities(t::Tape) -> utilde

Forward absolute sensitivities: utilde[k][a] >= abs(∂v_k/∂x_a), obtained by propagating utilde[k] = sum_{i in P_k} magnitude(∂φ_k/∂v_i) * utilde[i] from utilde[x_a] = e_a. The bound follows from the triangle inequality over paths, and is tight when no path cancellation occurs. The entries are in boundtype of the tape's element type.

Costs O(n * L) time and memory, so this is a diagnostic path: the sum merge is obtainable without it (see abs_reverse).

split_fallback(lo, hi, c12, z1, z2) -> rho1

The equalizing point balanced_split aims at, for blocks where sqrt(c12*z2/z1) misses it.

A pair with a ratio to read misses it only by overflow or underflow of the product under the root, and frexp recovers the point from the separated exponents: the exponent arithmetic is exact, and a rescaling by powers of two shifts e by an even amount, so this form is covariant to the bit exactly as the product form is. A pair with no ratio — a vanishing or non-finite reference — gets the geometric mean of the endpoints, the covariant, permutation-equivariant point of the feasible segment.

squared_seeds(t::Tape, rhos) -> theta2

Squared seed weights: theta2[j] sums rhos[k][i]^2 over every node k and parent slot i with t.parents[k][i] == j. This is the squared-channel counterpart of the linear seeds theta that local_covers returns, and it seeds the squared channel of dependency_sweep.

squared_surrogate(t::Tape, outidxs::AbstractVector{<:Integer}) -> U2

The squared channel U2 = R.^2 of quadrature_surrogate: U2[k, j] sums, over the paths from input j to output k, the squared product of the local partial magnitudes along the path, so U2 <= Utilde.^2 entrywise and, under independent equiprobable signs, U2[k, j] = E[J_kj^2]. Its range is that of a squared Jacobian entry: U2[k, j] loses precision once R[k, j] falls below sqrt(floatmin(T)), and underflows to zero below about the square root of the smallest subnormal, while R itself stays accurate to rounding. The sign-averaged pair is computed on R; this matrix is for inspection.

tape_potentials(t::Tape, outidxs::AbstractVector{<:Integer}; root=nothing, cgiters=0)
    -> (; theta, live, roots, nties)

Log-scale potentials on the live tape graph: one number per node, least-squares-fitting the slot equations θ_child - θ_parent = log|d|. exp(θ) at the output nodes is the start balance_cover uses for start = :covariant on a tape, because rescaling a node's value by s shifts its potential by log s and leaves every difference within a component covariant.

A node is live when it lies on some input-to-output path whose local derivatives are all nonzero; live is that mask over the tape's nodes, and theta is NaN outside it. The live nodes and slots form an undirected graph (parallel slots counted separately) whose components are rooted at their node of maximum undirected degree, lowest node index on ties. roots lists those node indices and nties[c] how many nodes of component c share the root's degree, so a caller can see when the root came down to the index rather than the graph. root overrides the rule with a node index or a vector of them, at most one per component and each live; components with no supplied root keep the default. A supplied root breaks no tie, so its component reports nties[c] = 1.

theta is the breadth-first layer average from each root: θ = 0 at the root, and a node in layer d takes the mean over the slots joining it to layer d-1 of θ_u ± log|d|, + when it is the slot's child. cgiters > 0 refines that with up to this many conjugate-gradient iterations on the normal equations, stopping earlier at a residual floor or at the Laplacian's rank, both of which are invariant under rescaling; past those the iteration only amplifies roundoff along the nullspace. The exact least-squares potentials are independent of the root up to one additive constant per component.

Costs O(L) for cgiters = 0 and O(L) per refinement iteration, with one logarithm per live slot and no derivative evaluations. The fit weights tape slots equally, so like the surrogate itself it depends on how the program is written.

Examples

julia> r(x) = [2 * x[1] * x[2], 4 * x[2] * x[3]];

julia> pot = DerivativeCovers.tape_potentials(tape_function(r, [1.0, 1.0, 1.0])...);

julia> pot.nties, all(isnan, pot.theta[.!pot.live])
([5], true)

value_bounds_split(c11, c12, c22, z1, z2, ratio) -> Bool

Whether the value pair (z1, z2) fixes the split of the block [c11 c12; c12 c22].

balanced_split aims at sqrt(c12*z2/z1) and clamps it to the feasible segment [sqrt(c11), c12/sqrt(c22)], so a pair carrying no scale still yields a covariant split wherever the endpoint it runs toward is positive and finite: the clamp supplies the answer, and both endpoints are covariant. What is left is the pair of cases where that endpoint is degenerate as well — c22 == 0 puts the upper endpoint at infinity, c11 == 0 puts the lower one at zero — together with two vanishing values, which leave the ratio z2/z1 undefined whatever the block. ratio is how lopsided a pair must be to count as carrying no scale.

The test compares z1 and z2 directly, which is not invariant when the two parents carry different units, so it is applied only where no covariant split exists. In particular a block whose bounding diagonal entry is proportional to the vanishing value — a quotient's numerator, a power's exponent — has its split fixed by the values until the value is exactly zero, and that split grows without bound on the approach; see balanced_split.

worst_severity(severity::AbstractArray) -> (s, idx)

Largest entry of severity and where it occurs, treating NaN as passing.

Checks over a whole array score every entry and reduce, rather than testing entries in turn and raising at the first failure: the score is an array computation and only its reduction is inspected, so the decision to raise is taken once, by the caller, on a single scalar. Each caller normalizes its severity so that entries <= 1 pass, which also makes the reported number the factor by which the worst entry misses. A NaN entry comes from a vacuous comparison — zero against zero — and scores as 0.

gramcover!(s, t::Tape, outidxs) -> s
gramcover!(s, t::Tape, outidxs, w::AbstractVector) -> s
gramcover!(s, t::Tape, outidxs, W::Diagonal) -> s
gramcover!(s, t::Tape, outidxs, W::AbstractMatrix; v=<symmetric cover of abs.(W)>) -> s

Write gramcover's cover of the Gram matrix of the tape outputs outidxs into s and return it. s must be indexed like the inputs (eachindex(s) == 1:t.ninputs).

gramcover(t::Tape, outidxs) -> s
gramcover(t::Tape, outidxs, w::AbstractVector) -> s
gramcover(t::Tape, outidxs, W::Diagonal) -> s
gramcover(t::Tape, outidxs, W::AbstractMatrix; v=<symmetric cover of abs.(W)>) -> s

Symmetric cover of the Gram matrix G = J'*W*J of the Jacobian J of the tape outputs outidxs: a nonnegative vector s, one entry per input (1:t.ninputs), with abs(G[i,j]) <= s[i]*s[j] for every pair. Neither J nor G is formed. W defaults to the identity (G = J'*J, the Gauss–Newton matrix), and a vector w means W = Diagonal(w); the weights may have either sign. w, v, and both axes of W are indexed like outidxs, one entry per output.

No cover of J is needed: rescaling the outputs by a diagonal C together with W -> inv(C)*W*inv(C) leaves G unchanged, so the weight alone fixes the output scale.

For any v >= 0 with abs(W[k,l]) <= v[k]*v[l], the result is s = Ũ'*v, where Ũ >= abs.(J) is the absolute-sensitivity surrogate (jacobian_surrogate), and abs(G[i,j]) <= ∑_{k,l} Ũ[k,i]*abs(W[k,l])*Ũ[l,j] <= s[i]*s[j]. For diagonal weights v = sqrt.(abs.(w)) (all ones when unweighted), so an output with zero weight contributes nothing. For a general W the default v is MatrixCovers.symcover of max.(abs.(W), abs.(transpose(W))), which also admits an asymmetric W; a supplied v must be finite, nonnegative, and satisfy abs(W[k,l]) <= v[k]*v[l]. A repeated entry of outidxs is a repeated row of J, and its seeds add.

s comes from one absolute-value reverse sweep seeded v[k] at output node outidxs[k], so it costs O(L) for a tape of length L and never forms Ũ. It sums the weighted rows of Ũ, and for w >= 0 (or no weight) it brackets the per-output quadrature h_gn[j] = sqrt(∑_k w[k]*Ũ[k,j]^2):

h_gn[j] <= s[j] <= sqrt(n_j) * h_gn[j],     n_j = count(>(0), Ũ[:, j]),

where n_j is the number of outputs that input j reaches. h_gn is the minimal symmetric cover of Ũ'*W*Ũ and costs one sweep channel per output, O(m*L). Unweighted, it is the h_gn of cover_nls with ngroups=nothing; with Ũ in hand from a WithDiagnostics result, a diagonal weight gives it as MatrixCovers.gramcover(Utilde, w).

On a tape whose sweeps round upward (see boundtype), such as one recorded over a box of intervals, the seeds and the sweep round upward, and v is raised by the few units in the last place needed for abs(W[k,l]) <= v[k]*v[l] to hold exactly, so s covers G at every point of the box. On any other tape, s carries the roundoff of the sweep without inflation, like the package's other sweep-based covers, so the bound holds to a relative error of a few eps per node on the longest path.

Examples

julia> r(x) = [x[1] * x[2], exp(x[1]) - x[3], x[2] + 2x[3]];

julia> t, outidxs = tape_function(r, [0.5, 2.0, -1.0]);

julia> round.(DerivativeCovers.gramcover(t, outidxs); digits=4)   # column sums of abs.(J)
3-element Vector{Float64}:
 3.6487
 1.5
 3.0

julia> round.(DerivativeCovers.gramcover(t, outidxs, [1.0, 0.0, 4.0]); digits=4)   # output 2 dropped
3-element Vector{Float64}:
 2.0
 2.5
 4.0

support_components(t::Tape, outidxs) -> sc::SupportComponents

Connected components of the bipartite support graph of the Jacobian of the outputs outidxs with respect to the tape's inputs: outputs and inputs are the vertices, and each nonzero Ũ_kj is an edge joining output k to input j. The result is a MatrixCovers.SupportComponents whose row axis runs over the outputs in the order of outidxs (1:length(outidxs)) and whose column axis runs over the inputs (1:t.ninputs); query it with rowcomponent, colcomponent, and ncomponents.

This is the tape counterpart of MatrixCovers.support_components(Ũ): it costs O(L) in the tape length and never forms Ũ. A tape edge transmits sensitivity exactly where its first partial is nonzero, so the surrogate's support is carried by those edges alone, and a node is live when it both depends on an input and reaches one of outidxs along them. Any undirected path between live tape nodes splits into ascending and descending runs whose nodes each reach some output, which forces the corresponding bipartite vertices into one component; conversely a support edge is a directed path of live nodes. The undirected components of the live tape graph therefore restrict to the bipartite components exactly.