---
title: Optimize-then-Discretize Framework
url: https://www.emergentmind.com/topics/optimize-then-discretize-framework
type: topic
---

# Optimize-then-Discretize Framework

Searching arXiv for the specified papers and closely related optimize-then-discretize work.
Searching arXiv: 2304.13117.
Searching arXiv: 2402.07318.
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 [2304.13117], derivation of optimization methods in Hilbert or function spaces before projection onto finite-dimensional parametrizations in scientific machine learning [2402.07318], 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 [2005.13420; 2510.04478]. 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 [2304.13117]. 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\) [2402.07318]. 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 [2005.13420]. In continuous-time linear-quadratic optimal control, Pontryagin’s Minimum Principle is first derived, and only the resulting Hamiltonian boundary-value subproblems are discretized [2510.04478].

| 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 \(r=1/\rho\) [2304.13117]. The component-wise plateau mapping is
\[
T_\rho(x)=x-(x \bmod \rho),
\]
and, to ensure that the known global optimum \(x^*\) remains feasible, the grid is shifted by its residue modulo \(\rho\):
\[
T_{x^*}(x_\rho)=x_\rho+(x^*\bmod \rho),
\]
followed by clamping to \([lb,ub]\). The resulting discretized objective is
\[
f_\rho(x)=f\bigl(T_{x^*}(T_\rho(x))\bigr).
\]

An equivalent integer view defines \(z\in\mathbb Z^n\) by
\[
z_i=\lceil x_i\cdot r\rceil,\qquad x_i=z_i/r,
\]
so that
\[
f_\rho(z)=f\bigl(T_{x^*}(z/r)\bigr).
\]
This equivalence permits continuous optimizers and integer-based solvers to evaluate exactly the same underlying real-valued function on the induced grid [2304.13117]. The source presents explicit pseudocode for evaluating any candidate \(x\in\mathbb R^n\) under plateau size \(\rho\), and for integer-based solvers it specifies that one maintains \(z\in\mathbb Z^n\), proposes mutations on \(z\), and then sets \(x=z\cdot\rho+(x^*\bmod\rho)\) before calling \(f(\cdot)\).

The landscape consequences are central. Each point \(x\) maps to a plateau of side length \(\rho\) on which \(f\) is constant. These plateaus introduce flat regions where continuous step-size adaptation may over-shrink \(\sigma\), 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 \(\rho\) grows [2304.13117]. The source also states an informal step-size stagnation theorem: if a Gaussian-mutation strategy adapts \(\sigma\) such that \(\sigma\ll\rho\), then the probability of leaving the current plateau tends to \(0\) as \(\sigma/\rho\to0\), 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 \(\rho\) can be lower-bounded by the inverse probability of a mutation exceeding \(\rho\) in any coordinate.

The empirical comparison in dimension \(n=10\), averaged over 100 runs, focuses on the canonical \((\mu_W,\lambda)\)-CMA-ES on BBOB F1 (Sphere) and F8 (Rosenbrock) [2304.13117]. On Sphere, the continuous case \((\rho\to0)\) attains median hitting time \(\approx 1\,000\) evaluations to \(f(x)-f(x^*)\le10^{-8}\) with \(100\%\) success; at \(\rho=0.1\), the median is \(\approx1\,200\) and success drops to \(\approx90\%\); at \(\rho=1.0\), about half the runs succeed quickly at \(\approx1\,500\) evaluations and half stagnate permanently, yielding overall success \(\approx50\%\). On Rosenbrock, the continuous case needs \(\sim3\,000\) evaluations with \(100\%\) success; \(\rho=0.1\) yields success \(\approx75\%\) with median \(\approx4\,000\); \(\rho=1.0\) yields success \(\approx90\%\) with median \(\approx2\,200\), because the largest plateau smooths the curved valley; and \(\rho=2.0\) introduces multiple spurious minima, reducing success to \(\approx60\%\) with median \(\approx3\,000\). Empirical CDF plots over a budget up to \(50\,000\) evaluations confirm that for intermediate \(\rho\) 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 \(\sigma\): once mutations become smaller than \(\rho\), the method remains trapped on one plateau [2304.13117]. 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 \(A\) so that the mutation distribution becomes
\[
N(m,\sigma^2 A C A^T)
\]
and is chosen to satisfy
\[
P\bigl(|\Delta x_i|\ge \rho\bigr)\ge \alpha
\]
for each axis \(i\), with \(\alpha\) typically set to \(1/(\lambda n)\). In the reported experiments, CMA-ESwM with \(\alpha=2/(\lambda n)\) restored nearly full success on all tested plateau sizes without degrading performance on the continuous problems [2304.13117]. The same source gives a practical rule of thumb for choosing \(r=1/\rho\): set \(r\) so that \(\rho\) is not smaller than the minimal meaningful step size one expects to explore, and ensure \(\sigma_{\text{initial}}\gtrsim \rho\).

