Computing and using covers

DerivativeCovers

Compute certified rank-1 bounds on Hessians and Jacobians without forming either matrix.

A Hessian cover is a nonnegative vector h satisfying abs(H[i,j]) <= h[i] * h[j] for every entry of H. cover_hessian records a scalar function on a Tape and computes such a cover, together with the gradient, in O(L) work for a tape of length L. The result is deterministic and certified.

A cover is covariant under rescaling: if x -> D*x, then h -> inv(D)*h; if f -> c*f, then h -> sqrt(c)*h. That holds exactly on any program whose blocks all split on parent values carrying a scale of their own. Where one does not, the split falls back on declared input scales and stays covariant when those are rescaled alongside, scales -> D*scales; guards_any_block says whether a given program is in that case.

using DerivativeCovers imports the main entry points and result types. Public utilities such as equilibrate are not exported; call them with the DerivativeCovers. prefix.

cover_hessian([mode], t::Tape, outidx::Integer; check=true, scales=nothing, ratio=nothing)
cover_hessian([mode], f, x0::AbstractVector; check=true, scales=nothing, ratio=nothing)
cover_hessian(CoverOnly(), t::Tape, vbar::AbstractVector; check=true, scales=nothing, ratio=nothing)

Compute a certified rank-1 cover of a scalar function's Hessian.

Pass either a tape and scalar output node, or a function and evaluation point. The default CoverOnly mode returns a HessianCover containing the gradient, the sum-merge cover, and the sign-averaged cover. WithDiagnostics additionally computes an explicit elementwise bound and a tighter quadrature cover.

The scales keyword declares the input scales the split falls back on where a block's parent values carry none, and ratio the threshold at which it does so; see split_block and split_ratio. Neither affects validity, only which of the valid covers is chosen, and guards_any_block says whether a program consults them at all. check=true verifies the local covering conditions.

Two merges are available:

• sum merge, h_sum = sum_k s_k: the cross terms of (sum_k s_k)(sum_k s_k)' are nonnegative slack;
• quadrature merge, h_quad[i] = sqrt(sum_atoms s_atom[i]^2): valid by Cauchy–Schwarz, and the minimal cover of the no-cancellation surrogate M = sum_atoms s_atom s_atom'. An atom (see atom_slots) is a whole node for unary and binary elementals and one parent pair for an inner-product node, so an inner-product node contributes q separate terms to the quadrature merge rather than one.

Both satisfy abs.(H) <= h * h'.

The gradient is in the tape's element type T and the covers in magtype(T). T may be a type whose values stand for sets, such as Interval{Float64} from IntervalArithmetic.jl: recording f over a box x0 = interval.(lo, hi) then gives covers (here Vector{Float64}) that hold at every point of the box, and a gradient whose entries enclose the gradient over the box. The sweeps over such a tape round every sum, product and square root upward (see boundtype), so the covers hold exactly; on a tape of points they hold up to the rounding of the sweeps. Control flow in f that the box does not decide, such as a comparison between overlapping intervals, raises an error, as does an elemental evaluated outside its domain somewhere on the box (see check_values).

CoverOnly()

This mode takes one forward pass and one reverse sweep, or O(L) time and memory. The reverse sweep carries the ordinary adjoint, the absolute cover channel, and the squared channel together:

vbar[p_i] += d_i * vbar[k]                              # signed channel: reverse-mode AD
abar[p_i] += magnitude(d_i) * abar[k] + rho_k[i]        # absolute channel: the cover
qbar[p_i] += magnitude(d_i)^2 * qbar[k] + rho_k[i]^2    # squared channel: h_var^2

with rho_k from local_cover_slots. Every channel of node k is final when the sweep reaches it: adjoint flow and theta seeds come only from k's consumers, which follow it on the tape and are processed first. That makes rho_k — which needs the finished vbar[k] — computable in place, and abar[1:ninputs] equal to abs_reverse(t, theta)[1:ninputs] up to roundoff.

abar sums, over every path from a parent slot of a node to an input, the slot's seed times the product of the absolute first partials along the path. qbar sums the squares of the same path contributions, and h_var = sqrt.(qbar[1:ninputs]) is the sign-averaged cover of HessianCover, which does not cover the Hessian.

A program that reaches the declared-scale fallback requires one additional forward sweep.

Passing the node adjoints vbar =adjoints(t, outidx) in place of outidx runs the same sweep without the signed channel, for a caller that already holds the adjoints of a reverse pass. It returns the same HessianCover bit for bit, with grad = vbar[1:ninputs].

WithDiagnostics()

This mode returns a HessianDiagnostics and costs O(n * L) time and O(n^2) memory. With check=true, it also verifies that the diagnostic and fused computations agree.

Examples

julia> f(x) = exp(x[1]) + x[1] * x[2]^2;

julia> c = cover_hessian(f, [0.5, 2.0]);

julia> round.(c.h_sum; digits=4)
2-element Vector{Float64}:
 4.1125
 2.4142

julia> round.(c.grad; digits=4)
2-element Vector{Float64}:
 5.6487
 2.0

julia> H = [exp(0.5) 4.0; 4.0 1.0];   # the Hessian of f at [0.5, 2.0]

