从天体物理时间序列研究历史看Continuous Delayed-Memory Stochastic Gradient Descent和Continuous-Time Reinforcement Learning

# Continuous Delayed-Memory Stochastic Gradient Descent and Continuous-Time Reinforcement Learning from History of Astrophysical Time Series Studies
Source: [https://arxiv.org/html/2609.20906](https://arxiv.org/html/2609.20906)
Debartha Paul Juncheng Yi debartha@iastate\.edu ccyi3@iastate\.edu May 8th 2026

###### Abstract

Quasars are luminous objects in the universe that exhibit stochastic brightness variations encoding information about the supermassive black holes powering them, and modeling these variations from ground\-based survey data time series, known as light curves, is a statistical challenge\[[7](https://arxiv.org/html/2609.20906#bib.bib18)\]\. This paper reviews how stochastic differential equations \(SDEs\) have been adapted with neural network parameterizations to overcome this challenge in history:\[[7](https://arxiv.org/html/2609.20906#bib.bib18)\]\[[6](https://arxiv.org/html/2609.20906#bib.bib4)\],\[[17](https://arxiv.org/html/2609.20906#bib.bib15)\]\. Following\[[17](https://arxiv.org/html/2609.20906#bib.bib15)\], we create the Continuous\-Delayed\-Memory Stochastic Gradient Descent which depend on the past state of the discrete iteration process\. We performed the simulation on some 2\-dimensional landscape and observed some wider\-exploration and more precise convergent behavior compared to Vanilla SGD by adjusting hyperparameters\. Besides, we proposed a reinforcement learning structure with continuous time policy gradients for exploratory policies without solving HJB PDE, and we show that its optimality conditions recover the Gibbs policy of\[[19](https://arxiv.org/html/2609.20906#bib.bib16)\]\.

Keywords:Stochastic process, Stochastic gradient descent, Continuous\-Delayed\-Memory Stochastic Gradient Descent, Stochastic Delay Differential Equation, Reinforcement Learning, Adjoint method

## 1Introduction of Pre\-Neural SDE Methods for Astrophysical Quasar Analysis

Three mathematical objects sit at the heart of this project:

1. 1\.the*Ornstein–Uhlenbeck SDE*, which has modeled quasar optical variability since\[[10](https://arxiv.org/html/2609.20906#bib.bib8)\]
2. 2\.*stochastic gradient descent*, the workhorse optimization algorithm of modern machine learning\[[1](https://arxiv.org/html/2609.20906#bib.bib1)\]
3. 3\.the*stochastic adjoint method*\(a backward SDE\) that makes training Neural SDEs practical\[[13](https://arxiv.org/html/2609.20906#bib.bib11)\]

At first glance, these live in different worlds\. The OU process is a physical model; SGD is a numerical algorithm; the adjoint is a tool for automatic differentiation\. The central observation that structures this project is that*all three are instances of the same idea*: a state \(or parameter\) vector evolving in time under a deterministic drift plus a noise term, with the noise either injected physically \(OU\), introduced by mini\-batch sampling \(SGD\), or inherited from the forward Brownian path \(the adjoint\)\.

Our literature review started with the modeling of quasar light curves from irregularly sampled photometric surveys\. Quasars are powered by accretion onto supermassive black holes and exhibit stochastic brightness variations whose statistics encode information about the central engine\. The result of\[[10](https://arxiv.org/html/2609.20906#bib.bib8)\]was that quasar optical variability is well described by a Damped Random Walk \(DRW\), mathematically, the OU SDE

\(1\.1\)d​X​\(t\)=−1τ​\(X⁡\(t\)−μ\)​d​t\+σ​d​W​\(t\)\\mathrm\{d\}X\(t\)\\;=\\;\-\\frac\{1\}\{\\tau\}\\bigl\(X\(t\)\-\\mu\\bigr\)\\,\\mathrm\{d\}t\\;\+\\;\\sigma\\,\\mathrm\{d\}W\(t\)
whereτ\\tauis a damping timescale,μ\\muis the long\-run mean magnitude,σ\\sigmais a short\-term volatility, andW⁡\(t\)W\(t\)is a standard Wiener process\. The DRW is Gaussian, Markov, and stationary, and it admits the exact transition density

\(1\.2\)X\(t\+Δt\)\|X\(t\)∼𝒩\(μ\+\(X\(t\)−μ\)e−Δt/τ,σ2​τ2\(1−e−2Δt/τ\)\)X\(t\+\\Delta t\)\\,\\big\|\\,X\(t\)\\;\\sim\\;\\mathcal\{N\}\\\!\\left\(\\mu\+\\bigl\(X\(t\)\-\\mu\\bigr\)e^\{\-\\Delta t/\\tau\},\\;\\frac\{\\sigma^\{2\}\\tau\}\{2\}\\bigl\(1\-e^\{\-2\\Delta t/\\tau\}\\bigr\)\\right\)
So the likelihood of an irregularly sampled light curve can be written exactly without any interpolation\. This characteristic made \([1\.1](https://arxiv.org/html/2609.20906#S1.E1)\) the de\-facto standard for time\-domain astrophysics for a decade\.

However, the OU/DRW model has several well\-documented limitations\. The following list documents some of those limitations, some extensions to other modelling methods and why they fail in this problem\.

PSD slope mismatchThe OU process has a Lorentzian power spectral densityP⁡\(f\)∝1/\(f02\+f2\)P\(f\)\\propto 1/\(f\_\{0\}^\{2\}\+f^\{2\}\), giving a spectral slope of exactly−2\-2at high frequencies\.\[[16](https://arxiv.org/html/2609.20906#bib.bib14)\]analysed four AGN observed at 30\-minute cadence by the*Kepler*space telescope and found PSD slopes ranging from−2\.6\-2\.6to−3\.3\-3\.3, significantly steeper than predicted by DRW\.\[[9](https://arxiv.org/html/2609.20906#bib.bib7)\]extended the analysis to 20*Kepler*AGN and found fewer than half consistent with DRW\. Ground\-based surveys with their sparser sampling had masked this mismatch\.

Linearity and single\-band limitationBy Doob’s theorem, the OU process is the unique process that is simultaneously Gaussian, Markov, and stationary\. Analytic tractability thus comes at the cost of linear dynamics and a single output\. Real quasar variability involves nonlinear physical processes \(accretion disk instabilities, corona–disk interactions\) and multi\-band correlations with inter\-band time lags that a single\-output OU process cannot capture\.

CARMA models\[[11](https://arxiv.org/html/2609.20906#bib.bib9)\]introduced continuous\-time autoregressive moving average \(CARMA\) models, which generalize DRW to higher\-order linear SDEs driven by a common Brownian motion\. CARMA\(2,1\)\(2,1\)\(the damped harmonic oscillator\) fits many AGN light curves better than DRW, but the family remains linear and parametric: its PSD is a rational function \(a sum of Lorentzians\)\.

Discrete\-time deep learning modelsRecurrent networks \(RNNs, LSTMs, GRUs\) operate on discrete, regularly spaced time steps\. Handling irregular sampling requires interpolation, binning, or imputation, all of which can introduce artifacts\. More fundamentally, these models are deterministic mappings from input sequences to outputs; they are not generative models for the underlying process and do not quantify uncertainty in the latent dynamics\.

Neural ODEs: continuous but deterministicNeural ODEs\[[4](https://arxiv.org/html/2609.20906#bib.bib3)\]model the latent state as the solution ofd​z/d​t=fθ​\(z⁡\(t\),t\)\\mathrm\{d\}z/\\mathrm\{d\}t=f\_\{\\theta\}\(z\(t\),t\)for a neural networkfθf\_\{\\theta\}, integrated by a numerical solver that can be queried at arbitrary times\. Neural ODEs are deterministic: given an initial condition, the trajectory is fully determined\. It makes them unsuitable for systems where different realizations from the same initial condition produce different trajectories—precisely the situation for quasar variability, where the stochasticity is physical \(turbulence in the accretion flow\) rather than merely observational\.

These models’ limitations collectively motivate Neural SDEs, which combine the continuous\-time, irregular\-sampling\-compatible framework of Neural ODEs with the stochastic dynamics that are physically appropriate for quasar variability, while using neural networks to move beyond linear parametric constraints\.

With the Rubin Observatory’s Legacy Survey of Space and Time \(LSST\), which delivers light curves for roughly10710^\{7\}quasars in six photometric bands, the DRW is no longer adequate: it is linear, single\-band, and has a power spectral slope that is systematically wrong at short timescales\[[16](https://arxiv.org/html/2609.20906#bib.bib14),[9](https://arxiv.org/html/2609.20906#bib.bib7)\]\. Neural SDEs replace the fixed parametric drift and diffusion of \([1\.1](https://arxiv.org/html/2609.20906#S1.E1)\) with neural networks while retaining the continuous\-time stochastic framework\.

OrganizationSections[2](https://arxiv.org/html/2609.20906#S2)reviews the astrophysical applications of latent neural SDE from\[[6](https://arxiv.org/html/2609.20906#bib.bib4)\], and Sections[3](https://arxiv.org/html/2609.20906#S3)reviews\[[17](https://arxiv.org/html/2609.20906#bib.bib15)\]and illustrate how it progresses from\[[6](https://arxiv.org/html/2609.20906#bib.bib4)\]by considering past state in stochastic process\. Section[4](https://arxiv.org/html/2609.20906#S4)derives Continuous\-Delayed\-Memory SGD and performs its simulations in 2D parameter space\. Section[5](https://arxiv.org/html/2609.20906#S5)develops the continuous\-time RL story, culminating in the Exploratory Backward Stratonovich SDE\. To conclude, Section[6](https://arxiv.org/html/2609.20906#S6)enumerates several compelling directions for future research\.

## 2Latent SDEs for Quasar Light Curves: Fagin et al\. \(2024\)

\[[6](https://arxiv.org/html/2609.20906#bib.bib4)\]applied the latent SDE framework to astrophysical time series for the first time\. Their model simultaneously \(i\) reconstructs multi\-band quasar light curves across seasonal gaps and \(ii\) infers physical properties of the accreting black hole\.

### 2\.1Problem Formulation

For a single quasar, the data consist of irregularly sampled magnitudes inB=6B=6bands \(u,g,r,i,z,yu,g,r,i,z,y\) over a 10\-year LSST baseline:\{\(ti,xti\(b\),σti\(b\),mti\(b\)\)\}\\\{\(t\_\{i\},x\_\{t\_\{i\}\}^\{\(b\)\},\\sigma\_\{t\_\{i\}\}^\{\(b\)\},m\_\{t\_\{i\}\}^\{\(b\)\}\)\\\}, wherexxis the measured magnitude,σ\\sigmais the photometric error, andm\(b\)∈\{0,1\}m^\{\(b\)\}\\in\\\{0,1\\\}is a band\-specific observation mask\. Spacings between successivetit\_\{i\}are irregular, with seasonal gaps of∼6\\sim 6months\.

The physical parameters to infer are the black hole masslog10⁡\(MBH/M⊙\)∈\[7,10\]\\log\_\{10\}\(M\_\{\\rm BH\}/M\_\{\\odot\}\)\\in\[7,10\], the disk inclinationii, the temperature profile exponentβ\\beta\(whereT⁡\(r\)∝r−βT\(r\)\\propto r^\{\-\\beta\}; standard thin\-disk theory predictsβ=3/4\\beta=3/4\), and the DRW parametersτ\\tauandSF∞\\mathrm\{SF\}\_\{\\infty\}\.

### 2\.2Model Architecture

The model has three main components \(see[1](https://arxiv.org/html/2609.20906#S2.F1)\), and generates about9×1059\\times 10^\{5\}total trainable parameters\.

![Refer to caption](https://arxiv.org/html/2609.20906v1/model.png)Figure 1:Model used for training the latent SDE by FaginEncoderA GRU\-D network\[[3](https://arxiv.org/html/2609.20906#bib.bib2)\]processes 12 features per time step \(6 magnitudes \+ 6 errors\)*backward*in time, producing a 64\-dimensional context vectorccthat conditions the posterior drifthξ​\(⋅,⋅,c\)h\_\{\\xi\}\(\\cdot,\\cdot;c\)\. GRU\-D handles missing values through a learned decay mechanism: when a band is unobserved, its hidden state decays exponentially toward a learnable mean, with the decay rate itself learned from data\. The backward encoding is important: the posterior drift at timettshould incorporate information from*future*observations \(this is inference, not causal prediction\)\.

Neural SDE decoderThe latent statez⁡\(t\)∈ℝ8z\(t\)\\in\\mathbb\{R\}^\{8\}evolves according to the prior and posterior SDEs as defined in\[[13](https://arxiv.org/html/2609.20906#bib.bib11)\], integrated forward with Euler\-Maruyama viatorchsde\. The step size is set to twice the minimum observation bin width\.

Parameter estimation headAn MLP takes the final encoder state and predicts Gaussian posteriors \(mean and variance\) over the physical parameters, trained with a supervised negative log\-likelihood loss alongside the ELBO\.

### 2\.3Training Data and Simulation

The training set consists of10510^\{5\}simulated 10\-year LSST quasar light curves with known ground\-truth physical parameters\. The simulations use a DRW driving signal mapped through general\-relativistic accretion disk transfer functions to six\-band UV/optical light curves, then sampled with realistic LSST cadences fromrubin\_simand noiseσtot2=σsys2\+σrand2\\sigma\_\{\\rm tot\}^\{2\}=\\sigma\_\{\\rm sys\}^\{2\}\+\\sigma\_\{\\rm rand\}^\{2\}withσsys=0\.005\\sigma\_\{\\rm sys\}=0\.005mag\.

### 2\.4Training Objective

The total loss is a weighted sum:

\(2\.1\)ℒ=ℒNLL\+λctx​ℒctx\+βk​ℒKL\+λparam​ℒparam\\mathcal\{L\}\\;=\\;\\mathcal\{L\}\_\{\\rm NLL\}\\;\+\\;\\lambda\_\{\\rm ctx\}\\,\\mathcal\{L\}\_\{\\rm ctx\}\\;\+\\;\\beta\_\{k\}\\,\\mathcal\{L\}\_\{\\rm KL\}\\;\+\\;\\lambda\_\{\\rm param\}\\,\\mathcal\{L\}\_\{\\rm param\}
withℒNLL\\mathcal\{L\}\_\{\\rm NLL\}the Gaussian negative log\-likelihood of the reconstructed light curve,ℒctx\\mathcal\{L\}\_\{\\rm ctx\}a weighted MSE at observed context points,ℒKL\\mathcal\{L\}\_\{\\rm KL\}the path\-space KL divergence, used in\[[13](https://arxiv.org/html/2609.20906#bib.bib11)\], plus the initial\-state KL \(with a cyclically annealed weightβk\\beta\_\{k\}ramping from00to11to prevent posterior collapse\), andℒparam\\mathcal\{L\}\_\{\\rm param\}a supervised NLL for the physical parameters\.

The important point is thatℒ\\mathcal\{L\}has*four*sources of stochasticity: the mini\-batch, the Brownian realization of the posterior SDE, the decoder sampling, and the parameter\-head noise\. The gradient estimator is unbiased but high\-variance, and is the main reason Fagin et al\. use Adam with a small learning rate and extensive gradient clipping\.

### 2\.5Results and Follow\-up

The latent SDE is compared to multi\-output Gaussian process regression \(GPR\)\. GPR is the standard baseline for single\-object light\-curve interpolation but requires per\-object fitting \(∼\\simminutes per object\) and assumes a fixed kernel family\. Fagin et al\. show that the latent SDE \(i\) produces better reconstructions in seasonal gaps \(where GPR reverts to its prior mean\), \(ii\) processes millions of light curves in seconds after training, and \(iii\) achieves good recovery of black hole mass and disk parameters directly from photometry, a task that would otherwise require spectroscopy\.

A follow\-up paper\[[5](https://arxiv.org/html/2609.20906#bib.bib5)\]makes the reconstruction and parameter inference*physically self\-consistent*by embedding an auto\-differentiable simulator of the accretion disk directly into the computational graph: a latent SDE generates the driving X\-ray variability, predicted disk parameters, determine the transfer functions, and the UV/optical light curves are obtained by convolution—all differentiable end\-to\-end, and all trained by SGD on the composite loss\.

## 3Neural SDDEs for Astronomical Time Series: Oh et al\. \(2025\)

### 3\.1Limitation of latent SDE

The Latent SDE framework is Markovian: the evolution ofz⁡\(t\)z\(t\)depends only on the statez⁡\(t\)z\(t\)at timett\. However, many astrophysical systems naturally violate this assumption\. In quasar accretion disk reprocessing problem, a X\-ray corona that itself produces variability illuminates the accretion disk, and the disk re\-emits UV or optical light with wavelength\-dependent time delaysτλ∝λ4/3\\tau\_\{\\lambda\}\\propto\\lambda^\{4/3\}\(from the thin\-disk temperature profile\)\. The optical variability at timettis a reprocessed signal of the X\-ray flux at timet−τλt\-\\tau\_\{\\lambda\}\. Consequently, the system exhibits historical dependencies that violate the Markov assumption inherent in standard Latent SDEs\. While it might attempt to implicitly encode past\-time features within its parametric latent spaceθ\\theta, the neural SDE model does not explicitly account for delay effects\. Therefore, latent SDE, or neural latent SDE, has limitations for modeling astronomical time series data like quasar light curve\.

Instead, the suggested model uses past observations in highlighted window \(see[2](https://arxiv.org/html/2609.20906#S3.F2)\), and learns the dynamics from incomplete data by capturing delayed & stochastic dependencies\.

![Refer to caption](https://arxiv.org/html/2609.20906v1/2.png)Figure 2:Neural SDDE learning window by Oh
### 3\.2From Neural SDE to Neural SDDE

Neural SDE is usually written in the form:

\(3\.1\)d​z​\(t\)=γ⁡\(z⁡\(t\),θγ\)​d​t\+σ⁡\(z⁡\(t\),θσ\)​d​W​\(t\),t≥0\\mathrm\{d\}z\(t\)\\;=\\;\\gamma\\bigl\(z\(t\);\\,\\theta\_\{\\gamma\}\)\\,\\mathrm\{d\}t\\;\+\\;\\sigma\(z\(t\);\\,\\theta\_\{\\sigma\}\)\\,\\mathrm\{d\}W\(t\),\\qquad t\\geq 0
where the driftγ\\gammagoverns the deterministic evolution, and diffusionσ\\sigmagoverns stochasticity\. Both ofγ\\gamma,σ\\sigmaare neural networks parameterized byθγ\\theta\_\{\\gamma\}andθσ\\theta\_\{\\sigma\}\. By explicitly considering delayed effects,\[[17](https://arxiv.org/html/2609.20906#bib.bib15)\]extend this system to define the Neural Stochastic Delay Differential Equations \(Neural SDDEs\):

\(3\.2\)d​z​\(t\)=γ⁡\(z⁡\(t\),z⁡\(t−τ\),θγ\)​d​t\+σ⁡\(t,θσ\)​d​W​\(t\),t≥0\\mathrm\{d\}z\(t\)\\;=\\;\\gamma\\bigl\(z\(t\),z\(t\-\\tau\);\\,\\theta\_\{\\gamma\}\)\\,\\mathrm\{d\}t\\;\+\\;\\sigma\(t;\\,\\theta\_\{\\sigma\}\)\\,\\mathrm\{d\}W\(t\),\\qquad t\\geq 0
with an initial segmentz⁡\(t\)=φ⁡\(t\)z\(t\)=\\varphi\(t\)fort∈\[−τ,0\]t\\in\[\-\\tau,0\], whereτ\>0\\tau\>0is a fixed delay hyperparameter\. Now, the driftγθ\\gamma\_\{\\theta\}depends on bothz⁡\(t\)z\(t\)and the delayed statez⁡\(t−τ\)z\(t\-\\tau\)\. This change allows the model to capture long\-range dependencies in the process\.

Function\-space Markov propertyAlthough \([3\.2](https://arxiv.org/html/2609.20906#S3.E2)\) is non\-Markovian as a process onℝd\\mathbb\{R\}^\{d\}, it becomes Markovian whenzt\(⋅\):=z\(t\+⋅\)z\_\{t\}\(\\cdot\):=z\(t\+\\cdot\)is considered as an element of the function spaceC:=C⁡\(\[−τ,0\],ℝd\)C:=C\(\[\-\\tau,0\];\\mathbb\{R\}^\{d\}\)\.zt​\(⋅\)z\_\{t\}\(\\cdot\)is called aCC\-valued Markov process\[[14](https://arxiv.org/html/2609.20906#bib.bib12)\]\.

Reconstruction PropertySDDEs exhibit a unique reconstruction property, which implies that the system’s initial history function can be recovered using only a future segment of the solution path, notably without requiring knowledge of the specific noise trajectory\. This characteristic enables the adjoint method for memory\-efficient backpropagation

Augmented state with controlled pathsTo inject information from irregular observations into the continuous\-time dynamics, Oh et al\. borrow from Neural CDEs\[[12](https://arxiv.org/html/2609.20906#bib.bib10)\]and interpolate the raw observations to a continuous pathX⁡\(t\)X\(t\)\. The state is then augmented asz¯​\(t\)=ζ⁡\(t,z⁡\(t\),X⁡\(t\),θζ\)\\bar\{z\}\(t\)=\\zeta\(t,z\(t\),X\(t\);\\theta\_\{\\zeta\}\), whereζ\\zetais a neural network,z⁡\(t\)z\(t\)is the current latent state, andX⁡\(t\)X\(t\)is a controlled path\.

### 3\.3Adjoint Method for SDDEs

Computing∂ℒ/∂θ\\partial\\mathcal\{L\}/\\partial\\thetathrough \([3\.2](https://arxiv.org/html/2609.20906#S3.E2)\) requires a delay\-aware version of the stochastic adjoint\. The key complication is that the delayτ\\tauin the forward pass becomes an*advance*in the backward pass: the adjointa⁡\(t\)=∂ℒ/∂z⁡\(t\)a\(t\)=\\partial\\mathcal\{L\}/\\partial z\(t\)receives contributions both from the current\-time derivative∂γ/∂z⁡\(t\)\\partial\\gamma/\\partial z\(t\)and from a future timet\+τt\+\\tau\(wherez⁡\(t\)z\(t\)appears as the delayed argument of the drift\)\. Oh et al\. implement this by partitioning the time interval into segments of length≤τ\\leq\\tauand solving the backward SDE on each segment, carrying forward the necessary future adjoint values\.

### 3\.4Experiments and Results

The experiments use the ELAsTiCC and PLAsTiCC datasets \(simulated LSST\-like light curves across many object classes\)\. Four scenarios are considered: standard supervised classification, classification with50%50\\%missing labels, novelty detection, and joint classification and novelty detection with missing labels\. Baselines include GRU\-D, Neural ODE, Neural SDE, Neural CDE, ODE\-RNN, and Neural LSDE\. The Neural SDDE consistently achieves the highest accuracy and weighted F1, with the largest margins in the missing\-label and novelty\-detection settings\. Sensitivity toτ\\tauis mild\.

## 4Delayed\-Memory Gradient Descent

### 4\.1Mathematical Definitions

#### 4\.1\.1Neural SDDEs in Langevin form

A special class of Neural SDDEs is the Neural Langevin\-type Stochastic Delay Differential Equation \(Neural LSDDE\), inspired from\[[17](https://arxiv.org/html/2609.20906#bib.bib15)\]and Langevin dynamics\. LetX⁡\(t\)∈ℝdXX\(t\)\\in\\mathbb\{R\}^\{d\_\{X\}\}be the dynamic state at current timett, then the Neural Stochastic Delay Differential Equation \(Neural SDDE\) can be written as:

\(4\.1\)d​X​\(t\)=γ⁡\(X⁡\(t\),X⁡\(t−τ\),θγ\)​d​t\+σ⁡\(t,X⁡\(t\),X⁡\(t−τ\),θσ\)​d​W​\(t\)\\mathrm\{d\}X\(t\)=\\gamma\(X\(t\),X\(t\-\\tau\);\\theta\_\{\\gamma\}\)\\mathrm\{d\}t\+\\sigma\(t,X\(t\),X\(t\-\\tau\);\\theta\_\{\\sigma\}\)\\mathrm\{d\}W\(t\)
whereγ:ℝdX×ℝdX→ℝdX\\gamma:\\mathbb\{R\}^\{d\_\{X\}\}\\times\\mathbb\{R\}^\{d\_\{X\}\}\\rightarrow\\mathbb\{R\}^\{d\_\{X\}\}, andσ:ℝ\+×ℝdX→ℝm\\sigma:\\mathbb\{R\}\_\{\+\}\\times\\mathbb\{R\}^\{d\_\{X\}\}\\rightarrow\\mathbb\{R\}^\{m\}are drift and diffusion functions, andγ\\gammais a latent neural network parameterized byθγ\\theta\_\{\\gamma\}, andθγ\\theta\_\{\\gamma\}is learned directly from available data \(Similar forθσ\\theta\_\{\\sigma\}\)\.

In \([4\.1](https://arxiv.org/html/2609.20906#S4.E1)\),τ\>0\\tau\>0is a fixed delay, andX⁡\(t−τ\)X\(t\-\\tau\)denote the dynamic state at past timet−τt\-\\tau\. Because of this delay, the system can not be initialized by a single pointX⁡\(0\)X\(0\)but with a full history segment of the trajectory:

\(4\.2\)X\(t\)=ϕ\(t\),∀t∈\[−τ,0\],and the functionϕ:\[−τ,0\]→ℝdXX\(t\)=\\phi\(t\),\\text\{ \}\\forall t\\in\[\-\\tau,0\],\\qquad\\text\{ and the function \}\\phi:\[\-\\tau,0\]\\rightarrow\\mathbb\{R\}^\{d\_\{X\}\}
Equation \([4\.1](https://arxiv.org/html/2609.20906#S4.E1)\) explicitly integrates past statesX⁡\(t−τ\)X\(t\-\\tau\), allowing the model to capture memory effects\. Also,τ\\tauserves as a lookback window for the model to learn temporal dependencies\.

#### 4\.1\.2Stochastic Delay Gradient Flow

The continuous time deterministic Delayed\-Memory Gradient Flow can be written as:

\(4\.3\)d​Xt=−\(∇f​\(Xt\)\+λ⁡\(Xt−Xt−τ\)\)​d​t\\mathrm\{d\}X\_\{t\}=\-\(\\nabla f\(X\_\{t\}\)\+\\lambda\(X\_\{t\}\-X\_\{t\-\\tau\}\)\)\\mathrm\{d\}t
with the gradient\-flow term−∇f​\(Xt\)\-\\nabla f\(X\_\{t\}\)and the delayed\-memory feedback term−λ⁡\(Xt−Xt−τ\)\-\\lambda\(X\_\{t\}\-X\_\{t\-\\tau\}\)\. The parameterλ\\lambdahere is a delay\-coupling strength parameter, determining how strong the delay\-memory effect is\.

We can also create an extension of \([4\.3](https://arxiv.org/html/2609.20906#S4.E3)\) by introducing stochasticity, capturing system uncertainty and robustness to noise:

\(4\.4\)d​Xt=−\(∇f​\(Xt\)\+λ⁡\(Xt−Xt−τ\)\)​d​t\+η​Σ​\(Xt,Xt−τ\)1/2​d​Wt\\mathrm\{d\}X\_\{t\}=\-\(\\nabla f\(X\_\{t\}\)\+\\lambda\(X\_\{t\}\-X\_\{t\-\\tau\}\)\)\\mathrm\{d\}t\+\\sqrt\{\\eta\}\\Sigma\(X\_\{t\},X\_\{t\-\\tau\}\)^\{1/2\}\\mathrm\{d\}W\_\{t\}
withWtW\_\{t\}add\-dimensional Brownian motion,Σ\\Sigmarepresenting the diffusion matrix, andη\\etarepresenting the learning rate\. Notice that \([4\.4](https://arxiv.org/html/2609.20906#S4.E4)\) is a special case of \([4\.1](https://arxiv.org/html/2609.20906#S4.E1)\), withγ=−\(∇f​\(Xt\)\+λ⁡\(Xt−Xt−τ\)\)\\gamma=\-\(\\nabla f\(X\_\{t\}\)\+\\lambda\(X\_\{t\}\-X\_\{t\-\\tau\}\)\), andσ=η​Σ​\(Xt,Xt−τ\)1/2\\sigma=\\sqrt\{\\eta\}\\Sigma\(X\_\{t\},X\_\{t\-\\tau\}\)^\{1/2\}\.

### 4\.2Discrete\-time Algorithms

#### 4\.2\.1Delayed\-Memory Stochastic Gradient Descent

We start with the basic discrete\-time gradient algorithm by considering non\-stochastic and non\-latent\-neural case, by transform \([4\.3](https://arxiv.org/html/2609.20906#S4.E3)\) into a iterated version\. Letm∈ℕm\\in\\mathbb\{N\}be a delay parameter representing how many steps furthermost does the algorithm consider, and letη\>0\\eta\>0be the learning rate,λ≥0\\lambda\\geq 0be the coupling parameter \(also called memory strength parameter\), we can define the discrete history segment as:

x−m,x−m\+1,…,x−1,x0∈ℝdX, then:x\_\{\-m\},x\_\{\-m\+1\},\\dots,x\_\{\-1\},x\_\{0\}\\in\\mathbb\{R\}^\{d\_\{X\}\}\\text\{, then:\}\(4\.5\)xk\+1=xk−η∇f\(xk\)\+ηλ\(xk−m−xk\)x\_\{k\+1\}=x\_\{k\}\-\\eta\\nabla f\(x\_\{k\}\)\+\\eta\\lambda\(x\_\{k\-m\}\-x\_\{k\}\)
is the basic Delayed\-Memory Gradient Descent \(DMGD\) iteration formula\.

In real optimization problems, the Gradient Descent algorithm is usually outperformed in practice by Stochastc Gradient Descent \(SGD\), which generally uses a mini\-batch subset of data instead of the entire dataset to determine the gradient at each iterationkk\. The advantages of SGD are faster updates, reduced memory capacity, and more likely to escape local minima or saddle points\. Therefore, we can design a counterpart of SGD for the basic DMGD:

\(4\.6\)xk\+1=xk−η∇f⁡\(xk\)^\+ηλ\(xk−m−xk\)x\_\{k\+1\}=x\_\{k\}\-\\eta\\nabla\\widehat\{f\(x\_\{k\}\)\}\+\\eta\\lambda\(x\_\{k\-m\}\-x\_\{k\}\)
where∇f⁡\(xk\)^=\(1/b\)​\[∑i∈Bk∇fi​\(xk\)\]\\nabla\\widehat\{f\(x\_\{k\}\)\}=\(1/b\)\[\\sum\\limits\_\{i\\in B\_\{k\}\}\\nabla f\_\{i\}\(x\_\{k\}\)\], andBkB\_\{k\}is a mini\-batch with\|Bk\|=b\|B\_\{k\}\|=brandomly selected from the entire dataset at eachkk\.

Now, focusing on the termλ⁡\(xk−m−xk\)\\lambda\(x\_\{k\-m\}\-x\_\{k\}\), which is called the delay feedback term, the memory strength parameterλ≥0\\lambda\\geq 0determines how strongly the past state affects the current dynamics\. Overall, this term’s role can be described as ”besides moving downhill, the system also feels a restoring/pulling force toward its past state, or it penalizes moving too far from the past\.”

#### 4\.2\.2Continuous\-Delayed\-Memory SGD

In many astronomical systems’ data, like irregularly sampled quasar light curves and quasar visibility, from different distances and different radii, bands, and reprocessing regions, the delays can not ideally be concentrated at a single lagτ\\tau\. Instead, the present response may depend on a distribution of past lags, reflecting propagation, scattering, and reprocessing across multiple spatial or physical regions\[[2](https://arxiv.org/html/2609.20906#bib.bib19)\]\. Thus, a single\-delay termXt−Xt−τX\_\{t\}\-X\_\{t\-\\tau\}may be too restrictive once the basic delayed\-memory mechanism has been understood\. The integral

∫0τΛ⁡\(s\)​\(Xt−Xt−s\)​𝑑s\\int\_\{0\}^\{\\tau\}\\Lambda\(s\)\\bigl\(X\_\{t\}\-X\_\{t\-s\}\\bigr\)\\,\\mathrm\{d\}sis used to model a continuum of delays instead of a single delayτ\\tau, andΛ⁡\(s\)\\Lambda\(s\)is a memory kernel function representing how much influence size\-ssdelay has, wheres∈\[0,τ\]s\\in\[0,\\tau\]\.Λ\\Lambdadescribes how strongly past states contributes to the present update through a distributed memory mechanism\. In this sense,Λ⁡\(s\)\\Lambda\(s\)serves as a distributed\-delay analogue of the coupling parameterλ⁡\(s\)\\lambda\(s\)in the Basic Delayed\-Memory model\. Now, if only focusing on the deterministic part of the continuous\-time flow:

\(4\.7\)d​Xt=−\(∇f​\(Xt\)\+∫0τΛ⁡\(s\)​\(Xt−Xt−s\)​𝑑s\)​𝑑t\\mathrm\{d\}X\_\{t\}=\-\(\\nabla f\(X\_\{t\}\)\+\\int\_\{0\}^\{\\tau\}\\Lambda\(s\)\\bigl\(X\_\{t\}\-X\_\{t\-s\}\\bigr\)\\,\\mathrm\{d\}s\)\\mathrm\{d\}t
Therefore, for a discrete\-time algorithm, let\{wj\}j=1m\\\{w\_\{j\}\\\}\_\{j=1\}^\{m\}be the set of nonnegative weights approximating the continuous kernel, then \([4\.7](https://arxiv.org/html/2609.20906#S4.E7)\) can be written as an iteration map:

\(4\.8\)xk\+1=xk−η∇f\(xk\)\+η∑j=1mwj\(xk−j−xk\)x\_\{k\+1\}=x\_\{k\}\-\\eta\\nabla f\(x\_\{k\}\)\+\\eta\\sum\_\{j=1\}^\{m\}w\_\{j\}\(x\_\{k\-j\}\-x\_\{k\}\)
Besides, the Basic Delayed\-Memory model can be viewed as a limiting\-special case of the Continuous\-Delayed\-Memory model withΛ⁡\(s\)\\Lambda\(s\):ifΛ⁡\(s\)=λ​δ​\(s−τ0\)\\Lambda\(s\)=\\lambda\\,\\delta\(s\-\\tau\_\{0\}\), whereδ\\deltais the Dirac\-delta function atτ0\\tau\_\{0\}\. Thus, \([4\.7](https://arxiv.org/html/2609.20906#S4.E7)\), \([4\.8](https://arxiv.org/html/2609.20906#S4.E8)\) can be simplified as

dXt=−\(∇f\(Xt\)\+λ\(Xt−Xt−τ0\)\)dt,xk\+1=xk−η∇f\(xk\)\+ηλ\(xk−m−xk\)\\mathrm\{d\}X\_\{t\}=\-\(\\nabla f\(X\_\{t\}\)\+\\lambda\(X\_\{t\}\-X\_\{t\-\\tau\_\{0\}\}\)\)\\mathrm\{d\}t,\\qquad\\text\{ \}x\_\{k\+1\}=x\_\{k\}\-\\eta\\nabla f\(x\_\{k\}\)\+\\eta\\lambda\(x\_\{k\-m\}\-x\_\{k\}\)Hence, the Continuous\-Delayed\-Memory model is a strictly generalized version of the previous model inSection 4\.2\.1\. Finally, by considering the stochastic part, we add the diffusion term to \([4\.7](https://arxiv.org/html/2609.20906#S4.E7)\):

\(4\.9\)dXt=−\(∇f\(Xt\)\+∫0τΛ\(s\)\(Xt−Xt−s\)ds\)dt\+ηΣ\(Xt,\{Xt−s\)1/2\}s∈\[0,τ\]dWt\\mathrm\{d\}X\_\{t\}=\-\(\\nabla f\(X\_\{t\}\)\+\\int\_\{0\}^\{\\tau\}\\Lambda\(s\)\\bigl\(X\_\{t\}\-X\_\{t\-s\}\\bigr\)\\,\\mathrm\{d\}s\)\\mathrm\{d\}t\+\\sqrt\{\\eta\}\\Sigma\(X\_\{t\},\\\{X\_\{t\-s\}\)^\{1/2\}\\\}\_\{s\\in\[0,\\tau\]\}\\mathrm\{d\}W\_\{t\}

### 4\.3Solution Analysis about Stochastic Delayed Differential Equation \(SDDE\)

In this section, the goal is to mathematically analyze the assumptions under which the SDDE that is a continuous\-time approximation to our designed discrete algorithms has existed, unique, and global solutions, specifically \([4\.4](https://arxiv.org/html/2609.20906#S4.E4)\), \([4\.9](https://arxiv.org/html/2609.20906#S4.E9)\)\. Because every SDDE is a special, often simpler, case of a Stochastic Functional Differential Equation \(SFDE\), and there have many existed theories about assumptions required for the existence and uniqueness of solutions to SFDE, like in\[[18](https://arxiv.org/html/2609.20906#bib.bib17)\], all findings are inherited from\[[18](https://arxiv.org/html/2609.20906#bib.bib17)\]but we create the variants only specific to SDDEs \([4\.4](https://arxiv.org/html/2609.20906#S4.E4)\)and \([4\.9](https://arxiv.org/html/2609.20906#S4.E9)\)\.

#### 4\.3\.1Necessary condition for SDDE \([4\.4](https://arxiv.org/html/2609.20906#S4.E4)\)

To prove that this exact model has a unique, global \(non\-exploding\) solution using the Khasminskii\-type approach \(bypassing strict linear growth\), the following two specific mathematical conditions must hold\.

Condition 1: Local Lipschitz Continuity \(Guarantees Uniqueness\)

LetXtX\_\{t\}be the current state andXt−τX\_\{t\-\\tau\}be the delayed state, for any local region bounded by radiusRR,∃KR∈ℝ\\exists\\text\{ \}K\_\{R\}\\in\\mathbb\{R\}, s\.t\.

∀two pairs of states,​\(Xt1,Xt1−τ1\),\(Xt2,Xt2−τ2\),\\displaystyle\\forall\\text\{ two pairs of states, \}\(X\_\{\\displaystyle t\_\{1\}\},X\_\{\\displaystyle t\_\{1\}\-\\tau\_\{1\}\}\),\(X\_\{\\displaystyle t\_\{2\}\},X\_\{\\displaystyle t\_\{2\}\-\\tau\_\{2\}\}\),\|−∇f​\(Xt1\)−λ⁡\(Xt1−Xt1−τ1\)−\(−∇f​\(Xt2\)−λ⁡\(Xt2−Xt2−τ2\)\)\|2\\displaystyle\\left\|\-\\nabla f\(X\_\{\\displaystyle t\_\{1\}\}\)\-\\lambda\(X\_\{\\displaystyle t\_\{1\}\}\-X\_\{\\displaystyle t\_\{1\}\-\\tau\_\{1\}\}\)\-\(\-\\nabla f\(X\_\{\\displaystyle t\_\{2\}\}\)\-\\lambda\(X\_\{\\displaystyle t\_\{2\}\}\-X\_\{\\displaystyle t\_\{2\}\-\\tau\_\{2\}\}\)\)\\right\|^\{2\}∨\\displaystyle\\vee\|η​Σ​\(Xt1,Xt1−τ1\)1/2−η​Σ​\(Xt2,Xt2−τ2\)1/2\|2\\displaystyle\\left\|\\sqrt\{\\eta\}\\Sigma\(X\_\{\\displaystyle t\_\{1\}\},X\_\{\\displaystyle t\_\{1\}\-\\tau\_\{1\}\}\)^\{1/2\}\-\\sqrt\{\\eta\}\\Sigma\(X\_\{\\displaystyle t\_\{2\}\},X\_\{\\displaystyle t\_\{2\}\-\\tau\_\{2\}\}\)^\{1/2\}\\right\|^\{2\}≤KR​\(\|Xt1−Xt2\|2\+\|Xt1−τ1−Xt2−τ2\|2\)\\displaystyle\\leq K\_\{R\}\(\\left\|X\_\{\\displaystyle t\_\{1\}\}\-X\_\{\\displaystyle t\_\{2\}\}\\right\|^\{2\}\+\\left\|X\_\{\\displaystyle t\_\{1\}\-\\tau\_\{1\}\}\-X\_\{\\displaystyle t\_\{2\}\-\\tau\_\{2\}\}\\right\|^\{2\}\)
Because the delay termλ⁡\(Xt−Xt−τ\)\\lambda\(X\_\{t\}\-X\_\{t\-\\tau\}\)is strictly linear, it is automatically globally Lipschitz\. Thus, whether Lipschitz Continuity satisfied or not is entirely determined on objective functionffand noiseΣ\\Sigma\.

Because the drift function, which contains the gradient∇f​\(Xt\)\\nabla f\(X\_\{t\}\), should be locally Lipschitz continuous, for any local bounded region \(a compact setEE\), there must exist a constantLE\>0L\_\{E\}\>0such that for anyXt1,Xt2∈EX\_\{t\_\{1\}\},X\_\{t\_\{2\}\}\\in E,‖∇f​\(Xt1\)−∇f​\(Xt2\)‖≤LE​‖Xt1−Xt2‖\\\|\\nabla f\(X\_\{t\_\{1\}\}\)\-\\nabla f\(X\_\{t\_\{2\}\}\)\\\|\\leq L\_\{E\}\\\|X\_\{t\_\{1\}\}\-X\_\{t\_\{2\}\}\\\|\.

Supposef⁡\(Xt\)∈C2f\(X\_\{t\}\)\\in C^\{2\}, then its second derivative, the Hessian matrix∇2f​\(Xt\)\\nabla^\{2\}f\(X\_\{t\}\), exists and is continuous everywhere\.

By Multidimensional Mean Value Theorem, the difference between the gradients at two stagesXt1X\_\{t\_\{1\}\}andXt2X\_\{t\_\{2\}\}can be expressed as an integral of the Hessian along the straight line path between them:

∇f​\(Xt1\)−∇f​\(Xt2\)=\(∫01∇2f​\(Xt2\+t⁡\(Xt1−Xt2\)\)​𝑑t\)​\(Xt1−Xt2\)\\nabla f\(X\_\{t\_\{1\}\}\)\-\\nabla f\(X\_\{t\_\{2\}\}\)=\\left\(\\int\_\{0\}^\{1\}\\nabla^\{2\}f\\big\(X\_\{t\_\{2\}\}\+t\(X\_\{t\_\{1\}\}\-X\_\{t\_\{2\}\}\)\\big\)dt\\right\)\(X\_\{t\_\{1\}\}\-X\_\{t\_\{2\}\}\)
Normalize both sides and by Submultiplicativity of Induced Matrix Norm, Triangle Inequality for Integrals:

\(4\.10\)‖∇f​\(Xt1\)−∇f​\(Xt2\)‖≤\(∫01‖∇2f​\(Xt2\+t⁡\(Xt1−Xt2\)\)‖​𝑑t\)​‖Xt1−Xt2‖\\\|\\nabla f\(X\_\{t\_\{1\}\}\)\-\\nabla f\(X\_\{t\_\{2\}\}\)\\\|\\leq\\left\(\\int\_\{0\}^\{1\}\\big\\\|\\nabla^\{2\}f\\big\(X\_\{t\_\{2\}\}\+t\(X\_\{t\_\{1\}\}\-X\_\{t\_\{2\}\}\)\\big\)\\big\\\|dt\\right\)\\\|X\_\{t\_\{1\}\}\-X\_\{t\_\{2\}\}\\\|
Becausef⁡\(Xt\)∈C2f\(X\_\{t\}\)\\in C^\{2\}, its Hessian∇2f​\(Xt\)\\nabla^\{2\}f\(X\_\{t\}\)is still a continuous function\. By Extreme value theorem,∇2f​\(Xt\)\\nabla^\{2\}f\(X\_\{t\}\)evaluated on a closed, bounded regionEEmust have a finite maximum value\. So, the normalized maximum bound valueLE=maxZt∈E⁡‖∇2f​\(Zt\)‖L\_\{E\}=\\max\_\{Z\_\{t\}\\in E\}\\\|\\nabla^\{2\}f\(Z\_\{t\}\)\\\|exists\.

Then, replace the integral in \([4\.10](https://arxiv.org/html/2609.20906#S4.E10)\) with this maximum bound value, we have

‖∇f​\(Xt1\)−∇f​\(Xt2\)‖≤\(∫01LE​𝑑t\)​‖Xt1−Xt2‖=LE​‖Xt1−Xt2‖\\\|\\nabla f\(X\_\{t\_\{1\}\}\)\-\\nabla f\(X\_\{t\_\{2\}\}\)\\\|\\leq\\left\(\\int\_\{0\}^\{1\}L\_\{E\}dt\\right\)\\\|X\_\{t\_\{1\}\}\-X\_\{t\_\{2\}\}\\\|=L\_\{E\}\\\|X\_\{t\_\{1\}\}\-X\_\{t\_\{2\}\}\\\|Similarly, the diffusion termσ⁡\(Xt,Xt−τ\)=η​Σ​\(Xt,Xt−τ\)1/2\\sigma\(X\_\{t\},X\_\{t\-\\tau\}\)=\\sqrt\{\\eta\}\\Sigma\(X\_\{t\},X\_\{t\-\\tau\}\)^\{1/2\}should also be locally Lipschitz:

\(4\.11\)‖σ⁡\(Xt1,Xt1−τ1\)−σ⁡\(Xt2,Xt2−τ2\)‖≤ME​\(‖Xt1−Xt2‖\+‖Xt1−τ1−Xt2−τ2‖\)\\\|\\sigma\(X\_\{t\_\{1\}\},X\_\{t\_\{1\}\-\\tau\_\{1\}\}\)\-\\sigma\(X\_\{t\_\{2\}\},X\_\{t\_\{2\}\-\\tau\_\{2\}\}\)\\\|\\leq M\_\{E\}\\big\(\\\|X\_\{t\_\{1\}\}\-X\_\{t\_\{2\}\}\\\|\+\\\|X\_\{t\_\{1\}\-\\tau\_\{1\}\}\-X\_\{t\_\{2\}\-\\tau\_\{2\}\}\\\|\\big\)
Now, suppose thatσ⁡\(Xt,Xt−τ\)∈C1\\sigma\(X\_\{t\},X\_\{t\-\\tau\}\)\\in C^\{1\}, then its Jacobian matrixJσJ\_\{\\sigma\}is continuous\. Then, on any bounded compact regionEE, this continuous Jacobian will attain a finite maximum normME=maxXt∈E⁡‖Jσ​\(Xt\)‖M\_\{E\}=\\max\_\{X\_\{t\}\\in E\}\\\|J\_\{\\sigma\}\(X\_\{t\}\)\\\|\.

By Extreme value theorem, this result guarantees \([4\.11](https://arxiv.org/html/2609.20906#S4.E11)\)\.

Therefore,∀Xt∈E,\(f⁡\(Xt\)∈C2\)∧\(σ⁡\(Xt,Xt−τ\)∈C1\)⟹\(LE=KRCLOSE\\forall X\_\{t\}\\in E,\\big\(f\(X\_\{t\}\)\\in C^\{2\}\\big\)\\land\\big\(\\sigma\(X\_\{t\},X\_\{t\-\\tau\}\)\\in C^\{1\}\\big\)\\implies\\big\(L\_\{E\}=K\_\{R\}is the Lipschitz constant\)\\big\)\.

\(Note:In stochastic analysis, authors often blanket\-assumeC2C^\{2\}smoothness for all coefficients—both drift and diffusion\. This is becauseC2C^\{2\}continuity is strictly required later to applyItô’s Lemma\. Itô’s Lemma fundamentally relies on a second\-order Taylor expansion to calculate the Infinitesimal Generatorℒ​V​\(x,y\)\\mathcal\{L\}V\(x,y\)that is used in the Khasminskii theorem\)\.

Condition 2: The Khasminskii Lyapunov Bound \(Guarantees Global Existence\)

ByTheorem 2\.6in\[[18](https://arxiv.org/html/2609.20906#bib.bib17)\],

we construct a non\-negative Lyapunov function \(also called ”containment” function\)V⁡\(Xt\)V\(X\_\{t\}\)\. Specifically, for our SDDE model, the Infinitesimal Generator operatorℒ\\mathcal\{L\}acting on a twice\-differentiable functionV⁡\(Xt\)V\(X\_\{t\}\)is defined using Itô’s Lemma as:

ℒV\(Xt,Xt−τ\)=∇V\(Xt\)⊤γ\(Xt,Xt−τ\)\+12Tr\[σ\(Xt,Xt−τ\)⊤∇2V\(Xt\)σ\(Xt,Xt−τ\)\]\\mathcal\{L\}V\(X\_\{t\},X\_\{t\-\\tau\}\)=\\nabla V\(X\_\{t\}\)^\{\\top\}\\gamma\(X\_\{t\},X\_\{t\-\\tau\}\)\+\\frac\{1\}\{2\}\\text\{Tr\}\\Big\[\\sigma\(X\_\{t\},X\_\{t\-\\tau\}\)^\{\\top\}\\nabla^\{2\}V\(X\_\{t\}\)\\sigma\(X\_\{t\},X\_\{t\-\\tau\}\)\\Big\]whereγ\\gammais the drift vector andσ\\sigmais the diffusion matrix\. Substituting our specific driftγ=−\(∇f​\(Xt\)\+λ⁡\(Xt−Xt−τ\)\)\\gamma=\-\\big\(\\nabla f\(X\_\{t\}\)\+\\lambda\(X\_\{t\}\-X\_\{t\-\\tau\}\)\\big\)and diffusionσ=η​Σ​\(Xt,Xt−τ\)1/2\\sigma=\\sqrt\{\\eta\}\\Sigma\(X\_\{t\},X\_\{t\-\\tau\}\)^\{1/2\}:

\(4\.12\)ℒ​V​\(Xt,Xt−τ\)=−∇V\(Xt\)⊤\(∇f\(Xt\)\+λ\(Xt−Xt−τ\)\)\+η2​Tr​\[Σ​\(Xt,Xt−τ\)1/2​∇2V​\(Xt\)​Σ​\(Xt,Xt−τ\)1/2\]\\begin\{split\}\\mathcal\{L\}V\(X\_\{t\},X\_\{t\-\\tau\}\)=&\-\\nabla V\(X\_\{t\}\)^\{\\top\}\\big\(\\nabla f\(X\_\{t\}\)\+\\lambda\(X\_\{t\}\-X\_\{t\-\\tau\}\)\\big\)\\\\ &\+\\frac\{\\eta\}\{2\}\\text\{Tr\}\\Big\[\\Sigma\(X\_\{t\},X\_\{t\-\\tau\}\)^\{1/2\}\\nabla^\{2\}V\(X\_\{t\}\)\\Sigma\(X\_\{t\},X\_\{t\-\\tau\}\)^\{1/2\}\\Big\]\\end\{split\}
To make this mathematically actionable, let us test the standard quadratic Lyapunov containment function:V⁡\(Xt\)=‖Xt‖2V\(X\_\{t\}\)=\\\|X\_\{t\}\\\|^\{2\}\. This implies its gradient is∇V​\(Xt\)=2​Xt\\nabla V\(X\_\{t\}\)=2X\_\{t\}and its Hessian is∇2V​\(Xt\)=2​I\\nabla^\{2\}V\(X\_\{t\}\)=2I\(whereIIis the identity matrix\)\.

Plugging these exact derivatives into our generator \([4\.12](https://arxiv.org/html/2609.20906#S4.E12)\) yields:

\(4\.13\)ℒ​V​\(Xt,Xt−τ\)=2​Xt⊤​\(−∇f​\(Xt\)−λ​Xt\+λ​Xt−τ\)\+η2​Tr​\[Σ​\(Xt,Xt−τ\)1/2​\(2​I\)​Σ​\(Xt,Xt−τ\)1/2\]\\begin\{split\}\\mathcal\{L\}V\(X\_\{t\},X\_\{t\-\\tau\}\)=&2X\_\{t\}^\{\\top\}\\Big\(\-\\nabla f\(X\_\{t\}\)\-\\lambda X\_\{t\}\+\\lambda X\_\{t\-\\tau\}\\Big\)\\\\ &\+\\frac\{\\eta\}\{2\}\\text\{Tr\}\\Big\[\\Sigma\(X\_\{t\},X\_\{t\-\\tau\}\)^\{1/2\}\(2I\)\\Sigma\(X\_\{t\},X\_\{t\-\\tau\}\)^\{1/2\}\\Big\]\\end\{split\}
Since matrix trace satisfies the cyclic propertyTr​\[A​B\]=Tr​\[B​A\]\\text\{Tr\}\[AB\]=\\text\{Tr\}\[BA\], the diffusion term simplifies perfectly toη​Tr​\[Σ⁡\(Xt,Xt−τ\)\]\\eta\\text\{Tr\}\[\\Sigma\(X\_\{t\},X\_\{t\-\\tau\}\)\]\. Expanding the inner product for the drift term, we obtain:

\(4\.14\)ℒV\(Xt,Xt−τ\)=−2Xt⊤∇f\(Xt\)−2λ∥Xt∥2\+2λ\(Xt⊤Xt−τ\)\+ηTr\[Σ\(Xt,Xt−τ\)\]\\mathcal\{L\}V\(X\_\{t\},X\_\{t\-\\tau\}\)=\-2X\_\{t\}^\{\\top\}\\nabla f\(X\_\{t\}\)\-2\\lambda\\\|X\_\{t\}\\\|^\{2\}\+2\\lambda\(X\_\{t\}^\{\\top\}X\_\{t\-\\tau\}\)\+\\eta\\text\{Tr\}\\Big\[\\Sigma\(X\_\{t\},X\_\{t\-\\tau\}\)\\Big\]
By applying Young’s inequality, which states that2​\(Xt⊤​Xt−τ\)≤‖Xt‖2\+‖Xt−τ‖22\(X\_\{t\}^\{\\top\}X\_\{t\-\\tau\}\)\\leq\\\|X\_\{t\}\\\|^\{2\}\+\\\|X\_\{t\-\\tau\}\\\|^\{2\}, we can establish a strict upper bound for the cross\-term\. Substituting this inequality into \([4\.14](https://arxiv.org/html/2609.20906#S4.E14)\):

ℒV\(Xt,Xt−τ\)≤−2Xt⊤∇f\(Xt\)−λ∥Xt∥2\+λ∥Xt−τ∥2\+ηTr\[Σ\(Xt,Xt−τ\)\]\\mathcal\{L\}V\(X\_\{t\},X\_\{t\-\\tau\}\)\\leq\-2X\_\{t\}^\{\\top\}\\nabla f\(X\_\{t\}\)\-\\lambda\\\|X\_\{t\}\\\|^\{2\}\+\\lambda\\\|X\_\{t\-\\tau\}\\\|^\{2\}\+\\eta\\text\{Tr\}\\Big\[\\Sigma\(X\_\{t\},X\_\{t\-\\tau\}\)\\Big\]
The Khasminskii condition dictates that for the system state to not explode to infinity, there must exist a bounding constantc\>0c\>0such that the generator is constrained by the current and delayed size ofVV:

ℒ​V​\(Xt,Xt−τ\)≤c⁡\(1\+‖Xt‖2\+‖Xt−τ‖2\)\\mathcal\{L\}V\(X\_\{t\},X\_\{t\-\\tau\}\)\\leq c\\big\(1\+\\\|X\_\{t\}\\\|^\{2\}\+\\\|X\_\{t\-\\tau\}\\\|^\{2\}\\big\)
Therefore, for the DMGD SDDE to admit a unique, global, non\-exploding solution, the objective functionf⁡\(Xt\)f\(X\_\{t\}\)and the diffusion matrixΣ\\Sigmamust rigorously satisfy the following inequality for some constantc\>0c\>0:

\(4\.15\)−2Xt⊤∇f\(Xt\)\+ηTr\[Σ\(Xt,Xt−τ\)\]≤c\(1\+∥Xt∥2\)\+\(c−λ\)∥Xt−τ∥2\+λ∥Xt∥2\-2X\_\{t\}^\{\\top\}\\nabla f\(X\_\{t\}\)\+\\eta\\text\{Tr\}\\Big\[\\Sigma\(X\_\{t\},X\_\{t\-\\tau\}\)\\Big\]\\leq c\(1\+\\\|X\_\{t\}\\\|^\{2\}\)\+\(c\-\\lambda\)\\\|X\_\{t\-\\tau\}\\\|^\{2\}\+\\lambda\\\|X\_\{t\}\\\|^\{2\}
In summary, \([4\.15](https://arxiv.org/html/2609.20906#S4.E15)\) and∀Xt∈E,\(f⁡\(Xt\)∈C2\)∧\(σ⁡\(Xt,Xt−τ\)∈C1\)\\forall X\_\{t\}\\in E,\\big\(f\(X\_\{t\}\)\\in C^\{2\}\\big\)\\land\\big\(\\sigma\(X\_\{t\},X\_\{t\-\\tau\}\)\\in C^\{1\}\\big\)should both satisfied to garantee unique, global solution for \([4\.4](https://arxiv.org/html/2609.20906#S4.E4)\)\.

#### 4\.3\.2Necessary condition for SDDE \([4\.9](https://arxiv.org/html/2609.20906#S4.E9)\)

Following the Khasminskii\-type approach as done inSection 4\.3\.1, and adapting it specifically for the integral delay functional, the model \([4\.9](https://arxiv.org/html/2609.20906#S4.E9)\) possesses a unique, global solution if the following conditions holds:

Condition 1: Local Lipschitz Continuity \(Guarantees Uniqueness\)

For any local compact region bounded by radiusRR, denoted asEE, let\(Xt1,\{Xt1−s\}s∈\[0,τ\]\)\(X\_\{t\_\{1\}\},\\\{X\_\{t\_\{1\}\-s\}\\\}\_\{s\\in\[0,\\tau\]\}\)and\(Xt2,\{Xt2−s\}s∈\[0,τ\]\)\(X\_\{t\_\{2\}\},\\\{X\_\{t\_\{2\}\-s\}\\\}\_\{s\\in\[0,\\tau\]\}\)be two state paths inEE\. To guarantee uniqueness,∃KR\>0\\exists K\_\{R\}\>0such that:

\(4\.16\)‖γ⁡\(Xt1,\{Xt1−s\}\)−γ⁡\(Xt2,\{Xt2−s\}\)‖2∨‖σ⁡\(Xt1,\{Xt1−s\}\)−σ⁡\(Xt2,\{Xt2−s\}\)‖2≤KR​\(‖Xt1−Xt2‖2\+sups∈\[0,τ\]‖Xt1−s−Xt2−s‖2\)\\begin\{split\}&\\left\\\|\\gamma\(X\_\{t\_\{1\}\},\\\{X\_\{t\_\{1\}\-s\}\\\}\)\-\\gamma\(X\_\{t\_\{2\}\},\\\{X\_\{t\_\{2\}\-s\}\\\}\)\\right\\\|^\{2\}\\vee\\left\\\|\\sigma\(X\_\{t\_\{1\}\},\\\{X\_\{t\_\{1\}\-s\}\\\}\)\-\\sigma\(X\_\{t\_\{2\}\},\\\{X\_\{t\_\{2\}\-s\}\\\}\)\\right\\\|^\{2\}\\\\ &\\leq K\_\{R\}\\Big\(\\left\\\|X\_\{t\_\{1\}\}\-X\_\{t\_\{2\}\}\\right\\\|^\{2\}\+\\sup\_\{s\\in\[0,\\tau\]\}\\left\\\|X\_\{t\_\{1\}\-s\}\-X\_\{t\_\{2\}\-s\}\\right\\\|^\{2\}\\Big\)\\end\{split\}
By linearity and the Triangle Inequality for Integrals, the distributed delay term satisfies:

‖∫0τΛ⁡\(s\)​\(Xt1−Xt1−s\)​ds−∫0τΛ⁡\(s\)​\(Xt2−Xt2−s\)​ds‖≤∫0τΛ⁡\(s\)​‖Xt1−Xt2‖​𝑑s\+∫0τΛ⁡\(s\)​‖Xt1−s−Xt2−s‖​𝑑s≤LΛ​\(‖Xt1−Xt2‖\+sups∈\[0,τ\]‖Xt1−s−Xt2−s‖\)\\begin\{split\}&\\left\\\|\\int\_\{0\}^\{\\tau\}\\Lambda\(s\)\(X\_\{t\_\{1\}\}\-X\_\{t\_\{1\}\-s\}\)ds\-\\int\_\{0\}^\{\\tau\}\\Lambda\(s\)\(X\_\{t\_\{2\}\}\-X\_\{t\_\{2\}\-s\}\)ds\\right\\\|\\\\ &\\leq\\int\_\{0\}^\{\\tau\}\\Lambda\(s\)\\left\\\|X\_\{t\_\{1\}\}\-X\_\{t\_\{2\}\}\\right\\\|ds\+\\int\_\{0\}^\{\\tau\}\\Lambda\(s\)\\left\\\|X\_\{t\_\{1\}\-s\}\-X\_\{t\_\{2\}\-s\}\\right\\\|ds\\\\ &\\leq L\_\{\\Lambda\}\\left\(\\left\\\|X\_\{t\_\{1\}\}\-X\_\{t\_\{2\}\}\\right\\\|\+\\sup\_\{s\\in\[0,\\tau\]\}\\left\\\|X\_\{t\_\{1\}\-s\}\-X\_\{t\_\{2\}\-s\}\\right\\\|\\right\)\\end\{split\}
whereLΛ=∫0τΛ⁡\(s\)​𝑑s<∞L\_\{\\Lambda\}=\\int\_\{0\}^\{\\tau\}\\Lambda\(s\)ds<\\infty\. Thus, the integral delay term is globally Lipschitz\.

For the objective functionffand diffusionσ\\sigma, assume\(f⁡\(Xt\)∈C2\)∧\(σ⁡\(Xt,\{Xt−s\}\)∈C1\)\\big\(f\(X\_\{t\}\)\\in C^\{2\}\\big\)\\land\\big\(\\sigma\(X\_\{t\},\\\{X\_\{t\-s\}\\\}\)\\in C^\{1\}\\big\)\. By the Extreme Value Theorem on the compact setEE, the continuous Hessian∇2f\\nabla^\{2\}fand JacobianJσJ\_\{\\sigma\}attain finite maximum norms:

LE=maxZ∈E⁡‖∇2f​\(Z\)‖,ME=maxZ∈E⁡‖Jσ​\(Z\)‖L\_\{E\}=\\max\_\{Z\\in E\}\\\|\\nabla^\{2\}f\(Z\)\\\|,\\quad M\_\{E\}=\\max\_\{Z\\in E\}\\\|J\_\{\\sigma\}\(Z\)\\\|
Applying the Multidimensional Mean Value Theorem yields:

‖∇f​\(Xt1\)−∇f​\(Xt2\)‖≤LE​‖Xt1−Xt2‖\\\|\\nabla f\(X\_\{t\_\{1\}\}\)\-\\nabla f\(X\_\{t\_\{2\}\}\)\\\|\\leq L\_\{E\}\\\|X\_\{t\_\{1\}\}\-X\_\{t\_\{2\}\}\\\|‖σ⁡\(Xt1,\{Xt1−s\}\)−σ⁡\(Xt2,\{Xt2−s\}\)‖≤ME​\(‖Xt1−Xt2‖\+sups∈\[0,τ\]‖Xt1−s−Xt2−s‖\)\\\|\\sigma\(X\_\{t\_\{1\}\},\\\{X\_\{t\_\{1\}\-s\}\\\}\)\-\\sigma\(X\_\{t\_\{2\}\},\\\{X\_\{t\_\{2\}\-s\}\\\}\)\\\|\\leq M\_\{E\}\\Big\(\\\|X\_\{t\_\{1\}\}\-X\_\{t\_\{2\}\}\\\|\+\\sup\_\{s\\in\[0,\\tau\]\}\\\|X\_\{t\_\{1\}\-s\}\-X\_\{t\_\{2\}\-s\}\\\|\\Big\)
Therefore, combining the bounded constantsLΛ,LEL\_\{\\Lambda\},L\_\{E\}, andMEM\_\{E\}algebraically guarantees the existence of a finiteKRK\_\{R\}, strictly satisfying \([4\.16](https://arxiv.org/html/2609.20906#S4.E16)\), and thus uniqueness of the solution\.

Condition 2: The Khasminskii Lyapunov Bound \(Guarantees Global Existence\)

FollowingTheorem 2\.6in\[[18](https://arxiv.org/html/2609.20906#bib.bib17)\]for Stochastic Functional Differential Equations, we apply the Infinitesimal Generatorℒ\\mathcal\{L\}to the quadratic containment functionV⁡\(Xt\)=‖Xt‖2V\(X\_\{t\}\)=\\\|X\_\{t\}\\\|^\{2\}\.

Applying Itô’s Lemma and substituting the derivatives∇V​\(Xt\)=2​Xt\\nabla V\(X\_\{t\}\)=2X\_\{t\}and∇2V​\(Xt\)=2​I\\nabla^\{2\}V\(X\_\{t\}\)=2I, the generator expands to:

ℒ​V​\(Xt,\{Xt−s\}\)=2​Xt⊤​\(−∇f​\(Xt\)−∫0τΛ⁡\(s\)​\(Xt−Xt−s\)​𝑑s\)\+η2​Tr​\[Σ​\(Xt,\{Xt−s\}\)1/2​\(2​I\)​Σ​\(Xt,\{Xt−s\}\)1/2\]\\begin\{split\}\\mathcal\{L\}V\(X\_\{t\},\\\{X\_\{t\-s\}\\\}\)=&2X\_\{t\}^\{\\top\}\\Big\(\-\\nabla f\(X\_\{t\}\)\-\\int\_\{0\}^\{\\tau\}\\Lambda\(s\)\(X\_\{t\}\-X\_\{t\-s\}\)ds\\Big\)\\\\ &\+\\frac\{\\eta\}\{2\}\\text\{Tr\}\\Big\[\\Sigma\(X\_\{t\},\\\{X\_\{t\-s\}\\\}\)^\{1/2\}\(2I\)\\Sigma\(X\_\{t\},\\\{X\_\{t\-s\}\\\}\)^\{1/2\}\\Big\]\\end\{split\}
Utilizing the cyclic property of the trace matrix and expanding the inner product into the integral, we obtain:

\(4\.17\)ℒ​V​\(Xt,\{Xt−s\}\)=−2Xt⊤∇f\(Xt\)−2∫0τΛ\(s\)∥Xt∥2ds\+2∫0τΛ\(s\)\(Xt⊤Xt−s\)ds\+ηTr\[Σ\(Xt,\{Xt−s\}\)\]\\begin\{split\}\\mathcal\{L\}V\(X\_\{t\},\\\{X\_\{t\-s\}\\\}\)=&\-2X\_\{t\}^\{\\top\}\\nabla f\(X\_\{t\}\)\-2\\int\_\{0\}^\{\\tau\}\\Lambda\(s\)\\\|X\_\{t\}\\\|^\{2\}ds\\\\ &\+2\\int\_\{0\}^\{\\tau\}\\Lambda\(s\)\(X\_\{t\}^\{\\top\}X\_\{t\-s\}\)ds\+\\eta\\text\{Tr\}\\Big\[\\Sigma\(X\_\{t\},\\\{X\_\{t\-s\}\\\}\)\\Big\]\\end\{split\}
By applying Young’s inequality,2​\(Xt⊤​Xt−s\)≤‖Xt‖2\+‖Xt−s‖22\(X\_\{t\}^\{\\top\}X\_\{t\-s\}\)\\leq\\\|X\_\{t\}\\\|^\{2\}\+\\\|X\_\{t\-s\}\\\|^\{2\}, we can establish a strict upper bound for the cross\-term inside the integral\. Assuming the weighting function is non\-negative \(Λ⁡\(s\)≥0\\Lambda\(s\)\\geq 0\), this yields:

ℒ​V​\(Xt,\{Xt−s\}\)≤−2Xt⊤∇f\(Xt\)−2∫0τΛ\(s\)∥Xt∥2ds\+∫0τΛ\(s\)∥Xt∥2ds\+∫0τΛ\(s\)∥Xt−s∥2ds\+ηTr\[Σ\(Xt,\{Xt−s\}\)\]\\begin\{split\}\\mathcal\{L\}V\(X\_\{t\},\\\{X\_\{t\-s\}\\\}\)\\leq&\-2X\_\{t\}^\{\\top\}\\nabla f\(X\_\{t\}\)\-2\\int\_\{0\}^\{\\tau\}\\Lambda\(s\)\\\|X\_\{t\}\\\|^\{2\}ds\\\\ &\+\\int\_\{0\}^\{\\tau\}\\Lambda\(s\)\\\|X\_\{t\}\\\|^\{2\}ds\+\\int\_\{0\}^\{\\tau\}\\Lambda\(s\)\\\|X\_\{t\-s\}\\\|^\{2\}ds\+\\eta\\text\{Tr\}\\Big\[\\Sigma\(X\_\{t\},\\\{X\_\{t\-s\}\\\}\)\\Big\]\\end\{split\}
Which simplifies perfectly to match the structure of our previous single\-delay bound:

ℒV\(Xt,\{Xt−s\}\)≤−2Xt⊤∇f\(Xt\)−∫0τΛ\(s\)∥Xt∥2ds\+∫0τΛ\(s\)∥Xt−s∥2ds\+ηTr\[Σ\(Xt,\{Xt−s\}\)\]\\mathcal\{L\}V\(X\_\{t\},\\\{X\_\{t\-s\}\\\}\)\\leq\-2X\_\{t\}^\{\\top\}\\nabla f\(X\_\{t\}\)\-\\int\_\{0\}^\{\\tau\}\\Lambda\(s\)\\\|X\_\{t\}\\\|^\{2\}ds\+\\int\_\{0\}^\{\\tau\}\\Lambda\(s\)\\\|X\_\{t\-s\}\\\|^\{2\}ds\+\\eta\\text\{Tr\}\\Big\[\\Sigma\(X\_\{t\},\\\{X\_\{t\-s\}\\\}\)\\Big\]
The general Khasminskii condition for SFDEs dictates that the generator must be bounded by a constantc\>0c\>0that scales with the current state and the supremum of the delayed state path:

ℒ​V​\(Xt,\{Xt−s\}\)≤c⁡\(1\+‖Xt‖2\+sups∈\[0,τ\]‖Xt−s‖2\)\\mathcal\{L\}V\(X\_\{t\},\\\{X\_\{t\-s\}\\\}\)\\leq c\\Big\(1\+\\\|X\_\{t\}\\\|^\{2\}\+\\sup\_\{s\\in\[0,\\tau\]\}\\\|X\_\{t\-s\}\\\|^\{2\}\\Big\)Therefore, for the distributed\-delay SDDE \([4\.9](https://arxiv.org/html/2609.20906#S4.E9)\) to guarantee a global solution, the objective functionf⁡\(Xt\)f\(X\_\{t\}\)and the diffusion matrixΣ\\Sigmamust satisfy the following inequality for some constantc\>0c\>0:

−2Xt⊤∇f\(Xt\)\+ηTr\[Σ\(Xt,\{Xt−s\}\)\]≤c⁡\(1\+‖Xt‖2\+sups∈\[0,τ\]‖Xt−s‖2\)\+∫0τΛ\(s\)\(∥Xt∥2−∥Xt−s∥2\)ds\\begin\{split\}\-2X\_\{t\}^\{\\top\}\\nabla f\(X\_\{t\}\)\+\\eta\\text\{Tr\}\\Big\[\\Sigma\(X\_\{t\},\\\{X\_\{t\-s\}\\\}\)\\Big\]\\leq&\\;c\\Big\(1\+\\\|X\_\{t\}\\\|^\{2\}\+\\sup\_\{s\\in\[0,\\tau\]\}\\\|X\_\{t\-s\}\\\|^\{2\}\\Big\)\\\\ &\+\\int\_\{0\}^\{\\tau\}\\Lambda\(s\)\\Big\(\\\|X\_\{t\}\\\|^\{2\}\-\\\|X\_\{t\-s\}\\\|^\{2\}\\Big\)ds\\end\{split\}

### 4\.42D Simulation Results

#### 4\.4\.1Simulation on 2D Convex Quadratic Loss Landscape

We first generate a 2\-dimensional synthetic dataset representing a loss landscape off⁡\(x,y\)=x2\+5​y2f\(x,y\)=x^\{2\}\+5y^\{2\}with 20000 random data points, which is a typical convex quadratic loss function\. Then, we control the random seed 1000, batch sizeb=4b=4, learning rateη=0\.05\\eta=0\.05delay\-coupling strength parameterλ=5\\lambda=5, number of stepsk=30k=30, and we gradually changes the delayed stepsmmin the set\{1,2,8\}\\\{1,2,8\\\}\. Then, we run both the Vanilla SGD and the Delayed\-Memory SGD onffstarting at initial point\(−9\.6,5\.6\)\(\-9\.6,5\.6\)\. The SGD trajectory is in red and the DMSGD trajectory is in green\. \(See[3](https://arxiv.org/html/2609.20906#S4.F3), where the loss map is restricted only to\[−10\.4,10\.4\]×\[−6,6\]\[\-10\.4,10\.4\]\\times\[\-6,6\]since the trajectories are all in this range\)

The simulation shows that when keeping every other parameter constant and with small batch size, the DMSGD is likely to have more fluctuated moves on the loss landscape as the delayed stepsmmincreases, which is caused by the ”long\-range” phase lag forcing the optimizer to react to the geometry of a more distant region of the landscape rather than the local curvature\.

Again, we use the same dataset, initial point, random seed 1000, and batch sizeb=4b=4, to perform both SGD and DMSGD\. In the 1st round, we set the DMSGD parametersm=5m=5,η=0\.004\\eta=0\.004,λ=245\\lambda=245,k=500k=500; in the 2nd round, we changed the DMSGD parametersλ=14\\lambda=14,k=12000k=12000\(changingkkto generate enough iteration that goes to the local minimum\)\.

See comparison plots in[3](https://arxiv.org/html/2609.20906#S4.F3), bigλ\\lambdawould lead DMSGD algorithm into divergent behavior, which helps the algorithm to explore a wider landscape in each single run\. However, asλ\\lambdabecomes smaller, the DMSGD trajectory is actually more smoothly\-convergent to the local minimum compared to SGD’s\.

![Refer to caption](https://arxiv.org/html/2609.20906v1/3.png)

![Refer to caption](https://arxiv.org/html/2609.20906v1/4.png)

![Refer to caption](https://arxiv.org/html/2609.20906v1/5.png)

Figure 3:SDG \(red\) & DMSGD \(green\) trajectories asm∈\{1,2,8\}m\\in\\\{1,2,8\\\}![Refer to caption](https://arxiv.org/html/2609.20906v1/6.png)

![Refer to caption](https://arxiv.org/html/2609.20906v1/8.png)

Figure 4:SDG \(red\) & Continuous\-DMSGD \(green\) trajectories asλ∈\{245,14\}\\lambda\\in\\\{245,14\\\}
#### 4\.4\.2Simulation on 2D Nonconvex Rastrigin Loss Landscape

The, we generate another 2\-dimensional synthetic dataset representing a loss landscape off2​\(𝐱\)=20\+∑i=12\[xi2−10​cos⁡\(2​π​xi\)\]f\_\{2\}\(\\mathbf\{x\}\)=20\+\\sum\_\{i=1\}^\{2\}\\left\[x\_\{i\}^\{2\}\-10\\cos\(2\\pi x\_\{i\}\)\\right\]with 20000 random data points and starting point\(−4\.5,3\.5\)\(\-4\.5,3\.5\)\. We use the same random seed 1000, and run Vanilla SGD and the Continuous\-Delayed\-Memory SGD onf2f\_\{2\}\. We set the batch sizeb=8b=8, learning rateη=0\.0075\\eta=0\.0075, number of stepsk=600k=600, delayed stepsm=5m=5with uniform initial delay\-coupling strength weight\[14,14,14,14,14\]⊤\[14,14,14,14,14\]^\{\\top\}\. \(See[5](https://arxiv.org/html/2609.20906#S4.F5)\)

![Refer to caption](https://arxiv.org/html/2609.20906v1/7.png)Figure 5:SDG \(red\) & Continuous\-DMSGD \(blue\) trajectories with uniform initial delay\-coupling strength weightIn this scenario, the vanilla SGD quickly descends into the nearest local minimum and becomes permanently trapped due to its lack of momentum to overcome the walls\. Conversely, the Continuous\-DMSGD algorithm introduces intentional fluctuations into the path due to its pulling\-back force to past locations stored in the history\. These ”restoring forces” allow the Continuous\-DMSGD optimizer to escape various local traps and successfully navigate across the rugged terrain to converge at the global minimum in this example\.

In another example, both SGD and Continuous\-DMSDG optimizers are navigating on the 2D Styblinski\-Tang Loss, represented byf⁡\(x,y\)=12​\[\(x4−16​x2\+5​x\)\+\(y4−16​y2\+5​y\)\]f\(x,y\)=\\frac\{1\}\{2\}\\left\[\(x^\{4\}\-16x^\{2\}\+5x\)\+\(y^\{4\}\-16y^\{2\}\+5y\)\\right\]\. \(See[6](https://arxiv.org/html/2609.20906#S4.F6)\) In this certain parameter set up, the 2 optimizers finally stuck into different local minimums, which additionally supports that SGD and Continuous\-DMSDG has distinct behaviors in low\-dimensional space\.

![Refer to caption](https://arxiv.org/html/2609.20906v1/9.png)Figure 6:SDG \(red\) & Continuous\-DMSGD \(green\) trajectories finalizes into different local minimums
#### 4\.4\.3Potential advantages of Continuous\-DMSGD over Vanilla SGD

The magnitude of theΛ\\Lambdadictates a critical trade\-off between convergence stability and exploration\. WhenΛ\\Lambdais small, the optimization steps becomes smaller in a local minimum region, ensuring faster convergence to the minimum\. WhenΛ\\Lambdais large, the step sizes progressively expand and allows the optimizer to explore larger regions, but it might lead to divergence\. So, finding the boundary ofΛ\\Lambda’s magnitude that determines the algorithm’s different behavior is a potential future direction\.

Observed Advantages from simulations:

1\. In Highly Non\-Convex Landscapes, especially those with alternating local maximums and local minimums, Continuous\-DMSGD is more likely to escape some local minimum compared to Vanilla SGD through the pulling force toward its past state\. This is because the pulling force from the memory term provides energy to jump over local barriers, which override the zero\-gradient traps of local minima\. \[Increase overall magnitude ofΛ\\Lambda\]

2\. In very sharp local minimum region, or even non\-smooth holes, Continuous\-DMSGD is more likely to finally converge into lower Loss regions than Vanilla SGD by penalizing moving too far from the past\. \[Use very smallΛ\\Lambdamagnitude\]

## 5From the Stochastic Adjoint to Continuous\-Time Reinforcement Learning

This section develops a backward Stratonovich SDE that carries the stochastic adjoint method described in\[[13](https://arxiv.org/html/2609.20906#bib.bib11)\], which is the close cousin of the backward SDE that carries the costate in stochastic optimal control, and hence of the Bellman/HJB machinery of continuous\-time reinforcement learning\. We proceed in three steps: \(i\) classical stochastic control and its BSDE adjoint; \(ii\) continuous\-time RL with deterministic policies, following\[[15](https://arxiv.org/html/2609.20906#bib.bib13)\]; and \(iii\) the exploratory \(entropy\-regularised\) formulation of\[[19](https://arxiv.org/html/2609.20906#bib.bib16)\], together with what we will call the*Exploratory Backward Stratonovich SDE*\.

### 5\.1Controlled SDEs and the Stochastic Control Problem

Let the stateX⁡\(t\)∈ℝnX\(t\)\\in\\mathbb\{R\}^\{n\}evolve according to the controlled SDE

\(5\.1\)d​X​\(t\)=b⁡\(X⁡\(t\),u⁡\(t\)\)​d​t\+σ⁡\(X⁡\(t\),u⁡\(t\)\)​d​W​\(t\),X⁡\(0\)=x0\\mathrm\{d\}X\(t\)\\;=\\;b\\bigl\(X\(t\),u\(t\)\\bigr\)\\,\\mathrm\{d\}t\\;\+\\;\\sigma\\bigl\(X\(t\),u\(t\)\\bigr\)\\,\\mathrm\{d\}W\(t\),\\qquad X\(0\)=x\_\{0\}whereu⁡\(t\)∈𝒰⊂ℝku\(t\)\\in\\mathcal\{U\}\\subset\\mathbb\{R\}^\{k\}is the control \(action\) at timett, chosen by the decision maker, andW⁡\(t\)W\(t\)is anmm\-dimensional Brownian motion\. The cost functional over a horizonTTis

\(5\.2\)J⁡\(u\)=𝔼⁡\[∫0TL⁡\(X⁡\(s\),u⁡\(s\)\)​𝑑s\+Φ⁡\(X⁡\(T\)\)\]J\(u\)\\;=\\;\\mathbb\{E\}\\\!\\left\[\\int\_\{0\}^\{T\}L\\bigl\(X\(s\),u\(s\)\\bigr\)\\,\\mathrm\{d\}s\+\\Phi\\bigl\(X\(T\)\\bigr\)\\right\]whereLLis a running cost andΦ\\Phiis a terminal cost\. The value function isV⁡\(t,x\):=infu𝔼⁡\[∫tTL​𝑑s\+Φ⁡\(X⁡\(T\)\)∣X⁡\(t\)=x\]V\(t,x\):=\\inf\_\{u\}\\mathbb\{E\}\[\\int\_\{t\}^\{T\}L\\,\\mathrm\{d\}s\+\\Phi\(X\(T\)\)\\mid X\(t\)=x\]and, under standard regularity, satisfies theHamilton–Jacobi–Bellman \(HJB\) equation

\(5\.3\)∂tV\+infu∈𝒰\{b​\(x,u\)⊤​∇xV\+12​tr​\(σ​σ⊤​\(x,u\)​∇x2V\)\+L⁡\(x,u\)\}=0,V⁡\(T,x\)=Φ⁡\(x\)\\partial\_\{t\}V\+\\inf\_\{u\\in\\mathcal\{U\}\}\\left\\\{b\(x,u\)^\{\\top\}\\nabla\_\{x\}V\+\\tfrac\{1\}\{2\}\\hbox\{\\rm tr\}\\bigl\(\\sigma\\sigma^\{\\top\}\(x,u\)\\,\\nabla\_\{x\}^\{2\}V\\bigr\)\+L\(x,u\)\\right\\\}=0,\\qquad V\(T,x\)=\\Phi\(x\)

### 5\.2The Adjoint BSDE from Pontryagin’s Principle

An alternative route to the optimal control is Pontryagin’s stochastic maximum principle, which expresses optimality through a pair\(Y⁡\(t\),Z⁡\(t\)\)\(Y\(t\),Z\(t\)\)of adapted processes satisfying a*backward stochastic differential equation*\(BSDE\) of the Pardoux–Peng type:

\(5\.4\)\{−d​Y​\(t\)=\[bx​\(X,u\)⊤​Y​\(t\)\+tr​\(σx​\(X,u\)⊤​Z​\(t\)\)\+Lx​\(X,u\)\]​d​t−Z⁡\(t\)​d​W​\(t\),Y⁡\(T\)=Φx​\(X​\(T\)\)\\left\\\{\\begin\{aligned\} \-\\mathrm\{d\}Y\(t\)&=\\bigl\[\\,b\_\{x\}\(X,u\)^\{\\top\}Y\(t\)\+\\hbox\{\\rm tr\}\\\!\\bigl\(\\sigma\_\{x\}\(X,u\)^\{\\top\}Z\(t\)\\bigr\)\+L\_\{x\}\(X,u\)\\,\\bigr\]\\,\\mathrm\{d\}t\-Z\(t\)\\,\\mathrm\{d\}W\(t\),\\\\ Y\(T\)&=\\Phi\_\{x\}\\bigl\(X\(T\)\\bigr\)\\end\{aligned\}\\right\.where subscripts denote partial derivatives inxx\. The processY⁡\(t\)Y\(t\)is the costate \(adjoint\), andZ⁡\(t\)Z\(t\)is the martingale representation term introduced to ensure the solution is adapted to the forward filtration\. The optimal control maximises the Hamiltonian

\(5\.5\)H⁡\(x,u,y,z\)=b​\(x,u\)⊤​y\+tr​\(σ​\(x,u\)⊤​z\)−L⁡\(x,u\)\.H\(x,u,y,z\)\\;=\\;b\(x,u\)^\{\\top\}y\+\\hbox\{\\rm tr\}\\\!\\bigl\(\\sigma\(x,u\)^\{\\top\}z\\bigr\)\-L\(x,u\)\.
The structural parallel with the stochastic adjoint of Li et al\.We find that the adjoint SDE of\[[13](https://arxiv.org/html/2609.20906#bib.bib11)\]is a pathwise backward SDE driven by the*forward*Brownian path, whereas \([5\.4](https://arxiv.org/html/2609.20906#S5.E4)\) is a Pardoux–Peng BSDE whose solution must be adapted and hence involves a new martingale termZ​d​WZ\\,\\mathrm\{d\}W\. But they are the same*kind*of object: in both, a costate is propagated backward in time, driven in part by a Brownian motion, and used to read off gradients \(or, in the control case, the Hamiltonian\-maximizing action\)\.

### 5\.3Continuous\-Time RL with Deterministic Policies \(Munos 2006\)

\[[15](https://arxiv.org/html/2609.20906#bib.bib13)\]studies policy gradient methods for deterministic feedback policiesu⁡\(t\)=πθ​\(X⁡\(t\)\)u\(t\)=\\pi\_\{\\theta\}\(X\(t\)\)in continuous time\. Substituting the policy into \([5\.1](https://arxiv.org/html/2609.20906#S5.E1)\) gives an SDE parameterised byθ\\theta:

\(5\.6\)d​Xθ​\(t\)=b⁡\(Xθ​\(t\),πθ​\(Xθ​\(t\)\)\)​d​t\+σ⁡\(Xθ​\(t\),πθ​\(Xθ​\(t\)\)\)​d​W​\(t\)\\mathrm\{d\}X^\{\\theta\}\(t\)\\;=\\;b\\bigl\(X^\{\\theta\}\(t\),\\pi\_\{\\theta\}\(X^\{\\theta\}\(t\)\)\\bigr\)\\,\\mathrm\{d\}t\\;\+\\;\\sigma\\bigl\(X^\{\\theta\}\(t\),\\pi\_\{\\theta\}\(X^\{\\theta\}\(t\)\)\\bigr\)\\,\\mathrm\{d\}W\(t\)Now, compare this to the neural SDE defined in\[[13](https://arxiv.org/html/2609.20906#bib.bib11)\]as

\(5\.7\)d​z​\(t\)=fθ​\(z⁡\(t\),t\)​d​t\+gϕ​\(z⁡\(t\),t\)​d​W​\(t\),t∈\[0,T\]\\mathrm\{d\}z\(t\)\\;=\\;f\_\{\\theta\}\\bigl\(z\(t\),t\\bigr\)\\,\\mathrm\{d\}t\\;\+\\;g\_\{\\phi\}\\bigl\(z\(t\),t\\bigr\)\\,\\mathrm\{d\}W\(t\),\\qquad t\\in\[0,T\]We find that these two are the similar, with parametersθ\\thetainside both the drift and the diffusion\. The policy gradient∇θJ​\(πθ\)\\nabla\_\{\\theta\}J\(\\pi\_\{\\theta\}\)can thus be computed in two equivalent ways:

1. 1\.via the HJB value function:∇θJ=𝔼⁡\[∫0T\(∂uH\)​\(X,u,∇V,∇2V\)​∇θπθ​\(X\)​𝑑t\]\\nabla\_\{\\theta\}J=\\mathbb\{E\}\\\!\\bigl\[\\int\_\{0\}^\{T\}\(\\partial\_\{u\}H\)\(X,u,\\nabla V,\\nabla^\{2\}V\)\\,\\nabla\_\{\\theta\}\\pi\_\{\\theta\}\(X\)\\,\\mathrm\{d\}t\\bigr\]
2. 2\.via the stochastic adjoint of\[[13](https://arxiv.org/html/2609.20906#bib.bib11)\]: solve the forward SDE \([5\.6](https://arxiv.org/html/2609.20906#S5.E6)\), then solve the backward Stratonovich SDE against the reward\-based lossℒ⁡\(Xθ​\(⋅\)\)=∫L​𝑑s\+Φ⁡\(Xθ​\(T\)\)\\mathcal\{L\}\(X^\{\\theta\}\(\\cdot\)\)=\\int L\\,\\mathrm\{d\}s\+\\Phi\(X^\{\\theta\}\(T\)\), and read off∇θJ\\nabla\_\{\\theta\}Jfrom the augmented adjoint\.

The second route makes the adjoint of\[[13](https://arxiv.org/html/2609.20906#bib.bib11)\]a tool for continuous\-time RL: policy gradients can be computed using the same Virtual\-Brownian\-Tree machinery used for Neural SDE VAEs\.

### 5\.4The Exploratory Formulation of Wang et al\. \(2020\)

The deterministic\-policy formulation above has a well\-known shortcoming from an RL perspective: it offers no systematic mechanism for exploration\.\[[19](https://arxiv.org/html/2609.20906#bib.bib16)\]proposed an*exploratory*relaxation of the stochastic control problem in which the control is, at each time, a*probability distribution*over actions rather than a single deterministic value, and the cost is augmented by an entropy term that rewards exploration\.

Letπt\(⋅∣x\)\\pi\_\{t\}\(\\cdot\\mid x\)be a \(measurable\) family of probability densities on𝒰\\mathcal\{U\}, one for each\(t,x\)\(t,x\)\. The*exploratory dynamics*are obtained by averaging the original drift and diffusion underπ\\pi:

\(5\.8\)d​X​\(t\)=b~​\(X⁡\(t\),πt\)​d​t\+σ~​\(X⁡\(t\),πt\)​d​W​\(t\)\\mathrm\{d\}X\(t\)\\;=\\;\\tilde\{b\}\\bigl\(X\(t\),\\pi\_\{t\}\\bigr\)\\,\\mathrm\{d\}t\\;\+\\;\\tilde\{\\sigma\}\\bigl\(X\(t\),\\pi\_\{t\}\\bigr\)\\,\\mathrm\{d\}W\(t\)where

\(5\.9\)b~​\(x,π\)=∫𝒰b⁡\(x,u\)​π​\(𝑑u∣x\),σ~​σ~⊤​\(x,π\)=∫𝒰σ​σ⊤​\(x,u\)​π​\(𝑑u∣x\)\\tilde\{b\}\(x,\\pi\)=\\int\_\{\\mathcal\{U\}\}b\(x,u\)\\,\\pi\(\\mathrm\{d\}u\\mid x\),\\qquad\\tilde\{\\sigma\}\\tilde\{\\sigma\}^\{\\top\}\(x,\\pi\)=\\int\_\{\\mathcal\{U\}\}\\sigma\\sigma^\{\\top\}\(x,u\)\\,\\pi\(\\mathrm\{d\}u\\mid x\)The matrix square rootσ~\\tilde\{\\sigma\}of the averaged diffusion matrix is chosen measurably\. The entropy\-regularised cost is

\(5\.10\)Jλ​\(π\)=𝔼⁡\[∫0T\(∫𝒰L⁡\(X,u\)​πs​\(𝑑u∣X\)\+λ​∫𝒰log⁡πs​\(u∣X\)​πs​\(𝑑u∣X\)\)​𝑑s\+Φ⁡\(X⁡\(T\)\)\]J^\{\\lambda\}\(\\pi\)\\;=\\;\\mathbb\{E\}\\\!\\left\[\\int\_\{0\}^\{T\}\\\!\\left\(\\int\_\{\\mathcal\{U\}\}L\(X,u\)\\,\\pi\_\{s\}\(\\mathrm\{d\}u\\mid X\)\+\\lambda\\int\_\{\\mathcal\{U\}\}\\log\\pi\_\{s\}\(u\\mid X\)\\,\\pi\_\{s\}\(\\mathrm\{d\}u\\mid X\)\\right\)\\mathrm\{d\}s\+\\Phi\(X\(T\)\)\\right\]whereλ\>0\\lambda\>0is the exploration temperature\. The classical Wang–Zhou result for linear–quadratic problems is that the optimalπt⋆\\pi^\{\\star\}\_\{t\}is*Gaussian*with mean determined by the classical LQ optimal control and variance determined byλ\\lambda; more generally the optimal exploratory policy is a Gibbs distribution proportional toexp\(−H¯/λ\)\\exp\(\-\\bar\{H\}/\\lambda\), whereH¯\\bar\{H\}is an appropriately defined Hamiltonian evaluated against the value function of the exploratory problem\.

### 5\.5The Exploratory Backward Stratonovich SDE

The natural question is:*what is the analog of the Li et al\. stochastic adjoint for the exploratory controlled SDE\([5\.8](https://arxiv.org/html/2609.20906#S5.E8)\)?*Answering it gives what we call anExploratory Backward Stratonovich SDE\(EB\-SSDE\)\.

Setup\.Fix an exploratory policyπθ\\pi\_\{\\theta\}with parametersθ\\theta\(e\.g\. a neural network outputting the parameters of a distribution on𝒰\\mathcal\{U\}\)\. Write the exploratory SDE \([5\.8](https://arxiv.org/html/2609.20906#S5.E8)\) in Stratonovich form:

\(5\.11\)d​X​\(t\)=b~θ​\(X⁡\(t\)\)​d​t\+σ~θ​\(X⁡\(t\)\)∘d​W​\(t\)\\mathrm\{d\}X\(t\)\\;=\\;\\tilde\{b\}\_\{\\theta\}\\bigl\(X\(t\)\\bigr\)\\,\\mathrm\{d\}t\\;\+\\;\\tilde\{\\sigma\}\_\{\\theta\}\\bigl\(X\(t\)\\bigr\)\\circ\\mathrm\{d\}W\(t\)whereb~θ\\tilde\{b\}\_\{\\theta\}includes the Itô–Stratonovich correction ofb~\\tilde\{b\}, and the drift and diffusion now depend onθ\\thetathrough the averaging operation\. Define the loss functional

\(5\.12\)ℒ⁡\(θ\)=∫0T\(L¯​\(X⁡\(s\),πθ,s\)\+λ​ℋ​\(πθ,s\)\)​𝑑s\+Φ⁡\(X⁡\(T\)\)\\mathcal\{L\}\(\\theta\)\\;=\\;\\int\_\{0\}^\{T\}\\\!\\left\(\\bar\{L\}\\bigl\(X\(s\),\\pi\_\{\\theta,s\}\\bigr\)\+\\lambda\\,\\mathcal\{H\}\\bigl\(\\pi\_\{\\theta,s\}\\bigr\)\\right\)\\mathrm\{d\}s\+\\Phi\\bigl\(X\(T\)\\bigr\)whereL¯​\(x,π\)=∫L⁡\(x,u\)​π​\(𝑑u∣x\)\\bar\{L\}\(x,\\pi\)=\\int L\(x,u\)\\,\\pi\(\\mathrm\{d\}u\\mid x\)andℋ⁡\(π\)=∫log⁡π​𝑑π\\mathcal\{H\}\(\\pi\)=\\int\\log\\pi\\,\\mathrm\{d\}\\pi\.

The EB\-SSDE\.Then \([5\.11](https://arxiv.org/html/2609.20906#S5.E11)\) with loss \([5\.12](https://arxiv.org/html/2609.20906#S5.E12)\) will yield an adjoint statea⁡\(t\)=∂ℒ/∂X⁡\(t\)a\(t\)=\\partial\\mathcal\{L\}/\\partial X\(t\)satisfying the backward Stratonovich SDE

\(5\.13\)da\(t\)=−a\(t\)⊤∂xb~θ\(X\(t\)\)dt−a\(t\)⊤∂xσ~θ\(X\(t\)\)∘dW\(t\)−∂xL¯\(X\(t\),πθ,t\)dt\\boxed\{\\;\\mathrm\{d\}a\(t\)\\;=\\;\-\\,a\(t\)^\{\\top\}\\,\\partial\_\{x\}\\tilde\{b\}\_\{\\theta\}\\bigl\(X\(t\)\\bigr\)\\,\\mathrm\{d\}t\\;\-\\;a\(t\)^\{\\top\}\\,\\partial\_\{x\}\\tilde\{\\sigma\}\_\{\\theta\}\\bigl\(X\(t\)\\bigr\)\\circ\\mathrm\{d\}W\(t\)\\;\-\\;\\partial\_\{x\}\\bar\{L\}\\bigl\(X\(t\),\\pi\_\{\\theta,t\}\\bigr\)\\,\\mathrm\{d\}t\\;\}with terminal conditiona⁡\(T\)=Φx​\(X⁡\(T\)\)a\(T\)=\\Phi\_\{x\}\(X\(T\)\), and solved backward along the*same*Brownian path used for the forward solve\. The parameter gradient then becomes

\(5\.14\)∇θℒ=−∫0T\[a\(t\)⊤∂θb~θ\(X\(t\)\)dt\+a\(t\)⊤∂θσ~θ\(X\(t\)\)∘dW\(t\)\+∂θ\(L¯\+λℋ\)\(πθ,t\)dt\]\\nabla\_\{\\theta\}\\mathcal\{L\}\\;=\\;\-\\int\_\{0\}^\{T\}\\\!\\left\[a\(t\)^\{\\top\}\\partial\_\{\\theta\}\\tilde\{b\}\_\{\\theta\}\\bigl\(X\(t\)\\bigr\)\\,\\mathrm\{d\}t\+a\(t\)^\{\\top\}\\partial\_\{\\theta\}\\tilde\{\\sigma\}\_\{\\theta\}\\bigl\(X\(t\)\\bigr\)\\circ\\mathrm\{d\}W\(t\)\+\\partial\_\{\\theta\}\\\!\\left\(\\bar\{L\}\+\\lambda\\mathcal\{H\}\\right\)\\\!\\bigl\(\\pi\_\{\\theta,t\}\\bigr\)\\mathrm\{d\}t\\right\]
We observe three characteristics:

1. 1\.EB\-SSDE is a pathwise backward SDE, not a Pardoux–Peng BSDE\.Like the Li et al\. adjoint, it reuses the forward Brownian path rather than introducing a new martingale\-representation term\. This is a computational advantage: we inherit the Virtual\-Brownian\-Tree trick and𝒪⁡\(log⁡L\)\\mathcal\{O\}\(\\log L\)memory\.
2. 2\.The entropy term appears as a directθ\\theta\-gradient, not as a backward\-SDE source\.Because the entropyℋ⁡\(πθ\)\\mathcal\{H\}\(\\pi\_\{\\theta\}\)depends onθ\\thetaonly through the policy \(not through the state\), it contributes a standard pathwise gradient and does not enter the drift of the adjoint \([5\.13](https://arxiv.org/html/2609.20906#S5.E13)\)\. This decoupling makes our EB\-SSDE a clean generalization of the Li et al\. adjoint: the classical adjoint is recovered in the limitλ→0\\lambda\\to 0with the policy pinched to a delta function at a deterministic action\.
3. 3\.The connection to Wang–Zhou\.The optimality condition∇θℒ=0\\nabla\_\{\\theta\}\\mathcal\{L\}=0, together with standard variational calculus, recovers the Gibbs formπ⋆\(u∣x\)∝exp\(−H¯\(x,u;a\)/λ\)\\pi^\{\\star\}\(u\\mid x\)\\propto\\exp\(\-\\bar\{H\}\(x,u;a\)/\\lambda\)for the optimal exploratory policy, whereH¯\\bar\{H\}is the Hamiltonian built from the learned adjointaa\. This makes the EB\-SSDE a constructive, gradient\-based route to exploratory policies that does not pass through solving an HJB PDE\.

## 6Extensions and Future Directions

Our immediate efforts will focus on three primary objectives: \(i\) solidifying the theoretical foundation by deriving \([5\.13](https://arxiv.org/html/2609.20906#S5.E13)\) for an exploratory linear–quadratic control problem, numerically verifying that Gibbs policies optimized viaθ\\thetaconverge to the closed\-form Gaussian solutions in\[[19](https://arxiv.org/html/2609.20906#bib.bib16)\]; \(ii\) providing a high\-performance implementation of EB\-SSDE within thetorchsdeecosystem, utilizing the Virtual Brownian Tree for memory\-efficient, reproducible path synthesis; and \(iii\) formalizing the astrophysical inference problem\[[6](https://arxiv.org/html/2609.20906#bib.bib4)\]as a reinforcement learning task, where the “control” represents physical perturbations to a latent SDE, optimized for light\-curve reconstruction fidelity\.

Beyond these immediate goals, several promising research trajectories remain:

1. 1\.Non\-Gaussianity and Jump\-Diffusion Processes\.To better capture the heavy\-tailed variability characteristic of certain AGN, we will extend the framework to include Lévy\-driven SDEs or jump\-diffusion models of\[[8](https://arxiv.org/html/2609.20906#bib.bib6)\]\. This necessitates a re\-derivation of the adjoint sensitivity equations, as the introduction of jumps breaks the standard continuity assumptions of the previously defined adjoint SDE in\[[13](https://arxiv.org/html/2609.20906#bib.bib11)\]\.
2. 2\.Theoretical Decision Boundaries for Delay\-Coupling\.The magnitude ofΛ\\Lambdadictates a critical trade\-off between convergence stability and exploration\. A highly promising future direction is to build rigorous mathematical theories that determine the exact decision boundary ofΛ\\Lambda\(whether characterized by its magnitude in norm or its specific distribution\)\. For instance, robust decision boundary in thr form of‖Λ‖<δ0\\\|\\Lambda\\\|<\\delta\_\{0\}, and‖Λ‖≥δ1\\\|\\Lambda\\\|\\geq\\delta\_\{1\}might be find with uncertainty intervals\. Establishing these boundaries helps define the conditions under which the Continuous\-DMSGD algorithm is guaranteed to behave convergently versus when it is prone to divergence\.
3. 3\.Reinforcement Learning Benchmarking\.To move beyond astrophysical toy cases, we will subject EB\-SSDE to standard continuous\-control benchmarks\. This will provide a rigorous evaluation of whether the pathwise backward Stratonovich adjoint can compete with established actor–critic methods in high\-dimensional state spaces\.
4. 4\.Bayesian Uncertainty Quantification\.We plan to treat the drift and diffusion networks \(θ\\thetaandϕ\\phiin \([5\.7](https://arxiv.org/html/2609.20906#S5.E7)\)\) as stochastic variables\. By employing Langevin\-style samplers, we can disentangle the aleatoric uncertainty inherent in the dynamical system from the epistemic uncertainty of the model itself, providing a more robust measure of confidence in physical parameter estimation\.

## Acknowledgments

The authors would like to thank the mathematicianDr\. Farzad Sabzikarfor his invaluable guidance during weekly meetings from2026/3to2026/5\. He helped finding relevant and useful literature and provided key insights into the use of gradient descent with a delayed term\.

## References

- \[1\]\(2018\)Optimization methods for large\-scale machine learning\.SIAM Review60\(2\),pp\. 223–311\.External Links:[Document](https://dx.doi.org/10.1137/16M1080173),[Link](https://doi.org/10.1137/16M1080173),https://doi\.org/10\.1137/16M1080173Cited by:[item 2](https://arxiv.org/html/2609.20906#S1.I1.i2.p1.1)\.
- \[2\]E\. M\. Cackett, K\. Horne, and H\. Winkler\(2007\)Testing thermal reprocessing in active galactic nuclei accretion discs\.Monthly Notices of the Royal Astronomical Society380\(2\),pp\. 669–682\.External Links:[Document](https://dx.doi.org/10.1111/j.1365-2966.2007.12098.x)Cited by:[§4\.2\.2](https://arxiv.org/html/2609.20906#S4.SS2.SSS2.p1.1)\.
- \[3\]Z\. Che, S\. Purushotham, K\. Cho, D\. Sontag, and Y\. Liu\(2018\)Recurrent neural networks for multivariate time series with missing values\.Vol\.8,pp\. 6085\.Cited by:[§2\.2](https://arxiv.org/html/2609.20906#S2.SS2.p2.1)\.
- \[4\]R\. T\. Q\. Chen, Y\. Rubanova, J\. Bettencourt, and D\. Duvenaud\(2018\)Neural ordinary differential equations\.InProceedings of the 32nd International Conference on Neural Information Processing Systems,NIPS’18,Red Hook, NY, USA,pp\. 6572–6583\.Cited by:[§1](https://arxiv.org/html/2609.20906#S1.p12.1)\.
- \[5\]J\. Fagin, J\. H\. Chan, H\. Best, M\. O’Dowd, K\. E\. S\. Ford, M\. J\. Graham, J\. W\. Park, and V\. A\. Villar\(2025\)Joint modeling of quasar variability and accretion disk reprocessing using latent stochastic differential equations\.The Astrophysical Journal988\(1\),pp\. 59\.External Links:[Document](https://dx.doi.org/10.3847/1538-4357/addabc),[Link](https://doi.org/10.3847/1538-4357/addabc)Cited by:[§2\.5](https://arxiv.org/html/2609.20906#S2.SS5.p2.1)\.
- \[6\]J\. Fagin, J\. W\. Park, H\. Best, J\. H\. H\. Chan, K\. E\. S\. Ford, M\. J\. Graham, V\. A\. Villar, S\. Ho, and M\. O’Dowd\(2024\)Latent stochastic differential equations for modeling quasar variability and inferring black hole properties\.The Astrophysical Journal965\(2\),pp\. 104\.External Links:[Document](https://dx.doi.org/10.3847/1538-4357/ad2988),[Link](https://doi.org/10.3847/1538-4357/ad2988)Cited by:[§1](https://arxiv.org/html/2609.20906#S1.p15.1),[§2](https://arxiv.org/html/2609.20906#S2.p1.1),[§6](https://arxiv.org/html/2609.20906#S6.p1.1),[Abstract](https://arxiv.org/html/2609.20906#abstract1.1)\.
- \[7\]Ž\. Ivezić and C\. L\. MacLeod\(2014\)Optical variability of quasars: a damped random walk\.InMultiwavelength AGN Surveys and Studies, Proceedings of the International Astronomical Union, IAU Symposium 304,A\. M\. Mickaelian, D\. B\. Sanders, and I\. S\. McLean \(Eds\.\),Vol\.304,Cambridge, UK,pp\. 131–136\.External Links:[Document](https://dx.doi.org/10.1017/S174392131400494X),[Link](https://arxiv.org/abs/1312.3966)Cited by:[Abstract](https://arxiv.org/html/2609.20906#abstract1.1)\.
- \[8\]J\. Jia and A\. R\. Benson\(2019\)Neural jump stochastic differential equations\.InAdvances in Neural Information Processing Systems,Vol\.32\.Cited by:[item 1](https://arxiv.org/html/2609.20906#S6.I1.i1.p1.1)\.
- \[9\]V\. P\. Kasliwal, M\. S\. Vogeley, and G\. T\. Richards\(2015\)Are the variability properties of the kepler agn light curves consistent with a damped random walk?\.Monthly Notices of the Royal Astronomical Society451\(4\),pp\. 4328–4345\.External Links:ISSN 0035\-8711,[Document](https://dx.doi.org/10.1093/mnras/stv1230),[Link](https://doi.org/10.1093/mnras/stv1230),https://academic\.oup\.com/mnras/article\-pdf/451/4/4328/3890355/stv1230\.pdfCited by:[§1](https://arxiv.org/html/2609.20906#S1.p14.1),[§1](https://arxiv.org/html/2609.20906#S1.p8.1)\.
- \[10\]B\. C\. Kelly, J\. Bechtold, and A\. Siemiginowska\(2009\)ARE the variations in quasar optical flux driven by thermal fluctuations?\.The Astrophysical Journal698\(1\),pp\. 895\.External Links:[Document](https://dx.doi.org/10.1088/0004-637X/698/1/895),[Link](https://doi.org/10.1088/0004-637X/698/1/895)Cited by:[item 1](https://arxiv.org/html/2609.20906#S1.I1.i1.p1.1),[§1](https://arxiv.org/html/2609.20906#S1.p2.1)\.
- \[11\]B\. C\. Kelly, A\. C\. Becker, M\. Sobolewska, A\. Siemiginowska, and P\. Uttley\(2014\)FLEXIBLE and scalable methods for quantifying stochastic variability in the era of massive time\-domain astronomical data sets\.The Astrophysical Journal788\(1\),pp\. 33\.External Links:[Document](https://dx.doi.org/10.1088/0004-637X/788/1/33),[Link](https://doi.org/10.1088/0004-637X/788/1/33)Cited by:[§1](https://arxiv.org/html/2609.20906#S1.p10.1)\.
- \[12\]P\. Kidger, J\. Morrill, J\. Foster, and T\. Lyons\(2020\)Neural controlled differential equations for irregular time series\.InProceedings of the 34th International Conference on Neural Information Processing Systems,NIPS ’20,Red Hook, NY, USA\.External Links:ISBN 9781713829546Cited by:[§3\.2](https://arxiv.org/html/2609.20906#S3.SS2.p8.1)\.
- \[13\]X\. Li, T\. L\. Wong, R\. T\. Q\. Chen, and D\. Duvenaud\(2020\)Scalable gradients for stochastic differential equations\.InProceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics,S\. Chiappa and R\. Calandra \(Eds\.\),Proceedings of Machine Learning Research, Vol\.108,pp\. 3870–3882\.External Links:[Link](https://proceedings.mlr.press/v108/li20i.html)Cited by:[item 3](https://arxiv.org/html/2609.20906#S1.I1.i3.p1.1),[§2\.2](https://arxiv.org/html/2609.20906#S2.SS2.p3.1),[§2\.4](https://arxiv.org/html/2609.20906#S2.SS4.p3.1),[item 2](https://arxiv.org/html/2609.20906#S5.I1.i2.p1.1),[§5\.2](https://arxiv.org/html/2609.20906#S5.SS2.p2.1),[§5\.3](https://arxiv.org/html/2609.20906#S5.SS3.p1.2),[§5\.3](https://arxiv.org/html/2609.20906#S5.SS3.p1.4),[§5](https://arxiv.org/html/2609.20906#S5.p1.1),[item 1](https://arxiv.org/html/2609.20906#S6.I1.i1.p1.1)\.
- \[14\]X\. Mao\(2007\)Stochastic differential equations and applications\.Woodhead Publishing\.External Links:ISBN 9780857099402,[Link](https://books.google.com/books?id=l5ejAgAAQBAJ)Cited by:[§3\.2](https://arxiv.org/html/2609.20906#S3.SS2.p6.1)\.
- \[15\]R\. Munos\(2006\)Policy gradient in continuous time\.Journal of Machine Learning Research7\(27\),pp\. 771–791\.External Links:[Link](http://jmlr.org/papers/v7/munos06b.html)Cited by:[§5\.3](https://arxiv.org/html/2609.20906#S5.SS3.p1.1),[§5](https://arxiv.org/html/2609.20906#S5.p1.1)\.
- \[16\]R\. F\. Mushotzky, R\. Edelson, W\. Baumgartner, and P\. Gandhi\(2011\)KEPLER observations of rapid optical variability in active galactic nuclei\.The Astrophysical Journal Letters743\(1\),pp\. L12\.External Links:[Document](https://dx.doi.org/10.1088/2041-8205/743/1/L12),[Link](https://doi.org/10.1088/2041-8205/743/1/L12)Cited by:[§1](https://arxiv.org/html/2609.20906#S1.p14.1),[§1](https://arxiv.org/html/2609.20906#S1.p8.1)\.
- \[17\]Y\. Oh, S\. Kam, D\. Lim, and S\. Kim\(2025\)Modeling irregular astronomical time series with neural stochastic delay differential equations\.InProceedings of the 34th ACM International Conference on Information and Knowledge Management,CIKM ’25,New York, NY, USA,pp\. 5068–5073\.External Links:ISBN 9798400720406,[Link](https://doi.org/10.1145/3746252.3760805),[Document](https://dx.doi.org/10.1145/3746252.3760805)Cited by:[§1](https://arxiv.org/html/2609.20906#S1.p15.1),[§3\.2](https://arxiv.org/html/2609.20906#S3.SS2.p3.1),[§4\.1\.1](https://arxiv.org/html/2609.20906#S4.SS1.SSS1.p1.1),[Abstract](https://arxiv.org/html/2609.20906#abstract1.1)\.
- \[18\]M\. Song, L\. Hu, X\. Mao, and L\. Zhang\(2013\)Khasminskii\-type theorems for stochastic functional differential equations\.Discrete and Continuous Dynamical Systems \- B18\(6\),pp\. 1697–1714\.External Links:[Document](https://dx.doi.org/10.3934/dcdsb.2013.18.1697)Cited by:[§4\.3\.1](https://arxiv.org/html/2609.20906#S4.SS3.SSS1.p18.1),[§4\.3\.2](https://arxiv.org/html/2609.20906#S4.SS3.SSS2.p12.1),[§4\.3](https://arxiv.org/html/2609.20906#S4.SS3.p1.1)\.
- \[19\]H\. Wang, T\. Zariphopoulou, and X\. Y\. Zhou\(2020\)Reinforcement learning in continuous time and space: a stochastic control approach\.Journal of Machine Learning Research21\(198\),pp\. 1–34\.External Links:[Link](http://jmlr.org/papers/v21/19-144.html)Cited by:[§5\.4](https://arxiv.org/html/2609.20906#S5.SS4.p1.1),[§5](https://arxiv.org/html/2609.20906#S5.p1.1),[§6](https://arxiv.org/html/2609.20906#S6.p1.1),[Abstract](https://arxiv.org/html/2609.20906#abstract1.1)\.