## 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 \(\mathcal H\), such as \(L^2(\Omega)\), \(H^1(\Omega)\), or a product thereof [2402.07318]. The objective is a Fréchet-differentiable functional
\[
J:\mathcal H\to\mathbb R,\qquad J(u)=E(u).
\]
Its derivative \(DE(u)\) lies in the dual space \(\mathcal H^*\), and the Riesz isomorphism \(\mathcal I:\mathcal H\to\mathcal H^*\), \(u\mapsto \langle u,\cdot\rangle_{\mathcal H}\), defines the gradient
\[
\nabla E(u):=\mathcal I^{-1}[DE(u)]\in\mathcal H,
\]
with dual norm \(\|DE(u)\|_{\mathcal H^*}=\|\nabla E(u)\|_{\mathcal H}\). More generally, the framework allows a state-dependent inner product induced by a symmetric positive-definite operator \(T_u:\mathcal H\to\mathcal H^*\),
\[
[v,w]_{T_u}:=\langle T_u v,w\rangle,\qquad \|v\|_{T_u}^2=[v,v]_{T_u},
\]
which defines a Riemannian metric on \(\mathcal H\) and acts as a function-space preconditioner.

The generic continuous update is
\[
u_{k+1}=u_k+\eta_k d_k,\qquad d_k=-T_{u_k}^{-1}[DE(u_k)]\in\mathcal H.
\]
Steepest descent corresponds to \(T_u=\mathcal I\), while Newton’s method corresponds to \(T_u=D^2E(u)\) [2402.07318]. The source notes that standard line-search or trust-region machinery can be carried out directly in \(\mathcal H\), including an Armijo-type rule and a trust-region subproblem.

Discretization enters only after the continuous direction has been defined. One approximates \(\mathcal H\) by a finite-dimensional parametric manifold
\[
\mathcal M=\{u_\theta:\theta\in\mathbb R^p\}\subset\mathcal H,
\]
with smooth parametrization \(P(\theta)=u_\theta\) and tangent space
\[
T_{u_\theta}\mathcal M=\operatorname{span}\{\partial_{\theta_i}u_\theta\}_{i=1}^p.
\]
The parameter increment \(w\in\mathbb R^p\) is chosen to best fit the continuous direction \(d_k\) in the \(T_u\)-norm via
\[
w_k=\arg\min_{w\in\mathbb R^p}\frac12\|DP(\theta_k)w-d_k\|_{T_{u_{\theta_k}}}^2.
\]
The normal equations produce the Gramian
\[
G(\theta_k)_{ij}=[\partial_{\theta_i}u_{\theta_k},\partial_{\theta_j}u_{\theta_k}]_{T_{u_{\theta_k}}},
\]
and the discrete update
\[
\theta_{k+1}=\theta_k-\eta_k G(\theta_k)^\dagger \nabla_\theta L(\theta_k),
\]
where \(L(\theta)=E(u_\theta)\) [2402.07318]. In practice, one often adds a damping term \(\epsilon I\) and uses \((G+\epsilon I)^{-1}\).