julia> all(abs.(H) .<= c.h_sum * c.h_sum')
true

cover_hybrid(t::Tape, outidx::Integer, ngroups::Integer; check::Bool=true, scales=nothing, ratio=nothing) -> h
cover_hybrid(t::Tape, outidx::Integer, groups::AbstractVector{<:Integer}; check::Bool=true, scales=nothing, ratio=nothing) -> h
cover_hybrid(f, x0::AbstractVector, spec; check::Bool=true, scales=nothing, ratio=nothing) -> h

Compute a Hessian cover by grouping cover contributions: contributions are summed within each group, then the group totals are combined in quadrature.

An integer ngroups groups the tape's cover atoms (see atom_slots), dealing them round-robin. ngroups=1 is the sum merge, and any ngroups at or above natoms(t) gives every atom a group of its own, which is exactly the quadrature merge h_quad of diagnostic cover_hessian. Because atoms rather than nodes are dealt, the whole range between the two merges stays reachable however the program was written: an elemental that carries several atoms — an inner product, a matrix-vector product — spreads them across groups like the scalar chain it stands for.

Alternatively, groups[k] assigns tape node k to a group, holding all of that node's atoms together; cyclic_partition and level_partition build such partitions.

Every grouping produces a valid cover, and one group is always the sum merge. The computation costs O(G * L) time and memory for G groups. scales and ratio govern the split, as in cover_hessian, and check=true verifies local covers.

Examples

The two endpoints, on a tape carrying three atoms:

julia> f(x) = exp(x[1]) + x[1] * x[2]^2;

julia> round.(DerivativeCovers.cover_hybrid(f, [0.5, 2.0], 1); digits=4)   # one group: the sum merge
2-element Vector{Float64}:
 4.1125
 2.4142

julia> round.(DerivativeCovers.cover_hybrid(f, [0.5, 2.0], 3); digits=4)   # a group per atom: the quadrature merge
2-element Vector{Float64}:
 3.1062
 1.7321

forward_cover(f, x0::AbstractVector; scales=nothing, ratio=nothing, check=true) -> ForwardCover
forward_cover(x::FVar) -> ForwardCover

Compute a certified rank-1 cover of a scalar function's Hessian by propagating it forward, in one pass and with no tape.

f must accept a vector of FVars and return one, using operations FVar supports. Per variable the propagation carries the value, the intermediate gradient ∇v, a squared cover of ∇²v, and the scalar declared scale sigma; the recursions are

∇v_k    = sum_i d_i * ∇v_i
sigma_k = sum_i magnitude(d_i) * sigma_i
G_k     = sum_i magnitude(d_i) * G_i + sum over atoms of (sum_{i in atom} rho_i * magnitude.(∇v_i))^2

with rho the same local cover a reverse sweep builds, except that the adjoint factor is 1: a forward pass cannot know how the output depends on the node it is at. The cover is sqrt.(G).

The gradient is in the element type T = float(eltype(x0)) and the cover in magtype(T). As for cover_hessian, x0 may be a box of intervals: the gradient then encloses the gradient over the box and the cover holds at every point of it, the propagation rounding its bounds upward (see boundtype). An elemental evaluated outside its domain somewhere on the box is an error (see check_values).

The second form returns the cover of any variable reached along the way, the propagation having maintained the invariant at every one of them, not only at the output.

Computational cost and tightness

O(n*L) time for n inputs and L operations, and O(n) memory per live variable — no tape, and nothing proportional to the length of the program. Its main advantage is therefore the ability to process programs that are too long to record.

The channels are stored in whatever x0 is made of, so a StaticVector point propagates in static arrays and no node allocates at all. That is worth doing for a handful of inputs and not for many: the propagation is then compiled for the exact input count, which for large n costs more in compilation than it saves at run time. A Vector point propagates in Vectors and compiles once for every size.

The cover keeps the signs of the intermediate gradients ∇v, which a forward pass carries exactly, and discards the signs of the adjoints, which it cannot see: each atom is weighted by the unsigned sum over paths from its node to the output rather than by the node's adjoint. That makes it incomparable with the quadrature merge h_quad of HessianDiagnostics, which discards the signs of the intermediate gradients and keeps those of the adjoints; which is tighter is a property of the program.

scales and ratio govern the split exactly as in cover_hessian: the propagation carries both of the references split_block reads, so it splits by the same rule a tape does and its cover is comparable with the tape's coordinate by coordinate. See forward_scale. check=true verifies every local cover.

Examples

julia> f(x) = exp(x[1]) + x[1] * x[2]^2;

julia> c = forward_cover(f, [0.5, 2.0]);

julia> round.(c.grad; digits=4)
2-element Vector{Float64}:
 5.6487
 2.0

julia> H = [exp(0.5) 4.0; 4.0 1.0];   # the Hessian of f at [0.5, 2.0]

julia> all(abs.(H) .<= c.h_quad * c.h_quad')
true

CoverOnly()

Return only the covers a sweep-tier computation yields: the gradient, sum-merge cover, and sign-averaged cover from cover_hessian at O(L), or the gradient, two covers, and two scales from cover_nls at O(G * L) for G residual channels, O(L) by default. This is the default mode of both.

WithDiagnostics()

Return the extended result, adding the quantities that must be materialized per input or per output: a HessianDiagnostics from cover_hessian, with the elementwise bound, quadrature merge, and an independent consistency check, at O(n * L) time and O(n^2) memory; an NLSDiagnostics from cover_nls, with the residual-sensitivity surrogate, at O(m * L) time and O(m * n) memory.

cover_nls([mode], r, x0; ngroups=1, check=true, scales=nothing, ratio=nothing)

Compute covers for the least-squares loss f(x) = ½∑ₖ rₖ(x)² at x0.

r must map a vector of TVars to a vector of residuals. The default CoverOnly mode returns an NLSCover holding a cover h of the full loss Hessian, a cover h_gn of the Gauss–Newton matrix JᵀJ, and the sign-averaged covers h_var <= h of the full Hessian and h_var_gn <= h_gn of JᵀJ (the latter from gramvariance), which cover the mean magnitudes of a sign ensemble rather than these matrices. WithDiagnostics returns an NLSDiagnostics, which adds the absolute-sensitivity surrogate Utilde for a further O(m * L).

ngroups sets the number of residual channels G, and with it the entire cost of the CoverOnly mode: O(G * L) time and O(G * n) memory. The residuals are spread cyclically over the G channels, and their sensitivities are summed within a channel before the quadrature merge. The default 1 sums them all, so h_gn is the gramcover of the residuals and both covers cost O(L). ngroups=nothing gives every residual a channel of its own (G = m), the tightest pair of covers available, at O(m * L); its h_gn[j] is smaller than the default's by a factor of at most sqrt(n_j), where n_j is the number of residuals input j reaches. That factor varies with j, so it changes the profile of the cover as well as its level. The two sign-averaged covers take one further O(L) sweep and do not depend on G.

scales and ratio govern the split, as in cover_hessian, and reach only the residual-curvature channel, and through it h and h_var: the residual channels carry unit atoms, which admit no split, so none of h_gn, h_var_gn and Utilde depends on them. check=true verifies the local covering conditions, and under WithDiagnostics also verifies h_gn against Utilde.

Examples

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

julia> nc = cover_nls(r, [0.5, 2.0]);

julia> round.(nc.h; digits=4)
2-element Vector{Float64}:
 3.7272
 0.5

julia> round.(nc.h_gn; digits=4)
2-element Vector{Float64}:
 3.6487
 0.5

julia> all(nc.h_gn .<= nc.h)
true

julia> round.(nc.h_var; digits=4)
2-element Vector{Float64}:
 2.7014
 0.5

julia> all(nc.h_var_gn .<= nc.h_var .<= nc.h)
true

gramvariance(t::Tape, outidxs) -> h
gramvariance(t::Tape, outidxs, w::AbstractVector) -> h
gramvariance(t::Tape, outidxs, W::Diagonal) -> h

Sign-averaged cover, the root-mean-square column scale, of the Gram matrix G = J'*W*J of the Jacobian J of the tape outputs outidxs: a nonnegative vector h, one entry per input (1:t.ninputs), with

h[j] = sqrt(∑_k abs(w[k]) * Ũ⁽²⁾[k,j]),

where Ũ⁽²⁾[k,j] sums, over the paths from input j to output k, the squared product of the local partial magnitudes along the path. W defaults to the identity, and a vector w means W = Diagonal(w), indexed like outidxs; only abs.(w) enters, so an output with zero weight contributes nothing, and a repeated entry of outidxs adds its weights.

Model each parent slot's local partial as carrying an independent, equiprobable sign (the sign ensemble). Distinct paths then contribute uncorrelated terms, so E[J[k,j]^2] = Ũ⁽²⁾[k,j] and h[j]^2 = E[(J'*abs.(W)*J)[j,j]]: h is the root-mean-square column norm of sqrt.(abs.(W))*J, the column (MINPACK) scaling of the Gram matrix in expectation. The entries of J for distinct outputs are uncorrelated too, so off-diagonal weights contribute nothing to this expectation; a general W::AbstractMatrix is an error. On an interval tape the magnitudes are the largest over the box, so h bounds the pointwise h at every point of the box.

