Uh oh!
There was an error while loading. Please reload this page.
Exploit symmetry in the Hessian, on top of #840 - #843
Open
devmotion wants to merge 7 commits into
Open
Conversation
instead of relying on the jacobian of gradient for the hessian explicitly seed dual numbers and only calculate the upper triangular part when chunking, gives ~2x speedup as input length becomes big
…#840 The symmetric sweep of #837 was written against `structural_eachindex`, which #840 replaced with `structural_indices`: linear indices of `x`, in seeding order, cached per work buffer on the config. Every conflict resolves to that API. Three consequences beyond the mechanical rename: * The result is `length(x)` by `length(x)`, both axes indexed by the linear indices of `x`, as #840 already established for the Jacobian and for the Hessian it computes from one. The rows and columns of the structurally zero entries are zeroed by `zero_unseeded_columns!`, and the gradient buffer by `zero_unseeded!`. * `extract_hessian_gradient_chunk!` can keep delegating to `extract_gradient_chunk!`, which now takes its positions from `x` rather than from the result -- that was #838, and it is what made `DiffResults.HessianResult`'s dense gradient buffer scatter its derivatives. * The entry points call `checkstructure` themselves. The sweep no longer goes through `gradient` and `jacobian`, so a config built for a differently structured input would otherwise no longer be caught (#842). The structured testset gains `n = 5`, the flat-vector result and the default chunk from #837's own, which it otherwise supersedes. `hessian!`'s dropped result shape check and the sweep's chunk-size dependence are left for the commits that follow. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
`HessianConfig` wrapped a `JacobianConfig` and a `GradientConfig` because the Hessian was `jacobian(gradient(f), x)`: the outer sweep seeded the Jacobian config's buffer and the inner `gradient` seeded the gradient config's. The symmetric sweep seeds both layers of the one nested buffer, so it reads four things -- the two seed tuples, the nested buffer and its positions -- and never touches the Jacobian config's buffer again after the constructor derives the nested element type from it. Holding those four directly drops the dead buffer, and `checkstructure(cfg, x)` now resolves through the generic `AbstractConfig` method like every other config's. At `length(x) == 1000` and a chunk size of 12: HessianConfig(f, x) 1.472 MB -> 1.368 MB HessianConfig(f, result, x) 1.576 MB -> 1.368 MB The result-aware constructor allocated two dead buffers rather than one, since it built the `f!(y, x)` `JacobianConfig`. Nothing about the work buffers depends on `result`, so it forwards to the plain constructor and the two now return the same type -- which is what the tests asserting the two configs interchangeable already implied. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The symmetric sweep writes the result entry by entry, which lost the validation that `reshape_jacobian` and the broadcast in `extract_jacobian!` used to provide. A vector result was still checked by `reshape`, a matrix one no longer was: hessian!(fill(NaN, 4, 4), f, rand(3)) # no error, row/col 4 left NaN hessian!(fill(NaN, 4, 4), f, SVector(1., 2., 3.)) `reshape_hessian` mirrors `reshape_jacobian`, down to its `DiffResult` method, so the `DiffResult` path is checked too -- it holds a buffer no entry point ever passes to `require_one_based_indexing`, hence the extra call here. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Diagonal blocks read `outer <= inner`, off-diagonal blocks read `outer > inner`, and the StaticArrays path reads `outer <= inner` throughout. Reading an entry with i in the outer layer rounds differently from reading it with j there, so the result was not reproducible: at `n = 16`, chunk sizes 1, 2, 3, 5, 7, 11 each differed from chunk 16 and from the `SVector` path by ~1e-16. Swapping which layer block q carries fixes it. Same number of evaluations and seed writes, and q is still seeded once outside the loop; afterwards every chunk size is bitwise identical to the `SVector` path. `log(sum(exp, z))` is the objective in the new test because its mixed partials actually round differently in the two orders -- `sum(z)^3`, `exp(sum(z))`, `prod(z)` and `sum(sin, z) * sum(cos, z)` all give bitwise equal results either way, so none of them would have caught this. It also makes `symmetric_static == hessian(symmetry_f, x)` robust rather than accidental: that only passed because `n = 9` is below `DEFAULT_CHUNK_THRESHOLD`, so the array path ran a single block. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
`seed_hessian_chunk!` materialised both zeros even when both seeds were supplied. Free
for an isbits value type, but for `BigFloat` at `N = 3` it cost 288 bytes per call,
identical for all four seed combinations. Each `=== nothing` is a compile-time constant,
so deciding per layer folds away:
iseeds oseeds before after
given given 288 0
given - 288 192
- given 288 96
- - 288 288
Also documents the result shape, in the wording `jacobian` uses for its own, and adds a
sentence on why the clearing calls cannot go through `seed_zero_partials!` and why the
`Partials{0}` method of `extract_hessian` is load-bearing -- for a constant `f` the
generic method builds a `0 × length(x)` result, not `length(x) × length(x)`.
New coverage: the extension's `reshape` branch and its two `HESSIAN_ERROR` throws, and
empty inputs on both paths.
Not added: the mixed seed forms in `SeedTest`. Which layer a seed lands in is enforced by
the types -- `oseeds` only fits the outer `Dual` -- and both mixed forms run in every
multi-block sweep, so the bitwise chunk-size test covers them with a failure mode that a
unit test would only relocate. This is unlike `seed_zero_partials!`, whose testset exists
because over-clearing is invisible through the public API.
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 @@## devmotion/structured-extraction-838-839 #843 +/- ##
===========================================================================
+ Coverage 91.54% 92.23% +0.69%
===========================================================================
Files 11 11 Lines 1076 1159 +83 ===========================================================================
+ Hits 985 1069 +84 + Misses 91 90 -1 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
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 freeto join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
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#253. Fixes#836.
Carries #837's two commits unmodified on top of #840, and answers its second review round.
Supersedes #841.
kc/symmetric_hessianis untouched.Based on #840 because the sweep needs its structural positions:
structural_indicescached onthe config, and an
extract_gradient_chunk!that takes its positions fromxrather than fromthe result (#838). On master the same review asks would need three new helpers that #840 already
has, all of which would then conflict with it.
Breaking
For structured inputs the result is now
length(x)bylength(x), both axes indexed by thelinear indices of
x, with zeros in the rows and columns of the structurally zero entries. Sohessian(f, UpperTriangular(rand(3, 3)))is 9x9 with rows and columns 2, 3, 6 zero. That is theJacobian's convention, and it makes
DiffResults.HessianResult(x)the right buffer, sohessian!(::DiffResult, f, x)works for structured inputs instead of erroring. It also makeshessian(f, ::Diagonal)quartic in the size of the diagonal.HessianConfigholdsiseeds,oseeds,dualsandindicesinstead of aJacobianConfiganda
GradientConfig. The sweep never reads the Jacobian config's buffer, which only existed becausethe Hessian used to be
jacobian(gradient(f), x). Dropping it takesHessianConfig(f, x)from1.472 MB to 1.368 MB at
length(x) == 1000, and the result-aware constructor (which built twodead buffers) from 1.576 MB. Separate commit, easy to drop.
Fixes
The result was not reproducible across chunk sizes: diagonal blocks read one triangle,
off-diagonal blocks the other, and the two round differently. At
n = 16every chunk size below16 differed from single-block by ~1e-16. Every chunk size is now bitwise equal to the StaticArrays
path.
hessian!had stopped checking the result shape when the sweep replacedextract_jacobian!, sohessian!(fill(NaN, 4, 4), f, rand(3))silently left row and column 4 untouched.reshape_hessianrestores it, mirroringreshape_jacobian.Smaller items:
seed_hessian_chunk!no longer builds a zero for a layer that was given seeds (288to 0 bytes per call for
BigFloatatN = 3), and the extension'sreshapebranch, its twoHESSIAN_ERRORthrows and empty inputs are now covered.Numbers
Suite 10420/10420. #837's table re-measured against the old algorithm rebuilt on this branch,
hessian!into a preallocated matrix:Three review items done differently
reshape_hessiancheckssizelikereshape_jacobianrather thanaxes, plus onerequire_one_based_indexingfor theDiffResultbuffer that no entry point checks.Diagonalcost is not in the docstring. Neithergradientnorjacobiandocumentsstructured inputs, and the new shape sentence (worded like
jacobian's) makes it derivable.SeedTestdoes not gain the mixed seed forms. Which layer a seed lands in is enforced by thetypes, and both mixed forms run in every multi-block sweep, so the bitwise chunk-size test
already covers them. That differs from
seed_zero_partials!, whose testset exists becauseover-clearing is invisible through the public API.