Neural Shooting Method: Theory & Practice

• Neural shooting method is a framework for training dynamical models by splitting long rollouts into short segments with trainable boundary variables and enforced continuity.
• It employs techniques such as quadratic penalties, augmented Lagrangians, and implicit root-finding to ensure inter-segment coherence and improve optimization geometry.
• This approach enhances stability and performance in tasks like oscillatory time-series fitting, latent state-space learning, and trajectory optimization.

Neural shooting method denotes a family of learning, inference, and control procedures that adapt direct shooting and, more prominently, multiple shooting from numerical analysis and optimal control to neural models. Its unifying idea is to replace a single long rollout of a dynamical or deep model by shorter segmentwise solves with trainable boundary variables, and then to recover global coherence through continuity equations, penalties, augmented Lagrangians, probabilistic priors, or implicit root-finding. In current usage, the term covers multiple-shooting training of neural differential equations on oscillatory time series (Turan et al., 2021), latent state-space learning from noisy sensor data and raw images (Jordana et al., 2021), differentiable multiple-shooting layers as implicit neural modules (Massaroli et al., 2021), shooting-based formulations of deep learning and continuous-depth networks (Vialard et al., 2020), and learning-aided shooting methods for trajectory optimization (Cheng et al., 2020).

1. Conceptual basis and optimal-control interpretation

The conceptual core of neural shooting is the distinction between single shooting and multiple shooting. In single shooting, one fixes an initial condition, integrates or iterates the learned dynamics over the full horizon, and optimizes model parameters from the resulting end-to-end loss. In multiple shooting, the horizon is partitioned, additional segment-initial variables are introduced, and the global trajectory is reconstructed by enforcing inter-segment continuity. This changes the optimization geometry: rather than forcing one parameter update to improve the fit everywhere along a long trajectory, the method allows local trajectory corrections before continuity is fully imposed (Turan et al., 2021, Jordana et al., 2021).

The optimal-control viewpoint makes this distinction precise. For a continuous problem with dynamics $\dot{x}(t)=f(x(t),u(t),t)$ and objective

$$J[x(\cdot),u(\cdot)] \;=\; \Phi\big(x(T)\big) \;+\; \int_0^T L\big(x(t),u(t),t\big)\,dt,$$

the Hamiltonian is

$$H(x,u,\lambda,t) \;=\; L(x,u,t) \;+\; \lambda^{\top} f(x,u,t).$$

A direct shooting parametrization $u(t;\theta)$ yields the adjoint gradient

$$\frac{\partial J}{\partial \theta} \;=\; \int_0^T \Big(\frac{\partial u}{\partial \theta}(t)\Big)^{\!\top}\frac{\partial H}{\partial u}(x,u,\lambda,t)\,dt.$$

In the discrete Euler setting studied in "A Direct Shooting Method is Equivalent to an Indirect Method" (Ross, 2020), the reduced direct-shooting gradient satisfies

$$\frac{\partial E^N}{\partial U^N} = h \begin{bmatrix} \partial_u H(\lambda_0,x_0,u_0)\ \vdots \ \partial_u H(\lambda_{N-1},x_{N-1},u_{N-1}) \end{bmatrix},$$

which establishes first-order equivalence between direct shooting and the indirect Pontryagin formulation.

This control-theoretic reading also recasts ordinary backpropagation. "Neural Network Training as an Optimal Control Problem: An Augmented Lagrangian Approach" explicitly identifies standard feedforward training with the sequential approach, i.e. single shooting, where state variables have been eliminated by forward propagation; the paper then motivates multiple shooting as the simultaneous alternative for mitigating ill-conditioning in deep compositions (Evens et al., 2021). "A Shooting Formulation of Deep Learning" extends the same logic to continuous-depth networks by replacing optimization over an entire time-varying weight trajectory $\theta(t)$ with a shooting parametrization in terms of initial conditions and compatibility equations (Vialard et al., 2020).