h covers the ensemble's mean magnitudes, E[abs(G[i,j])] <= h[i]*h[j], not G itself: where several paths reach one entry of J with the same sign it can fall below that column's norm. Per entry Ũ⁽²⁾ <= Ũ.^2 (a sum of squares is at most the square of the sum), where Ũ is the absolute-sensitivity surrogate (jacobian_surrogate), so

h[j] <= h_gn[j] = sqrt(∑_k abs(w[k]) * Ũ[k,j]^2) <= gramcover(t, outidxs, w)[j].

Unweighted, h_gn is the Gauss–Newton cover of cover_nls with ngroups=nothing. Where each entry of J has a single path, Ũ⁽²⁾ = Ũ.^2 = abs2.(J) on a point tape, and h equals h_gn and the column norms of sqrt.(abs.(W))*J. Rescaling the inputs, x -> D*x, gives h -> inv(D)*h; rescaling the outputs by a diagonal C together with w -> w ./ C.^2 leaves h unchanged.

h comes from one reverse sweep with squared local partials (the squared channel of dependency_sweep, without its linear channel) seeded abs(w[k]) at output node outidxs[k], so it costs O(L) for a tape of length L, where h_gn needs one channel per output. gramvariance! writes h into a supplied vector.

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.gramvariance(t, outidxs); digits=4)   # column norms of J
3-element Vector{Float64}:
 2.592
 1.118
 2.2361

