Repository navigation
Keep the BFGS inverse Hessian approximation exactly symmetric - #1291
Conversation
Replace the `invH isa Array` loop and the three `mul!` rank-1 updates with a single symmetric rank-2 broadcast. Both previous code paths left `invH` slightly asymmetric after every step, and non-`Array` CPU arrays such as `FixedSizeArray` took the `mul!` path meant for GPU arrays. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Benchmark Results (Julia v1)Time benchmarks
Memory benchmarks
A plot of the benchmark results has been uploaded as an artifact at . |
Benchmark Results (Julia vlts)Time benchmarks
Memory benchmarks
A plot of the benchmark results has been uploaded as an artifact at . |
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## master #1291 +/- ##
==========================================
+ Coverage 88.20% 88.24% +0.04%
==========================================
Files 52 52
Lines 5080 5072 -8
==========================================
- Hits 4481 4476 -5
+ Misses 599 596 -3 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
pkofod
left a comment
pkofod
left a comment
There was a problem hiding this comment.
The symmetric update is correct and I'd like to keep it. I checked exact symmetry on all 17 unconstrained OptimTestProblems and on 150 random complex problems: 0 asymmetric (or non-Hermitian) states with this PR, against about 99% on master, and every problem still converges.
The broadcast is a performance regression for Matrix, though. BFGS time per iteration on Extended Rosenbrock (analytic O(n) gradient, so the O(n²) update dominates), Ryzen 5 2600, one thread, relative to the base commit:
| Julia | n = 100 | n = 500 | n = 1000 |
|---|---|---|---|
| 1.10 | 1.76x | 1.94x | 1.48x |
| 1.11 | 1.03x | 1.05x | 1.02x |
| 1.12 | 1.76x | 1.95x | 1.54x |
| 1.13 | 1.82x | 1.95x | 1.52x |
An M-series Mac shows the same pattern. On 1.12 and 1.13 the update block alone is about 2x slower than the muladd loop from n = 32 on. On aarch64, LLVM emits 9 to 11 vectorized loop versions for this 2D broadcast on 1.12 and 1.13, against 38 to 46 on 1.10 and 1.11. On 1.10 the update block alone is only about 10% slower at n = 100 to 200, and I haven't found why it costs this much end to end. The CI benchmark misses all of it because BFGS only runs on the 2-parameter Gabor problem.
Suggestion: keep your formula, but restore the Array branch with a plain loop (inline suggestion below). It gives bitwise the same invH as the broadcast for real and complex inputs, so your argument and test still hold. Other arrays keep the broadcast, which for FixedSizeArray is still faster than master's three mul! calls (up to 12x at small n, 1.3 to 1.6x at n ≥ 50 on x86). End to end, the loop stays within 0.98 to 1.03 of master at n ≥ 100 on all four versions. It must not use muladd, since fusing rounds (i, j) and (j, i) differently.
The loop needs @inbounds: without it, it is up to 3x slower than master on 1.11 and 1.12 for n ≥ 16. If you'd rather avoid @inbounds, this column-wise version needs none and cost up to 8% per iteration in my runs:
for j in axes(invH, 2)
sj = dx[j]'
wj = u[j]'
@views invH[:, j] .+= vec(dx) .* wj .+ vec(u) .* sj
endIt would also be worth adding BFGS to the Rosenbrock benchmarks that scale with n, so CI catches this kind of change.
The benchmarks and this review were prepared with the assistance of Claude Code.
| u .= (c1 / 2) .* dx .- c2 .* u | ||
| invH .+= vec(dx) .* vec(u)' .+ vec(u) .* vec(dx)' |
There was a problem hiding this comment.
Gives bitwise the same invH as the broadcast, and matches master's speed for Array:
| u .= (c1 / 2) .* dx .- c2 .* u | |
| invH .+= vec(dx) .* vec(u)' .+ vec(u) .* vec(dx)' | |
| u .= (c1 / 2) .* dx .- c2 .* u | |
| if invH isa Array | |
| # for Array a loop is much faster than the broadcast below on Julia 1.10 and 1.12+ | |
| @inbounds for j = 1:length(dx) | |
| sj = dx[j]' | |
| wj = u[j]' | |
| for i = 1:length(dx) | |
| invH[i, j] += dx[i] * wj + u[i] * sj | |
| end | |
| end | |
| else | |
| invH .+= vec(dx) .* vec(u)' .+ vec(u) .* vec(dx)' | |
| end |
|
So performance is not impacted on Julia 1.11? Generally, given all the calculations involved in each iteration, is this arguably limited performance regression relevant in practice? Does the comparison need |
yes, as far as I understand it, 1.11 compiles some optimized version for your formula that are not there on v1.12 and v1.13.
That depends |
|
A note on the 1.11 answer. It comes from counting the vectorized loop bodies (
The loop on master, the I only counted them. I didn't check which ones were generated, so this doesn't show which version runs on the hot path, and I haven't found the compiler change behind the drop. It also doesn't explain 1.10: it has the vectorized versions but is still 1.48 to 1.94x slower per iteration at n ≥ 100 in the table above. Prepared with the assistance of Claude Code. |
|
I can try to see which ones were generated. |
|
I checked which versions run. I sampled instruction addresses with the profiler and mapped them onto the JIT'd machine code of
So the vectorized versions are generated on every version, but on 1.10, 1.12 and 1.13 they never run for this call. Why. On 1.11 the load and the store share one address and that check doesn't exist. In an isolated benchmark of the PR's two lines, lowering the limit with A second, smaller effect on 1.12 and 1.13. When the vector body does run, the loop vectorizer estimates its cost at 20 to 22 (18 to 19 on 1.11). That is at or above its Time per iteration, PR / base commit, n = 100 / 200 / 500, with the flags passed through
The flags leave the base commit unchanged within 1%. The n = 500 runs hit the 1000-iteration limit. Still open: why 1.12 hits the unswitch limit where 1.11 doesn't. The 1.12 IR computes the addresses with byte-offset GEPs and more instructions, but I haven't shown that this is what raises the cost. I only looked at the machine code on Apple Silicon. So to the 1.11 question: 1.11 is fast here because it happens to unswitch this loop. The same broadcast can fall back to the scalar loop on any version, depending on the code around it. Prepared with the assistance of Claude Code. |
|
How large is this compared to a whole Whole run, PR / base (Apple M2 Pro):
On 1.12 that is 5.6 to 9.7 µs per iteration at n = 100 and 24 to 41 µs at n = 200. A full run at n = 200 goes from 1.9 ms to 3.3 ms. The added time is a fixed amount per iteration and does not depend on Prepared with the assistance of Claude Code. |
|
we can condition the branch on the Julia version ;) |
Inlined into `accept_step!`, LLVM does not unswitch the broadcast loop on Julia 1.10 and 1.12+, so its vectorized body is never entered and BFGS iterations become up to 1.7x slower than master. In a separate function the loop is vectorized again, without an `Array` branch or `@inbounds`. Also benchmark BFGS on Rosenbrock with n = 2, 20, 100, where the O(n^2) update dominates for larger n. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
|
Thanks for tracking down the unswitching. That points to a fix without a branch: in 478212e the broadcast moved into a small function of its own, BFGS time per iteration on Extended Rosenbrock relative to master, n = 100 to 500, Apple M2 Pro:
This still depends on what LLVM decides, so I added BFGS on Rosenbrock with n = 2, 20 and 100 to the benchmarks. CI should catch it if this comes back. If we ever need a branch, I don't think Prepared with the assistance of Claude Code. |
The inverse Hessian update in
BFGShas two code paths, and neither keepsinvHsymmetric in floating point:invH isa Array, the fusedmuladdloop roundsinvH[i, j]andinvH[j, i]differently, so they drift apart by a few ulps after every step. For complex inputs, the diagonal also picks up nonzero imaginary parts.mul!branch, which was meant for GPU arrays (# i.e. not a CuArray). This includes CPU arrays such asFixedSizeArray(FixedSizeArrays.jl), whoseinvHis aFixedSizeMatrixbecause_init_identity_matrixbroadcasts overx. The branch does three separately rounded rank-1 updates. The first,mul!(invH, dx, dx', c1, 1), also changes its implementation depending on that asymmetry: withβ = 1, LinearAlgebra callsBLAS.syrk!only ifissymmetric(invH)holds bitwise, andBLAS.gemm!otherwise. SoMatrixandFixedSizeMatrixinputs follow different arithmetic, and the non-Arraypath depends on ulp-level noise.This PR replaces both branches with a single symmetric rank-2 update. It uses
c1 * s * s' - c2 * (u * s' + s * u') = s * w' + w * s'withw = c1 / 2 * s - c2 * u, written as one broadcast intoinvHthat reuses theubuffer. Element(i, j)and element(j, i)are computed from the same products in commuted order, soinvHstays exactly symmetric (Hermitian for complex inputs). One code path then servesArray,FixedSizeArray, and GPU arrays, without theisa Arraycheck and in a single pass.On a 16-dimensional SPD quadratic over 200 steps:
invHBigFloatmuladdloop (Array)mul!(other arrays, e.g.FixedSizeArray)The update is mathematically unchanged, but rounding differs, so BFGS iterates can change in the last bits.
A test checks that every
~inv(H)in the extended trace is exactly symmetric (real) or Hermitian (complex). It fails on master.This pull request was written with the assistance of generative AI (Claude Code).
🤖 Generated with Claude Code