2. Mathematical formulations across model classes

Despite shared terminology, neural shooting appears in several mathematically distinct forms. The common structure is a set of local solves plus a mechanism that ties neighboring segments together.

Formulation	Core variables	Continuity mechanism
Neural differential equations	$\theta$, shooting states $Z=\{z_k\}$	$c_k(\theta,Z)=x_k(t_{k+1};z_k,\theta)-z_{k+1}=0$
Discrete latent state-space models	$$J[x(\cdot),u(\cdot)] \;=\; \Phi\big(x(T)\big) \;+\; \int_0^T L\big(x(t),u(t),t\big)\,dt,$$0, shooting nodes $$J[x(\cdot),u(\cdot)] \;=\; \Phi\big(x(T)\big) \;+\; \int_0^T L\big(x(t),u(t),t\big)\,dt,$$1	$$J[x(\cdot),u(\cdot)] \;=\; \Phi\big(x(T)\big) \;+\; \int_0^T L\big(x(t),u(t),t\big)\,dt,$$2
Implicit multiple-shooting layers	$$J[x(\cdot),u(\cdot)] \;=\; \Phi\big(x(T)\big) \;+\; \int_0^T L\big(x(t),u(t),t\big)\,dt,$$3	$$J[x(\cdot),u(\cdot)] \;=\; \Phi\big(x(T)\big) \;+\; \int_0^T L\big(x(t),u(t),t\big)\,dt,$$4
Sparse Bayesian latent NODEs	block latents $$J[x(\cdot),u(\cdot)] \;=\; \Phi\big(x(T)\big) \;+\; \int_0^T L\big(x(t),u(t),t\big)\,dt,$$5	Gaussian continuity prior with variance $$J[x(\cdot),u(\cdot)] \;=\; \Phi\big(x(T)\big) \;+\; \int_0^T L\big(x(t),u(t),t\big)\,dt,$$6
Condensed MS-NODE	$$J[x(\cdot),u(\cdot)] \;=\; \Phi\big(x(T)\big) \;+\; \int_0^T L\big(x(t),u(t),t\big)\,dt,$$7, $$J[x(\cdot),u(\cdot)] \;=\; \Phi\big(x(T)\big) \;+\; \int_0^T L\big(x(t),u(t),t\big)\,dt,$$8, auxiliary $$J[x(\cdot),u(\cdot)] \;=\; \Phi\big(x(T)\big) \;+\; \int_0^T L\big(x(t),u(t),t\big)\,dt,$$9	linearized KKT system for $$H(x,u,\lambda,t) \;=\; L(x,u,t) \;+\; \lambda^{\top} f(x,u,t).$$0

In continuous-time neural ODE training, the segmentwise state on interval $$H(x,u,\lambda,t) \;=\; L(x,u,t) \;+\; \lambda^{\top} f(x,u,t).$$1 is defined by

$$H(x,u,\lambda,t) \;=\; L(x,u,t) \;+\; \lambda^{\top} f(x,u,t).$$2

with continuity constraints

$$H(x,u,\lambda,t) \;=\; L(x,u,t) \;+\; \lambda^{\top} f(x,u,t).$$3

This is the formulation used to fit oscillatory time series in (Turan et al., 2021). The same paper also considers more general neural differential equations of the form

$$H(x,u,\lambda,t) \;=\; L(x,u,t) \;+\; \lambda^{\top} f(x,u,t).$$4

and, for the cascading-tanks experiment,

$$H(x,u,\lambda,t) \;=\; L(x,u,t) \;+\; \lambda^{\top} f(x,u,t).$$5

In discrete deterministic state-space learning from noisy measurements, the dynamics are

$$H(x,u,\lambda,t) \;=\; L(x,u,t) \;+\; \lambda^{\top} f(x,u,t).$$6