julia> round.(DerivativeCovers.gramvariance(t, outidxs, [1.0, 0.0, 4.0]); digits=4)
3-element Vector{Float64}:
 2.0
 2.0616
 4.0

julia> t2, o2 = tape_function(x -> [x[1] * x[1]], [3.0]);   # two paths, J = [6.0;;]

julia> round.(DerivativeCovers.gramvariance(t2, o2); digits=4)   # sqrt(3^2 + 3^2)
1-element Vector{Float64}:
 4.2426

gramvariance!(h, t::Tape, outidxs) -> h
gramvariance!(h, t::Tape, outidxs, w::AbstractVector) -> h
gramvariance!(h, t::Tape, outidxs, W::Diagonal) -> h

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

nls_loss_tape(r, x0) -> (t, outidx, residx, sqidx)

Record the least-squares loss f = ½∑ₖ rₖ(x)² at x0, where r maps a vector of TVars to a vector of residual TVars. Returns the Tape, the loss output node index, the node indices residx of the residuals rₖ, and the node indices sqidx of their squares rₖ².

r may return residuals that are input or duplicated nodes; each is squared into its own fresh node, so sqidx has one distinct entry per residual regardless. A constant residual — a fixed offset of the loss, contributing ½c² to f and a zero row to J — gets a node of its own (constantnode); an all-constant residual vector depends on no parameter and is an error, as is a recording that check_values rejects.

gn_ab(nls::NLSDiagnostics) -> (a, b)

A valid rank-1 Jacobian cover (a, b) — a_k b_j ≥ |J_{kj}| — read off the residual sensitivities of nls: a = ones(m) and b[j] = maximum(Utilde[:, j]), so a_k b_j = b_j ≥ ũ_{rₖ,j} ≥ |J_{kj}|. b is the max merge for unit output weights: the tightest cover of J whose row scale the unweighted least-squares loss fixes at one. It serves to cross-check a general asymmetric Jacobian cover; the symmetric Gauss–Newton cover of JᵀJ is nls.h_gn.

cover_jacobian([mode,] t::Tape, outidxs; start=:covariant, rounds=4, a=nothing, check=true)
cover_jacobian([mode,] r, x0::AbstractVector; ...)

Asymmetric rank-1 cover of a vector-valued function's Jacobian: weights with a[k] * b[j] >= abs(J[k,j]). Pass either a tape and its output nodes, or a function and an evaluation point.

a is chosen by the balancing construction: start is :covariant, whose products on the support transform exactly under a rescaling of the program's inputs, outputs and intermediates, or :equal, which is permutation equivariant instead; rounds alternation rounds then refine it, and rounds = nothing runs balance_cover's alternation to convergence instead. Passing a yourself fixes the weights and rules out start and rounds. a must be strictly positive on every output some input reaches and may be zero on the others, which then contribute to no merge.

Modes

• CoverOnly() (the default) returns a JacobianCover from O(L) sweeps, never forming the surrogate. Checking a pair against the surrogate is what the other mode is for, so check has no effect here; a factor that vanishes on a row or column the program actually uses is an error either way.
• WithDiagnostics() returns a JacobianDiagnostics: the same a and b_sum, plus the surrogate and the two tighter merges, at O(m * L). check then verifies b_max against the surrogate.

Examples

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

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

julia> round.(jc.b_max; digits=4)
3-element Vector{Float64}:
 2.0
 1.0
 2.0

julia> all(jc.b_max .<= jc.b_quad .<= jc.b_sum)
true

cover_jacobian_hybrid(t::Tape, outidxs, groups::AbstractVector{<:Integer}; kwargs...) -> (; a, b)
cover_jacobian_hybrid(t::Tape, outidxs, ngroups::Integer; kwargs...) -> (; a, b)
cover_jacobian_hybrid(r, x0::AbstractVector, groups | ngroups; kwargs...) -> (; a, b)

Compute a Jacobian cover by grouping outputs.

groups[k] assigns output k to a group. Absolute sensitivities are summed within each group, then the group totals are combined with max. Alternatively, pass ngroups to assign outputs cyclically.

Every grouping is valid. One group reproduces b_sum; one group per output reproduces b_max. The cost is O(G * L) for G = maximum(groups), without forming the full surrogate.

The output weights are cover_jacobian's, chosen by the same keywords:

• start: :covariant (the default) or :equal.
• rounds: alternation rounds refining them, 4 by default.
• a: fixed weights of your own instead, strictly positive on every output some input reaches and optionally zero on the others. It rules out start and rounds.
• check = true: verify the returned b against the surrogate, at O(m * L).

Returns the weights used and the grouped merge for them: a[k] * b[j] >= abs(J[k,j]).

Examples

The two endpoints, on a program with three outputs:

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.cover_jacobian_hybrid(t, outidxs, 1; a = ones(3)).b; digits=4)   # one group: b_sum
3-element Vector{Float64}:
 3.6487
 1.5
 3.0

