Optimize-then-Discretize Framework

• Optimize-then-Discretize is a framework that defines optimization in continuous or infinite-dimensional spaces before introducing any discretization, preserving the inherent problem geometry.
• It is applied in diverse contexts such as black-box optimization, scientific machine learning, and neural ODE training, enabling resolution-controlled discretization and tailored numerical methods.
• The approach highlights trade-offs between preserving continuous structure and achieving computational efficiency, influencing convergence behavior and the emergence of spurious minima or stability issues.

Searching arXiv for the specified papers and closely related optimize-then-discretize work. Searching (Thomaser et al., 2023). Searching (Müller et al., 2024). Optimize-then-Discretize denotes a family of optimization paradigms in which the optimization problem, update rule, or search dynamics are formulated first in a continuous, infinite-dimensional, or nominally real-valued setting and only afterward discretized for computation or finite-resolution evaluation. In the cited literature, the term covers at least three distinct but structurally related constructions: explicit discretization of continuous black-box domains into grid-based or integer encodings (Thomaser et al., 2023), derivation of optimization methods in Hilbert or function spaces before projection onto finite-dimensional parametrizations in scientific machine learning (Müller et al., 2024), and continuous-time optimal-control formulations in neural ODEs or linear-quadratic control whose adjoint or Hamiltonian systems are discretized only at the numerical-solver stage (Onken et al., 2020, Zhao et al., 6 Oct 2025). The shared premise is that discretization is not treated as the primary mathematical object, but as a secondary approximation layer imposed on a problem whose geometry is specified first.

1. Core idea and scope

The literature uses the phrase in more than one technical sense. In black-box optimization, the framework makes finite precision explicit by starting from a continuous problem

$\min f(x)\quad \text{over } x\in[lb,ub]^n,$

and then introducing a user-controlled resolution so that each real variable is projected onto a finite grid (Thomaser et al., 2023). In scientific machine learning, the starting point is an infinite-dimensional optimization problem on a Hilbert space $\mathcal H$, with a Fréchet-differentiable functional $J:\mathcal H\to\mathbb R$, and the discretization step consists of fitting the continuous update into a finite-dimensional tangent space associated with a parametric manifold $$\mathcal M=\{u_\theta:\theta\in\mathbb R^p\}\subset\mathcal H$$ (Müller et al., 2024). In neural ODEs, training is posed as an optimal-control problem in continuous time, and the adjoint equations are derived before any time integrator is selected (Onken et al., 2020). In continuous-time linear-quadratic optimal control, Pontryagin’s Minimum Principle is first derived, and only the resulting Hamiltonian boundary-value subproblems are discretized (Zhao et al., 6 Oct 2025).

Setting	Continuous object specified first	Discretization imposed afterward
Black-box optimization	$f(x)$ on $x\in[lb,ub]^n$	Grid with plateau size $\rho$ or integer encoding
SciML in function space	$J:\mathcal H\to\mathbb R$, metric $T_u$	Tangent-space Gramian update in parameter space
Neural ODE training	Continuous optimal-control problem and adjoint	Time integrator for state and adjoint ODEs
Linear-quadratic OCP	PMP Hamiltonian system in continuous time	Numerical integration of local two-point BVPs

This suggests that Optimize-then-Discretize is less a single algorithm than a methodological orientation. The invariant feature is that the optimization principle is defined in the continuous model before being represented on a grid, in a basis, or in a time-stepping scheme.

2. Resolution-controlled discretization of continuous black-box problems

In the black-box setting, the framework introduces discretization through a plateau size $\rho>0$, equivalently a resolution parameter $\mathcal H$0 (Thomaser et al., 2023). The component-wise plateau mapping is

$\mathcal H$1

and, to ensure that the known global optimum $\mathcal H$2 remains feasible, the grid is shifted by its residue modulo $\mathcal H$3: $\mathcal H$4 followed by clamping to $\mathcal H$5. The resulting discretized objective is

$\mathcal H$6

An equivalent integer view defines $\mathcal H$7 by

$\mathcal H$8

so that

$\mathcal H$9

This equivalence permits continuous optimizers and integer-based solvers to evaluate exactly the same underlying real-valued function on the induced grid (Thomaser et al., 2023). The source presents explicit pseudocode for evaluating any candidate $J:\mathcal H\to\mathbb R$0 under plateau size $J:\mathcal H\to\mathbb R$1, and for integer-based solvers it specifies that one maintains $J:\mathcal H\to\mathbb R$2, proposes mutations on $J:\mathcal H\to\mathbb R$3, and then sets $J:\mathcal H\to\mathbb R$4 before calling $J:\mathcal H\to\mathbb R$5.