The special-case identifications are one of the central contributions of the framework. For variational Monte Carlo in \(L^2(\Omega)\), \(T_u=\mathcal I\) yields a Gramian
\[
G_{ij}=\int \partial_{\theta_i}\psi_{\text{norm}}\,\partial_{\theta_j}\psi_{\text{norm}},
\]
and the discretization is exactly the standard natural gradient descent in VMC, commonly solved via K-FAC [2402.07318]. For Deep Ritz with Newton in \(H^1(\Omega)\), the tangent-space Gramian
\[
G_{ij}=\int \bigl[\nabla(\partial_{\theta_i}u_\theta)\cdot\nabla(\partial_{\theta_j}u_\theta)+3u_\theta^2\,\partial_{\theta_i}u_\theta\,\partial_{\theta_j}u_\theta\bigr]\,dx
\]
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 \(u_\theta,\lambda_\psi,\mu_\xi\) leads to the block Gram matrix
\[
G=\begin{bmatrix}
0 & A & B\\
A^\top & 0 & 0\\
B^\top & 0 & 0
\end{bmatrix},
\]
which reproduces Competitive Gradient Descent and the CPINN formulation. For PINNs with residual \(R(u,p)\) and objective \(E(u,p)=\frac12\|R\|_{L^2}^2\), choosing \(T_u=\mathcal I\circ DR(u)^\top\circ DR(u)\) 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 [2402.07318].

The main theoretical statement, Theorem 2.3, asserts that when each \(T_{u_\theta}\) is symmetric positive-definite, the discretized dynamics recover the continuous-space update within \(\mathcal M\) up to an \(O(\eta^2)\) error [2402.07318]. For non-symmetric or indefinite \(T_u\), the source gives Céa-type quasi-optimality estimates of the form
\[
\|d_k-d_{\theta_k}\|_{\mathcal H}\le C\inf_{v\in T_u\mathcal M}\|d_k-v\|_{\mathcal H},
\]
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 \(T_u\), 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 [2402.07318].

## 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 [2005.13420]. The state \(y(t)\in\mathbb R^{n_f}\) evolves under a time-dependent control \(\theta(t)\in\mathbb R^{n_p}\) according to
\[
\partial_t y(t)=\ell(\theta(t),y(t),t),\qquad y(0)=y_0,
\]
and the reduced objective is
\[
J(\theta)=\frac1S\sum_{i=1}^S\int_0^T L(y^{(i)}(t),u^{(i)}(t))\,dt+R(\theta).
\]
Introducing an adjoint variable \(z(t)\in\mathbb R^{n_f}\) and forming the Lagrangian leads, after integration by parts and stationarity, to the continuous adjoint equations
\[
-\partial_t z(t)=\bigl[\partial_y\ell(\theta(t),y(t),t)\bigr]^\top z(t)+\nabla_yL(y(t),u(t)),\qquad z(T)=\nabla_yL(y(T),u(T)),
\]
and gradient formula
\[
\frac{dJ}{d\theta}(t)=\bigl[\partial_\theta\ell(\theta(t),y(t),t)\bigr]^\top z(t).
\]
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 [2005.13420].

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 [2005.13420]. Although the Opt-Disc and Disc-Opt gradients collapse to the same continuous limit as \(h\to0\), for finite \(h\) they differ unless \(h\) is very small. A Taylor-expansion test yields
\[
E_0(h):=|J(\theta+hv)-J(\theta)|=O(h),\qquad
E_1(h):=|J(\theta+hv)-J(\theta)-h\,v^\top\nabla J|=O(h^2),
\]
and the reported numerical tests show that Disc-Opt achieves the expected \(O(h^2)\) derivative-check behavior, whereas the Opt-Disc adjoint can be severely degraded if solver tolerances or step sizes are not sufficiently small [2005.13420].

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 \(39\%\) to \(97\%\), with an average approximately \(6\times\) speedup in CNFs and approximately \(20\times\) in the toy 2D time-series problem [2005.13420]. 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 \(h\) 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 \(y(0)\mapsto y(T)\) 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 [2005.13420]. The source notes that RK4 with \(h=0.25\) 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 \(h\) small enough, to re-discretize at inference using a finer integrator or smaller \(h\), 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 [2005.13420].