julia> round.(DerivativeCovers.cover_jacobian_hybrid(t, outidxs, 3; a = ones(3)).b; digits=4)   # per output: b_max
3-element Vector{Float64}:
 2.0
 1.0
 2.0

HessianCover{T, M}

The standard result returned by cover_hessian. T is the element type of the Tape it was computed from, and M = magtype(T) the type of the cover. HessianCover{T} stands for HessianCover{T, magtype(T)}.

• grad::Vector{T}: gradient of the scalar function at the evaluation point.
• h_sum::Vector{M}: certified Hessian cover computed with the sum merge.
• h_var::Vector{M}: the sign-averaged cover, the root-sum-square of the path contributions that h_sum adds. It covers the mean magnitude of the sign ensemble described below, not the Hessian.

h_sum[i] adds one contribution per pair of a seeded parent slot and a path from that slot to input i: the slot's seed times the product of the absolute first partials along the path. h_var[i]^2 sums the squares of the same contributions. That is the expected square of their sum when each parent slot carries an independent, equiprobable sign, on its first partial and on its seed: the sign ensemble, the model gramvariance uses. The estimates satisfy h_var <= h_quad <= h_sum coordinatewise, with h_quad the quadrature merge of HessianDiagnostics. A coordinate reached by a single contribution has h_sum == h_var; one built by summing many overlapping contributions has h_sum grow like their sum and h_var like their root-sum-square, so the ratio h_sum ./ h_var grows with the effective overlap.

A variable that enters one elemental through two slots is two paths whose contributions combine in quadrature, although in the program they share a sign: x[1] * x[1] seeds each of its two slots with 1, so h_sum = 2 adds them and h_var = √2, while x[1]^2, a unary elemental with a single slot, gives √2 for both. Where such paths reconverge further along, h_var falls below the exact scale: for x[1] * x[1] * x[1] at x = 3, h_var = √15 against sqrt(abs(f''(3))) = √18.

Examples

julia> cover_hessian(x -> x[1] * x[1], [3.0]).h_sum   # two slots, two paths
1-element Vector{Float64}:
 2.0

julia> cover_hessian(x -> x[1] * x[1], [3.0]).h_var ≈ [sqrt(2)]
true

julia> cover_hessian(x -> x[1] * x[1] * x[1], [3.0]).h_var ≈ [sqrt(15)]   # |f''| = 18
true

HessianDiagnostics{T, M}

The extended result returned by cover_hessian(WithDiagnostics(), ...). T is the element type of the Tape it was computed from, which grad holds, and M = magtype(T) the type of the bound and the covers. HessianDiagnostics{T} stands for HessianDiagnostics{T, magtype(T)}.

• grad: gradient of the scalar function.
• B: explicit elementwise bound satisfying B >= abs.(H).
• h_sum: cover computed with the sum merge.
• h_quad: tighter cover computed by combining cover atoms in quadrature.
• h_sum_rev: the absolute-value reverse sweep's own route to the sum merge, a consistency check.

Forming these diagnostics costs O(n * L) time and O(n^2) memory.

ForwardCover{T, M, VG, VH}

The result of forward_cover, with T as the propagation's element type, M = magtype(T) the type of the cover, and VG, VH the channel storage the propagation ran in (see ForwardContext). ForwardCover{T} stands for ForwardCover{T, magtype(T)}.

• grad::VG: the exact gradient of the output, in T.
• h_quad::VH: a certified Hessian cover formed by merging the cover atoms in quadrature, in M.

NLSCover{T, M}

Result of cover_nls for the least-squares loss f = ½∑ₖ rₖ². T is the element type of the Tape it was computed from, and M = magtype(T) the type of the covers. NLSCover{T} stands for NLSCover{T, magtype(T)}.

• grad::Vector{T}: gradient of the loss.
• h::Vector{M}: certified cover of the full Hessian JᵀJ + ∑ₖ rₖ ∇²rₖ.
• h_gn::Vector{M}: certified cover of the Gauss–Newton matrix JᵀJ; always h_gn <= h.
• h_var::Vector{M}: sign-averaged cover of the full Hessian, its root-mean-square column scale under random local signs, h_var² = h_var_gn² + c_var², where c_var is the sign-averaged form of the residual-curvature channel (one independent sign per parent slot, as for h_var of cover_hessian). It covers the sign ensemble's mean magnitude, not the Hessian; always h_var_gn <= h_var <= h.
• h_var_gn::Vector{M}: sign-averaged cover of JᵀJ, its root-mean-square column scale under random local signs, gramvariance of the residuals. It covers the sign ensemble's mean magnitude, not JᵀJ; always h_var_gn <= h_gn.
• m: the number of residuals.

