GSOC Goal: Nonlinear optimization

[msricher]

Explore nonlinear optimizers, both local/"quasi-Newton" and global/stochastic, for the exact solver. For instance, see https://github.com/jbrea/CMAEvolutionStrategy.jl, https://jump.dev/, https://github.com/SciML/Optimization.jl.

[PaulWAyers]

It may be possible to solve the equations analytically. One can always write them as a system of polynomial equations, then use a package to solve that. Which package is best (as in fastest; all of the below should work) is a question that is not obvious to me up front. My guess is that the homotopy methods will work best, especially if one accepts a numerical (not forcing a rational) solution. I'd start there. (However, because our polynomials are relatively low degree, it is possible, however, that the Groebner-basis-based-methods (first and third links) are faster for our problem.)

https://github.com/[algebraic-solving/AlgebraicSolving.jl](https://github.com/algebraic-solving/AlgebraicSolving.jl)
https://www.juliahomotopycontinuation.org/examples/nash/
https://docs.sciml.ai/Symbolics/stable/manual/solver/#Symbolics.symbolic_solve

I've referred to tutorials/examples where I found them helpful, but all three packages are on GitHub. One can also use sage but that is a bit more unpleasant to interface to (even if perhaps it has more available methods).

[sparsh-karna]

Hi @msricher, @PaulWAyers,

I’ve completed a comparison of quasi-Newton solvers from NonlinearSolve.jl for the exact solver in NICE.jl, focusing on "System 1" (3 species, 2 reactions). The system parameters are:

• Stoichiometry: [-0.5 -0.5; 1.0 -1.0; 0.0 1.0] (Reaction 1: 0.5A → B, Reaction 2: B → C)
• Equilibrium Constants: [1.0, 0.1]
• Initial Concentrations: [1.0, 0.2, 0.4]
• Expected Equilibrium Concentrations: [0.9261879203, 0.9623865752, 0.0926187920]
• Tolerances: atol=1.0e-9, maxiters=1000, reltol=0.0 for solve

I tested 13 solver configurations, varying algorithms, linear solvers, and line search methods. The configurations and their parameters are:

• TrustRegion (Simple): TrustRegion(radius_update_scheme=RadiusUpdateSchemes.Simple), direct linear solver (linsolve=nothing).
• TrustRegion (NocedalWright): TrustRegion(radius_update_scheme=RadiusUpdateSchemes.NocedalWright), direct linear solver (linsolve=nothing).
• TrustRegion (Simple, Krylov): TrustRegion(radius_update_scheme=RadiusUpdateSchemes.Simple), Krylov linear solver (linsolve=KrylovJL_GMRES()).
• NewtonRaphson (Basic): NewtonRaphson(), direct linear solver (linsolve=nothing), no line search.
• NewtonRaphson (BackTracking): NewtonRaphson(linesearch=BackTracking()), direct linear solver (linsolve=nothing).
• NewtonRaphson (HagerZhang): NewtonRaphson(linesearch=HagerZhang()), direct linear solver (linsolve=nothing).
• NewtonRaphson (Basic, Krylov): NewtonRaphson(), Krylov linear solver (linsolve=KrylovJL_GMRES()), no line search.
• Broyden (Identity, Good): Broyden(init_jacobian=Val(:identity), update_rule=Val(:good_broyden)), direct linear solver (linsolve=nothing).
• Broyden (True Jacobian, Good): Broyden(init_jacobian=Val(:true_jacobian), update_rule=Val(:good_broyden)), direct linear solver (linsolve=nothing).
• Broyden (Identity, Bad): Broyden(init_jacobian=Val(:identity), update_rule=Val(:bad_broyden)), direct linear solver (linsolve=nothing).
• GaussNewton: GaussNewton(), direct linear solver (linsolve=nothing).
• GaussNewton (Krylov): GaussNewton(), Krylov linear solver (linsolve=KrylovJL_GMRES()).
• RobustMultiNewton: RobustMultiNewton(), direct linear solver (linsolve=nothing).

The results are stored in solver_comparison_system1.csv:

Algorithm	Success	Time (ms)	Iterations	Residual	Accuracy Check
TrustRegion (Simple)	✅ Pass	0.016458	5	3.89e-15	Accurate
TrustRegion (NocedalWright)	✅ Pass	0.016583	5	3.89e-15	Accurate
TrustRegion (Simple, Krylov)	✅ Pass	0.016541	5	3.89e-15	Accurate
NewtonRaphson (Basic)	✅ Pass	0.012000	3	2.80e-14	Accurate
NewtonRaphson (BackTracking)	✅ Pass	0.012000	3	2.80e-14	Accurate
NewtonRaphson (HagerZhang)	✅ Pass	0.011959	3	2.80e-14	Accurate
NewtonRaphson (Basic, Krylov)	✅ Pass	0.012083	3	2.80e-14	Accurate
Broyden (Identity, Good)	✅ Pass	0.013875	9	2.10e-11	Accurate
Broyden (True Jacobian, Good)	✅ Pass	0.012166	4	8.42e-11	Accurate
Broyden (Identity, Bad)	✅ Pass	0.013416	8	4.93e-11	Accurate
GaussNewton	✅ Pass	0.011959	3	2.80e-14	Accurate
GaussNewton (Krylov)	✅ Pass	0.012000	3	2.80e-14	Accurate
RobustMultiNewton	✅ Pass	0.018042	5	3.89e-15	Accurate

Key Findings:

• All solvers converged accurately to the expected concentrations within atol=1.0e-9.
• Fastest Solvers: NewtonRaphson(linesearch=HagerZhang()) and GaussNewton(linsolve=nothing) (both 0.011959 ms, 3 iterations) are the fastest, followed by NewtonRaphson(linsolve=nothing), NewtonRaphson(linesearch=BackTracking()), and GaussNewton(linsolve=KrylovJL_GMRES()) (~0.012 ms, 3 iterations).
• Most Precise: TrustRegion variants and RobustMultiNewton achieved the smallest residuals (3.89e-15), but are slower (0.016–0.018 ms).
• Broyden Methods: Less efficient (4–9 iterations) due to approximate Jacobian updates.
• Recommendation: NewtonRaphson(linesearch=HagerZhang(), linsolve=nothing) or GaussNewton(linsolve=nothing) are optimal for small, non-stiff systems like "System 1" due to their speed and low iteration count.

Changes Made:

• Updated src/Exact.jl to manually specify solver configurations, fixing invalid keyword errors (stats, autodiff) that previously caused failures.
• Modified test/runtests.jl to test each solver individually and save results to solver_comparison_system1.csv using CSV.jl and DataFrames.jl.

Issue with simulate for System 2:
The simulate function in src/NEKMC.jl is not converging accurately for "System 2," even with maxiters=1e10 and atol=1.0e-3. This suggests potential stiffness or numerical instability, possibly due to large equilibrium constants. I’d appreciate guidance on debugging, such as scaling keq_vals or logging intermediate concentrations.

Next Steps:

1. Finalize Solver: Select NewtonRaphson(linesearch=HagerZhang(), linsolve=nothing) or GaussNewton(linsolve=nothing) and raise a PR to integrate into NICE.jl.
2. Test System 2: Apply the chosen solver with scaled equilibrium constants to handle stiffness:

   rxn_system = ReactionSystem(stoich_system2, initial_concs_system2, keq_vals_system2 ./ 1e9; φ=1.0)
   hybrid_solve(rxn_system, keq_vals_system2 ./ 1e9; maxiters=100, abstol=1e-6)
   solve(rxn_system, keq_vals_system2 ./ 1e9; alg=NewtonRaphson(linesearch=HagerZhang()), linsolve=nothing)

3. Symbolic Solvers: Explore AlgebraicSolving.jl or HomotopyContinuation.jl for analytical solutions, as suggested by @PaulWAyers:

   using HomotopyContinuation
   # Define polynomial system from rxn_system and K_eqs
   # Solve using HomotopyContinuation.solve

4. Optimize Parameters: Test maxiters=500 to reduce runtime:

   solve(rxn_system, keq_vals; alg=NewtonRaphson(linesearch=HagerZhang()), linsolve=nothing, maxiters=500, abstol=1.0e-9, reltol=0.0)

5. Enhance CSV Output: Add final concentrations to solver_comparison_system1.csv for verification:

   Final_Concs = stats.success ? stats.final_concs : nothing

Please advise on the solver choice, debugging simulate for "System 2" (e.g., sharing stoich_system2, keq_vals_system2), and priorities for symbolic solvers. I’m happy to start setting up AlgebraicSolving.jl or HomotopyContinuation.jl if that’s the next focus.