DerivativeCoverProblems

DerivativeCoverProblems

Parameterized families of scalar programs for studying cover quality.

The programs are written in generic Julia and depend on nothing here, so any operator-overloading tool can run them; DerivativeCovers tapes them by virtue of DerivativeCovers.TVar overloading the elementals they use.

Each generator returns a CoverProblem: a program f written from those elementals, the point x0 to evaluate it at, an optional idiomatic-Julia twin that computes the same function by an independent route, and an optional analytic hessian. The families differ in the property that drives cover looseness — how much the paths from an input to the output reconverge — and each has knobs for size, reconvergence, and dynamic range:

• quadratic: x'Ax/2 with A = W*W' of prescribed rank, written entrywise, as a matrix-vector product, or in factored form. The same A in three programs, so the three differ only in graph structure.
• deep_composition: layers of width channels around a shared trunk, so each layer fans a single node out to width consumers that immediately reconverge.
• random_dag: a random program whose sharing parameter interpolates from a tree (sharing = 0) to heavy subexpression reuse (sharing = 1).
• adjoint_cancellation: a chain of multiplication nodes the output is nearly stationary in, so the paths from each node to the output cancel in the adjoints rather than the intermediate gradients. converged_inner_solve reaches the same condition by the envelope theorem.

Wide dynamic range enters along two independent axes: sigma widens the log-magnitude spread of a family's coefficients, while rescale (with scales from logscales) spreads the units of the inputs.

CoverProblem

A scalar program together with the references needed to judge a cover of its Hessian.

• name: a short identifier, unique within a family sweep.
• f: the program. Built only from the elementals overloaded for DerivativeCovers.TVar, so f(x) records a tape when x is a vector of TVars, and evaluates normally on any other numeric vector (Float64, ForwardDiff.Dual, ...).
• x0: the point at which the Hessian is covered.
• twin: an idiomatic-Julia implementation of the same function, or nothing. It is written independently of f, by an independent route — array operations f cannot express (most families), or an explicit scalar loop when f is itself written in array notation (see mlp) — so agreement between the two is evidence that f is the function it claims to be.
• hessian: x -> H, the analytic Hessian, or nothing when no closed form is available. Use exact_hessian rather than the field.
• xstar: a global minimizer of f, or nothing when none is known. Problems built with a planted solution carry it.

A random program, stored as the recipe that generated it: ops[k] builds node ninputs + k from earlier nodes, and the output is the sum of outputs. Replaying the recipe on a vector of any numeric type evaluates the same function, so one recipe serves the tape, the reference derivatives, and the plain evaluation.

JacobianProblem

A vector-valued program together with the references needed to judge an asymmetric cover of its Jacobian — the vector-output counterpart of CoverProblem.

• name: a short identifier, unique within a family sweep.
• r: the program, mapping an input vector to an output vector. Built only from the elementals overloaded for DerivativeCovers.TVar, so r(x) records a tape when x is a vector of TVars and evaluates normally on any other numeric vector.
• x0: the point at which the Jacobian is covered.
• nout: the number of outputs of r.
• twin: an independent idiomatic-Julia implementation of the same map, or nothing.
• jacobian: x -> J, the analytic nout × n Jacobian, or nothing when no closed form is available. Use exact_jacobian rather than the field.

Coefficients with magnitude abs(randn()) * exp(sigma * randn()): unit scale at sigma = 0, and log-magnitudes spread by sigma beyond that.

The linear combination sum_i c[i] * z[i], recorded one term at a time.

The left-to-right chain w[1]*z[1] + w[2]*z[2] + ... + w[n]*z[n], the scalar order a dot node's parents follow.

Sum the entries of v using only the binary + of its element type. v must be non-empty.

Take a parent for a new node. With probability sharing (or when nothing is unconsumed) reuse an arbitrary existing node, giving it a second consumer and so creating reconvergence; otherwise consume a node from frontier, which leaves it with exactly one consumer.

adjoint_cancellation(n; delta=1e-2, sigma=0.0, rng=MersenneTwister(1)) -> CoverProblem

A chain of n - 1 multiplication nodes u_j = x_j * x_{j+1}, each reaching the output through u_j - sin(u_j), evaluated where every u_j is small.