The landscape consequences are central. Each point $J:\mathcal H\to\mathbb R$6 maps to a plateau of side length $J:\mathcal H\to\mathbb R$7 on which $J:\mathcal H\to\mathbb R$8 is constant. These plateaus introduce flat regions where continuous step-size adaptation may over-shrink $J:\mathcal H\to\mathbb R$9, preventing further progress. On non-convex or non-separable landscapes such as Rosenbrock, discretization can create spurious local minima, and one observes a transition from unimodality to multimodality as $$\mathcal M=\{u_\theta:\theta\in\mathbb R^p\}\subset\mathcal H$$0 grows (Thomaser et al., 2023). The source also states an informal step-size stagnation theorem: if a Gaussian-mutation strategy adapts $$\mathcal M=\{u_\theta:\theta\in\mathbb R^p\}\subset\mathcal H$$1 such that $$\mathcal M=\{u_\theta:\theta\in\mathbb R^p\}\subset\mathcal H$$2, then the probability of leaving the current plateau tends to $$\mathcal M=\{u_\theta:\theta\in\mathbb R^p\}\subset\mathcal H$$3 as $$\mathcal M=\{u_\theta:\theta\in\mathbb R^p\}\subset\mathcal H$$4, so the algorithm almost surely stagnates. No closed-form complexity bounds are available for general black-box functions, but the expected time to escape a plateau of radius $$\mathcal M=\{u_\theta:\theta\in\mathbb R^p\}\subset\mathcal H$$5 can be lower-bounded by the inverse probability of a mutation exceeding $$\mathcal M=\{u_\theta:\theta\in\mathbb R^p\}\subset\mathcal H$$6 in any coordinate.

The empirical comparison in dimension $$\mathcal M=\{u_\theta:\theta\in\mathbb R^p\}\subset\mathcal H$$7, averaged over 100 runs, focuses on the canonical $$\mathcal M=\{u_\theta:\theta\in\mathbb R^p\}\subset\mathcal H$$8-CMA-ES on BBOB F1 (Sphere) and F8 (Rosenbrock) (Thomaser et al., 2023). On Sphere, the continuous case $$\mathcal M=\{u_\theta:\theta\in\mathbb R^p\}\subset\mathcal H$$9 attains median hitting time $f(x)$0 evaluations to $f(x)$1 with $f(x)$2 success; at $f(x)$3, the median is $f(x)$4 and success drops to $f(x)$5; at $f(x)$6, about half the runs succeed quickly at $f(x)$7 evaluations and half stagnate permanently, yielding overall success $f(x)$8. On Rosenbrock, the continuous case needs $f(x)$9 evaluations with $x\in[lb,ub]^n$0 success; $x\in[lb,ub]^n$1 yields success $x\in[lb,ub]^n$2 with median $x\in[lb,ub]^n$3; $x\in[lb,ub]^n$4 yields success $x\in[lb,ub]^n$5 with median $x\in[lb,ub]^n$6, because the largest plateau smooths the curved valley; and $x\in[lb,ub]^n$7 introduces multiple spurious minima, reducing success to $x\in[lb,ub]^n$8 with median $x\in[lb,ub]^n$9. Empirical CDF plots over a budget up to $\rho$0 evaluations confirm that for intermediate $\rho$1 the difficulty peaks due to multimodality, whereas very coarse grids can sometimes simplify the search on highly non-convex ridges.

A key failure mode of standard CMA-ES under discretization is over-shrinkage of $\rho$2: once mutations become smaller than $\rho$3, the method remains trapped on one plateau (Thomaser et al., 2023). A second failure mode is covariance adaptation locking into a suboptimal region created by spurious local minima. The remedy discussed in the source is CMA-ES with Margin, which inserts a diagonal margin matrix $\rho$4 so that the mutation distribution becomes

$\rho$5

and is chosen to satisfy

$\rho$6

for each axis $\rho$7, with $\rho$8 typically set to $\rho$9. In the reported experiments, CMA-ESwM with $J:\mathcal H\to\mathbb R$0 restored nearly full success on all tested plateau sizes without degrading performance on the continuous problems (Thomaser et al., 2023). The same source gives a practical rule of thumb for choosing $J:\mathcal H\to\mathbb R$1: set $J:\mathcal H\to\mathbb R$2 so that $J:\mathcal H\to\mathbb R$3 is not smaller than the minimal meaningful step size one expects to explore, and ensure $J:\mathcal H\to\mathbb R$4.

