Affine constraints in non-reversible diffusions with degenerate noise

TL;DR

This work develops a framework for enforcing affine constraints on nonreversible diffusions with degenerate noise via stiff confinement, showing that the softly constrained OU dynamics converge as $\varepsilon\to0$ to a projected SDE on the constraint surface, explicitly given by $dY_t = P M Y_t dt + \sqrt{2} P C dW_t$ with projection $P$ determined by the confinement matrix $K$ through $P = I - \tilde{K}\tilde{K}^{\dag}$. It provides an ambient-space (projection-based) formulation that can realize orthogonal or oblique projections, and derives conditions under which the constrained process preserves conditional or invariant measures, including cases like $K = A-J$ or $K = \Sigma$. The results extend to nonlinear drifts and include a quantitative pathwise convergence rate, initial-layer analysis, and implications for sampling conditioned Gaussian measures and solving constrained stochastic optimization, complemented by numerical demonstrations on constrained OU and discretised elliptic operators. This work advances understanding of constrained nonreversible diffusions, offering practical guidelines for designing stiff confinement to achieve desired long-time behavior and measure-preservation properties in high-dimensional settings.

Abstract

This paper deals with the realisation of affine constraints on nonreversible stochastic differential equations (SDE) by strong confining forces. We prove that the confined dynamics converges pathwise and on bounded time intervals to the solution of a projected SDE in the limit of infinitely strong confinement, where the projection is explicitly given and depends on the choice of the confinement force. We present results for linear Ornstein-Uhlenbeck (OU) processes, but they straightforwardly generalise to nonlinear SDEs. Moreover, for linear OU processes that admit a unique invariant measure, we discuss conditions under which the limit also preserves the long-term properties of the SDE. More precisely, we discuss choices for the design of the confinement force which in the limit yield a projected dynamics with invariant measure that agrees with the conditional invariant measure of the unconstrained processes for the given constraint. The theoretical findings are illustrated with suitable numerical examples.

Theorems & Definitions (5)