Both covers are read off one grouped reverse sweep with G residual channels and one residual-curvature channel, so together they cost O(G * L) time and memory. With G = 1, the default, h_gn is the sum ∑ₖ ũ. With G = m, one residual per channel, h_gn is the per-residual quadrature √(∑ₖ ũ²), which equals sqrt.(diag(J'J)) when no signed paths cancel. The two scales share one further two-channel O(L) sweep through squared partials, whatever G is, so neither depends on G. h_var_gn equals h_gn (for G = m) where each entry of J has a single path, and h_var equals h_var_gn where the residuals have no curvature. h_var is the h_var of cover_hessian applied to the loss tape of nls_loss_tape.

Rescaling the inputs, x -> D*x with D diagonal, maps every channel h -> inv(D)*h on a program without guarded elementals; scales and ratio reach h and h_var only through the residual-curvature channel.

NLSDiagnostics{T, M}

The extended result returned by cover_nls(WithDiagnostics(), ...). T is the element type of the Tape it was computed from, which grad holds, and M = magtype(T) the type of the covers and of Utilde. NLSDiagnostics{T} stands for NLSDiagnostics{T, magtype(T)}. It carries the fields of NLSCover together with

• Utilde::Matrix{M}: m × n matrix of absolute residual sensitivities, satisfying Utilde[k,j] >= abs(J[k,j]).

Forming Utilde takes an m-channel reverse sweep of its own, O(m * L) time and O(m * n) memory, whatever number of residual channels the covers use.

JacobianCover{T}

Result of cover_jacobian: a[k] * b_sum[j] >= abs(J[k,j]). T is magtype of the element type of the Tape it was computed from, which for a real floating-point tape is that element type itself.

• a: output weights.
• b_sum: input weights from the sum merge, computed in one O(L) reverse sweep.
• roots, nties: the root data of the tape-potential start (see tape_potentials), or nothing when another start or an explicit a was used.
• m, n: numbers of outputs and inputs.

JacobianDiagnostics{T}

The extended result returned by cover_jacobian(WithDiagnostics(), ...). T is magtype of the element type of the Tape it was computed from. It carries the fields of JacobianCover together with

• b_quad: the quadrature merge b_j = √(∑_k (Ũ_kj/a_k)²).
• b_max: the smallest valid b for this a, b_j = max_k Ũ_kj/a_k.
• Utilde: m × n absolute-sensitivity surrogate satisfying Utilde >= abs.(J).

The merges satisfy b_max <= b_quad <= b_sum coordinatewise. Forming Utilde takes an m-channel reverse sweep, O(m * L) time and O(m * n) memory.

hvp(t::Tape, outidx::Integer, v::AbstractVector) -> Hv
hvp(f, x0::AbstractVector, v::AbstractVector) -> Hv
hvp!(Hv, t::Tape, outidx::Integer, v::AbstractVector) -> Hv

Hessian-vector product H*v of the taped program, formed without the Hessian by differentiating the adjoint recursion along v: one forward tangent sweep dv_k = sum_i (∂φ_k/∂v_i) dv_i seeded with v, then one reverse sweep carrying both the adjoints vbar and their tangents,

vbardot_i += (∂φ_k/∂v_i) * vbardot_k + vbar_k * sum_j (∂²φ_k/∂v_i∂v_j) * dv_j.

The curvature term of each node comes from curvature_apply!, which every Curvature supplies in the form its structure allows — a PairedCurvature node reads its 2q nonzeros off the pairs rather than scanning a 2q × 2q matrix, without which a program built from q-long inner products would cost O(q) gradients instead of two. A node whose second partials are only declared (has_exact_partials is false) makes an exact product impossible: hvp throws on any tape containing one, unconditionally, rather than letting success depend on whether the node's adjoint happens to vanish. The tape method of scaled_spectral_radius inherits this restriction; its dense method does not.

Costs about two gradients and O(L) memory.

hvp! writes into Hv, which shares the axes of v. It saves the result allocation only: the two sweeps still allocate their own O(L) working vectors, so it is not allocation-free.

Examples

julia> f(x) = exp(x[1]) + x[1] * x[2]^2;

julia> x0, v = [0.5, 2.0], [1.0, 0.0];

julia> round.(hvp(f, x0, v); digits=4)
2-element Vector{Float64}:
 1.6487
 4.0

julia> h = cover_hessian(f, x0).h_sum;   # the cover bounds the product

julia> all(abs.(hvp(f, x0, v)) .<= h .* (h' * abs.(v)))
true

hvp(t::Tape, outidx::Integer, v::AbstractVector) -> Hv
hvp(f, x0::AbstractVector, v::AbstractVector) -> Hv
hvp!(Hv, t::Tape, outidx::Integer, v::AbstractVector) -> Hv

Hessian-vector product H*v of the taped program, formed without the Hessian by differentiating the adjoint recursion along v: one forward tangent sweep dv_k = sum_i (∂φ_k/∂v_i) dv_i seeded with v, then one reverse sweep carrying both the adjoints vbar and their tangents,

vbardot_i += (∂φ_k/∂v_i) * vbardot_k + vbar_k * sum_j (∂²φ_k/∂v_i∂v_j) * dv_j.

The curvature term of each node comes from curvature_apply!, which every Curvature supplies in the form its structure allows — a PairedCurvature node reads its 2q nonzeros off the pairs rather than scanning a 2q × 2q matrix, without which a program built from q-long inner products would cost O(q) gradients instead of two. A node whose second partials are only declared (has_exact_partials is false) makes an exact product impossible: hvp throws on any tape containing one, unconditionally, rather than letting success depend on whether the node's adjoint happens to vanish. The tape method of scaled_spectral_radius inherits this restriction; its dense method does not.

Costs about two gradients and O(L) memory.

hvp! writes into Hv, which shares the axes of v. It saves the result allocation only: the two sweeps still allocate their own O(L) working vectors, so it is not allocation-free.

Examples

julia> f(x) = exp(x[1]) + x[1] * x[2]^2;

julia> x0, v = [0.5, 2.0], [1.0, 0.0];

julia> round.(hvp(f, x0, v); digits=4)
2-element Vector{Float64}:
 1.6487
 4.0

julia> h = cover_hessian(f, x0).h_sum;   # the cover bounds the product

julia> all(abs.(hvp(f, x0, v)) .<= h .* (h' * abs.(v)))
true

jvp(t::Tape, outidxs::AbstractVector{<:Integer}, v::AbstractVector) -> Jv
jvp(r, x0::AbstractVector, v::AbstractVector) -> Jv
jvp!(Jv, t::Tape, outidxs::AbstractVector{<:Integer}, v::AbstractVector) -> Jv

Jacobian-vector product J*v of the taped vector-valued program, formed without the Jacobian by a single forward tangent sweep dv_k = sum_i (∂φ_k/∂v_i) dv_i seeded with v at the input nodes and read off at outidxs. Costs one sweep and O(L) time and memory, independent of the output count m — unlike adjoints, which is O(m*L).

A repeated entry in outidxs yields an identical entry of Jv.

jvp! writes into Jv, which shares the axes of outidxs. It saves the result allocation only: the sweep still allocates its own O(L) working vector, so it is not allocation-free.

Examples

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

julia> round.(jvp(r, [0.5, 2.0, -1.0], [1.0, 0.0, 1.0]); digits=4)
3-element Vector{Float64}:
 2.0
 0.6487
 2.0

jvp(t::Tape, outidxs::AbstractVector{<:Integer}, v::AbstractVector) -> Jv
jvp(r, x0::AbstractVector, v::AbstractVector) -> Jv
jvp!(Jv, t::Tape, outidxs::AbstractVector{<:Integer}, v::AbstractVector) -> Jv

Jacobian-vector product J*v of the taped vector-valued program, formed without the Jacobian by a single forward tangent sweep dv_k = sum_i (∂φ_k/∂v_i) dv_i seeded with v at the input nodes and read off at outidxs. Costs one sweep and O(L) time and memory, independent of the output count m — unlike adjoints, which is O(m*L).

A repeated entry in outidxs yields an identical entry of Jv.

jvp! writes into Jv, which shares the axes of outidxs. It saves the result allocation only: the sweep still allocates its own O(L) working vector, so it is not allocation-free.

Examples

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

julia> round.(jvp(r, [0.5, 2.0, -1.0], [1.0, 0.0, 1.0]); digits=4)
3-element Vector{Float64}:
 2.0
 0.6487
 2.0

vjp(t::Tape, outidxs::AbstractVector{<:Integer}, u::AbstractVector) -> Jtu
vjp(r, x0::AbstractVector, u::AbstractVector) -> Jtu
vjp!(Jtu, t::Tape, outidxs::AbstractVector{<:Integer}, u::AbstractVector) -> Jtu

Vector-Jacobian product J'*u of the taped vector-valued program: one reverse adjoint sweep seeded with u[k] at node outidxs[k], read off at the input nodes. Costs one sweep and O(L) time and memory, independent of the output count m.

Outputs that name the same node accumulate their seeds, so a repeated entry in outidxs contributes the sum of the corresponding entries of u — the transpose of jvp's duplicated rows.

The seed u shares the axes of outidxs, the mirror of jvp's Jv.

vjp! writes into Jtu, which is indexed by the tape's inputs, Base.OneTo(t.ninputs). It saves the result allocation only: the sweep still allocates its own O(L) working vector, so it is not allocation-free.

Examples

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

julia> round.(vjp(r, [0.5, 2.0, -1.0], [1.0, 1.0, 1.0]); digits=4)
3-element Vector{Float64}:
 3.6487
 1.5
 1.0

vjp(t::Tape, outidxs::AbstractVector{<:Integer}, u::AbstractVector) -> Jtu
vjp(r, x0::AbstractVector, u::AbstractVector) -> Jtu
vjp!(Jtu, t::Tape, outidxs::AbstractVector{<:Integer}, u::AbstractVector) -> Jtu

Vector-Jacobian product J'*u of the taped vector-valued program: one reverse adjoint sweep seeded with u[k] at node outidxs[k], read off at the input nodes. Costs one sweep and O(L) time and memory, independent of the output count m.

Outputs that name the same node accumulate their seeds, so a repeated entry in outidxs contributes the sum of the corresponding entries of u — the transpose of jvp's duplicated rows.

The seed u shares the axes of outidxs, the mirror of jvp's Jv.

vjp! writes into Jtu, which is indexed by the tape's inputs, Base.OneTo(t.ninputs). It saves the result allocation only: the sweep still allocates its own O(L) working vector, so it is not allocation-free.

Examples

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

julia> round.(vjp(r, [0.5, 2.0, -1.0], [1.0, 1.0, 1.0]); digits=4)
3-element Vector{Float64}:
 3.6487
 1.5
 1.0

adjoints(t::Tape, outidx::Integer) -> vbar

Reverse-mode adjoints vbar[k] = ∂v_outidx/∂v_k for every node of t. The gradient of the program with respect to its inputs is vbar[1:t.ninputs], and cover_hessian(CoverOnly(), t, vbar) computes the Hessian cover from these adjoints without propagating them again.

adjoints(t::Tape, outidxs::AbstractVector{<:Integer}) -> Vbar

Vector-mode reverse-mode adjoints for several outputs at once: Vbar[k, j] = ∂v_{outidxs[k]}/∂v_j, every output k against every node j, all propagated in one traversal. The m × n Jacobian of the outputs with respect to the inputs is Vbar[:, 1:t.ninputs].

Each output owns a row seeded independently, so a repeated entry in outidxs yields identical rows and an input-node entry a unit row. Costs O(m*L) time and memory: a diagnostic path, like sensitivities, for the exact signed Jacobian; the production cover merges seed abs_reverse at the outputs instead and stay O(L)/O(G*L).

scaled_spectral_radius(H::AbstractMatrix, h::AbstractVector) -> lambda
scaled_spectral_radius(t::Tape, outidx::Integer, h::AbstractVector; maxiter=200, rtol=cbrt(eps(eltype(h))), rng=Xoshiro(0)) -> lambda

Spectral norm of the scaled Hessian Hhat = D^-1 H D^-1, D = diag(h), which measures how loose the cover h is in the spectral sense: h is valid iff every entry of Hhat is at most 1 in magnitude, and lambda >= max_ij abs(Hhat[i,j]).

Because every entry of the symmetric Hhat is bounded by its spectral norm, abs(H[i,j]) <= lambda * h[i] * h[j], so sqrt(lambda) * h is also a cover of H. It is tighter than h when lambda < 1. lambda must bound the scaled norm wherever the cover is meant to hold: for a cover computed over a set of points, such as a box of intervals, that is the supremum over the set.

The dense method is exact. Coordinates with h[i] == 0 are admissible only if row i of H vanishes, in which case they are dropped; otherwise h does not cover H and an error is thrown.

The tape method is matrix-free: power iteration on Hhat driven by hvp, stopped when the estimate changes by less than rtol relative. It requires h to be strictly positive, and a tape whose element type is its own magtype: power iteration evaluates at a point, so a set-valued tape (such as one recorded on intervals) is an error. Each iterate norm(Hhat*v) with norm(v) == 1 is a lower bound on lambda and the sequence converges to it from below, so rescaling a cover by sqrt(safety * lambda) is valid only for a safety > 1 large enough to make safety * lambda an upper bound. Failure to converge within maxiter is an error, not a silently returned partial result.

scaled_singular_value(J::AbstractMatrix, a::AbstractVector, b::AbstractVector) -> sigma
scaled_singular_value(t::Tape, outidxs, a::AbstractVector, b::AbstractVector; maxiter=200, rtol=cbrt(eps(promote_type(eltype(a), eltype(b)))), rng=Xoshiro(0)) -> sigma

Largest singular value of the equilibrated Jacobian Jhat = diag(a)^-1 J diag(b)^-1 (see equilibrate), which measures how loose the cover (a, b) is in the spectral sense: (a, b) is valid iff every entry of Jhat is at most 1 in magnitude, and sigma >= max_kj abs(Jhat[k,j]).

Because every entry of Jhat is bounded by its spectral norm, abs(J[k,j]) <= sigma * a[k] * b[j], so (sqrt(sigma) * a, sqrt(sigma) * b) is also a cover of J, tighter than (a, b) when sigma < 1. sigma must bound the scaled norm wherever the cover is meant to hold: for a cover computed over a set of points, such as a box of intervals, that is the supremum over the set.

The dense method is exact, and drops zero weights on the same terms as equilibrate.

The tape method is matrix-free: power iteration on Jhat'Jhat driven by jvp and vjp — about two gradients per iteration, O(L) memory — stopped when the estimate changes by less than rtol relative. It requires a and b to be strictly positive, and a tape whose element type is its own magtype: power iteration evaluates at a point, so a set-valued tape (such as one recorded on intervals) is an error. Each iterate norm(Jhat*v) with norm(v) == 1 is a lower bound on sigma and the sequence converges to it from below, so rescaling a cover by sqrt(safety * sigma) is valid only for a safety > 1 large enough to make safety * sigma an upper bound. Failure to converge within maxiter is an error, not a silently returned partial result.

equilibrate(J::AbstractMatrix, a::AbstractVector, b::AbstractVector) -> Jhat

The two-sided equilibrated Jacobian Jhat = diag(a)^-1 J diag(b)^-1, in which the cover pair (a, b) becomes the all-ones covering matrix: (a, b) covers J if and only if every entry of Jhat is at most 1 in magnitude. The equilibrated matrix is the scale-invariant object a downstream linear method should see, since it is unchanged by any rescaling of the outputs or inputs that the cover tracks.

A zero weight is admissible only if the row or column it scales vanishes, in which case that row or column of Jhat is zero; otherwise (a, b) does not cover J and an error is thrown.