3. Function-space geometry and tangent-space discretization in SciML

In scientific machine learning, Optimize-then-Discretize is formulated as an infinite-dimensional optimization procedure on a real Hilbert space $J:\mathcal H\to\mathbb R$5, such as $J:\mathcal H\to\mathbb R$6, $J:\mathcal H\to\mathbb R$7, or a product thereof (Müller et al., 2024). The objective is a Fréchet-differentiable functional

$J:\mathcal H\to\mathbb R$8

Its derivative $J:\mathcal H\to\mathbb R$9 lies in the dual space $T_u$0, and the Riesz isomorphism $T_u$1, $T_u$2, defines the gradient

$T_u$3

with dual norm $T_u$4. More generally, the framework allows a state-dependent inner product induced by a symmetric positive-definite operator $T_u$5,

$T_u$6

which defines a Riemannian metric on $T_u$7 and acts as a function-space preconditioner.

The generic continuous update is

$T_u$8

Steepest descent corresponds to $T_u$9, while Newton’s method corresponds to $\rho>0$0 (Müller et al., 2024). The source notes that standard line-search or trust-region machinery can be carried out directly in $\rho>0$1, including an Armijo-type rule and a trust-region subproblem.

Discretization enters only after the continuous direction has been defined. One approximates $\rho>0$2 by a finite-dimensional parametric manifold

$\rho>0$3

with smooth parametrization $\rho>0$4 and tangent space

$\rho>0$5

The parameter increment $\rho>0$6 is chosen to best fit the continuous direction $\rho>0$7 in the $\rho>0$8-norm via

$\rho>0$9

The normal equations produce the Gramian

$\mathcal H$00

and the discrete update

$\mathcal H$01

where $\mathcal H$02 (Müller et al., 2024). In practice, one often adds a damping term $\mathcal H$03 and uses $\mathcal H$04.

The special-case identifications are one of the central contributions of the framework. For variational Monte Carlo in $\mathcal H$05, $\mathcal H$06 yields a Gramian

$\mathcal H$07

and the discretization is exactly the standard natural gradient descent in VMC, commonly solved via K-FAC (Müller et al., 2024). For Deep Ritz with Newton in $\mathcal H$08, the tangent-space Gramian

$\mathcal H$09

coincides with the “energy natural gradient descent” of Müller–Zeinhofer and with a generalized Gauss-Newton in parameter space. For constrained PDE problems, discretization of the Lagrange-Newton system with networks $\mathcal H$10 leads to the block Gram matrix

$\mathcal H$11

which reproduces Competitive Gradient Descent and the CPINN formulation. For PINNs with residual $\mathcal H$12 and objective $\mathcal H$13, choosing $\mathcal H$14 yields a Gramian based on residual Jacobians, and with the same quadrature used to define the finite-point PINN loss, one recovers the ordinary Gauss-Newton step in parameter space (Müller et al., 2024).

The main theoretical statement, Theorem 2.3, asserts that when each $\mathcal H$15 is symmetric positive-definite, the discretized dynamics recover the continuous-space update within $\mathcal H$16 up to an $\mathcal H$17 error (Müller et al., 2024). For non-symmetric or indefinite $\mathcal H$18, the source gives Céa-type quasi-optimality estimates of the form

$\mathcal H$19

provided suitable inf-sup conditions hold, as in mixed finite-element theory. The stated advantages relative to discretize-then-optimize are discretization-invariance or mesh-independence, natural preconditioning through the choice of $\mathcal H$20, robustness against ad hoc scaling or weighting of loss terms in PINNs, and a unified view that places K-FAC, CPINNs, energy NG, and Gauss-Newton under a common infinite-dimensional origin (Müller et al., 2024).

4. Continuous adjoints, neural ODE training, and the Opt-Disc versus Disc-Opt distinction

For neural ODEs, Optimize-then-Discretize is presented as the direct continuous-time optimal-control approach (Onken et al., 2020). The state $\mathcal H$21 evolves under a time-dependent control $\mathcal H$22 according to

$\mathcal H$23

and the reduced objective is

$\mathcal H$24

Introducing an adjoint variable $\mathcal H$25 and forming the Lagrangian leads, after integration by parts and stationarity, to the continuous adjoint equations

$\mathcal H$26

and gradient formula

$\mathcal H$27