u - sin(u) has derivative 1 - cos(u), which vanishes to second order at u = 0, so the output is nearly stationary in each u_j: the signed adjoint of that node is ~u_j^2/2 while its two paths to the output contribute ~2 in absolute-value arithmetic. The adjoint index acheck/abs(vbar) at the node is therefore ~4/u_j^2 and the cover separation ~2/u_j, both set by delta alone.

x0 alternates delta, 1, delta, 1, ..., so every u_j is delta and every unit carries the same index. The units overlap — x_j feeds both u_{j-1} and u_j — so the Hessian is tridiagonal rather than block diagonal, and the separation is a property of each node rather than of the problem size: it is the same at every n, which is what makes the family usable in a size sweep.

Each input reaches the output along one path per unit, so no intermediate gradient cancels: its absolute-value bound is exact, and the family isolates the cancellation in the adjoints. sigma spreads the log-magnitudes of the odd coordinates, spreading the per-unit index without touching the mechanism.

cancellation(n; m=n, delta=1e-3, aligned=false, support=1:n, sigma=0.0, rng=MersenneTwister(1)) -> CoverProblem

An adversarial pair of programs built from m quartets of squared linear forms whose signed sum cancels exactly, riding on a delta-weighted rank-1 quadratic:

acc = (delta/2) * (q'x)^2
for j in 1:m
    acc += 0.5*(a_j'xS)^2 + 0.5*(b_j'xS)^2 + s*(c_j'xS)^2 + s*(d_j'xS)^2
end

where xS = x[support], c_j = (a_j+b_j)/sqrt(2), d_j = (a_j-b_j)/sqrt(2), and a_j, b_j are drawn orthogonal to x0[support]. At s = -0.5 (aligned = false, the default) the rotation identity c_j c_j' + d_j d_j' = a_j a_j' + b_j b_j' makes every quartet's contribution to the Hessian cancel exactly; at s = 0.5 (aligned = true) the same identity makes it add. The two variants differ only in the sign of the two constants multiplying the c_j and d_j squares: every node of the tape up to that sign carries the same value in both, because every quartet linear form vanishes at x0 (orthogonality to x0[support]), so f, its gradient, and the magnitude of every elemental partial and adjoint coincide between the two variants — only the curvature they carry differs.

Closed forms (the Hessian is constant in x):

H_cancel  = delta * q*q'
H_aligned = delta * q*q' + sum_j 2*(a_j*a_j' + b_j*b_j')   (scattered onto `support`)

The minimal cover of H_cancel is hstar = sqrt(delta) .* abs.(q) = sqrt.(diag(H_cancel)) (a rank-1 |H|, since H_cancel is itself rank 1). The quadrature merge of the cancel variant attains sqrt.(diag(H_aligned)) exactly — the sign flip that cancels the Hessian does not touch the elementwise bound B, which sees only magnitudes — so its looseness against hstar is sqrt.(kappa) with kappa_i = diag(H_aligned)_i / diag(H_cancel)_i, the per-coordinate cancellation index.

m sets the number of quartets and hence the size of kappa on the support (more quartets, more cancellation, larger kappa); delta sets the scale of the surviving rank-1 term relative to the quartets; support restricts the quartets to a coordinate subset, leaving off-support coordinates untouched by any cancellation; sigma widens the log-magnitude spread of the quartet coefficients and of q. Two calls that agree in n, m, delta, support, sigma, and rng state draw the same x0, q, and quartets regardless of aligned, so a cancel/aligned pair built from fresh, identically seeded rngs shares every random draw and differs only in the program's two constants.

converged_inner_solve(; steps=2) -> CoverProblem

A bilevel program: steps Newton iterations on the inner problem min_u (u - x1)^2/2 + x2 u^4/4, with the output that same inner objective at the final iterate.

By the envelope theorem the output is stationary in u at the exact inner minimizer, so the adjoint through the final iterate collapses as the solve converges while the unsigned sum over paths does not. Newton doubles the accuracy each step, so the adjoint index grows superexponentially in steps. At steps = 0 the inner variable keeps its initial value and nothing cancels, which is the control.

Unlike adjoint_cancellation this engages the intermediate-gradient channel too — the realistic case, since the same structure appears in implicit layers, self-consistent field iterations, and any objective differentiated through a converged inner solve. Its knob is steps rather than a size, so it takes no n and is not one of families.

Past four steps the inner residual is at roundoff and a further step measures the floor rather than the mechanism.

deep_composition(n; depth=4, width=4, sigma=0.0, planted=nothing, rng=MersenneTwister(1)) -> CoverProblem

A depth-depth composition of layers, each of which sums its width channels into a single trunk node and then feeds that node back into all width channels:

z = tanh.(W*x + b)
repeat depth times:
    t = sum(z)                       # one node, `width` consumers
    z = tanh.(a .* t .+ c .* z)
f = sum(z.^2) / 2

The trunk is the reconvergence: every channel of layer l+1 depends on every channel of layer l, so the paths from an input to the output multiply by width per layer. width and depth therefore tune reconvergence and program length independently of n.

sigma widens the input weights W and biases b only. The recurrent weights a and c keep unit scale (1/sqrt(width) per fan-in), because a wide spread there drives tanh into saturation and flattens the Hessian to nothing. Heterogeneous units belong on the inputs, which is where sigma puts them.

A real planted replaces the biases by b = -W*xstar with xstar = x0 .+ planted .* randn(n), so that the first layer, and with it every later layer, vanishes at xstar: f(xstar) = 0, the minimum of f. For width >= n and W of full rank, xstar is the only zero of the first layer. planted sets the distance of x0 from xstar in units of the inputs; W, a, c, and x0 are those of the unplanted problem drawn from the same rng state. The problem carries xstar as p.xstar.

No closed-form Hessian; differentiate p.f.

exact_hessian(p::CoverProblem, x=p.x0) -> H

The analytic Hessian of p at x. Errors when p has none; test with has_exact_hessian first.

exact_jacobian(p::JacobianProblem, x=p.x0) -> J

The analytic Jacobian of p at x. Errors when p has none; test with has_exact_jacobian first.

families(; sigma=0.0, seed=1) -> Vector{Pair{String,Any}}

The standard sweep: named generators n -> CoverProblem, one per family, differing in how severely the paths through the program reconverge. Pass the same n to each to compare families at fixed problem size, or sweep n within one to measure how a cover's looseness grows.

The three quadratic forms are seeded alike, so at a given n they carry the same Hessian and differ only as programs.

Every family but the last cancels in the intermediate gradients; adjcancel cancels in the adjoints, and is the only one whose looseness is set by a knob rather than by n.

has_exact_hessian(p::CoverProblem)

Whether p carries an analytic Hessian. When it does not, a reference Hessian must come from elsewhere (ForwardDiff.hessian(p.f, p.x0)).

has_exact_jacobian(p::JacobianProblem)

Whether p carries an analytic Jacobian. When it does not, a reference Jacobian must come from elsewhere (ForwardDiff.jacobian(p.r, p.x0)).

has_twin(p::CoverProblem)

Whether p carries an independent reference implementation of its program.

has_xstar(p::CoverProblem)

Whether p carries a known global minimizer p.xstar.

jacobian_cancellation(n; nout=4, m=n, delta=1e-3, aligned=false, support=1:n, sigma=0.0,
                      rng=MersenneTwister(1)) -> JacobianProblem

An adversarial pair of vector-valued programs built from m quartets of linear forms whose signed sum cancels exactly, riding on a delta-weighted rank-1 Jacobian and fanned out to nout outputs:

T = sum_i [ 0.5*(a_i'xS) + 0.5*(b_i'xS) + s*(c_i'xS) + s*(d_i'xS) ]
r_k = delta * q[k] * (p'x) + q[k] * T

where xS = x[support] and c_i .+ d_i == a_i .+ b_i entrywise (an unequal split c_i = t_i .* (a_i .+ b_i), t_i entrywise in (0, 1), so all four forms are distinct and none is a scalar multiple of another). At s = -0.5 (aligned = false, the default) each quartet's contribution to every Jacobian entry cancels exactly; at s = 0.5 (aligned = true) it adds. The two variants differ only in the sign of the constant multiplying the c_i and d_i forms: every quartet form vanishes at x0 (orthogonality, arranged by the draw) and within each support column all four forms share one sign, so node values and the magnitude of every elemental partial coincide between the two variants — the surrogate Ũ and every merge built from it are identical, and only the signed Jacobian differs.

Closed forms (the Jacobian is constant in x; ghat = sum_i (a_i .+ b_i) scattered onto support, g = abs.(ghat)):

J_cancel  = delta .* q * p'
J_aligned = q * (delta .* p .+ ghat)'
Utilde    = abs.(q) * (delta .* abs.(p) .+ g)'     (both variants)

The sign alignment makes Utilde == abs.(J_aligned) exactly, and both Utilde and the two Jacobians are rank 1, so the minimal cover of each is exact and analytic. The per-entry cancellation index is constant across outputs: kappa_j = 1 + g_j / (delta * abs(p_j)) on the support, 1 off it — more quartets or smaller delta, larger kappa; support restricts the cancellation to a column subset, making the kappa profile anisotropic.

sigma widens the log-magnitude spread of the quartet coefficients and of the output scales q, whose spread sets the sum/max merge gap sum(abs, q)/maximum(abs, q) — the effective output multiplicity, nout at sigma = 0 and approaching 1 as one output dominates. Two calls that agree in everything but aligned share every random draw and differ only in the program's one constant.

logscales(rng, n; sigma) -> d

Positive scale factors with log.(d) normal of standard deviation sigma. Pass to rescale to give a problem's inputs heterogeneous units.

lotka_volterra(m; delta=0.3, decades=2.0, ntraj=3, h=0.1, nsteps=50, every=5,
               rng=MersenneTwister(1)) -> CoverProblem

Least-squares estimation of the parameters of a generalized Lotka–Volterra system of m species,

dy_i/dt = y_i * (r_i + sum_j A[i,j] * y_j),

from its own simulated trajectories. The inputs are the parameters x = vcat(r, vec(A)), so n = m + m^2. The program integrates ntraj trajectories by nsteps classical Runge–Kutta steps of size h, observes every every-th step, and sums the squared misfits to the data, each divided by the equilibrium abundance of its species:

f(x) = sum over trajectories, observations o, species i of ((y_i(o; x) - data_i(o)) / ybar_i)^2 / 2

The data come from planted parameters xstar, built around a positive equilibrium ybar whose entries are spread log-uniformly over decades decades: A[i,j] = B[i,j] / ybar[j] with B the negative identity plus a random perturbation whose symmetric part is kept negative definite, and r = -A*ybar, so ybar is a stable equilibrium. The parameters of the abundant and the rare species therefore differ in units by decades decades. Initial conditions are ybar times log-normal factors of standard deviation 1/2. The data are simulated by f's own arithmetic, so f(xstar) = 0 exactly, the minimum of f; the problem carries xstar as p.xstar.

The start is x0 = xstar .+ delta .* abs.(xstar) .* randn(n): delta is the relative perturbation of each parameter. Far from xstar the dynamics can blow up within the horizon, and f then returns Inf or NaN.

Every step of the integrator reuses every species' state in every right-hand side, so paths from a parameter to the output reconverge at each of the 4 * nsteps stages.

No closed-form Hessian; differentiate p.f.

mlp(sizes; batch=4, sigma=0.0, activation=tanh, rng=MersenneTwister(1)) -> CoverProblem

The squared-error loss of a feedforward network, with the parameters as inputs. For sizes = (n0, n1, ..., nK), layer l = 1:K has a weight matrix W_l of size n_l x n_{l-1} and a bias vector b_l of length n_l, flattened (W1, b1, W2, b2, ..., WK, bK) into one vector of length p = sum_l n_l*(n_{l-1}+1):

z = data column                    # n0, fixed
for l = 1:K-1
    z = activation.(W_l*z .+ b_l)  # hidden layers
y = W_K*z .+ b_K                   # output layer, linear
f = sum_over_columns( sum(abs2, y - target) ) / 2

batch input/target column pairs are drawn once at construction and closed over; only the parameters vary thereafter. f is written in array notation: a matvec whose TVar factor meets a constant factor (the data column feeding layer 1, or feeding layer K directly when K = 1) records one linear tape node per row, while a matvec between two TVar factors (every hidden-to-hidden transition) records one inner-product node per row. Each column's error sums as one linear node, and the batch total sums as another, so the tape stays short despite the array notation.

twin computes the same loss by an explicit scalar loop: every row's dot product is the left-to-right chain w[1]*z[1] + w[2]*z[2] + ... in the order a dot node's parents follow, bias and activation are applied per scalar, and every sum (over output units, over the batch) is an explicit binary + chain. It therefore matches f node for node under a shared accumulation order, rather than taking an independent route through array operations as the other families' twins do.

sigma spreads the log-magnitude of the initial weights and biases, mirroring deep_composition's convention; the weights are additionally scaled by 1/sqrt(n_{l-1}) so a wider input layer does not by itself inflate the initial activations.

No closed-form Hessian; differentiate p.f.

noutputs(p::JacobianProblem)

Number of outputs of p.r.

nvars(p::CoverProblem)

Number of inputs of p.

quadratic(n; rank=n, form=:entrywise, sigma=0.0, rng=MersenneTwister(1)) -> CoverProblem

The quadratic f(x) = x'Ax/2 with A = W*W' positive semidefinite of the given rank, where W is n × rank and sigma sets the spread of its coefficients' log-magnitudes. The Hessian is A, independent of x.

form selects how the same function is written as a program, which is what the cover sees:

• :entrywise — sum_{i,j} A[i,j] * x[i] * x[j] / 2. Every entry of A gets its own nonlinear node, and all n^2 of them reconverge on the inputs.
• :matvec — y = A*x scalarized, then x'y/2. The n products x[i]*y[i] reconverge through the rows of A.
• :factored — s = W'x, then sum(s.^2)/2. Only rank nonlinear nodes, each carrying a whole rank-1 layer of the Hessian, so no two nonlinear contributions overlap.

For rank = 1 the factored form's minimal cover is known exactly: H = w*w' forces h >= abs.(w), and abs.(w) = sqrt.(diag(A)) attains it.

Two calls that agree in n, rank, sigma, and rng state produce the same A and the same x0 in every form, so a comparison across forms isolates the effect of program structure.

random_dag(n; nops=4n, sharing=0.5, sigma=0.0, rng=MersenneTwister(1)) -> CoverProblem

A random program of n inputs and at most nops elementary operations, drawn from +, -, *, scaling, tanh, squaring, and sin. The output sums every node that no other node consumes, so every input reaches the output. At sharing = 0 generation also stops as soon as a single unconsumed node remains (once any curvature-carrying op has been drawn): the only way to grow a tree past that point is one long unary chain, whose accumulated derivative factors wander into under- or overflow. A draw whose ops all happen to be linear gets one final squaring node, so every generated program has a nonzero Hessian.

sharing controls subexpression reuse, and with it the reconvergence the cover has to contend with. Each operand is either an unconsumed node (probability 1 - sharing) or an arbitrary existing node (probability sharing). At sharing = 0 no node is ever consumed twice, so the program's graph is a tree: each input reaches each node along at most one path. At sharing = 1 operands are drawn from the whole program and paths reconverge freely. Reaching the output is graph reachability, not influence: with reuse, an input can enter only linearly (a zero Hessian row), or even through linear paths whose unit weights cancel exactly (a zero gradient entry as well). Both are rare, correct, and must be tolerated by downstream code.

A node whose magnitude at x0 leaves [1e-3, 1e3] is followed by a normalizing scale node, so * and ^2 chains can neither overflow nor collapse to a numerically zero derivative. A tanh input outside [-2, 2] is likewise normalized to magnitude 1 first: 1 - tanh(v)^2 underflows to exactly zero for large v, and a derivative-dead node erases its whole subtree from the gradient and the Hessian. These guards are what keep the program nonlinear at every length, and what make the recipe safe to replay at points other than x0.

No closed-form Hessian; differentiate p.f.

rescale(p::CoverProblem, d::AbstractVector) -> CoverProblem

The same function in rescaled input coordinates: g(y) = f(y ./ d) at y0 = d .* x0, whose Hessian is Diagonal(inv.(d)) * H * Diagonal(inv.(d)).

This is the input-units axis of dynamic range, and it is exactly the transformation a scale-invariant cover must be covariant under: a cover h of p at x0 must reappear as h ./ d for the rescaled problem. Take d from logscales for a random spread.