and multiple shooting partitions a sequence of length $$H(x,u,\lambda,t) \;=\; L(x,u,t) \;+\; \lambda^{\top} f(x,u,t).$$7 into $$H(x,u,\lambda,t) \;=\; L(x,u,t) \;+\; \lambda^{\top} f(x,u,t).$$8 intervals of length $$H(x,u,\lambda,t) \;=\; L(x,u,t) \;+\; \lambda^{\top} f(x,u,t).$$9. The constrained problem is

$u(t;\theta)$0

which is then relaxed by a soft penalty (Jordana et al., 2021).

Differentiable Multiple Shooting Layers formulate the problem as an implicit layer. After partitioning an IVP

$u(t;\theta)$1

into subintervals and introducing shooting parameters $u(t;\theta)$2, the matching map is

$u(t;\theta)$3

and the layer output is defined through the exact condition

$u(t;\theta)$4

This is neither a penalty formulation nor an augmented-Lagrangian one; it is an implicit root-finding layer with structured Jacobian (Massaroli et al., 2021).

A further generalization is probabilistic. In sparse Bayesian multiple shooting for latent neural ODEs, continuity between blocks is encoded by the prior

$u(t;\theta)$5

which turns continuity control into a KL-regularized variational objective rather than a deterministic feasibility requirement (Iakovlev et al., 2022).

3. Constraint enforcement and optimization mechanisms

The dominant design choice in neural shooting concerns how continuity is enforced. The simplest approach is a quadratic penalty. For neural differential equations, the penalty objective is

$u(t;\theta)$6

The main drawback is the classical one: too small a $u(t;\theta)$7 leaves persistent shooting gaps, whereas too large a $u(t;\theta)$8 causes numerical ill-conditioning and training stalls (Turan et al., 2021).

A more robust alternative is the augmented Lagrangian,

$u(t;\theta)$9

with multiplier updates

$$\frac{\partial J}{\partial \theta} \;=\; \int_0^T \Big(\frac{\partial u}{\partial \theta}(t)\Big)^{\!\top}\frac{\partial H}{\partial u}(x,u,\lambda,t)\,dt.$$0

This is the mechanism used in oscillatory neural ODE fitting (Turan et al., 2021) and in the layerwise optimal-control formulation of DNN training, where an outer augmented-Lagrangian loop is combined with Gauss-Newton inner solves and a forward dynamic programming solver for the block-tridiagonal linear least-squares systems arising from the multiple-shooting constraints (Evens et al., 2021).

Discrete latent state-space learning from noisy images and sensors adopts a pure penalty form,

$$\frac{\partial J}{\partial \theta} \;=\; \int_0^T \Big(\frac{\partial u}{\partial \theta}(t)\Big)^{\!\top}\frac{\partial H}{\partial u}(x,u,\lambda,t)\,dt.$$1

optimized with Adam over both global parameters and trajectory-specific shooting nodes. The reported schedule initializes $$\frac{\partial J}{\partial \theta} \;=\; \int_0^T \Big(\frac{\partial u}{\partial \theta}(t)\Big)^{\!\top}\frac{\partial H}{\partial u}(x,u,\lambda,t)\,dt.$$2, increases it once after 200 epochs to $$\frac{\partial J}{\partial \theta} \;=\; \int_0^T \Big(\frac{\partial u}{\partial \theta}(t)\Big)^{\!\top}\frac{\partial H}{\partial u}(x,u,\lambda,t)\,dt.$$3–$$\frac{\partial J}{\partial \theta} \;=\; \int_0^T \Big(\frac{\partial u}{\partial \theta}(t)\Big)^{\!\top}\frac{\partial H}{\partial u}(x,u,\lambda,t)\,dt.$$4 depending on the task, reinitializes the optimizer, and reduces the learning rate by $$\frac{\partial J}{\partial \theta} \;=\; \int_0^T \Big(\frac{\partial u}{\partial \theta}(t)\Big)^{\!\top}\frac{\partial H}{\partial u}(x,u,\lambda,t)\,dt.$$5 at that point; another decay occurs at epoch 600 (Jordana et al., 2021).