These equations define the Opt-Disc procedure: at each training iteration one solves the forward ODE, solves the backward-in-time adjoint, and forms the gradient from the continuous formula (Onken et al., 2020).

The source contrasts this with Discretize-then-Optimize, where one first chooses a numerical integrator, writes out the discrete updates, and differentiates them directly, often by automatic differentiation. Forward Euler, classical Runge–Kutta 4, and adaptive Runge–Kutta methods such as dopri5 are explicitly discussed (Onken et al., 2020). Although the Opt-Disc and Disc-Opt gradients collapse to the same continuous limit as $\mathcal H$28, for finite $\mathcal H$29 they differ unless $\mathcal H$30 is very small. A Taylor-expansion test yields

$\mathcal H$31

and the reported numerical tests show that Disc-Opt achieves the expected $\mathcal H$32 derivative-check behavior, whereas the Opt-Disc adjoint can be severely degraded if solver tolerances or step sizes are not sufficiently small (Onken et al., 2020).

The computational trade-off is explicit. Across seven benchmark problems—one time-series regression problem and six continuous normalizing flows—Disc-Opt reduced training time by $\mathcal H$33 to $\mathcal H$34, with an average approximately $\mathcal H$35 speedup in CNFs and approximately $\mathcal H$36 in the toy 2D time-series problem (Onken et al., 2020). In one BSDS300 run, Opt-Disc required approximately nine days of GPU time while Disc-Opt finished in under 14 hours. Function-evaluation counts are predetermined by $\mathcal H$37 in Disc-Opt but adaptive and generally larger in Opt-Disc. On memory, Opt-Disc can avoid full state storage by recomputation or checkpointing, though recomputation is prone to instability, while Disc-Opt relies on standard backpropagation and is described as numerically stable and straightforward to implement.

Continuous-time structure nonetheless matters for task-specific constraints. In CNFs, the mapping $\mathcal H$38 is bijective in continuous time by construction, but a naive discrete solver in Disc-Opt can violate invertibility if the step size is too coarse (Onken et al., 2020). The source notes that RK4 with $\mathcal H$39 can produce a low training loss while yielding a discrete flow that is not one-to-one, with visual holes in the density and large inverse errors. The proposed remedies are to choose $\mathcal H$40 small enough, to re-discretize at inference using a finer integrator or smaller $\mathcal H$41, or to adopt multilevel training that starts on a coarse grid and switches to a finer grid mid-training. The practical recommendations distinguish regimes where Disc-Opt is preferable—fast training, data noisiness coarser than ODE-solver tolerances, availability of state storage or checkpointing—from regimes where Opt-Disc is preferable—strict adherence to continuous-time properties or highly sensitive dynamics where small discretization errors break invertibility (Onken et al., 2020).

5. Continuous-time optimal control and overlapping Schwarz decomposition

In continuous-time linear-quadratic optimal control, Optimize-then-Discretize is used to derive and decompose the Pontryagin system before any time discretization is chosen (Zhao et al., 6 Oct 2025). The Bolza-type problem on $\mathcal H$42 is

$\mathcal H$43

subject to

$\mathcal H$44

Pontryagin’s Minimum Principle introduces a costate $\mathcal H$45 and Hamiltonian $\mathcal H$46, yielding first-order necessary and sufficient conditions: the state equation, the adjoint equation

$\mathcal H$47

and the optimal control law

$\mathcal H$48

(Zhao et al., 6 Oct 2025).

The paper then introduces a modified overlapping Schwarz decomposition in time. Given a partition

$\mathcal H$49

with backward and forward overlaps $\mathcal H$50, extended subintervals $\mathcal H$51 are defined by

$\mathcal H$52

On each subinterval, a local parameterized subproblem is solved with boundary data $\mathcal H$53, using a terminal penalty

$\mathcal H$54

The Schwarz iteration updates

$\mathcal H$55

solves the local PMP boundary-value problems in parallel, and reassembles the global iterate from the unique subdomain covering each time point (Zhao et al., 6 Oct 2025).

The convergence analysis rests on an exponential decay of sensitivity property. Under Uniform Complete Controllability and coercivity, the shifted-system Riccati matrix $\mathcal H$56 satisfies uniform bounds

$\mathcal H$57

and the closed-loop generator

$\mathcal H$58

obeys the exponential stability estimate

$\mathcal H$59