• Theorem 3.1
• Remark 3.2
• Lemma 3.3: campbell1979singular
• proof : Proof of Theorem \ref{['thm:genOU-SC']}
• Remark 3.4: Initial conditions and convergence in $C([0,\infty))$ ] Theorem \ref{['thm:genOU-SC']} provides the convergence $X^\varepsilon_t \to Y_t$ as $\varepsilon \to 0$ for any $t>0$. This convergence statement is in fact wrong for $t=0$ since $\lim_{\varepsilon\to 0}X^{\varepsilon}_0 = Y_0$ if and only if $X_0 \in \xi^{-1}(0),$ where $Y_0 \coloneqq P ( X_0 - \left( b,\, 0 \right)^T ) + \left( b,\, 0 \right)^T$. Katzenberger Katzenberger91 addresses this issue by studying the convergence of a modified process $Y^\varepsilon_t \coloneqq X^\varepsilon_t - \psi(X^\varepsilon_0,\frac{t}{\varepsilon}) + \theta(X^\varepsilon_0)$ where $\psi$ defined in \ref{['eq:katz-ODE']} is the ODE flow corresponding to the stiff part of the SDE and $\theta$ is its long-time limit \ref{['eq:theta']}. By construction $\psi(X^\varepsilon_0,0)=X^\varepsilon_0$ and therefore at $t=0$, $Y^{\varepsilon}_0 = \theta(X^\varepsilon_0) \in \xi^{-1}(0)$. Using this modified process, uniform convergence on bounded time intervals in $C([0,\infty))$ follows, which leads to pathwise convergence in probability. We prefer to work with the original soft constrained process \ref{['eq:genOU-SC']} and study its limit. Of course, we could work with an analogue construction here and get the same kind of convergence as in Katzenberger91. Theorem \ref{['thm:pathwiseconvrate']} below implies the convergence result of Katzenberger for the unmodified process, if the limiting initial datum lies on the constraint manifold, i.e. $X_0 \in \xi^{-1}(0)$. The proof for Theorem \ref{['thm:genOU-SC']} can be easily generalised to the class of nonlinear SDEs where the drift can be written as a sum of a linear function and a smooth, bounded, nonlinear perturbation, see \ref{['eq:SC-Nonlin-gen']} below. This implies that the drift is Lipschitz which is typically required for well-posedness of strong solutions for SDEs. The major difficulty here is that we require compactness of the sequence $(X^\varepsilon)_{\varepsilon>0}\in C([0,\infty))$, which can be extracted via standard approaches (see Katzenberger91 for references). More precisely, let $X^\varepsilon_t$ solve dX^\varepsilon_t = MX^\varepsilon_t + f(X^\varepsilon_t) - \frac{1}{\varepsilon} \tilde{K} X^\varepsilon_t dt + \sqrt{2} C dW_t where $f\colon \mathbb{R}^d \to \mathbb{R}^d$ is smooth and uniformly bounded and the other coefficients are as before. Here the limiting dynamics $Y_t$ solves dY_t = PMY_t + Pf(Y_t)dt + \sqrt{2} PC dW_t, where $P=\left(00\alphaI \right)$ as before. By variation of constants (see Proposition \ref{['prop:varofconst']}), $X_t^\varepsilon$ admits the solution X^\varepsilon_t = e^{(M-\frac{1}{\varepsilon} \tilde{K})t}X_0^\varepsilon + \int_0^t e^{(M-\frac{1}{\varepsilon} \tilde{K})(t-s)} f(X^\varepsilon_s) ds + \sqrt{2}\int_0^t e^{(M-\frac{1}{\varepsilon} \tilde{K})(t-s)} C dW_s, and the convergence of the first and last term above to the corresponding terms in $Y_t$ follows as in the proof of Theorem \ref{['thm:genOU-SC']}. Assuming that $\left(X^\varepsilon_s\right)$ is compact in $C([0,t])$, i.e. converges up to subsequences, by dominated convergence we expect \lim\limits_{\varepsilon \to 0} \int_0^t e^{(M-\frac{1}{\varepsilon} \tilde{K})(t-s)} f(X^\varepsilon_s) ds = \int_0^t \lim\limits_{\varepsilon \to 0} e^{(M-\frac{1}{\varepsilon} \tilde{K})(t-s)} f(X^\varepsilon_s) ds = \int_0^t P f(Y_s) ds, and therefore $X^\varepsilon_t \to Y_t$ in probability as in the proof above. Consequently, Theorem \ref{['thm:genOU-SC']} generalises to a considerably larger class of SDEs. However, we do not provide a complete proof as this requires us to discuss technical results regarding compactness of the sequence $\left(X^\varepsilon\right)$ which are outside the scope of this work; see Katzenberger91 and KurtzProtter91 for details. We introduce a quantitative convergence result (see Theorem \ref{['thm:pathwiseconvrate']} below) which generalises Theorem \ref{['thm:genOU-SC']} by employing the following crucial lemma along with a Gronwall inequality argument. Consider the matrix $\tilde{K}=K_{11}0K_{21}0\in \mathbb{R}^{d\times d}$ where $K_{11} \in \mathbb{R}^{k \times k}$ is invertible and $K_{21} \in \mathbb{R}^{(d-k) \times k}$. Then e^{-\tilde{K}t} = e^{- K_{11} t}0-K_{21} K_{11}^{-1}+ K_{21} e^{- K_{11} t} K_{11}^{-1}I.If $-K_{11}$ is Hurwitz then e^{-\frac{1}{\varepsilon}\tilde{K}t}\xrightarrow{\varepsilon\to 0} P := 00- K_{21}K_{11}^{-1}I.Assume that $-K_{11}$ is Hurwitz. Let $V \in \mathbb{C}^{k \times k}$ be the matrix that transforms $K_{11}$ to its Jordan normal form, i.e. $V^{-1} K_{11} V = \Lambda + N$ where $\Lambda$ is a diagonal matrix with the eigenvalues of $K_{11}$ as entries and $N$ is an upper triangular nilpotent matrix of order $m \in \mathbb{N}$, i.e. $N^m=0$. Then for any $t>0$ \bigl\| e^{-\frac{1}{\varepsilon}\tilde{K} t} - P\bigr\|_F^2 \leq c_A e^{-\frac{2\lambda_1t}{\varepsilon} } \sum\limits_{j=0}^{m-1} t^j\frac{\|N^j\|_F^2}{\varepsilon^j } where $\lambda_1$ is the smallest real part of all eigenvalues of $K_{11}$ and $c_A= \Bigl( 1 + \bigl\| K_{21}\bigr\|_F^2 \bigl\| K^{-1}_{11}\bigr\|_F^2 \Bigr) \kappa(V) mk.$ Here $\kappa(V) \coloneqq \|V\|_F^2 \|V^{-1}\|_F^2$ is the condition number of the matrix $V$. Moreover, \int_0^t \| e^{-\frac{1}{\varepsilon}\tilde{K}(t-s)} - P \|_F^2 ds\leq \varepsilon c_P, where $c_P= \Bigl( 1 + \bigl\| K_{21}\bigr\|_F^2 \bigl\| K^{-1}_{11}\bigr\|_F^2 \Bigr) \kappa(V) k m \, \sum\limits_{j=0}^{m-1}\frac{\|N^j\|_F^2}{(2\lambda_1)^{j+1}}.$ These estimates simplify in two particular cases: If $K_{11}$ is non-defective then \ref{['eq:Frobeniusest_expK-P']} becomes \bigl\| e^{-\frac{1}{\varepsilon}\tilde{K} t} - P\bigr\|_F^2 \leq c_A e^{-\frac{2\lambda_1t}{\varepsilon} }, with $c_A=\Bigl( 1 + \bigl\| K_{21}\bigr\|_F^2 \bigl\| K^{-1}_{11}\bigr\|_F^2 \Bigr) \kappa(V) k ,$ and estimate \ref{['eq:matrixexp-int-estimate']} stays unchanged with $c_P$ is given by c_P = \left(1 + \bigl\| K_{21}\bigr\|_F^2 \bigl\| K^{-1}_{11}\bigr\|_F^2\right) \kappa(V) k \,.If $K_{11}$ is symmetric and non-defective, then the estimate above apply with $\kappa(V) \equiv 1$ in the constants. Let $\tilde{K}$ be defined as in Lemma \ref{['lem:matrixexp']}, $-K_{11} \in \mathbb{R}^{k \times k}$ be Hurwitz and defective. Then for any $\delta \in (0,\lambda_1)$ \bigl\|e^{-\frac{1}{\varepsilon}\tilde{K} t} - P\bigr\|_F^2 \leq c e^{-2\frac{\lambda_1 - \delta}{\varepsilon} t } where $c=c(\delta) = (m-1)^m \left(\frac{\|N\|^2_F}{2\delta e}\right)^{m-1}$ ,$\lambda_1$ is the smallest real part of all eigenvalues of $K_{11}$ and $\kappa(V)$, $m$ are as in Lemma \ref{['lem:matrixexp']} above. The proof of part \ref{['lem:matrixexpstatement']} above follows on the lines of SharmaZhang21. All proofs are presented in Appendix \ref{['app:MatrixExp']}. We now present the quantitative result of the pathwise convergence. Given $\varepsilon>0$, let $X^\varepsilon_t$ solve \ref{['eq:genOU-SC']} with initial datum $X^\varepsilon_0$ and $Y_t$ solve \ref{['eq:genOU-SC-limit']} with initial datum $Y_0=P (X_0 -\left( b,\, 0 \right)^T ) +\left( b,\, 0 \right)^T$. We have \mathbb{E}\biggl[\sup_{t\in [0,T]} \bigl|X_t^\varepsilon - Y_t\bigr|^2 \biggr] \leq c_1\biggl( \mathbb{E}\bigl[|X^{1,\varepsilon}_0-b\bigr|^2] + \mathbb{E}\bigl[|X^\varepsilon_0-X_0|^2\bigr] + \varepsilon\bigl[1+ e^{\lambda_{max}(PM)T} \bigr]\biggr) e^{c_2 T^2} where $c_1,c_2$ are independent of $\varepsilon$ and $T$, and $\lambda_{\max}(PM)\geq 0$ is the principal eigenvalue of $PM$. In particular, if the limiting initial datum is well-prepared, i.e. $X_0\in \xi^{-1}(0)$ \mathbb{E}\biggl[\sup_{t\in [0,T]} \bigl|X_t^\varepsilon - Y_t\bigr|^2 \biggr] \leq c_1\biggl(\mathbb{E}\bigl[|X_0^\varepsilon - X_0|^2\bigr] + \varepsilon\bigl[1+ e^{\lambda_{max}(PM)T} \bigr]\biggr) e^{c_2 T^2} and hence if $\mathbb{E}\bigl[|X_0^\varepsilon - X_0|^2\bigr] \to 0$ as $\varepsilon\to 0$ then $X^\varepsilon\xrightarrow{\varepsilon\to 0} Y \text{ in probability on } C([0,T];\mathbb{R}^d).$Let $K_{11}$ be non-defective. For any $t>0, \ X_0 \in \mathbb{R}^d$ we have \mathbb{E}\Bigl[ \bigl|X^\varepsilon_t - Y_t\bigr|^2 \Bigr] \leq c_1\biggl( e^{-\frac{2\lambda_1 t}{\varepsilon}} \mathbb{E}\bigl[ |X_0^{1} - b|^2\bigr] + 2\mathbb{E}\bigl[ | X_0^\varepsilon - X_0|^2 \bigr] + \varepsilon\bigl[1+ e^{\lambda_{max}(PM)t} \bigr] \biggr) e^{c_2 t^2}. Let $K_{11}$ be defective. For any $t>0, \ X_0 \in \mathbb{R}^d$ and any $\delta \in (0,\lambda_1)$ we have \mathbb{E}\Bigl[ \bigl|X^\varepsilon_t - Y_t\bigr|^2 \Bigr] \leq \tilde{c}_1\biggl( e^{-\frac{2(\lambda_1-\delta )t}{\varepsilon}} \mathbb{E}\bigl[ |X_0^{1} - b|^2\bigr] + 2\mathbb{E}\bigl[ | X_0^\varepsilon - X_0|^2 \bigr] + \varepsilon\bigl[1+ e^{\lambda_{max}(PM)t} \bigr] \biggr) e^{c_2 t^2}, where $c_1, \tilde{c}_1, c_2>0$ are independent of $\varepsilon$ and $t$, and $\lambda_1>0$ is the smallest real part of eigenvalues of $K_{11}$. In particular, with initial data $\mathbb{E}\bigl[|X_0^\varepsilon - X_0|^2\bigr] \to 0$, where $X_0 \in \mathbb{R}^d$ is arbitrary, for any $t>0$ we have $X_t^\varepsilon\xrightarrow{\varepsilon\to 0} Y_t \text{ in probability}.$ Before presenting the proof, a remark is in order to discuss: (a) role of initial conditions, (b) role of $\lambda_{\max}(PM)$, and (c) related estimates for Fokker-Planck equations. The key difference between the pathwise bound \ref{['eq:PathBound']} and the pointwise-in-time bound \ref{['eq:PointBound']} is the treatment of the initial condition $|X^{1,\varepsilon}_0-b|^2$ and the fact that while the pathwise bound compares the entire trajectory (including $t=0$) the pointwise bound only applies to $t>0$. Due to this key difference the pointwise bound does not see the initial boundary layer (characterized via the exponential decay in the first term) while the pathwise bound does. As a consequence, if $X^\varepsilon_0 \to X_0 \notin \xi^{-1}(0)$, i.e. $X^1_0\neq b$, then the pathwise estimate \ref{['eq:PathBound']} does not vanish, but the pointwise estimate \ref{['eq:PointBound']} does. Furthermore, if $X^{1,\varepsilon}_0\to b$ at a rate $f(\varepsilon)$ and $X^{\varepsilon}_0\to X_0$ at rate $g(\varepsilon)$, then $f(\varepsilon)$, $g(\varepsilon)$ determine the rate of convergence in \ref{['eq:PathBound']} as $\varepsilon\to 0$, i.e. we have $\mathbb{E}\biggl[\sup_{t\in [0,T]}\bigl|X_t^\varepsilon -Y_t\bigr|^2 \biggr] \leq C\min\bigl\{f(\varepsilon),g(\varepsilon),\varepsilon\bigr\},$ where $C$ is independent of $\varepsilon$. Next we discuss the role of $\lambda_{\max}(PM)$ that appears in the estimates above. In the proof of the quantitative estimate we need to control $\mathbb{E}[|Y_t|^2]$ (see discussion on $I_{2}$ in the proof below), where $Y_t$ is the limiting dynamics \ref{['eq:genOU-SC-limit']}. The kernel of the matrix product $PM \in \mathbb{R}^{d \times d}$ is at least a $k$-dimensional, by definition of the projection in \ref{['def:genOU-Proj']} and hence $PM$ has $k$ zero eigenvalues. The other eigenvalues of $PM$ agree with the eigenvalues of $(PM)_{22}$ since for any eigenvector $v \in \mathbb{R}^{d-k}$ of $(PM)_{22}$, i.e. $(PM)_{22} v = \lambda v$, we have $PM (0,v)^T = \lambda (0,v)^T$. Consequently, there are two possible cases. First, in the case where $(PM)_{22}$ is Hurwitz, which corresponds to $Y^2$ in the limiting dynamics \ref{['eq:genOU-SC-limit']} admitting an invariant measure (see Section \ref{['sec:SteadyState-K']} for a detailed study of this setting), it follows that $\lambda_{\max}(PM)=0$ and we have $e^{\lambda_{max}(PM)t}=1$ in the bounds above. Second, if $(PM)_{22}$ is not Hurwitz, we have an additional term that grows exponentially in time with rate given by $e^{\lambda_{max}(PM)t}$. It should be noted that similar exponential growth estimates also arise when using alternative methods for controlling $\mathbb{E}[|Y_s|^2]$, for instance see the approach via Gronwall's inequality employed in Remark \ref{['rem:QuantNonlin']}. Note that the defective or non-defective nature of $K_{11}$ only plays a role in the pointiwse estimates. In particular, this determines the exponential rate of decay when the initial conditions do not satisfy the constraint. As indicated by the estimates, the case of non-defective $K_{11}$ leads to better convergence rate of the initial data. Let us mention that these quantitative estimates generalise to the case of nonlinear SDEs with Lipschitz drift, see Remark \ref{['rem:QuantNonlin']} for a discussion. Finally, we point out that the quantitative pathwise estimates in Theorem \ref{['thm:pathwiseconvrate']} also provide estimates on the corresponding Fokker-Planck equations. In particular, using $\nu^\varepsilon_t = \mathrm{law}(X_t^\varepsilon)$ and $\rho_t = \mathrm{law}(Y_t)$, a quantitative estimate of the Wasserstein-2 distance $\mathcal{W}_2(\nu_t,\rho_t)$ follows since $\bigl(\mathcal{W}_2 (\nu^\varepsilon_t , \rho_t)\bigr)^2 \leq \mathbb{E}\bigl[|X_t^\varepsilon-Y_t|^2\bigr] \leq \mathbb{E}\biggl[\sup_{t\in [0,T]}|X_t^\varepsilon-Y_t|^2\biggr].$ \ref{['item:PathQuant']} First consider the case $\xi(x) = x^1$, i.e. $b=0$. The solutions to \ref{['eq:genOU-SC']} and \ref{['eq:genOU-SC-limit']} read X^\varepsilon_t= e^{ - \frac{1}{\varepsilon}\tilde{K} t} X^\varepsilon_0 + \int_0^t e^{-\frac{1}{\varepsilon} \tilde{K}(t-s)} MX^\varepsilon_s ds + \int_0^t e^{-\frac{1}{\varepsilon}\tilde{K}(t-s)} C dW_s \ \text{ and }Y_t= Y_0 + \int_0^t PMY_s ds + \int_0^t PC dW_s\,, where we have employed variation of constants \ref{['eq:VarofCon-sol']} to arrive at the solution of $X^\varepsilon_t$ and $Y_t$. Using Young's inequality we find $\mathbb{E}\biggl[ \sup\limits_{t \in [0,T]} \bigl|X^\varepsilon_t - Y_t\bigr|^2 \biggr]\leq 3\mathbb{E}\biggl[ \sup\limits_{t \in [0,T]} \bigl|e^{-\frac{1}{\varepsilon}\tilde{K} t} X^\varepsilon_0 - Y_0\bigr|^2 \biggr] + 3 \mathbb{E}\biggl[ \sup\limits_{t \in [0,T]} \biggl|\int_0^t \Bigl( e^{-\frac{1}{\varepsilon}\tilde{K}(t-s)}MX^\varepsilon_s - PMY_s\Bigr) ds \biggr|^2 \biggr]\ +3\mathbb{E}\biggl[ \sup\limits_{t \in [0,T]} \biggl|\int_0^t \Bigl( e^{-\frac{1}{\varepsilon} \tilde{K}(t-s)}C - PC \Bigr)dW_s\biggr|^2 \biggr] \eqqcolon 3 (I_1 + I_2 + I_3) \,.$ We treat each term separately and start with the last one. Applying Doob's inequality followed by the Itô isometry, then using the sub-multiplicativity of the Frobenius norm and lastly the bound \ref{['eq:matrixexp-int-estimate']} we find I_3\leq 4 \mathbb{E} \biggl[\biggl| \int_0^T \Bigl( e^{- \frac{1}{\varepsilon} \tilde{K} (T-s) } - P \Bigr) C dW_s \biggr|^2 \biggr] = 4 \int_0^T \Bigl\|\bigl(e^{- \frac{1}{\varepsilon} \tilde{K} (T-s) } - P\bigr)C \Bigr\|^2_F ds\leq 4 \int_0^T \bigl\|e^{- \frac{1}{\varepsilon} \tilde{K} (T-s) } - P \bigr\|_F^2 \|C \|^2_F ds \leq 4 \varepsilon c_P \|C\|_F^2. Next we consider $I_1$. Recall that for $b=0$ we have $Y_0 = PX_0$. Adding a zero and applying Young's inequality in the first step we have $\bigl|e^{- \frac{1}{\varepsilon} \tilde{K} t}X_0^\varepsilon - Y_0\bigr|^2 \leq 2\bigl|(e^{- \frac{1}{\varepsilon} \tilde{K} t} - P)X_0^\varepsilon \bigr|^2 + 2\bigl| PX_0^\varepsilon - Y_0\bigr|^2 \leq 2\bigl\|e^{- \frac{1}{\varepsilon} \tilde{K} t} - P\bigr\|_F^2 \bigl|X^{1,\varepsilon}_0\bigr|^2 + 2|PX_0^\varepsilon - PX_0|^2$ where the second inequality follows by the sub-multiplicativity of the norm and using \ref{['eq:matrixexp']} which gives (e^{- \frac{1}{\varepsilon} \tilde{K} t} - P)X_0^\varepsilon = e^{-\frac{1}{\varepsilon}K_{11}t}0K_{21}e^{-\frac{1}{\varepsilon}K_{11}t}K_{11}^{-1}0X^{1,\varepsilon}_0X^{2,\varepsilon}_0 = (e^{- \frac{1}{\varepsilon} \tilde{K} t} - P) X^{1,\varepsilon}_00. As a consequence, using $|PX_0^\varepsilon - PX_0|^2 \leq \|P\|_F^2 |X_0^\varepsilon - X_0|^2$ and Corollary \ref{['cor:expK-P']} or \ref{['eq:Frobeniusest_expK-Psymm']} depending on whether $K_{11}$ is defective or not, we can bound $I_1$ from above by I_1 \leq c \left( \mathbb{E}(|X_0^{1,\varepsilon}|^2 + \mathbb{E}(|X_0^\varepsilon - X_0|^2) \right) where $c>0$ is independent of $T$ and $\varepsilon$. If additionally $X_0 \in \xi^{-1}(0),$ i.e. $X_0^1=0$, so that $|X_0^{1,\varepsilon}|^2 = |X_0^{1,\varepsilon} - X_0^1|^2 \leq |X_0^{\varepsilon} - X_0|^2$ and we have the overall bound $I_1 \leq c \mathbb{E}\bigl[|X_0^\varepsilon - X_0|^2\bigr],$ where again $c>0$ is independent of $T$ and $\varepsilon$. For $I_2$ we first add a zero and calculate using Young's inequality and in the second step the Cauchy-Schwarz inequality, I_2= \mathbb{E}\biggl[ \sup\limits_{t \in [0,T]} \biggl|\int_0^t \Bigl( e^{-\frac{1}{\varepsilon}\tilde{K}(t-s)}M (X^\varepsilon_s - Y_s) + (e^{-\frac{1}{\varepsilon}\tilde{K} (t-s)} - P)MY_s \Bigr)ds \biggr|^2 \biggr]\leq 2\mathbb{E}\biggl[ \sup\limits_{t \in [0,T]} \biggl|\int_0^t e^{-\frac{1}{\varepsilon}\tilde{K}(t-s)}M (X^\varepsilon_s - Y_s) ds\biggr|^2\biggr] + 2\mathbb{E}\biggl[ \sup\limits_{t \in [0,T]} \biggl|\int_0^t (e^{-\frac{1}{\varepsilon}\tilde{K} (t-s)} - P)MY_s ds \biggr|^2 \biggr]\leq 2\mathbb{E}\biggl[ \sup\limits_{t \in [0,T]} t \int_0^t \biggl|e^{-\frac{1}{\varepsilon}\tilde{K}(t-s)}M (X^\varepsilon_s - Y_s) \biggr|^2 ds\biggr]\qquad\qquad\qquad+ 2\mathbb{E}\biggl[ \sup\limits_{t \in [0,T]} \left( \int_0^t \bigl\| (e^{-\frac{1}{\varepsilon}\tilde{K} (t-s)} - P)M \bigr\|_F^2 ds \right) \left(\int_0^t |Y_s|^2 ds \right)\biggr]\leq 2\|M\|_F^2 \mathbb{E}\biggl[ \sup\limits_{t \in [0,T]} t \int_0^t \bigl\|e^{-\frac{1}{\varepsilon}\tilde{K}(t-s)}\bigr\|^2_F \bigl|X^\varepsilon_s - Y_s\bigr|^2 ds \biggr]\qquad\qquad\qquad+ 2\|M\|_F^2 \mathbb{E}\biggl[ \sup\limits_{t \in [0,T]} \biggl(\int_0^t \bigl\|e^{-\frac{1}{\varepsilon}\tilde{K}(t-s)}-P\bigr\|_F^2 ds \biggr) \int_0^t |Y_s|^2ds \biggr]\eqqcolon 2\|M\|_F^2 \bigl(I_{2,1} + I_{2,2}\bigr)\,. Let us consider the terms separately. For $I_{2,1}$ we first add a zero and apply Young's inequality, using Fubini's theorem in the third step we find $I_{2,1}\leq \mathbb{E}\biggl[ \sup\limits_{t \in [0,T]} 2t \int_0^t \Bigl[ \bigl\| e^{-\frac{1}{\varepsilon}\tilde{K}(t-s)} - P\bigr\|^2_F + \|P\|^2_F \Bigr] |X^\varepsilon_s - Y_s|^2 ds \biggr]\leq 2T\mathbb{E}\biggl[ \int_0^T \bigl[ \bigl\| e^{-\frac{1}{\varepsilon}\tilde{K}(T-s)} - P\bigr\|^2_F + \|P\|^2_F \bigr] |X^\varepsilon_s - Y_s|^2 ds\biggr]= 2T \int_0^T \bigl[ \bigl\| e^{-\frac{1}{\varepsilon}\tilde{K}(T-s)} - P\bigr\|^2_F + \|P\|^2_F \bigr] \mathbb{E}\biggl[ |X^\varepsilon_s - Y_s|^2 \biggr] ds\leq 2 T\int_0^T \bigl[ \bigl\| e^{-\frac{1}{\varepsilon}\tilde{K}(T-s)} - P\bigr\|^2_F + \|P\|^2_F \bigr] \mathbb{E}\biggl[ \sup_{\tau\in [0,s]}|X^\varepsilon_\tau - Y_\tau|^2 \biggr] ds.$ For $I_{2,2}$ we also use \ref{['eq:matrixexp-int-estimate']} in the first step and Fubini's theorem to compute I_{2,2}\leq \varepsilon c_P \mathbb{E}\biggl[ \sup\limits_{t \in [0,T]} \int_0^t |Y_s|^2ds \biggr] = \varepsilon c_P \int_0^T \mathbb{E} \bigl[ |Y_s|^2 \bigr] ds\,. Since $Y_t$ is an OU process we know that $m_s,\Sigma_s$ are given in \ref{['linSDE-NormalSolution']}, and therefore \mathbb{E}\bigl[|Y_s|^2\bigr]= |m_s|^2 + \mathrm{Tr}(\Sigma_s)= \bigl|e^{PMs}\mathbb{E}\bigl[Y_0\bigr]\bigr|^2 + \mathrm{Tr}\bigl(e^{PMs} \Sigma_0 e^{M^TP^Ts}\bigr) + \int_0^s \mathrm{Tr}\bigl(e^{PM(s-r)} PCC^TP^T e^{M^TP^T(s-r)}\bigr) dr\leq \bigl\|e^{PMs}\bigr\|_F^2 |\mathbb{E}\bigl[Y_0\bigr]|^2 + \|e^{PMs}\|_F^2 \| \Sigma_0^{\frac{1}{2}} \|^2_F + \int_0^s \|e^{PM(s-r)}\|_F^2 \| PC\|_F^2 dr\,. Let $S$ be the matrix that transforms $PM$ into its Jordan normal form, i.e. $SPMS^{-1} = \Lambda + N,$ where $\Lambda$ is a diagonal matrix with the eigenvalues of $PM$ as entries and $N$ is an upper triangular nilpotent matrix. Then \|e^{PMs}\|_F^2= \| S^{-1}e^{\Lambda s} e^{N s}S\|_F^2 \leq \| S^{-1} \|_F^2 \|e^{\Lambda s}\|_F^2 \|e^{N s}\|_F^2 \|S\|_F^2 \leq d \kappa(S) \|e^{N s}\|_F^2 e^{2\lambda_{\max(PM)}s}\,, where $\kappa(S)$ is the condition number of $S$. Hence, by a similar calculation as \ref{['eq:calcintexptimespolynomial']}, we find \int_0^T \|e^{PMs}\|_F^2 ds \leq \hat{c} e^{\lambda_{\max}(PM)T}. for some constant $\hat{c}>0$ independent of $\varepsilon$ and $T$. This means, that overall $I_{2,2} \leq \varepsilon c e^{\lambda_{max}(PM) T}$, where $c>0$ is independent of $T$ and $\varepsilon$ but depends on $\mathbb{E}[Y_0]$, $\Sigma_0$ (and consequently $\mathbb{E}[X_0]$ by the definition of $Y_0$) through the first term in \ref{['eq:2ndmoment-Y']}. Substituting these bounds into \ref{['eq:PathEstSep']}, we find \mathbb{E}\biggl[\sup_{t\in [0,T]} \bigl|X_t^\varepsilon - Y_t\bigr|^2 \biggr]\leq \gamma(\varepsilon,T) + c_2 T \int_0^T\bigl[ \bigl\| e^{-\frac{1}{\varepsilon}\tilde{K}(T-s)} - P\bigr\|^2_F + \|P\|^2_F \bigr] \mathbb{E}\biggl[\sup_{t\in [0,s]} \bigl|X_t^\varepsilon - Y_t\bigr|^2 \biggr] ds , where $\gamma(\varepsilon,t) = c_1 (\mathbb{E}(|X_0^{1,\varepsilon}|^2 + \mathbb{E}(|X_0^{\varepsilon} - X_0|^2 + \varepsilon (1+e^{\lambda_{max}(PM) T}))$ for $X_0 \in \mathbb{R}^d$ and $\gamma(\varepsilon,t) = c_1 ( \mathbb{E}(|X_0^{\varepsilon} - X_0|^2) + \varepsilon (1+e^{\lambda_{max}(PM) T}))$ for $X_0 \in \xi^{-1}(0)$. Thus, by Gronwall's inequality and \ref{['eq:matrixexp-int-estimate']} \mathbb{E}\biggl[\sup_{t\in [0,T]} \bigl|X_t^\varepsilon - Y_t\bigr|^2 \biggr] \leq \gamma(\varepsilon,T) e^{c_3 T^2 + \varepsilon c_P c_2 T}. Next we discuss the case $b\neq 0$, i.e. $\xi(x) = x^1 - b$. With the coordinate-shifted variables $\bar{X}_t = X_t - \left(b,\,0 \right)^T$ and $\bar{Y}_t = Y_t + \left(b,\, 0\right)^T$, we are back to the previous case of a coordinate projection onto zero and note that $\mathbb{E}\biggl[\sup\limits_{t \in [0,T]} \bigl|\bar{X}_t - \bar{Y}_t\bigr|^2\biggr]= \mathbb{E}\biggl[\sup\limits_{t \in [0,T]} \bigl|X_t - Y_t\bigr|^2\biggr]\,.$ The proof goes through exactly as above, in particular the estimates for $I_1, I_2, I_3$ remain unchanged. Let us briefly discuss $I_1$, since this is the term where the shifted initial datum appears. First note that \ref{['eq:Quant-I1-0']} becomes (after adding a zero and using that $\bar{Y}_0 = P \bar{X}_0$) \bigl|e^{- \frac{1}{\varepsilon} \tilde{K} t}\bar{X}_0^\varepsilon - \bar{Y}_0\bigr|^2= \bigl|e^{- \frac{1}{\varepsilon} \tilde{K} t}\bar{X}_0^\varepsilon - P \bar{X}_0^\varepsilon + P \bar{X}_0^\varepsilon - \bar{Y}_0\bigr|^2 \leq 2\bigl| \left(e^{- \frac{1}{\varepsilon} \tilde{K} t} - P \right) \left(X_0^\varepsilon -\left(b0\right) \right) \bigr|^2 + \bigl| P\bar{X}_0^\varepsilon- P \bar{X}_0 \bigr|^2\leq \|e^{-\frac{}{\varepsilon}\tilde{K} t} - P\|_F^2 \bigl|X^{1,\varepsilon}_0 - b \bigr|^2 + \|P\|_F^2 \bigl| X_0^\varepsilon- X_0 \bigr|^2 Now, since $X_0^1=b$, we find essentially the same estimate for $I_1$, which is given by \ref{['eq:Quant-I1-X0inxi']}. The only change is an additional term, called $I_4$ below. It arises due to the shift in the initial conditions and is deterministic. We bound it as follows I_4:=\mathbb{E}\biggl[\sup\limits_{t \in [0,T]} \biggl|\int_0^t \bigl(e^{-\frac{1}{\varepsilon} \tilde{K} (t-s)} - P\bigr) M b0 ds \biggr|^2\biggr] \leq \sup\limits_{t \in [0,T]} t \int_0^t \bigl\|e^{-\frac{1}{\varepsilon} \tilde{K} (t-s)} - P\bigr\|^2_F \|M\|^2_F |b|^2 ds\leq \sup\limits_{t \in [0,T]} t \|M\|^2_F |b|^2 c_P \varepsilon = T \|M\|^2_F |b|^2 c_P \varepsilon where we have used \ref{['eq:matrixexp-int-estimate']} and which adds up to the constant. Note that also the constant factors change, since Young's inequality is applied to four summands in this case, so that \ref{['eq:PathEstSep']} will have a factor of 4 now. The convergence in probability follows from the bound above and Markov's inequality. This completes the proof of \ref{['item:PathQuant']} . Next we prove the pointwise in time estimate \ref{['item:PointQuant']}. We have \mathbb{E}\Bigl[ \bigl|X^\varepsilon_t - Y_t\bigr|^2 \Bigr]\leq 4\mathbb{E}\biggl[ \bigl|e^{-\frac{1}{\varepsilon}\tilde{K} t} X^\varepsilon_0 - Y_0\bigr|^2 \biggr] + 4 \mathbb{E}\biggl[\biggl|\int_0^t \Bigl( e^{-\frac{1}{\varepsilon}\tilde{K}(t-s)}MX^\varepsilon_s - PMY_s\Bigr) ds \biggr|^2 \biggr]\ +4\mathbb{E}\biggl[ \biggl|\int_0^t \Bigl( e^{-\frac{1}{\varepsilon} \tilde{K}(t-s)}C - PC \Bigr)dW_s\biggr|^2 \biggr] + 4\biggl|\int_0^t \bigl(e^{-\frac{1}{\varepsilon} \tilde{K} (t-s)} - P\bigr) M b0 ds \biggr|^2\eqqcolon 4 (\bar{I}_1 +\bar{I}_2 + \bar{I}_3 + \bar{I}_4)\,, and by repeating the calculation in (\ref{['eq:Quant-I1-0']}-\ref{['eq:Quant-I1']}), we obtain \bar{I}_1 \leq 2\bigl\|e^{- \frac{1}{\varepsilon} \tilde{K} t} - P\bigr\|_F^2 \mathbb{E}\bigl[ |X_0^{1} - b|^2\bigr] + 2 \|P\|_F^2\mathbb{E}\bigl[ | X_0^\varepsilon - X_0|^2 \bigr]\,. Depending on whether $K_{11}$ is defective or not, we apply the bound of Corollary \ref{['cor:expK-P']} or \ref{['eq:Frobeniusest_expK-Psymm']} to the term $\bigl\|e^{- \frac{1}{\varepsilon} \tilde{K} t} - P\bigr\|_F^2$. From Itô's isometry along with \ref{['eq:matrixexp-int-estimate']} we bound $\bar{I}_3$ by \bar{I}_3 \leq \varepsilon c_P \|C\|_F^2. Repeating the calculations for $I_4$ and $I_2$, we find \bar{I}_4 \leq \varepsilon t \|M\|^2_F |b|^2 c_P and \bar{I}_2\leq 2\|M\|_F^2 \left( 2t \int_0^t (\|e^{-\frac{1}{\varepsilon }\tilde{K}(t-s)} - P \|_F^2 + \|P\|_F^2 ) \mathbb{E}\bigl[\bigl|X^\varepsilon_s - Y_s\bigr|^2\bigr] ds + \varepsilon c e^{2\lambda_{\max}(PM)t} \right) . Finally, combining these bounds and applying Gronwall's inequality we arrive at (in the case that $K_{11}$ non-defective) \mathbb{E}\Bigl[ \bigl|X^\varepsilon_t - Y_t\bigr|^2 \Bigr] \leq c_1\biggl( e^{-\frac{2\lambda_1 t}{\varepsilon}} \mathbb{E}\bigl[ |X_0^{1} - b|^2\bigr] + \mathbb{E}\bigl[ | X_0^\varepsilon - X_0|^2 \bigr] + \varepsilon\bigl[1+ e^{\lambda_{max}(PM)t} \bigr] \biggr) e^{\varepsilon t + c_2 t^2}, where $c_1,c_2$ are independent of $\varepsilon,t$. The bound for $K_{11}$ being defective follows analogously replacing the exponential decay of the initial term by $e^{-\frac{\lambda_1 -\delta}{\varepsilon}t}$ as given by Corollary \ref{['cor:expK-P']}. The quantitative estimate in Theorem \ref{['thm:pathwiseconvrate']} also works for nonlinear SDEs where the drift $f\colon \mathbb{R}^d \to \mathbb{R}^d$ is smooth and Lipschitz continuous (see below for precise growth conditions). The SDE \ref{['eq:genOU-SC']} then reads $dX_t = f(X_t)dt -\frac{1}{2\varepsilon} K\nabla|\xi(X_t)|^2dt +\sqrt{2}CdW_t$ where the corresponding limit (cf. \ref{['eq:genOU-SC-limit']}) is given by $dY_t = P f(Y_t) dt + \sqrt{2 } PC dW_t,$ with $P$ as defined in \ref{['def:genOU-Proj']}. The only change in Theorem \ref{['thm:pathwiseconvrate']} in this case is in the $I_2$ term. In particular, adding and subtracting $e^{-\frac{1}{\varepsilon}\tilde{K} (t-s)}f(Y_s)$ we can estimate $I_2= \mathbb{E}\biggl[ \sup\limits_{t \in [0,T]} \biggl|\int_0^t \Bigl( e^{-\frac{1}{\varepsilon}\tilde{K}(t-s)}\bigl(f(X^\varepsilon_s) -f(Y_s)\bigr) + (e^{-\frac{1}{\varepsilon}\tilde{K} (t-s)} - P)f(Y_s) \Bigr)ds \biggr|^2 \biggr]\leq 2\mathbb{E}\biggl[ \sup_{t\in [0,T]} \biggl|\int_0^t e^{-\frac{1}{\varepsilon}\tilde{K}(t-s)}\bigl(f(X^\varepsilon_s) -f(Y_s)\bigr) ds\biggr|^2\biggr] + 2\mathbb{E}\biggl[\sup_{t\in [0,T]} \biggl|\int_0^t (e^{-\frac{1}{\varepsilon}\tilde{K} (t-s)} - P)f(Y_s) ds \biggr|^2 \biggr]\leq 2 L_f\mathbb{E}\biggl[ \sup_{t\in [0,T]} \biggl|\int_0^t e^{-\frac{1}{\varepsilon}\tilde{K}(t-s)}\bigl|X^\varepsilon_s -Y_s\bigr| ds\biggr|^2\biggr] + 2 \biggl(\int_0^T \bigl\|e^{-\frac{1}{\varepsilon}\tilde{K}(T-s)}-P\bigr\|_F^2 ds \biggr) \int_0^T \mathbb{E}\bigl[|f(Y_s)|^2]ds\leq 2 T L_f\int_0^T \bigl[ \|e^{-\frac{1}{\varepsilon }\tilde{K}(T-s)} - P \|_F^2 + \|P\|^2_F \bigr] \mathbb{E}\biggl[ \sup_{\tau\in [0,s]}|X^\varepsilon_\tau - Y_\tau|^2 \biggr] ds + 2\varepsilon c_P \int_0^T \mathbb{E}\bigl[|f(Y_s)|^2\bigr]ds$ where the first inequality follows by Young's inequality, the second inequality follows since $f$ is Lipschitz with constant $L_f$ and by using Cauchy-Schwarz inequality in the second integral, the third inequality follows from the bound \ref{['eq:I21-Bound']} above and by using \ref{['eq:matrixexp-int-estimate']} for the second integral. Next we provide a bound for $\mathbb{E}[|f(Y_s)|^2]$. A Lipschitz function $f$ has linear growth at infinity and therefore (under sufficient regularity) we can assume that there exists a constant $c_f>0$ such that $|y\cdot f(y)|\leq c_f(1+|y|^2)$ (or equivalently $|f(y)|\leq c(1+|y|)$ for some constant $c>0$). In the following we use $\rho_t=\mathrm{law}(Y_t)$ which solves $\partial_t\rho_t = -\nabla\cdot \bigl( Pf \rho_t \bigr)+ \nabla^2:\bigl(\bar{C} \bar{C}^T \rho_t \bigr)$ where $\bar{C}=PC$, $\nabla^2$ is the Hessian and $A:B=\mathrm{Tr}(A^TB)$. We have \frac{1}{2} \frac{d}{dt}\mathbb{E}\bigl[|Y_s|^2\bigr]= \frac{1}{2}\frac{d}{dt} \int_{\mathbb{R}^d} |y|^2 \rho_t(dy) = \int_{\mathbb{R}^d} \Bigl( \frac{1}{2} \nabla |y|^2\cdot Pf(y)+ \bar{C}\bar{C}^T : \frac{1}{2} \nabla^2 ( |y|^2) \Bigr) \rho_t(dy)= \int_{\mathbb{R}^d} \Bigl( y\cdot Pf(y)+ \mathrm{Tr}(\bar{C}\bar{C}^T ) \Bigr) \rho_t(dy) \leq \|P\|_F \int_{\mathbb{R}^d} c_f(1+|y|^2) \rho_t(dy) + \|\bar{C}\|_F^2 where the second equality follows by using the dynamics of $\rho_t$ and applying integration by parts. The inequality now follows by applying the growth bounds on $f$ and using the Young's inequality. Using Gronwall's inequality we find $\mathbb{E}\bigl[|f(Y_s)|^2\bigr]\leq c_f(1+\mathbb{E}[|Y_s|^2]) \leq \mathbb{E}\bigl[ |Y_0|^2\bigr] e^{\|P\|_F c_f t } + \frac{\|\bar{C}\|^2_{F}}{c_f\|P\|_F}\bigl(\|\bar{C}\|^2_F e^{\|P\|_F c_f s }-1\bigr) \leq m_1 e^{\|P\|_F c_f s} + m_2$ where $m_1,m_2>0$ are independent of $s$. Substituting this bound back into the bound for $I_2$ \ref{['eq:remI2']} and assuming well-prepared initial data we arrive at the following quantitative estimate for the nonlinear SDE: $\mathbb{E}\biggl[\sup_{t\in [0,T]} \bigl|X_t^\varepsilon - Y_t\bigr|^2 \biggr] \leq \varepsilon c_1e^{c_2 T^2}.$ We note that the results in this section directly apply to the general setting of affine constraint and non-zero mean OU processes by recasting the process to the form \ref{['eq:genOU-SC']} considered in this paper. For affine constraints, this can be done by applying an appropriate similarity transformation, which transforms the affine constraint to a coordinate projection. More precisely, let $\xi(x) = Bx - b, \ B \in \mathbb{R}^{k \times d}, b \in \mathbb{R}^k$, where $B$ has rank $k$. Then $\nabla |\xi|^2 =2 B^T(Bx -b)$ and the similarity transformation is given by $S = V^T00I \in \mathbb{R}^{d \times d},$ where $V$ is the orthonormal matrix that diagonalizes $B$, i.e. $V^TB^TBV = \mathrm{diag}(\lambda_1,\ldots,\lambda_k)$. For non-zero OU process one can consider a coordinate shift, which eventually leads to a shifted level set of the constraint map $\xi.$ In this section we compare our results to those of Katzenberger Katzenberger91, who considers the asymptotic problems of the type studied here for general semi-martingales, and in particular the setting of diffusion processes (see Katzenberger91). We now present the result in Katzenberger91 in the language of this article. To this end consider $dX_t = f(X_t) dt + \frac{1}{\varepsilon} F(X_t) dt + \sqrt{2}C dW_t,$ where $f\colon \mathbb{R}^d \to \mathbb{R}^d$ is locally Lipschitz, $C \in \mathbb{R}^{d \times d}$ and $W_t$ is a $d$-dimensional Brownian motion. The stiff drift term is characterised by the vector field $F\colon \mathbb{R}^d\to\mathbb{R}^d$. In the setting of this paper $f(x)=Mx$ and $F(x) = \frac{1}{2} K \nabla |\xi(x)|^2$. The main result in Katzenberger91 makes use of the following ordinary differential equation: $\dot \psi(z,t) = F(\psi(z,t)) \in \mathbb{R}^d \,, \ \psi(z,0) = z \in \mathbb{R}^d\,.$ Define $\Gamma = \left\{ x \in \mathbb{R}^d : F(x)= 0 \right\} \subset \mathbb{R}^{k}\,$ as the set of fixed points of the ODE or, put differently, the set of points satisfying the constraint (i.e. $\Gamma=\xi^{-1}(0)$ in the language of the previous paragraphs). Furthermore, for a given initial condition $z\in \mathbb{R}^d$, we define the long-time limit of the ODE above \theta(z) :=\lim\limits_{t \to \infty} \psi(z,t) \quad \text{ and } \quad U_{\Gamma} := \left\{x \in \mathbb{R}^d : \theta(x) \text{ exists and } \theta(x) \in \Gamma \right\}. Before stating the main result in Katzenberger91, we introduce the limit process $dY_t = \nabla\theta^T f(Y_t) + \sqrt{2} \nabla \theta^T \, C dW_t \,,$ and the stopping time $\lambda(K) = \left\{ \inf t \geq 0 : Y_t \notin \mathring{K}\right\}$ for a compact set $K \subset \Gamma$, and we write $(Y_t)_{\lambda(K)}$ for the process stopped at $\lambda(K)$. Here $\mathring{K}$ denotes the interior of $K$. Katzenberger91 Assume that $\forall \ y \in U_{\Gamma}: \ \nabla F(y) \in \mathbb{R}^{d \times d}$ has $k$ eigenvalues with negative real part$\theta \in C^2$ and $\nabla \theta, \nabla^2 \theta$ are locally Lipschitz$X^\varepsilon(0) \to Y(0) \in \Gamma$ in probability. Then the solution $X^\varepsilon$ to \ref{['eq:katz-ODE']} satisfies $(X^\varepsilon_t)_{\lambda(K)}\to (Y_t)_{\lambda(K)} \text{ as $\varepsilon\to 0$ }$ in probability in $C([0,\infty))$ uniformly on bounded time intervals. Let us now relate the assumptions of the above Theorem to our assumptions and give explicit expressions for the functions $\psi, \theta$ defining the limiting dynamics $Y$ in our setting. To this end note that for $\xi(x) = x^1 - b \in \mathbb{R}^k$ we have F(x) = -\frac{1}{2} K \nabla |x^1 -b|^2 = - K_{11}0K_{21}0 \left( x - b0 \right) \,, \quad \text{and } \ \nabla F \equiv -K_{11}0K_{21}0 \eqqcolon \tilde{K}. Assumption $(i)$ in Theorem \ref{['thm:katzenberger']} is equivalent to requiring that $\tilde{K}$ or equivalently (for an explanation see the proof of Theorem \ref{['thm:genOU-SC']}) $-K_{11}$ is Hurwitz (as in Theorem \ref{['thm:genOU-SC']}), and Assumption $(ii)$ always holds in our setting as will be illustrated below. Note that Katzenberger provides an implicit form for the soft-constrained limit defined via $\theta$. In contrast, Theorem \ref{['thm:genOU-SC']} directly states the explicit form of the limiting dynamics which includes a projection matrix $P$ (recall Remark \ref{['rem:genOU-proj']} for details). Such projections have recently been studied in related works projection_diffusionZhang20SharmaZhang21 which deal with non-degenerate diffusions (note that in our SDEs the diffusion matrix can be degenerate). Note, however, that the results in Theorem \ref{['thm:genOU-SC']} entail a pointiwse-in-time limit while Katzenberger91 provides convergence of the path over $[0,T]$. This difference is due to the treatment of initial conditions (see Remark \ref{['rem:genOU-IC']} for details). The results in Theorem \ref{['thm:pathwiseconvrate']} thus are a quantitative version of Katzenberger91. Note that these quantitative results easily generalise to nonlinear SDEs (see Remark \ref{['rem:QuantNonlin']}) without requiring any compactness arguments. Additionally, our quantitative pointwise-in-time estimate \ref{['eq:PointBound']} does not require any special treatment of the initial conditions (see Remark \ref{['rem:QuantDisc']} for details). We now examine the limiting dynamics above in the setting of coordinate-projection constraints. In order to explicitly calculate the limit SDE above we need to solve the ODE \ref{['eq:katz-ODE']}, which admits the solution \psi(z,t)= e^{-\tilde{K} t}\psi(z,0)+\int_0^t e^{-\tilde{K} (t-s)} \tilde{K} b0 ds = e^{-\tilde{K} t}z + \bigl(I-e^{-\tilde{K} t}\bigr)b0, and compute it's long-time limit $\theta(z)$. This is contained in Lemma \ref{['lem:matrixexp']} which states that $e^{-\tilde{K} t}\to P$ as $t\to\infty$ (this follows by rescaling time by $\varepsilon$), using which we find $\theta(z) = P \left(z - b0\right)+b0 \in \xi^{-1}(0) \ \forall \, z \in \mathbb{R}^d \quad \text{ and } \quad \nabla \theta^T=P.$ This shows that $\theta$ is indeed a projection onto the manifold $\Gamma = \xi^{-1}(0)$ and $P$ is a projection onto the tangent space of $\xi^{-1}$. So far we have treated soft constraint limits of OU processes with possibly degenerate noise, admitting a unique invariant measure. The underdamped Langevin dynamics with quadratic potential can be seen as a special case of such OU processes. More precisely in the same notation as above, the underdamped Langevin equation can be written asdX_t = (J-A) \nabla H(X_t) dt + \sqrt{2A} dW_t, where X= qp\in \mathbb{R}^{2d},\ A= 000\gamma I\in \mathbb{R}^{2d \times 2d}, \gamma \in \mathbb{R},\ J= 0I-I0 \in \mathbb{R}^{2d \times 2d}, and the quadratic Hamiltonian $H(q,p) = V(q) + \frac{1}{2} |p|^2$, where $V$ is a quadratic potential in $q$. One physically relevant constraint is given by the zero level set of the spatial CG map $\xi(q)=q^1$. A natural choice for the matrix $K$ is $K=(A-J)$; cf. Proposition \ref{['prop:genOU-choiceK']} below. The invariant measure of the corresponding soft constrained SDE dX_t = (J-A) \nabla H(Z_t) dt - (A-J)\frac{1}{2\varepsilon}\nabla |q^1|^2 +\sqrt{2A} dW_t, is given by $\mu^{\varepsilon} = \frac{1}{Z} e^{-V(q) - \frac{1}{2} |p|^2 - \frac{1}{2\varepsilon}|q^1|^2}$ which converges to the conditional Gaussian measure $\mu^{\varepsilon = 0} = \frac{1}{\tilde{Z}} e^{-V(0,q^2) - \frac{1}{2} |p|^2}$ as $\varepsilon \to 0$. Therefore we expect our limit result in Theorem \ref{['thm:genOU-SC']} with $K=A-J$ to hold here. However, we cannot apply Theorem \ref{['thm:genOU-SC']} here because the condition that $-K_{11}$ is Hurwitz is not satisfied as $K_{11} = (A-J)_{11} = 0_{k \times k}.$ Similarly, Theorem \ref{['thm:katzenberger']} of Katzenberger91 does not apply since \nabla F = 0_{d \times d}0_{d \times d}I_{k \times k}0_{k \times (2d-k)}0_{(d-k) \times k}0_{(d-k) \times (2d-k)} and thus all eigenvalues of $\nabla F$ have zero real parts. Hence, we have to resort to other methods in order to prove similar limit results for the underdamped Langevin equation. This is the topic of the companion paper HartmannNeureitherSharma25. The previous section dealt with soft-constrained limits of OU processes without any discussion of the long-time behaviour. However, as stated in the introduction, a goal of soft-constraining is to sample conditional measures on manifolds. In this section we discuss invariant measures in the context of soft-constraining and the role of the matrix $K$ which characterizes the constraint. We answer several questions: (a) does the limit of the soft-constrained yield the correct conditional measure, i.e. is the invariant measure of the projected dynamics \ref{['eq:genOU-SC-limit']} the same as the invariant measure of the unconstrained process conditional on the constraint $\xi$, and (b) what choice for $K$ leads to the correct conditional measure. In this section we will focus on OU processes of the type \ref{['eq:genOU-SC']} which admit a unique Gaussian invariant measure. Therefore, we make the following assumptions throughout this section. $M$ is Hurwitz;$\bigl(M,C\bigr)$ is controllable, i.e. $\mathrm{rank}[C, MC, M^2C,\ldots, M^{d-1}C]=d$. In case $C$ has full rank, the controllability is given. Under these assumptions, the unconstrained SDE \ref{['eq:genOU-SC']}, i.e. with $K\equiv 0$, admits the unique invariant measure (see Proposition \ref{['prop:invmeasOU']}) $\mu\in \mathcal{P}(\mathbb{R}^d), \ \ \mu = \mathcal{N}(0,\Sigma),$ where $\Sigma \in \mathbb{R}^{d\times d}$ is the unique symmetric positive definite solution to the Lyapunov equation $M\Sigma+\Sigma M^T=-2CC^T.$ It turns out that any OU process which admits $\mu$ as invariant measure can be rewritten in the form (see Proposition \ref{['prop:genOU-PHS']}) $dX_t = (J-A)\Sigma^{-1}X_t dt + \sqrt{2}C dW_t,$ where $A\coloneqq CC^T \geq 0 \in\mathbb{R}^{d\times d}$ is symmetric positive semi-definite and $J \coloneqq \frac{1}{2} \bigl(-\Sigma M^T+M\Sigma\bigr) = - J^T \in\mathbb{R}^{d\times d}$ is skew symmetric. The corresponding soft-constrained version reads $dX_t = (J-A)\Sigma^{-1}X_t dt - \frac{1}{2\varepsilon} K \nabla|\xi(X_t)|^2 dt + \sqrt{2}C dW_t,$ which corresponds to \ref{['eq:genOU-SC']} with the choice $M=(J-A)\Sigma^{-1}$. We use the particular form \ref{['eq:genOU-SC-Ham']} for two reasons. First, this form allows us to make an educated guess for the crucial matrix $K$ which encodes how the constraint submanifold $\xi^{-1}(0)$ is approached. Second, the softly-constrained (nonlinear) Langevin dynamics studied in the companion paper HartmannNeureitherSharma25 can also be written in this form with specific choices of $A$ and $J$ and $\Sigma^{-1}X$ will be replaced by the gradient of a given Hamiltonian. The effect of the $K$ matrix will be looked at in detail in Section \ref{['ssec:MatrixK']} and illustrated with a numerical example of a Langevin-type dynamics in Section \ref{['sec:numEx']}. Given a CG map $\xi(x)=x^1-b$ we denote the conditional probability measure of $\mu=\mathcal{N}(0,\Sigma)$ restricted to the level set $\xi^{-1}(0)$ by $\mu_c \in \mathcal{P}(\mathbb{R}^{d-k})$. It is explicitly given by (see for instance Eaton83) $\mu_c =\mathcal{N}(m_c,\Sigma_c) \ \text{ where } \ m_c = \Sigma_{21}\Sigma_{11}^{-1} b \ \text{ and } \Sigma_c = \Sigma_{22} - \Sigma_{21}\Sigma_{11}^{-1} \Sigma_{12}.$ Let us recall the projected dynamics \ref{['eq:genOU-SC-limit']} in this setup, which reads dY_t = P(J-A)\Sigma^{-1}Y_t + \sqrt{2} PC dW_t, \quad \text{ where } P = 00\alphaI \ \text{ and } \alpha= -K_{21} K_{11}^{-1}\,. In the following we discuss the projection $P$ that appears in the limiting dynamics above. We discuss the orthogonality of the projection $P$ for various choices of $K$. First, observe that the projection $P$ is orthogonal with respect to the standard inner product if and only if $\alpha =0$ , i.e. $K_{21}=0$. On the other hand, if $A$ is positive definite, i.e. $A>0$, we can consider the weighted inner product $\langle x,y\rangle_A \coloneqq x^T A^{-1} y$ for $x,y \in \mathbb{R}^d$. In the weighted inner product space, the choice $K=A$ will result in an orthogonal projection, as we now show. The map $P$ is an orthogonal projection with respect to the inner product weighted by $A^{-1}$ if $\langle Pv - v, Pu \rangle_A = 0$ for any $u,v\in\mathbb{R}^d$. Writing $Pu = \left(0,w^2\right)^T$ and using $\alpha = - A_{11}^{-1}A_{12}$, we calculate \langle Pv - v, Pu \rangle_A= \left(v^1\right)^T \left[ -(A_{11})^{-1} A_{12} (A^{-1})_{22} - (A^{-1})_{12} \right] w^2 \,. By the expressions for block-matrix inversion (similar to \ref{['eq:inversesigma']}) we have $-(A_{11})^{-1} A_{12} (A^{-1})_{22} - (A^{-1})_{12}=0$, i.e. $P$ is indeed orthogonal with respect to the inner product weighted by $A^{-1}$ if $K=A>0.$ For a similar discussion see SharmaZhang21. A careful look at the asymptotic result in Theorem \ref{['thm:genOU-SC']} reveals that the limit dynamics $Y_t$ \ref{['eq:lim-HamForm']} is in fact the same as $PX_t$ with $K=0$, i.e the projection of the original unconstrained dynamics. The soft-constrained OU process \ref{['eq:genOU-SC-Ham']} can explicitly be written as dX^1_tdX^2_t = (J-A)\Sigma^{-1} X^1_tX^2_t dt -\frac{1}{\varepsilon} K_{11} X^1_tK_{21} X^1_t dt + \sqrt{2}CdW_t, where we have used the explicit expression \ref{['eq:genOU-ExplConst']} for $\nabla|\xi|^2$. Choosing $K_{21}=0$ implies that the $X^2$ dynamics has no stiff terms (containing $\varepsilon$) and follows the original unconstrained dynamics (i.e. with $K=0$). On the other hand, if $K_{21} \neq 0$ the $X^2$ dynamics is shifted by $\varepsilon^{-1} K_{21} X^1_t$ and in the limit as $\varepsilon \to 0$ this will result in an oblique projection with respect to the standard inner product. The following result identifies general conditions under which the limit dynamics \ref{['eq:lim-HamForm']}, or equivalently \ref{['eq:genOU-SC-limit-Expl']}, admits $\mu_c$ as the correct invariant measure.Let $-K_{11} \in \mathbb{R}^{k \times k}$ be Hurwitz. Define $\hat{M} \coloneqq \alpha M_{12} + M_{22} \in \mathbb{R}^{(d-k) \times (d-k)}$ where $\alpha=-K_{21}K_{11}^{-1}$ and $M=(J-A)\Sigma^{-1}$ (see Theorem \ref{['thm:genOU-SC']}). Assume that $\hat{M}$ is Hurwitz and that $(\hat{M},\hat{C})$ is controllable, where $\hat{C}$ is defined in Theorem \ref{['thm:genOU-SC']}. Then the (limiting) dynamics $(Y^1_t,Y^2_t)$ \ref{['eq:lim-HamForm']} admits the unique invariant measure $\mu^{\varepsilon=0}\in\mathcal{P}(\mathbb{R}^d)$ given by $\mu^{\varepsilon=0}(dy^1,dy^2) = \delta_{Y^1_0}(dy^1) \ \hat{\mu}^{\varepsilon=0}(dy^2), \ \ \ \text{where } \ \ \mathcal{P}(\mathbb{R}^{d-k}) \ni \hat{\mu}^{\varepsilon=0} = \mathcal{N}(\hat{m},\hat{\Sigma}),$ and $\delta$ is the Dirac-delta measure. Here the mean $\hat{m} \in \mathbb{R}^{d-k}$ given by $\hat{m} = - \hat{M}^{-1}(\alpha M_{11}+ M_{21}) Y^1_0$ and the variance $\hat{\Sigma} \in \mathbb{R}^{(d-k)\times (d-k)}$ is the unique positive definite solution to $\hat{M}\hat{\Sigma} + \hat{\Sigma} \hat{M}^T = -2\hat{C} \hat{C}^T.$ Moreover, $\hat{\Sigma} = \Sigma_c$ according to \ref{['eq:genOU-condSteSta']} if and only if the matrix $(\alpha+\Sigma_{21}\Sigma_{11}^{-1}) \left( (J_{11} + A_{11})\alpha^T + (J_{12} + A_{12}) \right) \in \mathbb{R}^{(d-k)\times (d-k)}$ is skew-symmetric. Assuming well-prepared initial datum for $Y^1_0$, i.e. $Y^1_0=b$, we have $\hat{m} = m_c$ if and only if $b\in \mathbb{R}^k$ is in the kernel of the matrix $\Sigma_{21}\Sigma_{11}^{-1} + \hat{M}^{-1}(\alpha M_{11}+M_{21}) \in \mathbb{R}^{(d-k)\times k}.$ In particular, $\hat{m} =m_c$ for $b=0$. Since $Y^1_t \equiv Y^1_0$ for the limiting dynamics \ref{['eq:genOU-SC-limit-Expl']}, the delta measure in $y^1$ follows. The invariant measure for $Y^2_t$ (for fixed value of $Y^1$) follows by using Proposition \ref{['prop:invmeasOU']}. Next we want to discuss conditions under which $\hat{\Sigma}=\Sigma_c$ for various choices of $K$. First note that, using $A = CC^T$ as in \ref{['eq:genOU-SC-Ham']}, we have $\hat{C} \hat{C}^T= (\alpha C_{11} + C_{21}) (\alpha C_{11} + C_{21})^T + (\alpha C_{21} + C_{22}) (\alpha C_{21} + C_{22})^T= \alpha A_{11} \alpha^T + \alpha A_{12} + A_{21} \alpha^T + A_{22} = (PAP^T)_{22}.$ This requires a study of the Lyapunov equation \ref{['eq:Steady-Lyap']} which we now compute explicitly. Using block-matrix inversion, $\Sigma^{-1}=\Sigma_{11}\Sigma_{12}\Sigma_{21}\Sigma_{22}^{-1} = \Sigma_{11}^{-1} - \Sigma_{11}^{-1} \Sigma_{12} \Sigma_c^{-1} \Sigma_{21} \Sigma_{11}^{-1}-\Sigma_{11}^{-1}\Sigma_{12}\Sigma_c^{-1}-\Sigma_c^{-1}\Sigma_{21}\Sigma_{11}^{-1}\Sigma_c^{-1}$ where $\Sigma_c=\Sigma_{22}-\Sigma_{21}\Sigma_{11}^{-1}\Sigma_{12}$ and we have used $\Sigma^T_{12}=\Sigma_{21}$. Using the definition of $M$, it follows that M_{12}=(J_{11}-A_{11})(\Sigma^{-1})_{12} + (J_{12}-A_{12})(\Sigma^{-1})_{22},M_{22}=(J_{21}-A_{21})(\Sigma^{-1})_{12} + (J_{22}-A_{22})(\Sigma^{-1})_{22}, and therefore, using the explicit form of $\Sigma^{-1}$ in \ref{['eq:inversesigma']}, we can expand $\hat{M}=\alpha M_{12} + M_{22}$ to arrive at $\hat{M} = \Bigl[ -\alpha (J_{11}-A_{11}) \Sigma_{11}^{-1}\Sigma_{12} + \alpha(J_{12}-A_{12}) - (J_{21}-A_{21})\Sigma_{11}^{-1}\Sigma_{12} + (J_{22}-A_{22}) \Bigr] \Sigma_c^{-1}.$ Since we are interested in showing that $\hat{\Sigma}=\Sigma_c$, we need to show that the Lyapunov equation \ref{['eq:Steady-Lyap']} holds, i.e. $2\hat{C}\hat{C}^T + \hat{M}\Sigma_c+ \Sigma_c \hat{M}^T = 0\,.$ Using the explicit formulae above and \ref{['eq:CC^T-form']} we have $2\hat{C}\hat{C}^T +\hat{M}\Sigma_c+ \Sigma_c \hat{M}^T= 2\alpha A_{11} \alpha^T + 2\alpha A_{12} + 2 A_{21} \alpha^T + 2 A_{22}\ \ -\alpha (J_{11}-A_{11})\Sigma_{11}^{-1}\Sigma_{12} + \alpha(J_{12}-A_{12}) - (J_{21}-A_{21}) \Sigma^{-1}_{11}\Sigma_{12} + (J_{22}-A_{22})\ \ - \Sigma_{21}\Sigma_{11}^{-1} (-J_{11}-A_{11})\alpha^T + (-J_{21}-A_{21})\alpha^T -\Sigma_{21}\Sigma_{11}^{-1} (-J_{12}-A_{12}) + (-J_{22}-A_{22})= 2\alpha A_{11} \alpha^T - \alpha(J_{11}-A_{11})\Sigma_{11}^{-1}\Sigma_{12} + \alpha(J_{12}+ A_{12}) + \Sigma_{21}\Sigma_{11}^{-1} (J_{11}+A_{11})\alpha^T - (J_{21}-A_{21})\alpha^T\ \ -(J_{21}-A_{21})\Sigma_{11}^{-1}\Sigma_{12} + \Sigma_{21}\Sigma_{11}^{-1}(J_{12}+A_{12})= (\alpha+\Sigma_{21}\Sigma_{11}^{-1})(J_{11}+A_{11})\alpha^T + \Bigl( [\alpha+\Sigma_{21}\Sigma_{11}^{-1}](J_{11}+A_{11})\alpha^T \Bigr)^T\quad +(\alpha+\Sigma_{21}\Sigma_{11}^{-1})(J_{12}+A_{12}) + \Bigl( (\alpha+\Sigma_{21}\Sigma_{11}^{-1})(J_{12}+A_{12}) \Bigr)^T= \underbrace{(\alpha+\Sigma_{21}\Sigma_{11}^{-1}) \left( (J_{11} + A_{11})\alpha^T + (J_{12} + A_{12}) \right)}_{= R} + \underbrace{\left((\alpha+\Sigma_{21}\Sigma_{11}^{-1}) \left( (J_{11} + A_{11})\alpha^T + (J_{12} + A_{12}) \right) \right)^T}_{=R^T}\,,$ where the second equality follows since $A=A^T$ and $J=-J^T$ which implies that $J_{11}^T=-J_{11}$, $J_{12}^T=-J_{21}$, $J_{22}^T=-J_{22}$, and 2\alpha A_{11} \alpha^T= \alpha (A_{11}+J_{11})\alpha^T + \alpha(A_{11}-J_{11})\alpha^T . The Lyapunov equation \ref{['eq:mod-Lyap']} is thus satisfied if and only if the last line in \ref{['eq:mod-Lyap-LHS']} equates to zero. In case $\alpha\neq 0$, this requires $R \in \mathbb{R}^{(d-k) \times (d-k)}$ to be skew symmetric. The condition for $\hat{m}=m_c$ follows directly from the definitions of $\hat{m}$ and $m_c$. The general conditions outlined in the result above do not provide intuition about the structure of $K$, which is required to ensure that the limiting invariant measure $\hat{\mu}^{\varepsilon=0}$\ref{['eq:genOU-limSteSta']} matches the conditional distribution $\mu_c$ \ref{['eq:genOU-condSteSta']}. However, there are some natural choices for $K$, such as $K=A,A-J$, which appear in the literature. Let us first consider the choice $K=A$. Simple choices of the matrices $A,J,\,\Sigma$ lead to $\hat{\mu}^{\varepsilon=0} \neq \mu_c$. For instance, if $A=\Sigma=I\quad\textrm{and} \quad J=0-II0\,,$ i.e. the invariant measure for the unconstrained OU process is $\mu=\mathcal{N}(0,I)$, then $b =\hat{m} \neq m_c=0$ for any $b\neq 0$. Alternatively, for $\beta>1$ and $A=I\,,\quad \Sigma=III\beta I\,,\quad J=0J_{12}-J_{12}^T0\,,$ with $J_{12}\neq -J^T_{12}$, we have that $\hat{\Sigma}\neq \Sigma_c$. The result below discusses other admissible choices for $K$. Note that the choice $K=-(J-A)$ is expected to work since in this special case the softly constrained OU process \ref{['eq:genOU-SC-Ham']} admits the unique invariant measure $\mu^\varepsilon(dx) = Z^{-1}\exp\biggl(-\frac{1}{2} x^T \Sigma^{-1}x - \frac{1}{2\varepsilon}|x^1-b|^2 \biggr),$ where $Z=Z^\varepsilon$ is a normalisation constant, which is expected to converge to 0 as $\varepsilon\to 0$ since the (unnormalised) Gaussian density becomes singular. This is made precise in the following result. Under the assumptions in Proposition \ref{['prop:genOU-SteSt']} the following holds. If $K=A-J$ then $\hat{\mu}^{\varepsilon=0}=\mu_c$. In particular if $J=0$ and $K=A$ then $\hat{\mu}^{\varepsilon=0}=\mu_c$.If $K=\Sigma$ then $\hat{\Sigma}=\Sigma_c$. In general, $\hat{m}\neq m_c$ for $b\neq 0$. Let us mention that in Proposition \ref{['prop:genOU-choiceK']} it is enough to assume that either $\hat{M}$ is Hurwitz or $(\hat{C}, \hat{M})$ is controllable, as the other assumption is implied. This follows because we have a positive definite solution to the Lyapunov equation given by $\Sigma_c$. Recall from the proof of Proposition \ref{['prop:genOU-SteSt']} that $\hat{\Sigma}=\Sigma_c$ if and only if $R=-R^T,$ where $R \in \mathbb{R}^{(d-k)\times (d-k)}$ is defined in \ref{['eq:mod-Lyap-LHS']}. The simplest choice would be $R=0$. Since $R$ is given by the product of two factors, we can choose $\alpha$ such that one factor vanishes to yield $R=0$; this corresponds to two the following choices $K=A-J$ which results in $\alpha=-(J_{21}-A_{21})(J_{11}-A_{11})^{-1}$,$K=\Sigma$ which results in $\alpha=-\Sigma_{21}\Sigma_{11}^{-1}$. Note that the previous asymptotic results only hold if $-K_{11}$ is Hurwitz, which guarantees that $J_{11}-A_{11}$ is invertible when $K=A-J$. Recalling Proposition \ref{['prop:genOU-SteSt']}, we shall now compare the means $m_c=\Sigma_{21} \Sigma_{11}^{-1}b \ \text{ and } \ \hat{m} = -\hat{M}^{-1}(\alpha M_{11} + M_{21})b.$ For $K=A-J$, using the explicit expressions for $\Sigma^{-1}$ in \ref{['eq:inversesigma']} we have \hat{M}= \left[(A_{21} - J_{21})(A_{11} - J_{11})^{-1} (A_{12} - J_{12}) - (A_{22} - J_{22}) \right] \Sigma^{-1}_c\,,\alpha M_{11} + M_{21}= -\left[(A_{21} - J_{21})(A_{11} - J_{11})^{-1} (A_{12} - J_{12}) - (A_{22} - J_{22}) \right] \Sigma^{-1}_c\Sigma_{21}\Sigma_{11}^{-1} and hence $\hat{m}= \Sigma_{21} \Sigma_{11}^{-1} b =m_c$. In summary, $\hat{\mu}^{\varepsilon=0}=\mu_c$ if $K=A-J$. For $K=\Sigma$, choosing $k,d$ even with $2k\leq d$, $b\neq 0$, $A=\Sigma=I$ and $J=-J^T$ the canonical symplectic matrix, we have $b =\hat{m} \neq m_c=0$. Therefore, in general, for the case $K=\Sigma$, $\hat{\mu}^{\xi=0}$ and $\mu_c$ need not agree. We briefly discuss the reversible case as an example, i.e. $A=A^T>0$, $J=0$, and a soft constraint with $K=I$. We show that the projected dynamics does not necessarily leave $\mu_c$ invariant, which means that it is the structure of the noise rather than the reversibility or irreversibility of the dynamics that determines whether the invariant measure is robust under constraining or not. Consider the case $J=0$ and $K=I$, such that $\alpha=-K_{21}K_{11}^{-1}=0$. For $\hat{\Sigma}=\Sigma_c$, Proposition \ref{['prop:genOU-SteSt']} states that the matrix $\Sigma_{21}\Sigma_{11}^{-1}A_{12}$ needs to be skew-symmetric, i.e. $\Sigma_{21}\Sigma_{11}^{-1}A_{12} + A_{21}\Sigma_{11}^{-1}\Sigma_{12} = 0.$ This is true if $A_{12}=0=A_{21}$ or $\Sigma_{12}=0=\Sigma_{21}$, in which case even $\mu^{\xi=0} = \mu_c$. An example, for which $\mu^{\xi=0} \neq \mu_c$ is given by the matrices $A=III\beta I\quad\textrm{and}\quad \Sigma=III\theta I\,,$ with $\beta,\theta>1$. The matrices $A$ and $\Sigma$ are symmetric positive definite, and in the language of Markov chain Monte Carlo algorithms (e.g. girolami2011riemann), the associated OU process is a preconditoned version of $d\tilde{X}_t = -\Sigma^{-1} \tilde{X}_t\,dt + \sqrt{2}\,dW_t\,,$ with uncorrelated noise. By construction, $\tilde{X}$ and the preconditioned system, $dX_t = -A\Sigma^{-1} X_t\,dt + \sqrt{2A}\,dW_t\,,$ with correlated noise, have the same invariant measure $\mu=\mathcal{N}(0,\Sigma)$, but the expression above gives $\Sigma_{21}\Sigma_{11}^{-1}A_{12} + A_{21}\Sigma_{11}^{-1}\Sigma_{12} = 2I \neq 0,$ which entails $\mu^{\xi=0} \neq \mu_c$. Note that if instead we had chosen $A=I=K$ and $J=0$ we would have been back to the case \ref{['item:K=A-J']} and $\mu^{\xi=0} = \mu_c$. The proof of Proposition \ref{['prop:genOU-choiceK']} suggests that there are infinitely many choices for $K$ which lead to $\hat{\mu}^{\varepsilon=0}=\mu_{c}$ which is made precise in the following result. Under the assumptions in Proposition \ref{['prop:genOU-SteSt']} the following holds. If $K=K_{11}*(A-J)_{21}(A-J)_{11}^{-1}K_{11}*$ with $(A-J)_{11}$ invertible, then $\hat{\mu}^{\varepsilon=0}=\mu_c$.If $K=K_{11}*\Sigma_{21}\Sigma_{11}^{-1}K_{11}*$ then $\hat{\Sigma}=\Sigma_c$. In general, $\hat{m}\neq m_c$ for $b\neq 0$.If $\Sigma_{12}=0$ then $K=K_{11}*(A_{21}-J_{21}) (A_{11} + \bar{J} )^{-1}K_{11}*$ where $\bar{J} \in \mathbb{R}^{k \times k}$ is any skew symmetric matrix such that $(A_{11} + \bar{J})$ is invertible, then $\hat{\Sigma}=\Sigma_c$. Before we prove this proposition let us give some intuition about its meaning. It provides us with infinitely many choices for the matrix $K$ in the soft-constrained problem \ref{['eq:genOU-SC-Ham']}, so that we sample the correct target conditional measure. To be more precise, the choice of $-K_{11}$ is free as long as it is Hurwitz. Vividly speaking this means that $X^1$ can approach the constraint manifold $\xi^{-1}(0)$ in any way, e.g. via a gradient flow ($K_{11}=I$) or spiralling towards it (e.g. $K_{11} = I + J$, where $J=-J^T\neq 0$). On the other hand the choice for $K_{21}$ is not free at all and is dictated by the first $k$ columns of the matrices $A, J$ (case \ref{['case:infKA-J']})and $\Sigma$ (case \ref{['case:infKSigma']}) as well as $K_{11}$. This special form of $K_{21}$ guarantees that $X^1$ enters the $X^2$ dynamics in the exact same way as in Proposition \ref{['prop:genOU-choiceK']}. Case \ref{['case:infKSpec']} tells us that if the invariant measure of the unconstrained process \ref{['eq:genOU-SC-Ham']} with $K=0$ is uncorrelated in the first $k$ and remaining $d-k$ dimensions, we can freely choose $K_{21}$ via the skew symmetric matrix $\bar{J} \in \mathbb{R}^{k \times k}.$ The proof of the first two cases is as in the proof of Proposition \ref{['prop:genOU-choiceK']}, realizing that $\alpha$ is unchanged. If $\Sigma_{12}=0$, then for any $\bar{J} = -\bar{J}^T \in \mathbb{R}^{d\times d}$ the Lyapunov equation \ref{['eq:mod-Lyap']} can be written as (cf. \ref{['eq:mod-Lyap-LHS']}) 2\alpha A_{11} \alpha^T + \alpha(J_{12} + A_{12}) + (A_{21} - J_{21})\alpha^T= \alpha \bigl[(\bar{J} + A_{11}) \alpha^T + (J_{12} + A_{12})\bigr] + \left(\alpha \bigl[(\bar{J} + A_{11}) \alpha^T + (J_{12} + A_{12})\bigr] \right)^T =0 \,, where we have used $\bar{J} = - \bar{J}^T \in \mathbb{R}^{d \times d}$, which gives $2\alpha A_{11}\alpha^T = \alpha (A_{11} + \bar{J})\alpha^T + \alpha (A_{11} - \bar{J})\alpha^T.$ Hence, it is easy to check that for any skew symmetric $\bar{J}$ such that $A_{11} + \bar{J}$ is invertible, $\alpha = -(A_{21}- J_{21})(A_{11} + \bar{J})^{-1}$ satisfies \ref{['eq:mod-Lyap']} and hence leads to the correct covariance $\hat{\Sigma} = \Sigma_c$. Equivalently, for a given Hurwitz matrix $K_{11} \in \mathbb{R}^{d \times d}$, choose $K=K_{11}*(A-J)_{21} (A-\bar{J})_{11}^{-1}K_{11}*\,.$ We consider two illustrative examples that demonstrate (a) the uniform pathwise convergence of the softly contrained process on any time interval $[\delta, T]$ for some small $\delta>0$ and (b) the preservation of the invariant measure under hard and soft constraints for a suitably chosen matrix $K$. Consider the $(1+2n)$-dimensional system $dX_t = (J - A)X_t\,dt + \sqrt{2}C\,dW_t\,,$ with $W=(V,U)$ denoting $(1+n)$-dimensional Brownian motion and $J-A = -L0\lambda^T00I-\lambda-I-\gamma\in\mathbb{R}^{(1+2n)\times (1+2n)}\,,\quad C = \sqrt{L}0000\sqrt{\gamma}\in\mathbb{R}^{(1+2n)\times (1+n)}\,.$ where $L>0$, $\gamma\in\mathbb{R}^{n\times n}$ symmetric positive definite, and $\lambda\in\mathbb{R}^{n}$. It can be readily seen that under these assumptions, the matrix $J-A$ is stable, and the pair $(J-A,C)$ is completely controllable. In what follows we use the shorthands $X^1=\zeta\in\mathbb{R}$ and $X^2=(q,p)\in\mathbb{R}^{2n}$, and we consider the codimension 1 constraint $\zeta = 0\,.$ The softly constrained system now reads $dX_t = (J - A)X_t\,dt - \frac{1}{\varepsilon}K X_t\,dt + \sqrt{2}C\,dW_t\,,$ with $K = 100000000\in\mathbb{R}^{(1+2n)\times (1+2n)}\,.$ Following Theorem \ref{['thm:genOU-SC']}, the limit system can be recast as $dY^2_t = (\bar{J} - \bar{A})Y^2_t dt + \sqrt{2}\bar{C}dU_t$ where $Y^2$ denotes the limit of the unconstrained component $X^2=(q,p)$, and $U$ denotes standard $n$-dimensional Brownian motion. The coefficients of the constrained system read $\bar{J} - \bar{A} = 0I-I-\gamma\in\mathbb{R}^{2n\times 2n}\,, \quad \bar{C} = 0\sqrt{\gamma}\in\mathbb{R}^{2n\times n}\,.$ Typical realisations of the softly constrained heat bath model (\ref{['heatbath_eps']}) for $n=1$ and its limit (\ref{['heatbath_lim']}) for the constrained thermostat variable $X^1=0$. The left panel shows the position variable $q$, the right panel shows the momentum variable $p$. Figure \ref{['fig:heatbath1']} shows typical realisations of the softly constrained 3-dimensional heat bath model (\ref{['heatbath_eps']}) for $\varepsilon=1$ and $\varepsilon=0.005$. In the simulation, the initial condition of the constrained variable $X^1=\zeta$ has been set to $X^1_0\neq 0$, which is inconsistent with the requirement $\zeta=0$. Both panels, for the position and variable $q$ (left) and the momentum variable $p$ (right) display the transient initial layer phenomenon for inconsistent initial conditions that is in accordance with Remark \ref{['rem:genOU-IC']}. For $\gamma>0$, the softly constrained dynamics has a unique Gaussian invariant measure with positive definite covariance matrix $\Sigma_\varepsilon$ for all $\varepsilon>0$. Even for $d=3$, the covariance matrix has a complicated expression for $\varepsilon>0$, with correlations between the thermostat variable $X^1=\zeta$ and position and momentum variables. Yet, the limit $\varepsilon\to 0$ is consistent with the constrained dynamics, in that $\lim_{\varepsilon\to 0}\Sigma_\varepsilon = 000010001\,.$ (Here and in what follows, we confine our attention the 3-dimensional case.) As a consequence, the approximation of the time marginal $\mathrm{law}(q_t,p_t)$ by the constrained dynamics, $Y_t$ is uniform in time. The situation changes when $\gamma=0$, where we confine the following discussion to the case $n=1$. The matrix $J-A$ is Hurwitz and so is $J-A- \frac{1}{\varepsilon}K$ for every $\varepsilon>0$ (cf. maddocks1995stability), moreover the pair $(J-A,C)$ is completely controllable. The softly constrained dynamics has a unique zero-mean Gaussian invariant measure $\rho^\varepsilon$ with $3\times 3$ covariance matrix $\tilde{\Sigma}_\varepsilon=\frac{\varepsilon}{1+\varepsilon} I$ that converges to the zero matrix as $\varepsilon\to 0$. As a consequence, the finite time marginal of the process converges in distribution to a Dirac centred at $X=0$, in other words: $\rho^\varepsilon \stackrel{*}{\rightharpoonup} \delta_0\quad\textrm{as}\quad \varepsilon\to 0\,.$ On the other hand, the dynamics (\ref{['heatbath']}) for $\gamma=0$ and subject to the constraint $X^1=0$ is the deterministic Hamiltonian system for $X^2=(q,p)$: $\dot{q}= p\dot{p}= -q\,.$ The dynamics has infinitely many invariant measures, and since (\ref{['hamiltonODE']}) conserves the energy $H(q,p) = \frac{1}{2}\left(q^2 + p^2\right)\,,$ none of these invariant measures of (\ref{['hamiltonODE']}) agrees with the weak--* limit of $\rho^\varepsilon$, unless $q_0=p_0=0$. Nevertheless, assuming consistent initial conditions for $X^1$, the softly constrained dynamics converges to (\ref{['hamiltonODE']}) in a pathwise sense, uniformly on any compact time interval $[0,T]$. While this is in accordance with Theorem \ref{['thm:genOU-SC']}, it implies that the approximation cannot be uniform in time if $\gamma=0$. Indeed, as the left panel of Figure \ref{['fig:heatbath3']} shows, the soft constraint dynamics for initial conditions $X^1=0$ and $X^2=(1,1)$ converges uniformly on $[0,T]$ to the solution of the Hamiltonian system (\ref{['hamiltonODE']}), which is a periodic motion on a circle of radius $\sqrt{2}$. Conversely, it is demonstrated in the right panel of the figure that the long term dynamics for small $\varepsilon$ is dissipative and spirals into the unique fixed point at $(q,p)=(0,0)$. The numerical integration of the stochastic dynamics has been carried out with an Euler-Maruyama scheme with step size $h=10^{-4}$, whereas the deterministic dynamics has been discretised using a symplectic Euler method that guarantees long-term stability and approximate energy conservation. On any bounded time interval, the softly constrained dynamics for $\gamma=0$ converges pathwise to the the limit dynamics (left panel), while the long term dynamics for small $\varepsilon$ departs from the deterministic Hamiltonian limit dynamics and spirals inwards towards the limit invariant measure $\delta_0$ (right panel). It is worth mentioning, that if one makes the correct choice for $K$ according to Proposition \ref{['prop:genOU-choiceK']}, e.g. K = L00000\lambda00, then the softly constrained dynamics admits the unique invariant measure $\rho^\varepsilon=\mathcal{N}(0,\tilde{\Sigma}_\varepsilon)$, with covariance $\tilde{\Sigma}_\varepsilon = \frac{\varepsilon}{1+\varepsilon}00010001\in\mathbb{R}^{3\times 3}\,.$ For every $\varepsilon>0$, the invariant measure $\rho^\varepsilon$ is approached at an exponential rate that is determined by the real part of the principal eigenvalue of the system matrix (J-A)I - \frac{1}{\varepsilon} K = -L(1+\frac{1}{\varepsilon})0\lambda^T00I-\lambda(1+\frac{1}{\varepsilon})-I0, (For the sake of transferrability of the heat bath model, we adopt the matrix notation from the multidimensional case $n>1$, even though $\lambda^T=\lambda\in\mathbb{R}$ for $n=1$.) While one eigenvalue that corresponds to the constrained direction diverges as $\mathcal{O}(-\varepsilon^{-1})$, the matrix has another pair of eigenvalues given by $\nu_{1,2}\simeq -\frac{\lambda^2}{2}\pm\sqrt{\frac{\lambda^2}{4}-L}\quad\textrm{as}\quad \varepsilon\to 0\,.$ For example, for $\lambda=1$ and $L=1$, the spectral abscissa converges to $-1/2$. As a consequence, the (degenerate) Gaussian limit measure $\rho^0 = \delta_{z=0}\otimes \mathcal{N}(0,I_{2\times 2})\,,$ which is the weak--* limit of $\rho^\varepsilon$ as $\varepsilon\to 0$, is approached at a finite exponential rate. This is consistent with the limit dynamics that is given by the projected equation $dY^2_t = (\bar{J}-\bar{A})Y^2_t dt + \sqrt{2}\bar{C} dW_t$ for $Y^2=(q,p)$ that is derived from the oblique projection $P = 0000I0- \lambda L^{-1}0I,$ with the resulting drift and diffusion coefficients (neglecting zero rows and columns) (\bar{J}-\bar{A}) = (P(J-A))_{2 2} = 0I-I- \lambda L^{-1}\lambda^T\,, \quad \bar{C} = (PC)_2 = 0- \lambda L^{-1/2}. As a consequence we recognise (\ref{['eq:proj-GLE']}) as an underdamped Langevin equation with unique invariant measure, given by the $(q,p)$-marginal of $\rho^0$. The negative sign in front of the diffusion coefficient is irrelevant for convergence of the law of paths, but it guarantees a pathwise approximation that is uniform in time. The convergence of the softly constrained dynamics to (\ref{['eq:proj-GLE']}) for all $t>0$ is confirmed by the numerical simulations shown in Figure \ref{['fig:heatbath4']}; the left panel shows a typical realisation for different values of $\varepsilon$. For $\varepsilon\to 0$ the realisation converges almost surely and uniformly on $[0,T]$. On the other hand, the long term dynamics correctly reproduces the (unnormalised) marginal density $\exp(-H)$ as the right panel of the figure shows. Left panel: typical realisations for $\varepsilon\in\{0.1,\,0.01,\,0.001\}$ and the limiting underdamped Langevin dynamics (\ref{['eq:proj-GLE']}), with parameters $\lambda=1$ and $L=1$. Right panel: the long term dynamics correctly reproduces the (unnormalised) standard Gaussian marginal density $\exp(-H(q,p))$. The simulations have been carried out using the gle-BAOAB scheme leimkuhler2022efficient with step size $h=10^{-5}$. As a second example, we consider the computation of the Green's function of a discretised elliptic differential operator on an bounded open domain, specifically we consider the symmetric operator $\mathcal{L} = \Delta - k^2$ on the unit interval $\Omega=(0,1)$, equipped with homogeneous Dirichlet boundary conditions. In the following, we do not distinguish between the infinite-dimensional operator $\mathcal{L}$ as an operator on the domain $\mathcal{D}=H^1_0(0,1)\cap H^2(0,1)$ and its finite-difference approximation as a $d\times d$ matrix. We denote both, the linear operator and its finite-dimensional approximation by the same symbol $\mathcal{L}$. Our aim is to obtain a Monte Carlo approximation of the Green's function, $G$, i.e. the solution to the linear elliptic boundary value problem $\mathcal{L} G(\cdot,a) = \delta_a\,,$ where $\delta_a$ denotes the Dirac delta function centered at $a\in(0,1)$. To this end we adopt ideas put forward in white2009green and consider the OU process $dX_t = (\mathcal{L} X_t - g)\,dt + \sqrt{2}dW_t\,.$ Since the noise term is non-degenerate and $\mathcal{L}$ is symmetric negative definite, the process converges to a unique Gaussian invariant measure with mean $\mu=\mathcal{L}^{-1}g$ and covariance $-\mathcal{L}^{-1}$. Setting $g=\delta_a$, we see that the mean formally agrees with the solution of the linear boundary value problem (\ref{['greens']}). For practical computations, we now discretise the unit interval $[0,1]$ into $d-1$ equally spaced intervals of length $\Delta u=(d-1)^{-1}$ and discretise the operator $\mathcal{L}$ by finite differences, which yields $\mathcal{L} = \frac{1}{(\Delta u)^2}-210\ldots01-21\ldots00\ddots\ddots\ddots0\vdots\vdots\ddots-210\ldots01-2 - k^2 100\ldots0010\ldots00\ddots\ddots\ddots0\vdots\vdots\ddots100\ldots001\in\mathbb{R}^{d\times d}\,.$ By construction, the matrix is Hurwitz, yet the solution does not satisfy the desired boundary conditions. Therefore we impose the constraint $X_1 = X_d = 0$ to implement homogeneous Dirichlet boundary conditions on the solution. Note that $X_1$ and $X_d$ constitute the constrained variable $X^1$, while the interior nodes $X_2,\ldots,X_{d-1}$ form the vector $X^2$. Yet, we refrain from permuting the states in order for the matrix $\mathcal{L}$ to have the standard form of a discrete differential operator. The Green's function is then approximated by computing the stationary (ergodic) mean under hard and soft constraints where the softly constrained dynamics is governed by the SDE $dX_t = (\mathcal{L} X_t - g)\,dt - \frac{1}{\varepsilon}K X_t + \sqrt{2}dW_t\,,$ with $K = 100\ldots0000\ldots00\ddots\ddots\ddots0\vdots\vdots\ddots000\ldots001\in\mathbb{R}^{d\times d}\,.$ Since the resulting projection is orthogonal and the overall dynamics is reversible, the invariant measure of the constrained dynamics agrees with the soft constraint limit. Note that, since the first and the last component, $X_1$ and $X_d$, represent the boundary values of the sampled Green's function, i.e. the constrained variables, the upper left and lower right entries of the matrix $K$ together form the invertible $K_{11}$ block of Theorem \ref{['thm:genOU-SC']}. Numerically computed Green's functions of the elliptic operator $\mathcal{L}=\Delta-k^2$ with homogeneous Dirichlet boundary conditions (left panel) and the resulting relative errors (right panel). Figure \ref{['fig:bvp']} shows the running average of the hard and soft constraint solutions (i.e. the ergodic mean) for $\Delta u=10^{-3}$ together with the exact Green's function $G(u,a) = \exp(-k) \frac{2k \, \exp(ka) - \exp(k(2-a))}{\exp(k) - \exp(-k)} (\exp(ku) - \exp(-ku))\text{if } u \leq a,\exp(-k) \frac{2k \, \exp(ka) - \exp(-ka)}{\exp(k) - \exp(-k)} (\exp(ku) - \exp(k(2-u)))\text{if } u > a\,.$ The Dirac delta function in (\ref{['elliptic']}) and (\ref{['elliptic_eps']}) is approximated by $g = (\Delta u)^{-1}e_{\lceil ad\rceil} \,,$ where $e_k\in\mathbb{R}^d$ denotes the $k$-th unit vector, and $\lceil s\rceil$ is the smallest integer greater than $s\in\mathbb{R}$. The numerical solutions have been computed by averaging over trajectories of length $T=10$ using an Euler-Maruyama discretisation with step size $h=10^{-7}$. The trajectory is subsampled using a macro time step $\Delta t=10^{-4}$ (The small step size $h$ is required since the largest eigenvalue of the matrix $-\mathcal{L}$ is approximately $4\cdot 10^6$.) The confinement parameter $\varepsilon$ has been set to $10^{-3}$. The right panel of Figure \ref{['fig:bvp']} shows the relative error $\frac{|G_{\rm num}(u,a) - |G(u,a)|}{|G(u,a)|}\,,$ which demonstrates good agreement of both hard and soft constraint averages with the exact solution, except at the boundaries $u\in\{0,1\}$ where the relative error is inevitably large, since $G(0,a)=0=G(1,a)$. We have studied the realisation of constraints on linear stochastic differential equations (SDE) by strong confining forces. Specifically, we have considered affine constraints and proved that the dynamics converges pathwise to the solution of a linear SDE that lives solely on the constraint subspace in the limit of infinitely strong confinement. The limit SDE can be recast as a linear projection of the original SDE onto the constraint subspace where the projection can be orthogonal or oblique, depending on the choice of the confining force. Under certain assumptions, the original SDE has a unique Gaussian invariant measure, and we have given necessary and sufficient conditions for the confining force and the resulting projection matrix that guarantee that the constrained dynamics preserves the invariant measure (on the constraint subspace). We have illustrated the theoretical findings with two numerical examples: a linear Hamiltonian system that is coupled to a stochastic heat bath, with a constraint imposed on the heat bath variables, and a linear SDE with a drift matrix that is a discretised elliptic differential operator on a one-dimensional domain, with the constraint representing boundary conditions. Future work ought to address the partially or fully nonlinear case, i.e. nonlinear SDEs with linear or nonlinear constraints. One case that is of particular relevance, both from a theoretical point of view and for practical purposes, is the underdamped Langevin system with different confinement mechanisms (e.g. stiff springs, high friction, or large masses) that act on configuration and/or momentum variables. Langevin systems will be dealt with in a forthcoming paper HartmannNeureitherSharma25. The authors thank Scott Hottovy and Mark Peletier for helpful discussions on the compactness argument related to the work Katzenberger91. This research has been partially funded by the German Federal Government, the Federal Ministry of Education and Research and the State of Brandenburg within the framework of the joint project EIZ: Energy Innovation Center (project numbers 85056897 and 03SF0693A) and by the Deutsche Forschungsgemeinschaft (DFG) through the grant CRC 1114: Scaling Cascades in Complex Systems (project no. 235221301). Throughout the proofs in this article we will make use of the following simple result which summarises an explicit solution for a class of SDEs with linear and nonlinear drift terms. For $t>0$, let $(x_t,y_t)\in\mathbb{R}^\ell\times\mathbb{R}^m$ be a strong unique solution to a coupled SDE system, where $x_t$ evolves according to dx_t = A x_t dt+ f(x_t,y_t,t) dt + C dW_t. Here $A \in \mathbb{R}^{\ell \times \ell}$ and $C \in \mathbb{R}^{\ell \times \ell}$ are constant matrices, $W_t$ is a standard Brownian motion in $\mathbb{R}^\ell$ and $x_0\in\mathbb{R}^\ell$ is the initial condition. Furthermore assume that $f\colon \mathbb{R}^{\ell+m} \times \mathbb{R}_{\geq 0} \to \mathbb{R}^\ell$ satisfies uniform Lipschitz continuity in space, i.e. there exists $L_f>0$ such that for any $z,z'\in\mathbb{R}^{\ell+m}$ we have $|f(z,t)-f(z',t)|\leq L_f|z-z'|$;sublinear growth, i.e. there exists a constant $c>0$ such that for any $z\in\mathbb{R}^{\ell+m}$ and $t>0$ we have $|f(z,t)|\leq c(1+|z|)$. Then $x_t$ can be explicitly written as $x_t = e^{At}x_0 + \int_0^t e^{A(t-s)} f(x_s,y_s,s) ds + \int_0^t e^{A(t-s)} C dW_s\,.$ The existence and uniqueness of the strong solution is standard (see for instance Klenke13). The integral form of the solution follows by using variations of constants along with integration by parts for Itô integrals. Using the stochastic integration by parts formula (see for instance Kallenberg21) and $W_0=0$ almost surely we find \int_0^t e^{A(t-s)} C dW_s = C W_t + \int_0^t A e^{A(t-s)} C W_s ds almost surely. Hence $x_t - CW_t \eqqcolon g_t$, where $x_t$ is given by \ref{['eq:VarofCon-sol']}, is differentiable and satisfies $\dot g_t = A \biggl( e^{At}x_0 + \int_0^t e^{A(t-s)} f(x_s,y_s,s) ds + \int_0^t A e^{A(t-s)} C W_s ds \biggr) + f(x_t,y_t,t) + ACW_t = A x_t + f(x_t,y_t,t).$ Therefore using the definition of $g_t$ we have x_t = g_0 + \int_0^t \frac{dg_s}{ds} ds + CW_t = x_0 + \int_0^t \bigl[Ax_s + f(x_s,y_s,s)\bigr]ds + \int_0^t C\, dW_s, i.e. $x_t$ solves \ref{['eq:SDEmatrixnot']}. In this section we present fundamental results on pathwise solutions of OU processes and the existence of a unique invariant measure. Consider a general OU process d X_t = M(X_t + v) dt + \sqrt{2} C dB_t where $X_t \in \mathbb{R}^d, \, M \in \mathbb{R}^{d \times d}, \, \mathrm{rank}(M)=d, \, C \in \mathbb{R}^{d \times m}, \, \mathrm{rank}(C) \leq d, \, v\in\mathbb{R}^d$ and $B_t$ is $d$-dimensional Brownian motion. Using variation of constants (see Proposition \ref{['prop:varofconst']}) it follows that the solution to \ref{['eq:OUgeneral']} is X_t = e^{Mt} X_0 + (e^{Mt} - I)v +\sqrt{2} \int_0^t e^{M(t-s)} C dB_s\,. Using $m_t=\mathbb{E}[X_t] \in\mathbb{R}^d$ and $\Sigma_t=\mathbb{E}[(X_t-m_t)((X_t-m_t))^T]\in\mathbb{R}^{d\times d}$ for the mean and variance respectively and assuming that $X_0 \sim \mathcal{N}(m_0, \Sigma_0)$ we find that $X_t \sim \mathcal{N}(m_t, \Sigma_t)$ with m_t = e^{Mt} m_0 + e^{Mt}v - v, \ \Sigma_t = e^{Mt} \Sigma_0 e^{M^T t} + 2 \int_0^t e^{M(t-s)} CC^T e^{M^T (t-s)} ds. Here $\mathcal{N}(\cdot,\cdot)$ is the normal distribution. The next proposition states the assumptions under which \ref{['eq:OUgeneral']} admits a unique invariant measure. The OU process defined by \ref{['eq:OUgeneral']} admits a unique invariant measure if and only if the following two conditions are satisfied $M$ is Hurwitz, i.e. its spectrum lies in the open left-plane;$(M,C)$ is controllable, i.e., $\mathrm{rank}[C, MC, M^2C,\ldots, M^{d-1}C]=d$. The invariant measure $\mu$ is given by \mu = \mathcal{N}(-v,\Sigma), where $\Sigma=\Sigma^T >0$ is the unique solution to the Lyapunov equation M \Sigma + \Sigma M^T = -2 CC^T. Note that $e^{Mt} \to 0$ as $t \to \infty$ since $M$ is Hurwitz. Consequently $\lim\limits_{t \to \infty} m_t = -v$ and $\Sigma \coloneqq \lim\limits_{t \to \infty} \Sigma_t = 2 \int_0^\infty e^{Ms} CC^T e^{M^Ts} ds$ if and only of $M$ is Hurwitz. Moreover $\Sigma >0$ if and only if $(M,C)$ is controllable zabczyk2020mathematical. As $\Sigma_t$ \ref{['linSDE-NormalSolution']} solves the differential equation $\frac{d}{dt} \Sigma_t = M \Sigma_t + \Sigma_t M^T + 2 CC^T,$ the limiting covariance $\Sigma$ is the steady state of this equation, i.e. it satisfies \ref{['eq:lyapunovgen']}. The following result shows that the OU process \ref{['eq:linSDE-v0']} can be transformed into the form \ref{['eq:OU_PHS']} used throughout this paper. Consider the linear SDE $dX_t = M(X_t+v) dt + \sqrt{2} C dB_t$ where $v\in\mathbb{R}^d$, $M\in\mathbb{R}^{d\times d}$ is Hurwitz, $C\in \mathbb{R}^{d\times d}$ such that the matrix pair $(M,C)$ is controllable and $W_t$ is the standard $d$-dimensional Brownian motion. Equation \ref{['eq:linSDE-v0']} can be rewritten as $dX_t = (J-A)\Sigma^{-1} (X_t+v) dt + \sqrt{2}CdB_t,$ where $CC^T\eqqcolon A\in\mathbb{R}^{d\times d}$ is symmetric positive semi-definite, $J\in\mathbb{R}^{d\times d}$ is skew symmetric and $\Sigma \in \mathbb{R}^{d \times d}$ is symmetric positive definite, i.e. $A=A^T\geq 0, \ J=-J^T$ and $\Sigma=\Sigma^T>0$. Conversely, for given $J=-J^T, \, \Sigma=\Sigma^T>0$, and $C \in \mathbb{R}^{d \times d}$ which defines $A=CC^T$, \ref{['eq:OU_PHS']} can be reformulated as \ref{['eq:linSDE-v0']} using $M = -(CC^T - J)\Sigma^{-1}\,.$ Furthermore, if $(M,C)$ is controllable, then $M$ is Hurwitz. Proposition \ref{['prop:invmeasOU']} implies the existence of a symmetric positive definite matrix $\Sigma$ (covariance for the invariant measure) which satisfies the Lyapunov equation \ref{['eq:lyapunovgen']}. Multiplying the Lyapunov equation from the right by $\Sigma^{-1}$ we find $M= -\bigl(\Sigma M^T + 2CC^T\bigr)\Sigma^{-1}.$ Next we introduce the skew-symmetric matrix $J$ $J\coloneqq \frac{1}{2} \bigl(-\Sigma M^T+M\Sigma\bigr) = -\Sigma M^T +\frac{1}{2}\bigl(\Sigma M^T+M\Sigma\bigr) = -\Sigma M^T - CC^T,$ using which we can write the drift in \ref{['eq:linSDE-v0']} as $MX= -\bigl(\Sigma M^T + 2CC^T\bigr)\Sigma^{-1}X = \bigl(J-CC^T \bigr)\Sigma^{-1}X.$ This leads to the reformulated SDE \ref{['eq:OU_PHS']} with $A\coloneqq CC^T$. For the converse statement, note that $M$ defined in \ref{['eq:driftdecomp']} satisfies the the Lyapunov equation \ref{['eq:lyapunovgen']}. Since $\Sigma$ is symmetric positive definite, we can define $\bar{M} = \Sigma^{-\frac{1}{2}}M\Sigma^{\frac{1}{2}}$ with positive definite $\Sigma^{\frac{1}{2}},\Sigma^{-\frac{1}{2}}$. Note that $M$ and $\bar{M}$ have the same spectrum with the property that if $(v,\lambda)$ is an eigenpair for $M^T$ then $(\bar{v},\lambda)$, with $\bar{v}=\Sigma^{\frac{1}{2}}v$, is an eigenpair for $\bar{M}^T$. Therefore using \ref{['eq:driftdecomp']} we find (with $v^*$ as the complex conjugate of $v$) \lambda |\bar{v}|^2 = \lambda \bar{v}^* \bar{v} = \bar{v}^* \bar{M}^T \bar{v} = - \bar{v}^* \Sigma^{-\frac{1}{2}} CC^T \Sigma^{-\frac{1}{2}} \bar{v} + \bar{v}^*\Sigma^{-\frac{1}{2}} J^T \Sigma^{-\frac{1}{2}} \bar{v} = - v^* CC^T v - v^* J v. Note that the controllability of $(M,C)$ is equivalent to the statement that no eigenvector of $M^T$ is in the kernel of $C^T$ (see zabczyk2020mathematical and arnold2014sharp), i.e. $v^* CC^T v>0$. Using $v^* J v=0$ since $J= - J^T$ and taking the real part of the equation above we find that $\mathrm{Real}(\lambda)<0$, i.e. $M$ is Hurwitz. \ref{['lem:matrixexpstatement']} Let $\xi\colon \mathbb{R}^d \to \mathbb{R}^k$, $\xi(x) =x^1 - b \in \mathbb{R}^k$, so that $\nabla \xi = I_{k \times k}0_{(d-k) \times k}.$ Also introduce $V \coloneqq 0_{k \times k}I_{(d-k) \times k}\,$ and note that $\nabla \xi^T V = 0 \in \mathbb{R}^{k\times k}$. Computing e^{-\tilde{K} t}=\exp\left(- K_{11}0K_{21}0 t \right) = \exp\left(- K \nabla \xi \nabla \xi^T t \right) = I_{d \times d} + \sum\limits_{j \geq 1} \left(- t\right)^j (K \nabla \xi \nabla \xi^T)^j, it follows that $e^{-\tilde{K} t}V=V$. Similarly, noting that $\nabla \xi^T K \nabla \xi = K_{11}$, we have $e^{-\tilde{K} t} K \nabla \xi = K \nabla \xi e^{-K_{11}t}.$ Writing $[A,B] = a_{11}\cdotsa_{1n}b_{11}\cdotsb_{1m}\vdots\vdotsa_{r1}\cdotsa_{rn}b_{r1}\cdotsb_{rm} \in \mathbb{R}^{r\times (n+m)}$ for two matrices $A \in \mathbb{R}^{r \times n} , \ B \in \mathbb{R}^{r \times m}$ with entries $a_{ij}, b_{il}$ for $1\leq i \leq r, 1\leq j \leq n, 1\leq l \leq m,$ respectively, the last two equations can be combined to give $e^{-\tilde{K} t}\left[ V , K \nabla \xi \right] = \left[ V , K \nabla \xi e^{-K_{11}t} \right]$ or equivalently e^{-\tilde{K} t}= \left[ V , K \nabla \xi e^{-K_{11}t} \right] \left[ V , K \nabla \xi \right]^{-1}. Observe that $\left[ V , K \nabla \xi \right] = 0_{k \times k}K_{11}I_{(d-k) \times k)}K_{21}$ and its inverse is given by $\left[ V , K \nabla \xi \right]^{-1} = -K_{21}K_{11}^{-1}I_{k \times (d-k)}K_{11}^{-1}0$ as can be checked. Therefore e^{-\tilde{K} t}= \left[ V , K \nabla \xi e^{-K_{11}t} \right] \left[ V , K \nabla \xi \right]^{-1} = 0K_{11}e^{- K_{11} t}IK_{21} e^{- K_{11} t} -K_{21}K_{11}^{-1}I_{k \times (d-k))}K_{11}^{-1}0= K_{11} e^{-K_{11}t} K_{11}^{-1}0-K_{21} K_{11}^{-1} + K_{21} e^{-K_{11}t} K_{11}^{-1}I = e^{-K_{11}t}0-K_{21} K_{11}^{-1} + K_{21} e^{-K_{11}t} K_{11}^{-1}I\,, where the last equality follows noting that $K_{11} e^{-K_{11}t} K_{11}^{-1} = e^{-K_{11}t}$. \ref{['lem:matrixexpconvrate']} If $K_{11}$ has eigenvalues with positive real parts (i.e. $-K_{11}$ is Hurwitz) then $e^{-K_{11}t} \to 0_{k \times k}$ as $t \to \infty$ and therefore e^{-\frac{1}{\varepsilon} \tilde{K} t} \xrightarrow{\varepsilon\to 0} P =00-K_{21} K_{11}^{-1}I. \ref{['lem:integratedexpmatrix-projection']} Using the explicit firm of the matrix exponential we have \bigl\|e^{- \frac{1}{\varepsilon} \tilde{K} (t-s) } - P \bigr\|_F^2= \mathrm{Tr}\bigl(e^{-\frac{1}{\varepsilon}K^T_{11}(t-s)} e^{-\frac{1}{\varepsilon}K_{11}(t-s)} + K_{11}^{-T} e^{-\frac{1}{\varepsilon}K_{11}^T(t-s)} K_{21}^T K_{21} e^{-\frac{1}{\varepsilon}K_{11}(t-s)} K_{11}^{-1}\bigr)= \mathrm{Tr}\bigl(e^{-\frac{1}{\varepsilon}K^T_{11}(t-s)} e^{-\frac{1}{\varepsilon}K_{11}(t-s)} \bigr)+ \mathrm{Tr}\bigl(K_{11}^{-T} e^{-\frac{1}{\varepsilon}K_{11}^T(t-s)} K_{21}^T K_{21} e^{-\frac{1}{\varepsilon}K_{11}(t-s)} K_{11}^{-1}\bigr)= \bigl\| e^{-\frac{1}{\varepsilon}K_{11}(t-s)}\bigr\|_F^2+ \bigl\|K_{21}e^{-\frac{1}{\varepsilon}K_{11}(t-s)} K_{11}^{-1}\bigr\|_F^2\leq \bigl\| e^{-\frac{1}{\varepsilon}K_{11}(t-s)}\bigr\|_F^2 \Bigl( 1 + \bigl\| K_{21}\bigr\|_F^2 \bigl\| K^{-1}_{11}\bigr\|_F^2 \Bigr), where we used the sub-multiplicativity of the Frobenius norm to arrive at the inequality. We can write K_{11} = V \left( \Lambda + N \right) V^{-1} in its Jordan normal form, where $\Lambda$ is a diagonal matrix with the eigenvalues of $K_{11}$, $N$ is an upper triangular nilpotent matrix of order $m$, i.e. $N^m = 0$ for some $m \in \mathbb{N}$, and $V$ consists of the (generalized) eigenvectors. Using $e^{-K_{11}} = V e^{-\Lambda} e^{-N} V^{-1}$ together with the sub-multiplicativity of the Frobenius norm, we have \| e^{-\frac{1}{\varepsilon}K_{11}(t-s)} \|_F^2= \| V e^{-\frac{1}{\varepsilon}\Lambda(t-s)} e^{-\frac{1}{\varepsilon}N(t-s)}V^{-1} \|_F^2 \leq \| V \|_F^2 \| e^{-\frac{1}{\varepsilon}\Lambda(t-s)} \|_F^2 \| e^{-\frac{1}{\varepsilon}N(t-s)} \|_F^2 \| V^{-1} \|_F^2\,. Denote the real parts of the eigenvalues of $K_{11}$ by $\lambda_i>0$ with $i=1,\ldots,k$, and assume them to be ordered, i.e. $0<\lambda_1 \leq \ldots \leq \lambda_k$. Using this we find \bigl\| e^{-\frac{1}{\varepsilon}\Lambda(t-s)} \bigr\|_F^2 = \mathrm{Tr}( e^{-\frac{1}{\varepsilon}(\Lambda + \Lambda^*) (t-s)} ) = \sum\limits_{i=1}^k e^{-\frac{2}{\varepsilon}\lambda_i (t-s)} \leq k e^{-\frac{2}{\varepsilon} \lambda_{1}(t-s)}, where $\Lambda^*$ is the complex conjugate of $\Lambda$. Next, observe that since $N$ is a nilpotent matrix of order $m$, i.e. $N^m=0$, we have $e^{-\frac{1}{\varepsilon}N} = \sum\limits_{j=0}^{m-1} \left(\frac{-N}{\varepsilon}\right)^j\frac{1}{j!}$, which leads to \| e^{-\frac{1}{\varepsilon}N(t-s)}\|_F^2 \leq m\sum_{j=0}^{m-1} \left(\frac{(t-s)}{\varepsilon}\right)^j\frac{1}{j!} \|N^j\|_F^2, where the constant $m$ arises due to Young's inequality. Using the condition number $\kappa(V)=\|V\|_F^2\|V^{-1}\|_F^2$ we arrive at the overall estimate \| e^{-\frac{1}{\varepsilon}K_{11}(t-s)} \|_F^2 \leq \kappa(V) km e^{-\frac{2}{\varepsilon} \lambda_1(t-s)} \sum\limits_{j=0}^{m-1}(t-s)^j \frac{\|N^j \|_F^2}{j! \varepsilon^j}. Substituting into \ref{['eq:app-Exp-pre-explicit']} we arrive at \| e^{-\frac{1}{\varepsilon}K_{11}(t-s)} -P \|_F^2 \leq \Bigl( 1 + \bigl\| K_{21}\bigr\|_F^2 \bigl\| K^{-1}_{11}\bigr\|_F^2 \Bigr) \kappa(V) k m e^{-\frac{2}{\varepsilon} \lambda_1(t-s)} \sum\limits_{j=0}^{m-1}(t-s)^j \frac{\|N^j \|_F^2}{j! \varepsilon^j}, which completes the proof of \ref{['eq:Frobeniusest_expK-P']}. Next we consider the special case that $K_{11}$ is symmetric (consequently non-defective), and therefore all its eigenvalues are real. Then $K_{11} = O \Lambda O^T$ for some orthonormal matrix $O$ (i.e. $OO^T=I$) and we arrive at \|e^{-\frac{1}{\varepsilon}K_{11}(t-s)}\|_F^2= \mathrm{Tr}(O e^{-\frac{1}{\varepsilon}\Lambda(t-s)}O^T O e^{-\frac{1}{\varepsilon}\Lambda(t-s)}O^T) = \| e^{-\frac{1}{\varepsilon}\Lambda(t-s)} \|_F^2 \leq k e^{-\frac{2\lambda_1}{\varepsilon}(t-s)} , where $\lambda_1>0$ is the smallest eigenvalue of $K_{11}$ as before. In this case $N\equiv 0$ and using $\kappa(V)=1$ we arrive at part \ref{['item:sym-nonDef']}. Note that the same arguments as above also apply for any non-defective matrix (now $\kappa(V)\neq 1$ but $N\equiv 0$) which leads to part \ref{['item:nonDef']} (note that we have used the time-integral bound derived below as well). Finally, we prove the bound \ref{['eq:matrixexp-int-estimate']} on the time integral. To achieve this we use the bound ($j\in \{0,\ldots, m-1\}$) $\int_0^t e^{-\frac{2\lambda_1}{\varepsilon}(t-s)} \frac{(t-s)^j}{\varepsilon^j} ds= \int_0^t e^{-\frac{2\lambda_1}{\varepsilon}s} \frac{s^j}{\varepsilon^j} ds = - \frac{t^j}{2\lambda_1 \varepsilon^{j-1}} e^{-\frac{2\lambda_1}{\varepsilon}t} + \frac{j}{2\lambda_1}\int_0^t e^{-\frac{2\lambda_1}{\varepsilon}s} \frac{s^{j-1}}{\varepsilon^{j-1}} ds= \ldots = - e^{-\frac{2\lambda_1}{\varepsilon}t} \sum\limits_{\ell=1}^{j+1}\frac{ t^{j-\ell+1}}{\varepsilon^{j-\ell}\lambda_1^\ell } \frac{j!}{(j-\ell+1)!} + \frac{j!}{(2\lambda_1)^{j+1}}\varepsilon \leq \frac{j!}{(2\lambda_1)^{j+1}}\varepsilon,$ where the first equality follows via the substitution $r=t-s$ and the rest of the equalities follow by performing a series of integration by parts. The bound \ref{['eq:matrixexp-int-estimate']} then follows immediately by using \ref{['eq:app-Exp-pre-explicit']}-\ref{['eq:Frobeniusest_expK-P-Explicit']} along with the bound above, since \int_0^t \| e^{-\frac{1}{\varepsilon}K_{11}(t-s)} - P \|_F^2 ds\leq \Bigl( 1 + \bigl\| K_{21}\bigr\|_F^2 \bigl\| K^{-1}_{11}\bigr\|_F^2 \Bigr) \kappa(V) k m\int_0^t e^{-\frac{2}{\varepsilon}\lambda_1(t-s)}\sum_{j=0}^{m-1} \frac{(t-s)^j}{\varepsilon^j}\frac{1}{j!}\|N^j\|_F^2 ds\leq \varepsilon \Bigl( 1 + \bigl\| K_{21}\bigr\|_F^2 \bigl\| K^{-1}_{11}\bigr\|_F^2 \Bigr) \kappa(V) k m\sum\limits_{j=0}^{m-1}\frac{\|N^j\|_F^2}{(2\lambda_1)^{j+1}}. Let $\delta,\varepsilon>0$ and introduce f_j(t) \coloneqq e^{-\delta t} \beta_j \left(\frac{t}{\varepsilon}\right)^j\,, \ \beta_j > 0, where $t\in [0,\infty)$ and $j\in\mathbb{N}$. Then $f_j(0) = 0$ and $\lim\limits_{t \to \infty} f_j(t) = 0$ and $f_j'(0) = \beta_j > 0$ for any $j$. Hence, there exists $t_j^* > 0$ such that $f_j(t^*_j) = \max_{t \geq 0} f_j(t) > 0$ and computing $f'_j(t) = \left(- \delta \beta_j t^j \varepsilon^{-j} + \beta_j j t^{j-1}\varepsilon^{-j}\right)e^{-\delta t}$ we find $t^*_j = \frac{j}{\delta}$. Hence $f_j(t^*_j) = \beta_j \left(\frac{j}{\delta \varepsilon e}\right)^j$. With this e^{-2\frac{\lambda}{\varepsilon} t } \sum\limits_{j=0}^{m-1} \beta_j \left(\frac{t}{\varepsilon}\right)^j = e^{-2\frac{\lambda}{\varepsilon} t + \delta t} \sum\limits_{j=0}^{m-1} f_j(t) \leq e^{-2\frac{\lambda}{\varepsilon} t + \delta t} \sum\limits_{j=0}^{m-1} \beta_j \left(\frac{j}{\delta\varepsilon e}\right)^j\,. Now let $\delta = 2\frac{\tilde{\delta}}{\varepsilon}$ for $\tilde{\delta}\in (0,\lambda)$. Inserting $\beta_j = \|N^j\|_F^2 \leq \|N\|_F^{2j}$ which holds by sub-multiplicativity of the norm, we find that for any $\tilde{\delta} \in (0,\lambda)$ e^{-2\frac{\lambda}{\varepsilon} t } \sum\limits_{j=0}^{m-1} \|N^j\|_F^{2} \left(\frac{t}{\varepsilon}\right)^j \leq e^{-2\frac{\lambda -\tilde{\delta}}{\varepsilon} t } \sum\limits_{j=0}^{m-1} \|N^j\|_F^{2} \left(\frac{j}{2\tilde{\delta} e}\right)^j\, \leq e^{-2\frac{\lambda -\tilde{\delta}}{\varepsilon} t } (m-1) \left(\frac{{m-1}}{2\tilde{\delta} e}\right)^{m-1} \|N\|_F^{2(m-1)}. Anton Arnold and Jan Erb. Sharp entropy decay for hypocoercive and non-symmetric fokker-planck equations with linear drift. arXiv preprint arXiv:1409.5425, 2014.Joris Bierkens. Non-reversible Metropolis-Hastings. Statistics and Computing, 26(6):1213--1228, 2016.Folkmar A Bornemann and Christof Schütte. Homogenization of hamiltonian systems with a strong constraining potential. Physica D: Nonlinear Phenomena, 102(1-2):57--77, 1997.M. Brubaker, M. Salzmann, and R. Urtasun. A family of MCMC methods on implicitly defined manifolds. In N. D. Lawrence and M. Girolami, editors, Proceedings of the Fifteenth International Conference on Artificial Intelligence and Statistics, volume 22 of Proceedings of Machine Learning Research, pages 161--172. PMLR, 2012.G. Ciccotti, R. Kapral, and E. Vanden-Eijnden. Blue moon sampling, vectorial reaction coordinates, and unbiased constrained dynamics. ChemPhysChem, 6(9):1809--1814, 2005.G. Ciccotti, T. Lelièvre, and E. Vanden-Eijnden. Projection of diffusions on submanifolds: Application to mean force computation. Commun. Pure Appl. Math., 61(3):371--408, 2008.Stephen L Campbell and Nicholas J Rose. Singular perturbation of autonomous linear systems. SIAM Journal on Mathematical Analysis, 10(3):542--551, 1979.Morris L. Eaton. Multivariate statistics. Wiley Series in Probability and Mathematical Statistics: Probability and Mathematical Statistics. John Wiley & Sons, Inc., New York, 1983.Jaap Eldering. Realizing nonholonomic dynamics as limit of friction forces. Regular and Chaotic Dynamics, 21:390--409, 2016.Mark Girolami and Ben Calderhead. Riemann manifold langevin and hamiltonian monte carlo methods. Journal of the Royal Statistical Society Series B: Statistical Methodology, 73(2):123--214, 2011.Carsten Hartmann. An ergodic sampling scheme for constrained hamiltonian systems with applications to molecular dynamics. Journal of Statistical Physics, 130:687--711, 2008.Chii-Ruey Hwang, Shu-Yin Hwang-Ma, and Shuenn-Jyi Sheu. Accelerating Gaussian diffusions. Ann. Appl. Probab., 3(3):897--913, 1993.EJ Hinch. Brownian motion with stiff bonds and rigid constraints. Journal of Fluid Mechanics, 271:219--234, 1994.Carsten Hartmann, Lara Neureither, and Upanshu Sharma. Soft constraints in Langevin dynamics. in preparation, 2025.Olav Kallenberg. Foundations of modern probability. Springer, 2021.GS Katzenberger. Solutions of a stochastic differential equation forced onto a manifold by a large drift. The Annals of Probability, 19(4):1587--1628, 1991.Achim Klenke. Probability theory: a comprehensive course. Springer Science & Business Media, 2013.Matthias Klar, Karsten Matthies, Celia Reina, and Johannes Zimmer. Second-order fast--slow dynamics of non-ergodic hamiltonian systems: Thermodynamic interpretation and simulation. Physica D: Nonlinear Phenomena, 428:133036, 2021.VV Kozlov and AI Neishtadt. On the realization of holonomic constraints. Journal of Applied Mathematics and Mechanics, 54(5):705--708, 1990.Thomas G Kurtz and Philip Protter. Weak limit theorems for stochastic integrals and stochastic differential equations. The Annals of Probability, pages 1035--1070, 1991.Ian Langmore, Michael Dikovsky, Scott Geraedts, Peter Norgaard, and Rob Von Behren. Hamiltonian Monte Carlo in inverse problems. Ill-Conditioning and multimodality. International Journal for Uncertainty Quantification, 13(1), 2023.Tony Lelièvre, Francis Nier, and Grigorios A Pavliotis. Optimal non-reversible linear drift for the convergence to equilibrium of a diffusion. J. Stat. Phys., 152(2):237--274, 2013.Tony Lelievre, Francis Nier, and Grigorios A Pavliotis. Optimal non-reversible linear drift for the convergence to equilibrium of a diffusion. Journal of Statistical Physics, 152(2):237--274, 2013.Benedict Leimkuhler, Timothée Pouchon, Tiffany Vlaar, and Amos Storkey. Constraint-based regularization of neural networks. arXiv e-prints, arXiv:2006.10114, 2020.Tony Lelièvre, Mathias Rousset, and Gabriel Stoltz. Langevin dynamics with constraints and computation of free energy differences. Math Comput., 81(280):2071 -- 2125, 2012.Benedict Leimkuhler and Matthias Sachs. Efficient numerical algorithms for the generalized Langevin equation. SIAM Journal on Scientific Computing, 44(1):A364--A388, 2022.John H. Maddocks, Michael L. Overton, and Maddocks. Stability theory for dissipatively perturbed Hamiltonian systems. Communications on Pure and Applied Mathematics, 48(6):583--610, 1995.Michela Ottobre, Natesh S. Pillai, and Konstantinos Spiliopoulos. Optimal scaling of the MALA algorithm with irreversible proposals for gaussian targets. Stochastics and Partial Differential Equations: Analysis and Computations, 8:311--361, 2020.Sebastian Reich. Smoothed langevin dynamics of highly oscillatory systems. Physica D: Nonlinear Phenomena, 138(3-4):210--224, 2000.Soumava Kumar Roy and Mehrtash Harandi. Constrained stochastic gradient descent: The good practice. In 2017 International Conference on Digital Image Computing: Techniques and Applications (DICTA), pages 1--8. IEEE, 2017.Hanan Rubin and Peter Ungar. Motion under a strong constraining force. Communications on pure and applied mathematics, 10(1):65--87, 1957.Zexi Song and Zhiqiang Tan. On irreversible Metropolis sampling related to Langevin dynamics. SIAM Journal on Scientific Computing, 44(4):A2089--A2120, 2022.Upanshu Sharma and Wei Zhang. Nonreversible sampling schemes on submanifolds. SIAM Journal on Numerical Analysis, 59(6):2989--3031, 2021.Jessika Walter, Carsten Hartmann, and John H Maddocks. Ambient space formulations and statistical mechanics of holonomically constrained langevin systems. The European Physical Journal Special Topics, 200(1):153--181, 2011.David White and Andrew Stuart. Green’s functions by Monte Carlo. In Monte Carlo and Quasi-Monte Carlo Methods 2008, pages 627--636. Springer, 2009.Kerun Xu and Miranda Holmes-Cerfon. Monte carlo on manifolds in high dimensions. Journal of Computational Physics, 506:112939, 2024.Jerzy Zabczyk. Mathematical control theory. Springer, 2020.Wei Zhang. Ergodic SDEs on submanifolds and related numerical sampling schemes. ESAIM: Mathematical Modelling and Numerical Analysis, 54(2):391--430, 2020.Benjamin Jiahong Zhang. Efficient Sampling Methods of, by, and for Stochastic Dynamical Systems. PhD thesis, Massachusetts Institute of Technology, 2022.