Fix LBFGS/BFGS callback receiving Dual numbers instead of scalar loss values

Fixes #1073

When using LBFGS or BFGS with bounds, Optim.jl wraps the optimizer in Fminbox,
which may use ForwardDiff internally for gradient computation. This resulted in
the callback receiving ForwardDiff.Dual numbers instead of scalar loss values,
causing incorrect (sometimes negative) values to be reported.

Changes:
- Added ForwardDiff as a dependency in OptimizationOptimJL
- Added _scalar_value() utility function to extract scalar values from Dual numbers
- Updated all three _cb callback functions to extract scalar values before passing to user callbacks
- Added comprehensive test case verifying callbacks receive correct scalar non-negative values

🤖 Generated with [Claude Code](https://claude.com/claude-code)

Co-Authored-By: Claude <[email protected]>

Author: ChrisRackauckas

lib/OptimizationOptimJL/Project.toml

@@ -3,6 +3,7 @@ uuid = "36348300-93cb-4f02-beb5-3c3902f8871e"
 authors = ["Vaibhav Dixit <[email protected]> and contributors"]
 version = "0.4.6"
 [deps]
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 OptimizationBase = "bca83a33-5cc9-4baa-983d-23429ab6bcbb"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
@@ -11,7 +12,6 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 
 [extras]
-ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
@@ -28,4 +28,4 @@ Reexport = "1.2"
 SciMLBase = "2.58"
 
 [targets]
-test = ["ForwardDiff", "ModelingToolkit", "Random", "ReverseDiff", "Test", "Zygote"]
+test = ["ModelingToolkit", "Random", "ReverseDiff", "Test", "Zygote"]

lib/OptimizationOptimJL/src/OptimizationOptimJL.jl

@@ -3,9 +3,19 @@ module OptimizationOptimJL
 using Reexport
 @reexport using Optim, OptimizationBase
 using SciMLBase, SparseArrays
+
+# Import ForwardDiff to handle Dual numbers in callbacks
+import ForwardDiff
+
 decompose_trace(trace::Optim.OptimizationTrace) = last(trace)
 decompose_trace(trace::Optim.OptimizationState) = trace
 
+# Extract scalar value from potentially Dual-valued trace values
+# This is needed because Optim.jl may use ForwardDiff internally for gradient computation,
+# resulting in Dual numbers in the trace, but callbacks should receive scalar values
+_scalar_value(x) = x  # Default case for regular numbers
+_scalar_value(x::ForwardDiff.Dual) = ForwardDiff.value(x)  # Extract value from Dual
+
 SciMLBase.allowsconstraints(::IPNewton) = true
 SciMLBase.allowsbounds(opt::Optim.AbstractOptimizer) = true
 SciMLBase.allowsbounds(opt::Optim.SimulatedAnnealing) = false
@@ -149,14 +159,16 @@ function SciMLBase.__solve(cache::OptimizationBase.OptimizationCache{
         trace_state = decompose_trace(trace)
         metadata = trace_state.metadata
         θ = metadata[cache.opt isa Optim.NelderMead ? "centroid" : "x"]
+        # Extract scalar value from potentially Dual-valued trace (issue #1073)
+        loss_val = _scalar_value(trace_state.value)
         opt_state = OptimizationBase.OptimizationState(iter = trace_state.iteration,
             u = θ,
             p = cache.p,
-            objective = trace_state.value,
+            objective = loss_val,
             grad = get(metadata, "g(x)", nothing),
             hess = get(metadata, "h(x)", nothing),
             original = trace)
-        cb_call = cache.callback(opt_state, trace_state.value)
+        cb_call = cache.callback(opt_state, loss_val)
         if !(cb_call isa Bool)
             error("The callback should return a boolean `halt` for whether to stop the optimization process.")
         end
@@ -270,14 +282,16 @@ function SciMLBase.__solve(cache::OptimizationBase.OptimizationCache{
         θ = !(cache.opt isa Optim.SAMIN) && cache.opt.method == Optim.NelderMead() ?
             metadata["centroid"] :
             metadata["x"]
+        # Extract scalar value from potentially Dual-valued trace (issue #1073)
+        loss_val = _scalar_value(trace_state.value)
         opt_state = OptimizationBase.OptimizationState(iter = trace_state.iteration,
             u = θ,
             p = cache.p,
-            objective = trace_state.value,
+            objective = loss_val,
             grad = get(metadata, "g(x)", nothing),
             hess = get(metadata, "h(x)", nothing),
             original = trace)
-        cb_call = cache.callback(opt_state, trace_state.value)
+        cb_call = cache.callback(opt_state, loss_val)
         if !(cb_call isa Bool)
             error("The callback should return a boolean `halt` for whether to stop the optimization process.")
         end
@@ -357,14 +371,16 @@ function SciMLBase.__solve(cache::OptimizationBase.OptimizationCache{
 
     function _cb(trace)
         metadata = decompose_trace(trace).metadata
+        # Extract scalar value from potentially Dual-valued trace (issue #1073)
+        loss_val = _scalar_value(trace.value)
         opt_state = OptimizationBase.OptimizationState(iter = trace.iteration,
             u = metadata["x"],
             p = cache.p,
             grad = get(metadata, "g(x)", nothing),
             hess = get(metadata, "h(x)", nothing),
-            objective = trace.value,
+            objective = loss_val,
             original = trace)
-        cb_call = cache.callback(opt_state, trace.value)
+        cb_call = cache.callback(opt_state, loss_val)
         if !(cb_call isa Bool)
             error("The callback should return a boolean `halt` for whether to stop the optimization process.")
         end

lib/OptimizationOptimJL/test/runtests.jl

@@ -199,6 +199,35 @@ end
     @test_throws ArgumentError (sol = solve(prob, Optim.BFGS())) isa Any  # test exception is thrown
     @test 10 * sol.objective < l1
 
+    # Test for issue #1073: callbacks should receive scalar non-negative loss values
+    # when using (L)BFGS with bounds and automatic differentiation
+    @testset "Issue #1073: LBFGS/BFGS callback receives correct scalar loss with bounds" begin
+        # Create a non-negative loss function (sum of squares)
+        loss_vals = Float64[]
+        function test_callback(state, loss_val)
+            # Verify loss_val is a scalar, not a Dual number
+            @test loss_val isa Real
+            @test !(loss_val isa ForwardDiff.Dual)
+            # For a sum-of-squares loss, values should be non-negative
+            push!(loss_vals, loss_val)
+            return false
+        end
+
+        # Test with LBFGS + bounds (triggers Fminbox wrapping)
+        optprob = OptimizationFunction(rosenbrock, OptimizationBase.AutoForwardDiff())
+        prob = OptimizationProblem(optprob, x0, _p; lb = [-1.0, -1.0], ub = [0.8, 0.8])
+        empty!(loss_vals)
+        sol = solve(prob, Optim.LBFGS(); callback = test_callback, maxiters = 10)
+        @test all(>=(0), loss_vals)  # All loss values should be non-negative
+        @test length(loss_vals) > 0  # Callback should have been called
+
+        # Test with BFGS + bounds
+        empty!(loss_vals)
+        sol = solve(prob, Optim.BFGS(); callback = test_callback, maxiters = 10)
+        @test all(>=(0), loss_vals)  # All loss values should be non-negative
+        @test length(loss_vals) > 0  # Callback should have been called
+    end
+
     @testset "cache" begin
         objective(x, p) = (p[1] - x[1])^2
         x0 = zeros(1)