Skip to content

Keep the BFGS inverse Hessian approximation exactly symmetric #1291

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
pkofod merged 4 commits into from
Oct 28, 2026
Merged

Keep the BFGS inverse Hessian approximation exactly symmetric #1291

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
97d8cd3
Keep the BFGS inverse Hessian approximation exactly symmetric
devmotion Oct 13, 2026
478212e
Move the symmetric rank-2 update into a separate function
devmotion Oct 28, 2026
9ee7283
Merge branch 'master' into dmw/bfgs-symmetric-update
devmotion Oct 28, 2026
82d75df
Import `ishermitian` in tests
devmotion Oct 28, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
24 changes: 4 additions & 20 deletions src/multivariate/solvers/first_order/bfgs.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -204,26 +204,10 @@ function accept_step!(d, state::BFGSState, method::BFGS, options)
c1 = (dx_dg + real(dot(dg, u))) / abs2(dx_dg)
c2 = 1 / dx_dg

# invH = invH + c1 * (s * s') - c2 * (u * s' + s * u')
if (invH isa Array) # i.e. not a CuArray
n = length(dx)
@inbounds for j = 1:n
c1dxj = c1 * dx[j]'
c2dxj = c2 * dx[j]'
c2uj = c2 * u[j]'
for i = 1:n
invH[i, j] = muladd(
dx[i],
c1dxj,
muladd(-u[i], c2dxj, muladd(c2uj, -dx[i], invH[i, j])),
)
end
end
else
mul!(invH, vec(dx), vec(dx)', c1, 1)
mul!(invH, vec(u), vec(dx)', -c2, 1)
mul!(invH, vec(dx), vec(u)', -c2, 1)
end
# invH = invH + c1 * (s * s') - c2 * (u * s' + s * u') = invH + s * w' + w * s'
# with w = c1 / 2 * s - c2 * u, which keeps invH exactly symmetric
u .= (c1 / 2) .* dx .- c2 .* u
invH .+= vec(dx) .* vec(u)' .+ vec(u) .* vec(dx)'

@pkofod pkofod Oct 27, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Gives bitwise the same invH as the broadcast, and matches master's speed for Array:

Suggested change
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

end
return true
end
Expand Down
26 changes: 26 additions & 0 deletions test/multivariate/solvers/first_order/bfgs.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,4 +6,30 @@
skip,
show_name = debug_printing,
)

@testset "Symmetric inverse Hessian approximation" begin
options = Optim.Options(store_trace = true, extended_trace = true)

prob = MVP.UnconstrainedProblems.examples["Extended Rosenbrock"]
res = optimize(
MVP.objective(prob),
MVP.gradient(prob),
prob.initial_x,
BFGS(),
options,
)
@test Optim.converged(res)
@test all(s -> issymmetric(s.metadata["~inv(H)"]), Optim.trace(res))

Random.seed!(0)
n = 4
A = randn(ComplexF64, n, n)
A = A'A + I
b = randn(ComplexF64, n)
fcomplex(x) = real(dot(x, A * x) / 2 - dot(b, x)) + sum(abs.(x) .^ 4)
gcomplex!(g, x) = copyto!(g, A * x - b + 4 * (abs.(x) .^ 2) .* x)
res = optimize(fcomplex, gcomplex!, randn(ComplexF64, n), BFGS(), options)
@test Optim.converged(res)
@test all(s -> ishermitian(s.metadata["~inv(H)"]), Optim.trace(res))
end
end
2 changes: 1 addition & 1 deletion test/runtests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@ using Random
import LineSearches
import NLSolversBase
import NLSolversBase: clear!
import LinearAlgebra: norm, diag, I, Diagonal, dot, eigen, issymmetric, mul!
import LinearAlgebra: norm, diag, I, Diagonal, dot, eigen, ishermitian, issymmetric, mul!
import SparseArrays: normalize!, spdiagm

import ForwardDiff
Expand Down
Loading