## 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 [2510.04478]. The Bolza-type problem on \([0,T]\) is
\[
\min_{u(\cdot),x(\cdot)} \;
\frac12\int_0^T
\begin{pmatrix}x(t)\\u(t)\end{pmatrix}^\top
\begin{pmatrix}Q(t)&H(t)^\top\\H(t)&R(t)\end{pmatrix}
\begin{pmatrix}x(t)\\u(t)\end{pmatrix}\,dt
+\frac12 x(T)^\top Q_T x(T),
\]
subject to
\[
\dot x(t)=A(t)x(t)+B(t)u(t),\qquad x(0)=x_0.
\]
Pontryagin’s Minimum Principle introduces a costate \(\lambda(t)\) and Hamiltonian \(\mathcal H(t,x,u,\lambda)\), yielding first-order necessary and sufficient conditions: the state equation, the adjoint equation
\[
\dot\lambda^*(t)=-(Q(t)x^*(t)+H(t)^\top u^*(t))-A(t)^\top\lambda^*(t),
\qquad \lambda^*(T)=Q_T x^*(T),
\]
and the optimal control law
\[
u^*(t)=-R(t)^{-1}\bigl(H(t)x^*(t)+B(t)^\top\lambda^*(t)\bigr)
\]
[2510.04478].

The paper then introduces a modified overlapping Schwarz decomposition in time. Given a partition
\[
0=t_0<t_1<\cdots<t_m=T
\]
with backward and forward overlaps \(\tau_j^0,\tau_j^1>0\), extended subintervals \([t_j^0,t_j^1]\) are defined by
\[
t_j^0=\max\{t_{j-1}-\tau_j^0,0\},\qquad
t_j^1=\min\{t_j+\tau_j^1,T\}.
\]
On each subinterval, a local parameterized subproblem is solved with boundary data \((p_j,q_j)\in\mathbb R^{n_x}\times\mathbb R^{n_x}\), using a terminal penalty
\[
L_j(x,q)=
\begin{cases}
\frac12 x^\top Q(t_j^1)x-x^\top Q(t_j^1)q,& j<m,\\[4pt]
\frac12 x^\top Q_T x,& j=m.
\end{cases}
\]
The Schwarz iteration updates
\[
p_j^{(k+1)}=x^{(k)}(t_j^0),\qquad
q_j^{(k+1)}=x^{(k)}(t_j^1)-Q^{-1}(t_j^1)\lambda^{(k)}(t_j^1),
\]
solves the local PMP boundary-value problems in parallel, and reassembles the global iterate from the unique subdomain covering each time point [2510.04478].

The convergence analysis rests on an exponential decay of sensitivity property. Under Uniform Complete Controllability and coercivity, the shifted-system Riccati matrix \(S(t)\) satisfies uniform bounds
\[
0<c_0 I\preceq S(t)\preceq c_1 I,
\]
and the closed-loop generator
\[
Z(t)=A-BR^{-1}H-BR^{-1}B^\top S(t)
\]
obeys the exponential stability estimate
\[
\|\Phi_Z(t,s)\|\le c_Z e^{-\rho_Z(t-s)},\qquad \rho_Z>0
\]
[2510.04478]. The paper proves EDS estimates both for point perturbations and for boundary perturbations, and uses them to derive the main Schwarz contraction theorem:
\[
\omega^{(k+1)}\le \gamma e^{-\rho_Z\tau}\,\omega^{(k)},
\]
a linear contraction in the sup-norm when the minimal overlap \(\tau\) 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 [2510.04478]. 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 \(m=3\) subdomains and overlaps \(\tau\in\{1\%,5\%,10\%,20\%,30\%,60\%\}\) of subinterval length [2510.04478]. 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 \(c\,\exp(-\rho_Z\tau)\); higher-order integrators reduce subproblem error per iteration; and in the stiff regime \((\xi\gg1)\), 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 [2402.07318]. 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 [2005.13420]. In black-box search, large plateau sizes can make classical discrete optimizers such as an integer-EA or GA outperform continuous ES, whereas as \(\rho\to0\) continuous methods regain their advantage because int-EA or GA suffer from high branching factors [2304.13117].

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 [2304.13117]. 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 [2402.07318]. 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 [2510.04478].

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 \(\rho\), an integer encoding, a tangent-space Gramian, a time integrator, or a domain decomposition [2304.13117; 2402.07318; 2005.13420; 2510.04478]. 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 \(\rho\) or hybrid continuous-discrete search strategies, and extension to multi-objective and constrained black-box settings [2304.13117]. 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 [2510.04478]. 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.

Source: https://www.emergentmind.com/topics/optimize-then-discretize-framework