Divide instead of multiplying by cached reciprocals (fixes #239) - #298
Merged
Merged
Conversation
The reverse rules for `mean` and elementwise `/` and `\` multiplied the adjoint by a precomputed reciprocal. That reciprocal was formed in the type of the denominator, so an integer or `Float64` denominator capped the gradient at `Float64` accuracy regardless of the deriv type, and `inv(d)` overflowed to `Inf` for subnormal `d`. Both rules now divide. `mean` divides by `length(input)` directly, and the `/` and `\` instructions no longer cache the numerator partial at all: `broadcast_increment_div_deriv!` divides by `value(b)` during propagation, mirroring how `broadcast_mul` already reads `value(b)` in its reverse pass. This drops one cache array per instruction. Scalar accumulations in `propagation.jl` now go through `sum`, whose pairwise reduction is both more accurate and several times faster than the sequential loops they replace. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## master #298 +/- ##
==========================================
+ Coverage 86.08% 86.69% +0.60%
==========================================
Files 19 19
Lines 1947 1916 -31
==========================================
- Hits 1676 1661 -15
+ Misses 271 255 -16 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
`denom_partials_kernel` computed `-n / (d * d)`, but `d * d` overflows to
`Inf` (and underflows to zero) while `-n / d^2` is still finite and normal,
so the denominator adjoint was wrong at the extremes of the range:
ReverseDiff.gradient(v -> sum([1e-200] ./ v), [1e-200]) # [-Inf]
ReverseDiff.gradient(v -> sum([1e200] ./ v), [1e200]) # [0.0]
Dividing twice never forms the square and is the same two roundings, so
accuracy in the normal range is unchanged: over 200k random `(n, d)` pairs
the mean relative error against a `BigFloat` reference is 4.007e-17 before
and 4.017e-17 after. The kernel backs both `./` and `.\`, so one change
covers the numerator adjoint's counterpart in each.
`PrecisionTests` also moves after `LinAlgTests`, which asserts that
`StatisticsExt` is still unloaded when it starts; `PrecisionTests` loads
Statistics and so has to run later.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
`numer_partials(d::Real) = Ref(inv(d))` was the last producer of a `Ref` partial, and dividing by `value(b)` directly retired it. Nothing constructs a `Ref` any more: the remaining `Ref` methods of `denom_partials`, `base_partials` and `exp_partials` need both operand values to be `<:Real`, and none of the nine operand combinations generated at `elementwise.jl:419` (or the matching `materialize` overloads in `broadcast.jl`) is scalar/scalar — at least one side is always array-valued. So the `partial::Ref` forwarder is unreachable, as is `_broadcast_increment_deriv!` for a `TrackedReal` input with a scalar partial: when the input is a tracked scalar the other operand is an array, so its partial is an `AbstractArray`. Codecov flagged exactly these two methods on #298. Dropping them is fail-loud rather than silently wrong: no less specific method survives for either signature, so a mistaken deletion is a `MethodError`, never a wrong gradient. `ElementwiseTests` covers all nine combinations. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
devmotion
added a commit
that referenced
this pull request
Sep 22, 2026
Almost all of ReverseDiff's broadcasting hung off entry points Julia never
documented as extension points, and the documented ones were unused: 2191
`Base.broadcast` overloads, 1777 `Base.broadcast!`, and a 2001-method
`materialize` table keyed on *any* package's style, against five `copy`
methods and no `copyto!`. Because `broadcast` was overloaded directly it
short-circuited the pipeline before a `Broadcasted` was ever built, so `f.(x)`
and `broadcast(f, x)` could diverge, which `Base` guarantees they cannot.
`TrackedStyle` becomes `TrackedStyle{N} <: AbstractArrayStyle{N}`, so it tracks
dimensionality and picks up `Base`'s scalar, tuple and `Unknown` precedence
rules for free. The `Base.broadcast` overloads and the `materialize` table are
deleted, making `f.(x)` and `broadcast(f, x)` the same computation by
construction. `copyto!` now throws a descriptive error instead of failing after
polluting the tape; `.=` with a tracked right-hand side has never worked.
`copy` no longer asks `Core.Compiler.return_type` what a broadcast will produce
in order to pick an implementation. The only question the argument types cannot
answer is whether a tracked value reaches `f` other than through the arguments,
such as a closure capturing one, and `mayhavetracked` already reports that.
What the output turns out to be is settled after the forward pass instead:
anything `f` hands back that is not a `Dual` does not depend on the tracked
arguments, so it carries no partials and is left off the tape. That covers a
`Bool` from a predicate, an `Int` from `exponent`, a non-`Real`, and a function
constant in its argument, with one test. The `:tracker` implementation goes
with it, along with three `@generated` functions.
`∇broadcast` seeds `Dual`s only for the tracked arguments. `DiffRules` gives
`besselj` a literal `NaN` partial for its order, and `NaN * 0.0` is `NaN`, so
dualizing an untracked argument poisoned every partial of the same `Dual`.
Narrowing the seeding also narrows the `Dual` to the number of tracked
arguments and removes the `SVector`/`DiffResults` machinery. The tag is now
derived once for the whole broadcast from the arguments that are actually
seeded, rather than per element from the element types.
Where `f`'s partials are known in closed form, `knownpartials` dispatches on the
operation and its arguments and the cache keeps a `Contract(op, args)` per
tracked argument instead of a `Dual` per element. The entry says that argument
collects `op(seed, args...)` from every element of the output's derivative,
with `args` naming broadcast arguments by position:
`+`, `-`, `identity` `identity` or `-`, applied to the seed alone
`*` `*` with the other argument
`/`, `\` `/` by an untracked denominator
Nothing is stored. A `Contract` is a zero-size type, and an argument a partial
names is read back through `instruction.input`, where `record!` has already
captured it, so the partial cannot go stale on a replay and an untracked array
is not copied a second time. Every `knownpartials` method annotates all of
`f`'s arguments, which is what makes naming one by position safe: those are
exactly the arguments `splitargs` keeps, in order, so an `f` that also receives
arguments held back from the tape matches no method. `op` meets the seed rather
than the partial, so `/` forms no reciprocal, which is how `sum(y ./ 1e-310)`
keeps the finite gradient that #298 gave it.
A tracked denominator keeps a `Dual` per element, its partial `-x/y^2` being no
argument of the broadcast, as does an argument whose shape is not the output's,
which could not be read at the output's own index.
Arrays of `TrackedReal` are differentiated like a `TrackedArray` rather than
traced scalar-wise: 64x faster on `x .* y` and 166x on a fused three-argument
expression, with the tape going from ~3n entries to 2.
Against master, with `gradient` and Chairmarks, best of five alternating runs:
case master.us branch.us t.ratio bytes.ratio
mulscalar_1k 2.41 2.19 0.91 1.00
mulscalar_10k 15.75 14.33 0.91 1.00
mulscalar_100k 159.58 144.62 0.91 1.00
addscalar_10k 15.62 14.13 0.90 1.00
addarray_10k 17.33 15.75 0.91 1.00
mularray_10k 21.58 18.83 0.87 1.00
twoscalars_10k 24.67 23.00 0.93 0.87
sinscalar_10k 207.21 88.46 0.43 1.00
exp_10k 51.75 53.25 1.03 1.00
exp_100k 530.96 535.50 1.01 1.00
`exp` keeps a `Dual` per element either way and is at parity; the spread across
runs is about 3%. Precompilation drops from 10.45s to 6.78s and load time from
2.61s to 2.39s, and `src` loses 417 lines.
Fixes #265: a 0-d broadcast collapsed to a single `ImmutableDiffResult`, whose
`length` was then taken.
Fixes #214: a scalar and a one-element-vector argument now agree, where the
scalar form raised a `ForwardOptimize` ambiguity.
Addresses #226: the counted methods on `broadcast`, `broadcast!` and
`materialize` go from 5969 to 0, and `Test.detect_ambiguities(ReverseDiff)`
from 877 to 581, with the 296 broadcast-related ones all gone. `map`, `hcat`
and `vcat` are untouched and still account for 4567.
`seed!`'s catch-all is narrowed to `::Real`: it silently did nothing for the
arrays and tuples that reach it from a broadcast whose output element type is
not `Real`, and `jacobian` now covers those outputs.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Fixes #239.
The reverse rules for
meanand elementwise/and\multiplied the adjoint by a precomputed reciprocal formed in the type of the denominator, so an integer orFloat64denominator capped the gradient atFloat64accuracy regardless of the deriv type. Both now divide instead, and the/and\instructions no longer cache the numerator partial at all —broadcast_increment_div_deriv!divides byvalue(b)during propagation, mirroringbroadcast_mul.This also fixes two overflows.
inv(d)isInffor subnormald, so1e-10 * sum(y ./ 1e-310)had an infinite gradient where the exact1e300is representable. Anddenom_partials_kernelformed-n / (d * d), butd * doverflows toInf(and underflows to zero) while-n / d^2is still finite and normal:It now divides twice, which never forms the square and is the same two roundings, so accuracy in the normal range is unchanged: over 200k random
(n, d)pairs the mean relative error against aBigFloatreference is 4.007e-17 before and 4.017e-17 after. The kernel backs both./and.\.Separately, scalar accumulations in
propagation.jlnow go throughsum, whose pairwise reduction is both more accurate and faster than the sequential loops it replaces. On a 100k-element broadcast into a scalar adjoint this cuts the relative error from 6.9e-15 to 3.6e-17 — the failure mode the removed docstring predicted, and the divide-and-conquer fix it proposed.test/PrecisionTests.jlcompares against closed-form gradients at 256-bit precision. It runs afterLinAlgTests, which asserts thatStatisticsExtis still unloaded when it starts. The one@test_brokenis a fused broadcast, which needs an upstream fix inForwardDiff._div_partials.🤖 Generated with Claude Code