ExaModels.jl: Sparse NLP & GPU Modeling

• ExaModels.jl is an open-source Julia tool for modeling nonlinear programs and automatically generating sparse derivative kernels optimized for CPU and GPU execution.
• It preserves repeated algebraic structures and sparsity patterns, compiling objectives, constraints, gradients, Jacobians, and Hessians into specialized code.
• The system acts as a bridge between high-level model formulation and low-level solver interoperability, enhancing performance in optimal control, AC OPF, and physics-constrained modeling.

ExaModels.jl is an open-source Julia-based algebraic modeling system and sparse automatic differentiation engine for large-scale nonlinear programming. In the literature, it functions as the modeling and derivative-evaluation layer in GPU-accelerated optimization workflows: it represents optimization variables, constraints, and objectives in sparse form, preserves repeated algebraic structure, and compiles objectives, constraints, gradients, Jacobians, and Hessians into specialized CPU or GPU code for downstream solvers such as MadNLP.jl and sparse linear algebra backends such as cuDSS (Montoison et al., 4 Oct 2025, Shin et al., 2024, Pacaud et al., 2024, Kolli et al., 30 Jun 2026).

1. Core characterization

ExaModels.jl is described as an algebraic modeling system and sparse automatic differentiation engine, and also as a GPU-capable modeling and derivative-evaluation layer for nonlinear programming (Shin et al., 2024, Pacaud et al., 2024). In the 2025 optimal-control workflow, it is the component that takes a discretized optimal control problem expressed in Julia and turns it into a GPU-executable sparse nonlinear program whose expensive pieces can be evaluated in parallel on NVIDIA hardware (Montoison et al., 4 Oct 2025).

Its technical profile is defined by several recurrent features in the source literature. It is built in Julia using multiple dispatch; it uses operator-overloading-based AD and pattern-based reverse-mode AD; it evaluates sparse Jacobians and Hessians directly in COO format, with optional compression to CSC; and it exploits pre-analyzed sparsity patterns together with GPU/SIMD execution of derivative kernels (Shin et al., 2024). The literature therefore presents ExaModels.jl not merely as a front-end syntax for optimization, but as a compilation-oriented modeling layer whose output is already structured for sparse second-order nonlinear optimization.

This role is significant because the cost of large-scale nonlinear programming is not limited to sparse linear algebra. Repeated evaluation of residuals, gradients, Jacobians, and Hessians is itself a dominant expense in many engineering and scientific models. ExaModels.jl is used precisely at that interface between high-level model description and low-level derivative execution (Pacaud et al., 2024).

2. Modeling abstraction and code generation

A central concept attached to ExaModels.jl is a “SIMD abstraction” of nonlinear programs. One formulation used in the AC optimal power flow literature is

$$\min_{x^\flat\le x\le x^\sharp} \sum_{l\in[L]}\sum_{i\in[I_l]} f^{(l)}(x; p_i^{(l)}) \quad \text{s.t.} \quad \sum_{n\in[N_m]}\sum_{k\in[K_n]} h^{(n)}(x; s_k^{(n)}) = 0,\quad \forall m\in[M].$$

In that description, the model is stored as repeated patterns, derivatives for each pattern are compiled once, the compiled kernels are applied across many data instances, and the work can be carried out in parallel on GPUs (Shin et al., 2024).

The optimal-control literature gives a more concrete account of how this abstraction is exposed in code. After direct transcription, a continuous-time control problem becomes a sparse nonlinear program on a time grid. With dynamics

$\dot{x}(t)=f(x(t),u(t)),$

a forward-Euler discretization yields

$X_{i+1}-X_i-h_i f(X_i,U_i)=0,\quad i=0,\dots,N-1,$

while a Bolza objective

$g(x(0),x(t_f))+\int_0^{t_f} f^0(x(t),u(t))\,dt$

becomes

$g(X_0,X_N)+\sum_{i=0}^{N-1} h_i f^0(X_i,U_i).$

Boundary and path constraints discretize as

$$b(X_0,X_N)\le 0,\qquad c(X_i,U_i)\le 0,\quad i=0,\dots,N-1.$$

The key structural observation is that the resulting NLP consists of repeated evaluations of a small set of functions applied across many grid points, and ExaModels.jl exploits exactly this repetition (Montoison et al., 4 Oct 2025).

The same source states that ExaModels.jl “uses generators in the form of for loop like statements to model the SIMD abstraction,” with the function at the core of each statement becoming a GPU kernel via KernelAbstractions.jl (Montoison et al., 4 Oct 2025). In effect, model expressions are not only stored symbolically; they are compiled into kernels that preserve the indexed structure of the original transcription. This suggests that ExaModels.jl is best understood as a structural compiler for sparse NLP families rather than only as a declarative modeling syntax.

3. Sparse derivatives and solver interoperability