More recent work has emphasized first-order-compatible exact enforcement. In condensed MS-NODE training, the stacked equality constraints $$\frac{\partial J}{\partial \theta} \;=\; \int_0^T \Big(\frac{\partial u}{\partial \theta}(t)\Big)^{\!\top}\frac{\partial H}{\partial u}(x,u,\lambda,t)\,dt.$$6 are linearized and the reduced system

$$\frac{\partial J}{\partial \theta} \;=\; \int_0^T \Big(\frac{\partial u}{\partial \theta}(t)\Big)^{\!\top}\frac{\partial H}{\partial u}(x,u,\lambda,t)\,dt.$$7

is solved for an auxiliary increment $$\frac{\partial J}{\partial \theta} \;=\; \int_0^T \Big(\frac{\partial u}{\partial \theta}(t)\Big)^{\!\top}\frac{\partial H}{\partial u}(x,u,\lambda,t)\,dt.$$8, after which $$\frac{\partial J}{\partial \theta} \;=\; \int_0^T \Big(\frac{\partial u}{\partial \theta}(t)\Big)^{\!\top}\frac{\partial H}{\partial u}(x,u,\lambda,t)\,dt.$$9 and $$\frac{\partial E^N}{\partial U^N} = h \begin{bmatrix} \partial_u H(\lambda_0,x_0,u_0)\ \vdots \ \partial_u H(\lambda_{N-1},x_{N-1},u_{N-1}) \end{bmatrix},$$0 are recovered and passed to Adam. The stated purpose is to incorporate shooting equality constraints while training a multiple-shooting neural ordinary differential equation using first-order optimization methods such as Adam, without tuning penalty weights (Prabhu et al., 31 May 2025).

Gradient computation likewise varies by formulation. Continuous-time methods may differentiate through the ODE solver segmentwise, use forward sensitivities, or integrate adjoint equations. Whole-brain dynamic causal modeling combines multiple shooting with an adjoint implementation called Adaptive Checkpoint Adjoint, chosen to preserve numerical consistency between forward and backward passes under noisy observations (Zhuang et al., 2021). Differentiable Multiple Shooting Layers instead exploit the nilpotent block-lower structure of $$\frac{\partial E^N}{\partial U^N} = h \begin{bmatrix} \partial_u H(\lambda_0,x_0,u_0)\ \vdots \ \partial_u H(\lambda_{N-1},x_{N-1},u_{N-1}) \end{bmatrix},$$1 to compute implicit gradients efficiently, and also provide a forward-sensitivity formulation in which $$\frac{\partial E^N}{\partial U^N} = h \begin{bmatrix} \partial_u H(\lambda_0,x_0,u_0)\ \vdots \ \partial_u H(\lambda_{N-1},x_{N-1},u_{N-1}) \end{bmatrix},$$2 satisfies

$$\frac{\partial E^N}{\partial U^N} = h \begin{bmatrix} \partial_u H(\lambda_0,x_0,u_0)\ \vdots \ \partial_u H(\lambda_{N-1},x_{N-1},u_{N-1}) \end{bmatrix},$$3

allowing segment sensitivities to be obtained in one solver call (Massaroli et al., 2021).

4. Empirical domains and representative results

