Precision loss in Hessian computation
Les mainteneurs répondent en général sous 1 jour
Évaluation
- Difficulté
- 5/5
- Temps estimé
- Plus d'une semaine
- Accessibilité débutants
- 15/100
- Type d'issue
- Bug
- Clarté
- À clarifier
- Activité
- Active
- Domaine
- backend
Piste de recherche
Start by running the included MWE, which compares the sparse Hessian from MOI.eval_hessian_lagrangian on MOI.Nonlinear.SparseReverseMode() with ForwardDiff.hessian, and check the entry H[4, 3] (about 1.2e-11 versus 1.2e-11 with a larger absolute error). The cause is the coloring-based Hessian recovery, not the derivative formulas. The maintainers have not chosen a fix among the three options (document, Double64 intermediates, or SparseMatrixColorings star coloring), so decide on the approach with them before writing code.
Rédigé par le modèle d'indexation à partir du texte de l'issue.
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-16
You 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-11
And 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.
- Langage dominant
- Julia
- Étoiles
- 506
- Forks
- 103
- Merge moyen
- 1 j 1 h
- PR mergées (30 j)
- 9
Préparer son environnement
Ce projet ne fournit ni conteneur de développement, ni Dockerfile, ni guide de contribution : l'installation est à votre charge. Commencez par son README, et consultez notre guide de la première contribution pour les étapes générales.
Par où commencer
- Lisez l'issue en entier, puis le guide de contribution du projet.
- Signalez en commentaire que vous la prenez — cela évite que deux personnes fassent le même travail.
- Forkez le dépôt et travaillez sur une branche.
- Ouvrez une pull request qui référence le numéro de l'issue.
Autres issues de jump-dev/MathOptInterface.jl
-
Difficulté 2/5 1-3 heures Accessibilité débutants 78/100
jump-dev/MathOptInterface.jl#3113 ·
Les mainteneurs répondent en général sous 1 jour
-
[Bridges] MultirowChange in set_map.jl drops zero coefficientsPeut-être pris @odow l’a pris il y a 6 jours. Ouverte
jump-dev/MathOptInterface.jl#3106 · 2 commentaires · 1 personne assignée ·
Les mainteneurs répondent en général sous 1 jour
-
Layer instantiatorOuverte
Difficulté 5/5 Plus d'une semaine Accessibilité débutants 35/100
jump-dev/MathOptInterface.jl#2986 · 3 commentaires ·
Les mainteneurs répondent en général sous 1 jour
-
Differentiable Optimization APIOuverte
Difficulté 5/5 Plus d'une semaine Accessibilité débutants 25/100
jump-dev/MathOptInterface.jl#2969 · 2 commentaires · 1 réaction ·
Les mainteneurs répondent en général sous 1 jour
-
Add Batched{S} setOuverteType: Set Request
Difficulté 5/5 Plus d'une semaine Accessibilité débutants 35/100
jump-dev/MathOptInterface.jl#2904 · 15 commentaires · 7 réactions ·
Les mainteneurs répondent en général sous 1 jour
Toutes les issues de jump-dev/MathOptInterface.jl
Issues similaires
-
Difficulté 1/5 Moins d'une heure Accessibilité débutants 83/100
JuliaMath/Polynomials.jl#647 ·
-
Difficulté 2/5 1-3 heures Accessibilité débutants 74/100
OpenSourceAWE/SymbolicAWEModels.jl#431 ·
Les mainteneurs répondent en général sous 1 jour
-
Make PDAQP an optional extraOuverte
Difficulté 2/5 1-3 heures Accessibilité débutants 66/100
-
Difficulté 2/5 1-3 heures Accessibilité débutants 75/100
JuliaArrays/StaticArrays.jl#1365 ·
-
`@fastmath` reductions are 2–5× slower: no `neutral_element` for `Base.FastMath.max_fast` and co.Ouverte
Difficulté 2/5 1-3 heures Accessibilité débutants 86/100
JuliaGPU/GPUArrays.jl#833 ·
Les mainteneurs répondent en général sous 1 jour