ExaModels.jl’s derivative machinery is central to its use with second-order nonlinear solvers. The optimal-control workflow explicitly states that it automatically differentiates model expressions to obtain objective gradients, constraint Jacobians, and Hessians of the Lagrangian, and that this is essential for modern second-order nonlinear solvers such as IPOPT and MadNLP.jl (Montoison et al., 4 Oct 2025). In the condensed-space GPU literature, ExaModels evaluates the nonlinear objective, constraints, gradients, Jacobians, and second derivatives using thread-parallel kernels tailored to sparsity and repetition patterns (Pacaud et al., 2024).

In the solver stack most often discussed with ExaModels.jl, the division of labor is explicit. ExaModels.jl provides model compilation and derivative kernels; MadNLP.jl serves as the nonlinear programming engine or “central orchestrator” of the optimization loop; and cuDSS, through CUDSS.jl, provides sparse factorizations and triangular solves for the Newton or condensed KKT systems (Montoison et al., 4 Oct 2025). In the GPU implementations of condensed-space interior-point methods, ExaModels is used to encode the nonlinear objective and constraints and generate custom derivative kernels, while MadNLP assembles KKT matrices and right-hand sides on the GPU and dispatches the sparse solves to cuDSS (Pacaud et al., 2024).

This interface is particularly important because the condensed-space papers emphasize that classical indefinite sparse KKT factorizations require numerical pivoting, which is difficult to parallelize efficiently on GPUs. Their alternative is to solve a symmetric positive-definite condensed system with pivot-free factorizations, while ExaModels.jl supplies the model-side computations needed at each interior-point iteration (Pacaud et al., 2024). ExaModels.jl is therefore not the optimizer itself; it is the model and derivative substrate on which these GPU-oriented optimization methods operate.

4. Representative application domains

The published use cases for ExaModels.jl span optimal control, power systems, and physics-constrained generative modeling.

In large-scale sparse nonlinear optimal control, the workflow is explicitly layered as OptimalControl.jl → ExaModels.jl → MadNLP.jl → cuDSS. Continuous-time dynamics are modeled in OptimalControl.jl, discretized by direct transcription, compiled into GPU kernels with ExaModels.jl, and then solved entirely on GPU using MadNLP.jl and cuDSS (Montoison et al., 4 Oct 2025). In this setting, ExaModels.jl is the “bridge” package that converts a discretized control problem into GPU-friendly sparse NLP code.

In multi-period AC optimal power flow, ExaModels.jl appears as one of the two open-source GPU-accelerated nonlinear programming frameworks used to solve extreme-scale OPF instances on high-memory GPUs (Shin et al., 2024). The same paper introduces ExaModelsPower.jl as an open-source modeling tool for static AC OPF, multi-period AC OPF, and GPU-oriented formulations (Shin et al., 2024). A later paper states that ExaModelsPower.jl is built on ExaModels.jl, provides a high-level interface for creating GPU-compatible nonlinear AC optimal power flow models, and automatically generates all necessary callback functions for GPU solvers; it is designed for large-scale problem instances that may include multiple time periods and security constraints (Johnson et al., 14 Oct 2025). The available excerpt does not provide further internal architectural detail, but it clearly places ExaModels.jl as the underlying modeling framework.

In physics-constrained generative modeling, ExaModels.jl is used inside SNAP-FM to accelerate the repeated projection step in Physics-Constrained Flow Matching. The projection subproblem is posed as

$$u_1=\arg\min_u \tfrac12\|u-\hat u_1\|^2 \quad \text{s.t.} \quad h(u)=0,$$

and the paper argues that batching across samples and local PDE couplings induce block-sparse Jacobian and KKT systems. ExaModels.jl exposes this structure and compiles the projection NLP into sparsity-aware kernels; MadNLP.jl then solves the equality-constrained NLP, with cuDSS handling GPU sparse factorization (Kolli et al., 30 Jun 2026). In that paper, ExaModels.jl is contrasted with JuMP as a baseline: ExaModels preserves sparsity and SIMD structure in the computational graph and compiles it through ExaCore into parallel kernels (Kolli et al., 30 Jun 2026).

5. Empirical performance characteristics

Reported performance results consistently present ExaModels.jl as part of an end-to-end stack rather than as an isolated component, but several studies identify its modeling and derivative compilation as a major source of runtime reduction.

In optimal control, the reported crossover depends on problem class and hardware. For the Goddard problem, with problem size roughly $10N$, GPU execution becomes faster than CPU after $N=2000$ on A100 and after $N=5000$ on H100, with speedup about $\dot{x}(t)=f(x(t),u(t)),$0; the largest H100 run is about 2 million total size and about 15 seconds. For the Quadrotor problem, with problem size roughly $\dot{x}(t)=f(x(t),u(t)),$1, GPU execution becomes faster than CPU after $\dot{x}(t)=f(x(t),u(t)),$2 on A100 and after $\dot{x}(t)=f(x(t),u(t)),$3 on H100, with speedup about $\dot{x}(t)=f(x(t),u(t)),$4; the largest H100 run is about $\dot{x}(t)=f(x(t),u(t)),$5 size and about 13 seconds. The same study notes comparable convergence behavior on CPU and GPU, with similar iteration counts and essentially the same solutions (Montoison et al., 4 Oct 2025).