A central empirical driver for neural shooting has been failure of single shooting on oscillatory or long-horizon dynamics. In neural ODE time-series fitting, standard single-shooting training can converge to a flattened, low-frequency trajectory that cuts through oscillations rather than tracking them; the paper attributes this to spectral bias, long-horizon sensitivity, and poor local minima. Multiple shooting with penalty or augmented Lagrangian resolves this behavior on both a synthetic spiral and a real cascading-tanks dataset. In the spiral example, single shooting with Nadam at learning rate $$\frac{\partial E^N}{\partial U^N} = h \begin{bmatrix} \partial_u H(\lambda_0,x_0,u_0)\ \vdots \ \partial_u H(\lambda_{N-1},x_{N-1},u_{N-1}) \end{bmatrix},$$4 fails to capture oscillations, whereas multiple shooting with $$\frac{\partial E^N}{\partial U^N} = h \begin{bmatrix} \partial_u H(\lambda_0,x_0,u_0)\ \vdots \ \partial_u H(\lambda_{N-1},x_{N-1},u_{N-1}) \end{bmatrix},$$5 intervals and an augmented-Lagrangian method with LBFGS successfully fits the data and generalizes to $$\frac{\partial E^N}{\partial U^N} = h \begin{bmatrix} \partial_u H(\lambda_0,x_0,u_0)\ \vdots \ \partial_u H(\lambda_{N-1},x_{N-1},u_{N-1}) \end{bmatrix},$$6 despite training only to $$\frac{\partial E^N}{\partial U^N} = h \begin{bmatrix} \partial_u H(\lambda_0,x_0,u_0)\ \vdots \ \partial_u H(\lambda_{N-1},x_{N-1},u_{N-1}) \end{bmatrix},$$7. In the cascading-tanks experiment, multiple shooting with Bayesian-optimized $$\frac{\partial E^N}{\partial U^N} = h \begin{bmatrix} \partial_u H(\lambda_0,x_0,u_0)\ \vdots \ \partial_u H(\lambda_{N-1},x_{N-1},u_{N-1}) \end{bmatrix},$$8, $$\frac{\partial E^N}{\partial U^N} = h \begin{bmatrix} \partial_u H(\lambda_0,x_0,u_0)\ \vdots \ \partial_u H(\lambda_{N-1},x_{N-1},u_{N-1}) \end{bmatrix},$$9, and $\theta(t)$0 achieves average square-root error $\theta(t)$1 on training, $\theta(t)$2 on validation, and $\theta(t)$3 on test (Turan et al., 2021).

The same stability advantage appears under indirect and noisy observations. In the discrete multiple-shooting state-space framework, the method is reported as robust to noisy measurements and able to handle chaotic systems. On pendulum images with Gaussian pixel noise $\theta(t)$4, test performance averaged over 1000 trajectories reached MSE $\theta(t)$5 for SSM-FC and $\theta(t)$6 for SSM-LL, with BCE $\theta(t)$7 and $\theta(t)$8 respectively. On Moving MNIST, the method trained stably on $\theta(t)$9 million images and $\theta$0 shooting nodes, obtaining BCE $\theta$1 for the reported SSM-LL configuration. On partially observed noisy Lorenz trajectories of length $\theta$2, ten-step prediction MSE for SSM-RNN improves from $\theta$3 at $\theta$4 to $\theta$5 at $\theta$6, and the learned dynamics reproduce the Lorenz strange attractor under long iteration of $\theta$7 (Jordana et al., 2021).

In computational neuroscience, the Multiple-Shooting Adjoint method extends the paradigm to whole-brain dynamic causal modeling. On large-scale linear systems with $\theta$8 compartments, the reported mean squared error $\theta$9 between estimated and ground-truth parameters is $Z=\{z_k\}$0, $Z=\{z_k\}$1, $Z=\{z_k\}$2, and $Z=\{z_k\}$3 for MSA, whereas EM is reported as $Z=\{z_k\}$4 and $Z=\{z_k\}$5 for 10 and 20 nodes and out-of-memory for 50 and 100 nodes. The same framework is applied to fMRI from 82 children with ASD and 48 controls, where time-resolved effective-connectome features estimated by MSA yield better accuracy, F1, and ROC-AUC than functional connectome features for both Random Forest and InvNet classifiers (Zhuang et al., 2021).

