abstol and or reltol work differently with ForwardDiff vs FiniteDiff

[elbert5770]

Using very small values for abstol and reltol (1e-14) with DynamicSS on a SteadyStateProblem, FiniteDiff solves out to t = 70,000, while ForwardDiff stops at t = 228. The true gradient is [0.0, -296.29629633013286], and FiniteDiff gets that with [-2.6518679901607456e-7, -296.2962963356087], while ForwardDiff is [0.0015469937539040605, -296.29629629628164], good but not great. The gradients are calculated at the bottom of the code with FiniteDiff.finite_difference_gradient and ForwardDiff.gradient

Are the tolerances actually being met with ForwardDiff? Is there a maxiter in the termination conditions in DiffEqBase that is getting hit? Reducing the tolerances to 1e-16 crashes FiniteDiff but ForwardDiff improves to [1.2214211430632144e-5, -296.2962962962962] and goes to t = 1219.

While DynamicSS is not needed on this MWE, it very often is the best choice.

using FiniteDiff
using ForwardDiff
using Optimization
using OptimizationOptimJL
using NonlinearSolve
using DifferentialEquations
using SteadyStateDiffEq
using Plots

function sys!(du, u, (p,p_fix,kf,plateau,kd,start_time,end_time), t)
    #@show t,forcing_function(t,plateau,kd,start_time,end_time),p,p_fix,kf
    du[1] = kf[1]*(1.0 - forcing_function(t,plateau,kd,start_time,end_time)) - p[1]*u[1]
    du[2] = p[1]*u[1] - p_fix[1]*u[2]
    du[3] = p_fix[1]*u[2] - p[2]*u[3]
    if t > 50
     #   @show u
      #  @show du
      @show t
    end
    return nothing
end

function forcing_function(t,plateau,kd,start_time,end_time)
    if t < start_time
        return plateau*t
    else
        if t <= end_time
            return plateau
        else
            if t > end_time && t < Inf
                return plateau*(exp(-kd*(t-end_time)))
            else
                return 0.0
            end
        end
    end
end



function set_kf_outer(probss,u03)
    function set_kf_inner(kf_guess,(u0,p_guess,p_fix,plateau,kd,start_time,end_time))
        # Solve for steady state with current p_guess and an initial guess of kf and u0
        plateau = 0.0
        sol = solve(remake(probss,u0=eltype(kf_guess).(u0),p=(p_guess,p_fix,kf_guess,plateau,kd,start_time,end_time)),DynamicSS(AutoTsit5(Rosenbrock23()),abstol = 1e-14,reltol = 1e-14))
        # Compare kf at steady state to measured u0[3]
        # @show sol.u
        # @show sol.resid
        # @show sol.prob
        @show (sol[3]-u03)^2
        return (sol[3]-u03)^2
    end
end

function main()
    #### Parameters and state variables
    kf0= [10.0] #true value of zero order production rate
    p0 = [1.0,0.25,0.5] #true values of parameters
    p_fix = [p0[2]] # a parameter that is known, note that the system is otherwise unidentifiable
    p_var = [p0[1],p0[3]] # Unknown parameters
    u0_guess = [1000.0,750.0,20.0] #initial guess at steady state of state variables, except u0[3] is known to be 20.0
    u03 = 20.0 # measured value
    plateau = 50.0 # Zero order production rate changes at time = 0, reaches a fractional change given by plateau at start_time
    kd = 1.0 # Rate constant for exponential decay back to original value starting at end_time
    start_time = 1.0
    end_time = 10.0

    #### Set up ODE problem
    t_end = 36.0 # Length of simulation
    tspan = (0.0,t_end) # timespan of simulation
    probODE = ODEProblem(sys!,u0_guess,tspan,(p_var,p_fix,kf0,plateau,kd,start_time,end_time))
    
    #### Set up steady state problem
    probss = SteadyStateProblem(probODE)
    sol_ss = solve(probss)
    u0_ss = sol_ss.u
    @show "u0_ss true",sol_ss.u
    
    #### Initial guesses
    p_guess = [0.1,3.0] #inital guess at values in p_val
    kf_guess = [100.0] #initial guess at kf0
    

    #### Set up optimization problem to determine kf for current value of p_guess
    set_kf = set_kf_outer(probss,u03)
    optkf = OptimizationFunction(set_kf,Optimization.AutoForwardDiff())
    prob_set_kf = OptimizationProblem(optkf,kf_guess,(u0_guess,p_guess,p_fix,plateau,kd,start_time,end_time))
    
    @show FiniteDiff.finite_difference_gradient(x -> set_kf(kf_guess,(u0_guess,x,p_fix,plateau,kd,start_time,end_time)),p_guess)
    @show ForwardDiff.gradient(x -> set_kf(kf_guess,(u0_guess,x,p_fix,plateau,kd,start_time,end_time)),p_guess )

    return nothing
end