In multi-period AC OPF, experiments on an NVIDIA GH200 show that a multi-period OPF instance with more than 10 million variables can be solved up to $\dot{x}(t)=f(x(t),u(t)),$6 precision in less than 10 minutes, and benchmarked cases reach 15 million variables and 24 million constraints (Shin et al., 2024). The paper attributes feasibility to the GPU-accelerated NLP frameworks ExaModels.jl and MadNLP.jl together with the GH200’s 480 GB unified memory (Shin et al., 2024).

In condensed-space GPU interior-point methods, experiments on PGLIB and COPS show that GPUs can attain up to a tenfold speed increase compared to CPUs when solving large-scale instances (Pacaud et al., 2024). On the large PGLIB OPF instance 78484epigrids, once model evaluation is pushed to the GPU with ExaModels.jl, KKT solving dominates runtime; matrix assembly time drops from about 8 seconds on CPU to about 0.2 seconds on GPU, and cuDSS factorization of the condensed matrix is about $\dot{x}(t)=f(x(t),u(t)),$7 faster than CHOLMOD in the cited first-iteration comparison (Pacaud et al., 2024).

In SNAP-FM, ExaModels + MadNLP is reported as the fastest method overall across six PDE-constrained projection problems, and the GPU version is usually the best or among the best. For a heat-equation benchmark with simpler constraints, ExaModels + MadNLP on GPU finishes in $\dot{x}(t)=f(x(t),u(t)),$8 s versus $\dot{x}(t)=f(x(t),u(t)),$9 s for JuMP + MadNLP and $X_{i+1}-X_i-h_i f(X_i,U_i)=0,\quad i=0,\dots,N-1,$0 s for JuMP + Ipopt, while preserving comparable infeasibility. For the heat equation with an energy-evolution constraint, the GPU version runs in $X_{i+1}-X_i-h_i f(X_i,U_i)=0,\quad i=0,\dots,N-1,$1 s while JuMP-based alternatives do not converge within the runtime budget (Kolli et al., 30 Jun 2026).

6. Limitations, caveats, and interpretation

The literature also places clear limits on what ExaModels.jl should be taken to imply. GPU execution is not uniformly superior. In the optimal-control benchmarks, the GPU becomes advantageous only after sufficiently large discretizations, indicating that overheads must be amortized by problem scale (Montoison et al., 4 Oct 2025). In SNAP-FM, the CPU ExaModels + MadNLP baseline can beat the GPU version on some simpler or less demanding projections; for the reaction-diffusion problem and one Burgers regime, the CPU ExaModels backend is reported as faster than the GPU (Kolli et al., 30 Jun 2026). The condensed-space NLP paper likewise notes that on small or super-sparse cases, CPU solvers can remain competitive (Pacaud et al., 2024).

A second caveat is numerical tolerance. The power-systems and generative-modeling papers repeatedly report results at $X_{i+1}-X_i-h_i f(X_i,U_i)=0,\quad i=0,\dots,N-1,$2 precision or feasibility tolerance, and tighter tolerances may change runtime or robustness (Shin et al., 2024, Kolli et al., 30 Jun 2026, Johnson et al., 14 Oct 2025). The SNAP-FM study is explicit that ExaModels does not eliminate the intrinsic difficulty of nonlinear projection: infeasibility values remain governed by solver tolerances, and the hardest cases can retain nontrivial residuals even when the optimizer terminates within tolerance (Kolli et al., 30 Jun 2026).

A common misconception would be to regard ExaModels.jl as interchangeable with a generic algebraic front end. The papers instead attribute much of its utility to compile-time preservation of sparsity and repetition structure, direct evaluation of sparse derivatives, and the ability to generate kernels for objectives, constraints, gradients, Jacobians, and Hessians (Shin et al., 2024, Montoison et al., 4 Oct 2025, Kolli et al., 30 Jun 2026). Another misconception would be to identify it with a complete solver. The published workflows are consistent in separating ExaModels.jl from MadNLP.jl and cuDSS: ExaModels models and differentiates, MadNLP optimizes, and cuDSS performs sparse direct linear algebra (Montoison et al., 4 Oct 2025, Pacaud et al., 2024).

Taken together, the arXiv literature presents ExaModels.jl as a structural modeling and sparse AD substrate for nonlinear optimization on CPUs and GPUs. Its distinctive role is to preserve problem structure—especially repeated indexed expressions and sparse derivative patterns—so that large-scale NLP workflows in optimal control, AC OPF, and physics-constrained inference can be executed with device-resident derivative evaluation and solver-compatible sparse linear algebra (Montoison et al., 4 Oct 2025, Shin et al., 2024, Pacaud et al., 2024, Kolli et al., 30 Jun 2026).