Neural shooting has also been used as a computational primitive rather than merely a training trick. Multiple Shooting Layers are presented as drop-in replacements for Neural ODEs and Neural CDEs. In latent sequence generation they are reported to require 60% less NFEs during both training and inference while maintaining equivalent sample quality to Latent Neural ODEs; in Neural CDE training on PhysioNet 2019 sepsis, total wall-clock time is reported as one order of magnitude smaller while preserving AUROC trends; and for Timoshenko beam boundary control they are reported as three times faster than the sequential rk4 with same accuracy (step size) (Massaroli et al., 2021).

In trajectory optimization, Neural-iLQR uses shooting over the unconstrained control space but replaces an inaccurate model with a learned local neural model. In numerical simulation, conventional model-based iLQR achieves objective values $Z=\{z_k\}$6 for vehicle tracking and $Z=\{z_k\}$7 for cartpole, while Neural-iLQR achieves $Z=\{z_k\}$8 and $Z=\{z_k\}$9. In MuJoCo cartpole with model inaccuracy, iLQR fails at 40% and 60% mismatch, whereas Neural-iLQR remains successful with objective values $c_k(\theta,Z)=x_k(t_{k+1};z_k,\theta)-z_{k+1}=0$0 and $c_k(\theta,Z)=x_k(t_{k+1};z_k,\theta)-z_{k+1}=0$1 (reported as Obj.Val $c_k(\theta,Z)=x_k(t_{k+1};z_k,\theta)-z_{k+1}=0$2) (Cheng et al., 2020).

5. Relation to adjacent methods and recurrent misconceptions

A recurrent misconception is that neural shooting is merely heuristic trajectory splitting. Several formulations are explicitly principled. The equivalence result for direct and indirect shooting shows that the direct reduced gradient and Hamiltonian stationarity encode the same first-order conditions (Ross, 2020). Sparse Bayesian multiple shooting replaces ad hoc stitching by a well-defined probabilistic continuity prior and an ELBO with closed-form Gaussian KL terms (Iakovlev et al., 2022). The augmented-Lagrangian and condensed-KKT formulations likewise treat continuity as a constrained optimization problem rather than an informal regularizer (Evens et al., 2021, Prabhu et al., 31 May 2025).

A second misconception is that neural shooting is synonymous with continuous-time neural ODE training. The literature is broader. The discrete state-space model of (Jordana et al., 2021) contains no numerical ODE solve; per-segment “integration” is $c_k(\theta,Z)=x_k(t_{k+1};z_k,\theta)-z_{k+1}=0$3 repeated applications of $c_k(\theta,Z)=x_k(t_{k+1};z_k,\theta)-z_{k+1}=0$4. Multiple shooting has also been used for feedforward neural networks interpreted as finite-horizon optimal-control systems (Evens et al., 2021), for implicit layers (Massaroli et al., 2021), and for shooting-based parametrization of continuous-depth or even discrete-time networks via particle ensembles (Vialard et al., 2020).

A third misconception is that continuity must always be enforced by large quadratic penalties. The literature contains at least five distinct mechanisms: quadratic penalties and augmented Lagrangians for neural differential equations (Turan et al., 2021); soft penalties in discrete latent state-space models (Jordana et al., 2021); exact root-finding in implicit multiple-shooting layers (Massaroli et al., 2021); Gaussian continuity priors in Bayesian latent neural ODEs (Iakovlev et al., 2022); and condensed linearized KKT solves designed for first-order optimizers such as Adam (Prabhu et al., 31 May 2025).

Neural shooting is also closely related to, but not identical with, collocation and two-stage smoothing. The oscillatory NODE paper contrasts multiple shooting with two-stage approaches that first fit a spline and estimate derivatives, noting that such methods require all states to be measured and degrade with noise or sparsity (Turan et al., 2021). The same paper observes that collocation enforces both dynamics and continuity algebraically at grid points, while multiple shooting solves IVPs on segments and enforces only segment matching; this often yields better conditioning and simpler parallelization for oscillatory dynamics (Turan et al., 2021). In latent neural ODE training, the Bayesian multiple-shooting paper argues that heuristic sub-trajectory training or progressive lengthening lacks principled continuity control and underperforms its blockwise probabilistic formulation, especially on irregular grids (Iakovlev et al., 2022).