main() ```

[ChrisRackauckas]

@avik-pal could this be due to the new termination conditions?

[avik-pal]

That is not entirely surprising to me. The tolerances for time-stepping specified via solve but tolerances for termination is specified via abstol and reltol in DynamicSS

SteadyStateDiffEq.jl/src/algorithms.jl

Line 47 in 063e958

termination_condition = NLSolveTerminationCondition(NLSolveTerminationMode.SteadyStateDefault; abstol, reltol))

[elbert5770]

The abstol and reltol are inside DynamicSS, not attached to solve. Moving them outside of DynamicSS into solve makes the termination much worse, as expected. But these are definitely inside DynamicSS and affecting termination.

[elbert5770]

Corrected in edit: xxx[Okay, I think the root issues are the aggressive time step growth and the reltol = 1e-6 is too large. The delta t is growing fast and the residuals are bouncing around. Setting dtmax in solve and lowering reltol to say 1e-10 leads to a well-behaved algorithm for both the FiniteDiff and ForwardDiff gradient. sol = solve(remake(probss,u0=eltype(kf_guess).(u0),p=(p_guess,p_fix,kf_guess,plateau,kd,start_time,end_time)),DynamicSS(AutoTsit5(Rosenbrock23()),reltol=1e-10),dtmax=1.0)] xxx

Nope, this only worked because I had changed p_guess to [10.0,3.0]. With the original p_guess, it terminates too quickly with reltol=1e-10.

As Avik mentioned, the time stepping can be affected by the abstol and reltol to solve, so this might be better than dtmax:

sol = solve(remake(probss,u0=eltype(kf_guess).(u0),p=(p_guess,p_fix,kf_guess,plateau,kd,start_time,end_time)),DynamicSS(AutoTsit5(Rosenbrock23()),abstol=1e-16,reltol=1e-16,tspan=1e6),abstol=1e-11,reltol=1e-11) 

Here, unfortunately, the solver must run to t=1e6 to get good results and avoid early termination consistently, or hit these very low values of reltol and abstol within DynamicSS.

[ChrisRackauckas]

Is this fixed with the new termination conditions? https://docs.sciml.ai/NonlinearSolve/stable/basics/TerminationCondition/

[elbert5770]

It appears to be exactly the same, but I'm only 99% sure of my implementation of TerminationCondition:

using FiniteDiff
using ForwardDiff
using Optimization
using OptimizationOptimJL
using NonlinearSolve
using DifferentialEquations
using SteadyStateDiffEq
using Plots

function sys!(du, u, (p,p_fix,kf,plateau,kd,start_time,end_time), t)
    #@show t,forcing_function(t,plateau,kd,start_time,end_time),p,p_fix,kf
    du[1] = kf[1]*(1.0 - forcing_function(t,plateau,kd,start_time,end_time)) - p[1]*u[1]
    du[2] = p[1]*u[1] - p_fix[1]*u[2]
    du[3] = p_fix[1]*u[2] - p[2]*u[3]
    if t > 50
     #   @show u
      #  @show du
      @show t
    end
    return nothing
end

function forcing_function(t,plateau,kd,start_time,end_time)
    if t < start_time
        return plateau*t
    else
        if t <= end_time
            return plateau
        else
            if t > end_time && t < Inf
                return plateau*(exp(-kd*(t-end_time)))
            else
                return 0.0
            end
        end
    end
end



function set_kf_outer(probss,u03,termination_condition)
    function set_kf_inner(kf_guess,(u0,p_guess,p_fix,plateau,kd,start_time,end_time))
        # Solve for steady state with current p_guess and an initial guess of kf and u0
        plateau = 0.0
        #sol = solve(remake(probss,u0=eltype(kf_guess).(u0),p=(p_guess,p_fix,kf_guess,plateau,kd,start_time,end_time)),DynamicSS(AutoTsit5(Rosenbrock23()),abstol = 1e-14,reltol = 1e-14))
        
        sol = solve(remake(probss,u0=eltype(kf_guess).(u0),p=(p_guess,p_fix,kf_guess,plateau,kd,start_time,end_time)),DynamicSS(AutoTsit5(Rosenbrock23()),termination_condition=termination_condition))
        # Compare kf at steady state to measured u0[3]
        # @show sol.u
        # @show sol.resid
        # @show sol.prob
        @show (sol[3]-u03)^2
        return (sol[3]-u03)^2
    end
end

function main()
    #### Parameters and state variables
    kf0= [10.0] #true value of zero order production rate
    p0 = [1.0,0.25,0.5] #true values of parameters
    p_fix = [p0[2]] # a parameter that is known, note that the system is otherwise unidentifiable
    p_var = [p0[1],p0[3]] # Unknown parameters
    u0_guess = [1000.0,750.0,20.0] #initial guess at steady state of state variables, except u0[3] is known to be 20.0
    u03 = 20.0 # measured value
    plateau = 50.0 # Zero order production rate changes at time = 0, reaches a fractional change given by plateau at start_time
    kd = 1.0 # Rate constant for exponential decay back to original value starting at end_time
    start_time = 1.0
    end_time = 10.0

    #### Set up ODE problem
    t_end = 36.0 # Length of simulation
    tspan = (0.0,t_end) # timespan of simulation
    probODE = ODEProblem(sys!,u0_guess,tspan,(p_var,p_fix,kf0,plateau,kd,start_time,end_time))
    
    #### Set up steady state problem
    probss = SteadyStateProblem(probODE)
    sol_ss = solve(probss)
    u0_ss = sol_ss.u
    @show "u0_ss true",sol_ss.u
    
    #### Initial guesses
    p_guess = [0.1,3.0] #inital guess at values in p_val
    kf_guess = [100.0] #initial guess at kf0
    
    termination_condition = NLSolveTerminationCondition(NLSolveTerminationMode.NLSolveDefault;
        abstol = 1e-14,
        reltol = 1e-14)
    #### Set up optimization problem to determine kf for current value of p_guess
    set_kf = set_kf_outer(probss,u03,termination_condition)
    optkf = OptimizationFunction(set_kf,Optimization.AutoForwardDiff())
    prob_set_kf = OptimizationProblem(optkf,kf_guess,(u0_guess,p_guess,p_fix,plateau,kd,start_time,end_time))
    
    @show FiniteDiff.finite_difference_gradient(x -> set_kf(kf_guess,(u0_guess,x,p_fix,plateau,kd,start_time,end_time)),p_guess)
    @show ForwardDiff.gradient(x -> set_kf(kf_guess,(u0_guess,x,p_fix,plateau,kd,start_time,end_time)),p_guess )

    return nothing
end

main()