(Zhao et al., 6 Oct 2025). The paper proves EDS estimates both for point perturbations and for boundary perturbations, and uses them to derive the main Schwarz contraction theorem: $\mathcal H$60 a linear contraction in the sup-norm when the minimal overlap $\mathcal H$61 exceeds a threshold. The stated interpretation is that the convergence rate improves exponentially with overlap and is independent of the full horizon.

The contrast with discretize-then-optimize is explicit. In the latter, one first discretizes the ODE or PDE, producing a finite-dimensional nonlinear program to be solved by methods such as SQP or interior-point algorithms. In the optimize-then-discretize formulation, one derives the continuous-time first-order conditions first, decomposes the Hamiltonian boundary-value problem in time, and discretizes only during subproblem solves (Zhao et al., 6 Oct 2025). This permits explicit, implicit, symplectic, or adaptive time-stepping per subproblem, preserves continuous-time structure such as symplecticity and turnpike behavior, and yields a naturally parallel-in-time method.

The numerical experiment considers a linearization of a constrained nonlinear OCP with $\mathcal H$62 subdomains and overlaps $\mathcal H$63 of subinterval length (Zhao et al., 6 Oct 2025). Gradient descent is applied to the Hamiltonian boundary-value subproblems using Forward Euler, Backward Euler, RK45, and adaptive MATLAB ode23. All overlap sizes yield linear convergence in Schwarz iterations; larger overlap produces a faster drop; empirical contraction rates fit $\mathcal H$64; higher-order integrators reduce subproblem error per iteration; and in the stiff regime $\mathcal H$65, fixed-step FE may become unstable whereas adaptive ode23 remains robust and efficient.

6. Comparative principles, misconceptions, and design implications

A recurring misconception is that Optimize-then-Discretize admits a uniform ranking over discretize-then-optimize or over native discrete methods. The cited literature does not support such a ranking. In function-space SciML, the optimize-then-discretize viewpoint is advocated because the quality of the update is governed by the continuous geometry, leading to mesh-independence, natural preconditioning, and a direct route to line-search, trust-region globalizations, and error estimates (Müller et al., 2024). In neural ODE training, however, discretize-then-optimize can achieve similar inference performance with drastically reduced training costs, provided that numerical treatment is careful and that continuous-time constraints such as invertibility are restored by finer inference discretization or multilevel training (Onken et al., 2020). In black-box search, large plateau sizes can make classical discrete optimizers such as an integer-EA or GA outperform continuous ES, whereas as $\mathcal H$66 continuous methods regain their advantage because int-EA or GA suffer from high branching factors (Thomaser et al., 2023).

Another misconception is that discretization is merely an implementation detail. The black-box framework explicitly argues that finite precision is a fundamental problem feature that can and often should influence algorithm choice and design (Thomaser et al., 2023). The SciML position paper makes a parallel argument at the level of geometry: if optimization is formulated directly in parameter space without reference to the ambient function space, one may lose the problem-specific metric structure encoded by the Riesz map, Hessian, Fisher operator, or residual Jacobian (Müller et al., 2024). The linear-quadratic control work extends the same logic to decomposition and solver design: deriving the PMP system first permits adaptive-time integrators and overlap-based parallelism that are not naturally visible in a fixed finite-dimensional NLP formulation (Zhao et al., 6 Oct 2025).

Across the surveyed domains, a common design pattern emerges. One first identifies the continuous optimization principle—gradient descent, Newton, Gauss-Newton, natural gradient, Lagrange-Newton, or PMP-based Hamiltonian dynamics—then chooses the operator or metric that defines the appropriate geometry, and only afterward selects the discretization mechanism: a plateau size $\mathcal H$67, an integer encoding, a tangent-space Gramian, a time integrator, or a domain decomposition (Thomaser et al., 2023, Müller et al., 2024, Onken et al., 2020, Zhao et al., 6 Oct 2025). This suggests that the central question is not whether to discretize, but at what stage discretization should enter the modeling and algorithmic pipeline.

The open questions listed in the black-box work reinforce this interpretation. They include systematic landscape characterization of discretized problems, online adaptation of $\mathcal H$68 or hybrid continuous-discrete search strategies, and extension to multi-objective and constrained black-box settings (Thomaser et al., 2023). The continuous-time Schwarz paper similarly points to extensions to nonlinear and PDE-constrained OCPs, second-order subproblem solvers using second-order adjoints, and multilevel or space-time Schwarz variants (Zhao et al., 6 Oct 2025). Taken together, these directions indicate that Optimize-then-Discretize functions as a general research program for preserving continuous structure while exposing, rather than hiding, the consequences of finite representation.