6. Limitations, tuning parameters, and current research directions

The main practical degrees of freedom are segment length, number of segments, initialization of shooting variables, and continuity strength. Across formulations, shorter segments improve stability by reducing gradient path length, integration drift, or repeated composition of an untrained transition model, but increase the number of auxiliary variables and the burden of satisfying continuity. The oscillatory neural ODE study uses $c_k(\theta,Z)=x_k(t_{k+1};z_k,\theta)-z_{k+1}=0$5 intervals on $c_k(\theta,Z)=x_k(t_{k+1};z_k,\theta)-z_{k+1}=0$6 and explicitly notes interval-length effects as future work (Turan et al., 2021). The discrete latent state-space model recommends segment lengths $c_k(\theta,Z)=x_k(t_{k+1};z_k,\theta)-z_{k+1}=0$7 as a stability–efficiency compromise (Jordana et al., 2021). Sparse Bayesian multiple shooting reports a “sweet spot” for block size and shows that overly strong continuity, i.e. very small $c_k(\theta,Z)=x_k(t_{k+1};z_k,\theta)-z_{k+1}=0$8, can cause collapse on complex data (Iakovlev et al., 2022).

Another persistent issue is the trade-off between optimization robustness and computational overhead. Multiple shooting introduces auxiliary variables such as $c_k(\theta,Z)=x_k(t_{k+1};z_k,\theta)-z_{k+1}=0$9, $$J[x(\cdot),u(\cdot)] \;=\; \Phi\big(x(T)\big) \;+\; \int_0^T L\big(x(t),u(t),t\big)\,dt,$$00, multipliers, or block latents, and each iteration typically solves many local IVPs. This cost is often acceptable because segments are independent and highly parallelizable, but it can still be memory-intensive. The discrete state-space model notes that storing many shooting nodes can be expensive and suggests streaming datasets and on-the-fly node management as future improvements (Jordana et al., 2021). Multiple Shooting Layers achieve large speedups in full-batch settings through one-step fixed-point tracking, but the paper explicitly states that mini-batch SGD breaks the tracking assumption (Massaroli et al., 2021).

Noise, partial observability, and identifiability remain central challenges. The DCM work emphasizes that deconvolution, fixed hemodynamic response functions, and high-dimensional effective-connectivity estimation may require structural priors or regularization (Zhuang et al., 2021). The latent state-space formulation separates dynamics learning from deployment-time inference and therefore requires either probabilistic filtering such as UKF or a learned encoder such as CNN-LSTM to infer latent states from observations (Jordana et al., 2021). Bayesian multiple shooting addresses part of this issue with amortized inference through a transformer-based recognition network with temporal attention and relative positional encoding, but its performance remains sensitive to encoder design and posterior variance flooring (Iakovlev et al., 2022).

Current research directions therefore cluster around better constraint handling, better inference, and better scalability. Proposed directions include augmented Lagrangian or SQP alternatives to manually tuned penalties in latent state-space learning (Jordana et al., 2021), adaptive segmentation and principled continuity schedules (Jordana et al., 2021, Prabhu et al., 31 May 2025), stronger first-order-compatible equality handling through condensing (Prabhu et al., 31 May 2025), joint training of inference modules with shooting-based latent dynamics (Jordana et al., 2021), and incorporation of controls, heteroscedastic or non-Gaussian observation models, Hamiltonian or symmetry structure, and sensor-fusion constraints (Jordana et al., 2021, Iakovlev et al., 2022). A broader implication is that neural shooting has evolved from a stabilization technique for long-horizon fitting into a general design pattern for differentiable constrained dynamical learning, spanning parameter estimation, latent representation learning, optimal control, and implicit neural architectures.