Repository navigation
Precision loss in Hessian computation #3104
Description
I'm not sure this is actionable, but it also seems worth doing something about: maybe just documentation, but there are a couple of other options I'll detail below.
This came up while analyzing (with Claude's help) a standard benchmark problem, the CUTEst problem nondquar at n = 5:
(x1 - x2)^2 + (x4 - x5)^2 + Σ_{i=1}^{3} (x_i + x_{i+1} + x_5)^4
Here's the MWE to demonstrate the precision loss:
import MathOptInterface as MOI
import ForwardDiff
# MOI implementation
n = 5
x = [1.0, -1.0, 1.0, -1.0, 1e-6]
v = MOI.VariableIndex.(1:n)
obj = Expr(:call, :+, :(($(v[1]) - $(v[2]))^2), :(($(v[n-1]) - $(v[n]))^2),
(:(($(v[i]) + $(v[i+1]) + $(v[n]))^4) for i in 1:n-2)...)
model = MOI.Nonlinear.Model()
MOI.Nonlinear.set_objective(model, obj)
ev = MOI.Nonlinear.Evaluator(model, MOI.Nonlinear.SparseReverseMode(), v)
MOI.initialize(ev, [:Hess, :HessVec])
st = MOI.hessian_lagrangian_structure(ev)
vals = zeros(length(st))
MOI.eval_hessian_lagrangian(ev, vals, x, 1.0, Float64[])
H = zeros(n, n)
for ((i, j), w) in zip(st, vals)
H[i, j] += w
i == j || (H[j, i] += w)
end
# ForwardDiff implementation
f(z) = (z[1] - z[2])^2 + (z[n-1] - z[n])^2 + sum((z[i] + z[i+1] + z[n])^4 for i in 1:n-2)
Hf = ForwardDiff.hessian(f, x)Now compare the two Hessians:
julia> H - Hf
5×5 Matrix{Float64}:
0.0 0.0 0.0 0.0 0.0
0.0 4.44089e-16 -4.3416e-17 0.0 0.0
0.0 -4.3416e-17 0.0 4.3416e-17 0.0
0.0 0.0 4.3416e-17 0.0 0.0
0.0 0.0 0.0 0.0 4.44089e-16You might think everything is fine: errors are ~eps. But they are against a baseline of 1e-11:
julia> H[4, 3]
1.2000043416036534e-11
julia> Hf[4, 3]
1.2e-11And that's a very large error indeed. At x = [1, -1, 1, -1, 1e-6] every inner sum equals 1e-6 exactly in floating point, so the cross partial ∂²f/∂x3∂x4 = 12·(1e-6)² = 1.2e-11; ForwardDiff's answer is the correct one.
It seems pretty clear this comes from the coloring, and because your coloring is both good and important the paths to fix this are narrow. I currently see three possible options:
- don't fix, just document
- support
Double64, if the user requests it, for intermediate storage - support SparseMatrixColorings (xref #3806), which would allow one to use star coloring
Any thoughts on what you prefer? The documentation or Double64 might be something I could try; I'd probably defer to others for SparseMatrixColorings if they'd be willing to pitch in.
Activity
odow commentedon Oct 15, 2026
I don't think we need to do anything about this. Solvers like Ipopt use tolerances that are much larger anyway. Does the error in the Hessian actually cause any numerical problem?
timholy commentedon Oct 15, 2026
I found it because it did. But this was a special case: I was working on a certified descent method, and when you have a numerical result that violates a theorem, you check it. But it doesn't mean it would cause trouble in general.
mlubin commentedon Oct 15, 2026
I'd be supportive of integrating SparseMatrixColorings.
timholy commentedon Oct 15, 2026
Pinging an expert on that option: @gdalle
timholy commentedon Oct 15, 2026
I should also emphasize that this is low-stakes for me: the easy workaround for my certified descent method is just to switch to other ways of computing the sparse Hessian, so no one will offend me if you decide not to do it.
gdalle commentedon Oct 15, 2026
I can probably help integrate SMC's star coloring to JuMP but I'm not sure it will be faster than edge pushing in general.
timholy commentedon Oct 15, 2026
I agree, it's probably not faster. I was wondering more about making the coloring choice an option, but I am entirely naive about the details.
gdalle commentedon Oct 15, 2026
Pinging @blegat since he's been working with SMC in ArrayDiff.jl already
amontoison commentedon Oct 16, 2026
We discussed about this issue this morning with @gdalle.
I did a quick test with star and acyclic coloring:
I don't have the same error for the acyclic coloring.
The BFS traversal of the trees is probably different.
One option for MOI.jl is to perform an iterative refinement when the small triangular systems (defined by the two-colored trees) are solved during the decompression.
timholy commentedon Oct 16, 2026
We'd also want some benchmarks. I haven't tested thoroughly, but so far nothing I've tried has been as fast as MOI's edge-pushing approach. But this wasn't a focus, so I wouldn't take it too seriously. Let me know if I can help with that aspect.
blegat commentedon Oct 16, 2026
Star coloring should be usable from JuMP (using SparseMatrixColorings) using #3105