\section{Introduction}
In many areas of physics and engineering such as data science, machine-learning, molecular dynamics, astronomy, robotics, etc., we find that the state space of the dynamical system is a manifold. When these systems are subject to noise and operation/performance related constraints (e.g. no-slip velocity constraint in physics, data manifold as constraint in the subject of data science, trajectory constraint in robotics, etc.), it is natural to model the dynamical system using constrained Stochastic Differential Equations (SDEs) on manifolds. In this article, we focus on SDEs on manifolds that are constrained to evolve on a level set of a map between two manifolds. We call these equations as \textit{explicit Stochastic Differential-Algebraic Equations (SDAEs)}.

In our earlier work \cite{suthar2021explicit}, we have defined a system of explicit Stochastic Differential Algebraic Equations on Euclidean spaces as
\begin{subequations}
\label{eq:SDAE1}
\begin{equation}
\label{eq:SDAE-SDE}
dx = f(x,u)dt + \sigma(x,u)dW_t,  
\end{equation}
\begin{equation}
\label{eq:SDAE-AE}
    h(x,u) = 0
\end{equation}
\end{subequations}
where $W_t = (W^1_t, W^2_t, ..., W^d_t)$ is a d-dimensional Wiener process, $f:\mathbb{R}^n\times \mathbb{R}^m\to \mathbb{R}^n$, and $\sigma:\mathbb{R}^n\times \mathbb{R}^m\to L(\mathbb{R}^d,\mathbb{R}^n)$, $h:\mathbb{R}^n\times \mathbb{R}^m \to \mathbb{R}^q$, $x(t,\omega) \in \mathbb{R}^n$, and $u(t,\omega) \in \mathbb{R}^m$. In case the stochastic processes are evolving on a manifold, then this definition can be directly extended by considering a generic form of explicit SDAEs on manifolds as
\begin{subequations}
\label{equation:SDAE on manifold strato}
\begin{equation}
\label{equation:SDAE on manifold - eq strato}
\delta X_t = V(X_t,U_t) dt + \sum_{l = 1}^d\sigma_l(X_t,U_t)\circ dW^l_t,
\end{equation}
\begin{equation}
\label{equation:SDAE on manifold - constraint strato}
h(X_t,U_t) = p;
\end{equation}
\end{subequations}
where the SDE in equation \ref{equation:SDAE on manifold - eq strato} is in Stratonovich representation, $X_t\in M$, $U_t\in N$, $p\in P$, $h:M\times N \to P$, $V:M\times N\to TM$ such that $\tau_M(V(\cdot, u)) = Id_M$ for all $u\in N$, $\sigma_l:M\times N\to TM$ such that $\tau_M(\sigma_l(\cdot, u)) = Id_M$ for all $u\in N$ and $l\in \{1, 2, ..., d\}$; and $M$, $N$, and $P$ are smooth manifolds of dimension $n$, $m$, and $q$ respectively. However, on manifolds, Stratonovich representation of SDEs make the stochastic analysis in the mean sense difficult. This difficulty is avoided if the SDE is represented in the standard Ito form. In \cite{gliklikh2011global}, readers will find that the standard Ito SDE representation requires the manifold to be equipped with a connection. If the manifolds is not equipped with a connection then the only way to avoid the difficulty in stochastic analysis in mean sense, is by considering the Intrinsic representation. The Intrinsic form of representing SDEs depends on Schwartz morphism that morphs a semi-martingale from one manifold to the other. In our earlier work \cite{suthar2022Intrinsic}, we have demonstrated a few techniques of constructing Schwartz morphism from $\mathbb{R}^{p+1}\to M$ such that it morphs the stochastic process $(t,W_t)\in \mathbb{R}^{p+1}$ into a semi-martingale on $M$, where $W_t$ is a $p$-dimensional Wiener process. These constructions relied on a special type of fiber preserving maps from the tangent bundle to the diffusion bundle. We have called these maps as  diffusion generators. This allowed us to write the SDE in terms of the differential $dt$ and the Ito differential $dW_t$, hence making the stochastic analysis in mean sense easier. In this article, we combine the ideas on explicitly constrained SDEs on Euclidean spaces from \cite{suthar2021explicit} with that on Intrinsic SDEs via diffusion generator approach from \cite{suthar2022Intrinsic}.

The main results of this manuscript can be found in section \ref{section:explicit SDAEs on mani} and section \ref{section:approximate solution}. In section \ref{section:explicit SDAEs on mani} we have introduced the notion of index for SDAEs on manifolds. We have proven the existence of a class of SDAEs that have no solution. We call these equations as \textit{ill-posed} SDAEs. The existence of these type of equations make it difficult to have a general method to solve high index SDAEs. However, if we relax the constraint, we find that the stochastic process may hit the cut-locus where the distance function loses its differentiability. We have discussed such approximate relaxed solution in section \ref{section:approximate solution} and have also shown that it is possible to reduce the probability of hitting the cut-locus. We have presented an algorithm for numerical computations in section \ref{section:numerical}.

\section{Basic definitions and notations}
In this section we summarize the key ideas of Stochastic Differential Equations on Manifolds that have been discussed in our earlier work \cite{suthar2022Intrinsic}. We will be following the notational convention from the same work \cite{suthar2022Intrinsic}. To recall the notational conventions and the basic definitions,
\begin{itemize}
\item Sections of any fiber bundle $F$ are denoted by $\Gamma(F)$ and the set of all smooth vector fields on manifold $M$ by $\mathfrak{X}(M)$. The set of all smooth functions is denoted by $\mathfrak{F}(M)$.
\item We call the {\textit{second order tangent bundle}} of an n-manifold $M$ as \textbf{\textit{diffusion bundle}} to avoid confusing it with second tangent bundle TTM.  We denote the diffusion bundle by $\mathfrak{D}M$. The fiber of the diffusion bundle at point $x\in M$ is called diffusion space and is denoted by $\mathfrak{D}_xM$. The elements of the diffusion space are called \textbf{\textit{diffusors}}. A \textbf{\textit{diffusor field}} is defined as the smooth sections of the diffusion bundle $\mathfrak{D}M$. Following our notational convention for section of a fiber bundle, the set of all smooth diffusor field is denoted as $\Gamma(\mathfrak{D}M)$.
\item If $\phi:M\to N$, the fiber preserving map over $\phi$, $\mathfrak{D}\phi:\mathfrak{D}M\to \mathfrak{D}N$ is called the \textit{\textbf{diffusion map}}.
\item Locally in charts $(U,\Upsilon)$ on $M$ and $(V,\chi)$ on $N$, for all $L\in \mathfrak{D}M$ such that $L|_U = a^i\partial_i + b^{ij}\partial^2_{ij}$,
\[\mathfrak{D}\phi \left(L\right)|_V = \left[ a^i\partial_i\phi^k + b^{ij}\partial^2_{ij}\phi^k\right] \partial_k + \left[b^{ij}\partial_i\phi^k \partial_j\phi^l\right]\partial^2_{kl}.\]
\item Given $L\in \mathfrak{D}_xM$, there exists a symmetric contravariant tensor $\hat{L}\in T^2_0M$ such that
\[\hat{L}(df, dg) = \dfrac{1}{2}(L[fg] - fL[g] - gL[f]).\] Moreover, if $L|_U = a^i\partial_i + b^{ij}\partial^2_{ij}$ then locally \[\hat{L}(df,dg) = b^{ij}\partial_if\partial_jg.\]
\noindent In other words, $\hat{\cdot} : \mathfrak{D}M\to T_0^2M$ such that $\hat{L}$ is the  symmetric part of the diffusor $L$.
\item Consider two manifolds $M$ and $N$ with $x\in M$ and $y\in N$. If there exists a linear map $J(x,y):\mathfrak{D}_xM\to \mathfrak{D}_yN$ such that $Img(J|_{T_xM})\subset T_yN$ and $\widehat{JL} = (J|_{T_xM}\otimes J|_{T_xM}) \hat{L},$ then such a map $J$ is called a \textbf{\textit{Schwartz morphism}}.
\item In Schwartz's second order stochastic differential geometry, the key idea is the infinitesimal representation of a semi-martingale. For a semi-martingale $X_t$ on a manifold $M$, an infinitesimal increment $\mathbf{d}X_t$ is defined as an element of the diffusion space at $X_t\in M$.
\item This infinitesimal increment is called the Intrinsic differential or the Schwartz differential. Using this differential, and the Schwartz morphism $J$ from $M$ to $N$, it is possible to morph the semi-martingale $X_t\in M$ into a semi-martingale $Y_t\in N$, such that locally we have
\[\mathbf{d}Y_t = J(X_t,Y_t)\mathbf{d}X_t.\]
These equation have been called intrinsic SDEs in \cite{emery2012stochastic}.
\end{itemize}
In \cite{suthar2022Intrinsic}, we have presented techniques of constructing Schwartz morphism from $\mathbb{R}^{p+1}\to M$. This means that this Schwartz morphism can transform the stochastic process $(t,W_t)\in \mathbb{R}^{p+1}$ into a semi-martingale on $M$, where $W_t$ is a $p$-dimensional Wiener process. The key idea of the constructions presented in \cite{suthar2022Intrinsic} is that of \textit{diffusion generator}. A fiber preserving map $G:TM\to \mathfrak{D}M$ is called a \textbf{\textit{diffusion generator}} if $\forall$ $Y\in TM$, $\widehat{G(Y)} = Y\otimes Y$. As per the diffusion generator approach, a semi-martingale $X_t\in M$ is given by the Intrinsic SDE, 
\begin{equation}
\label{equation:Intrinsic_SDE_diffusion_generator}
\mathbf{d}X_t = \left[V(X_t) + \dfrac{1}{2} \sum_{l = 1}^p G(\sigma_l(X_t))\right]dt + \sum_{l = 1}^p\sigma_l(X_t) dW^l_t,
\end{equation}
where $V\in\mathfrak{X}(M)$ and $\sigma_l\in \mathfrak{X}(M)$ for all $l\in \{1, 2, ..., p\}$.
In \cite{suthar2022Intrinsic}, we have also shown that the equivalent Stratonovich representation of the intrinsic SDE \ref{equation:Intrinsic_SDE_diffusion_generator} is given by
\begin{equation}
\label{equation:Stratonovich for Intrinsic}
\delta X_t = \left[V(X_t) + \dfrac{1}{2} \sum_{l = 1}^p \nabla^G_S(\sigma_l(X_t))\right]dt + \sum_{l = 1}^p\sigma_l(X_t)\circ dW^l_t,
\end{equation}
where $\nabla^G_S = (G - G_S)\in \Omega^1(M;TM)$ is a tangent bundle valued 1-form, and $G_S$ is a diffusion generator called the \textbf{\textit{Stratonovich diffusion generator}}; and the equivalent  standard Ito representation is given by
\begin{equation}
\label{equation: standard Ito for Intrinsic}
dX_t = \left[V(X_t) + \dfrac{1}{2} \sum_{l = 1}^p \nabla^G_I(\sigma_l(X_t))\right]dt + \sum_{l = 1}^p\sigma_l(X_t) dW^l_t,
\end{equation}
where $\nabla^G_I = (G - G_I)\in \Omega^1(M;TM)$ is a tangent bundle valued 1-form, and $G_I$ is a diffusion generator called the \textbf{\textit{Ito diffusion generator}}. However, in order to define $G_I$, the manifold $M$ must be equipped with a connection.
In short hand notation we can write this conversion as,
\begin{equation}
\label{equation:conversion_shorthand_stratonovich}
\delta X_t = \mathbf{d}X_t  - \dfrac{1}{2}\sum_{l = 1}^pG_S(\sigma_l(X_t)) dt,
\end{equation}
and
\begin{equation}
\label{equation:conversion_shorthand_stratonovich}
d X_t = \mathbf{d}X_t  - \dfrac{1}{2}\sum_{l = 1}^pG_I(\sigma_l(X_t)) dt
\end{equation}
respectively.

\section{Explicit Stochastic Differential-Algebraic Equations on manifolds}
\label{section:explicit SDAEs on mani}
In  \cite{suthar2021explicit}, we have defined a system of explicit Stochastic Differential Algebraic Equations in Euclidean spaces as
\begin{subequations}
\label{eq:SDAE1}
\begin{equation}
\label{eq:SDAE-SDE}
dx = f(x,u)dt + \sigma(x,u)dW_t,  
\end{equation}
\begin{equation}
\label{eq:SDAE-AE}
    h(x,u) = - \int_0^t\Gamma(x(s),u(s))dW_s;
\end{equation}
\end{subequations}
where $W_t = (W^1_t, W^2_t, ..., W^d_t)$ is a d-dimensional Wiener process, $f:\mathbb{R}^n\times \mathbb{R}^m\to \mathbb{R}^n$, and $\sigma:\mathbb{R}^n\times \mathbb{R}^m\to L(\mathbb{R}^d,\mathbb{R}^n)$, $h:\mathbb{R}^n\times \mathbb{R}^m \to \mathbb{R}^q$, $\Gamma:\mathbb{R}^n\times \mathbb{R}^m \to L(\mathbb{R}^d,\mathbb{R}^q)$, $x(t,\omega) \in \mathbb{R}^n$, and $u(t,\omega) \in \mathbb{R}^m$. In the same article we have also demonstrated that it is possible to convert a given SDAE with noisy constraint into an SDAE with noiseless constraint. Hence, for SDAEs on manifolds we can simply consider constraints without noise. A generic form of explicit SDAE on manifolds can be given as
\begin{subequations}
\label{equation:SDAE on manifold}
\begin{equation}
\label{equation:SDAE on manifold - eq}
\mathbf{d}X_t = \left[V(X_t,U_t) + \dfrac{1}{2} \sum_{l = 1}^d  {}^MG(\sigma_l(X_t,U_t))\right]dt + \sum_{l = 1}^d\sigma_l(X_t,U_t) dW^l_t,
\end{equation}
\begin{equation}
\label{equation:SDAE on manifold - constraint}
h(X_t,U_t) = p;
\end{equation}
\end{subequations}
where $X_t\in M$, $U_t\in N$, $p\in P$, $h:M\times N \to P$, $V:M\times N\to TM$ such that $\tau_M(V(\cdot, u)) = Id_M$ for all $u\in N$, $\sigma_l:M\times N\to TM$ such that $\tau_M(\sigma_l(\cdot, u)) = Id_M$ for all $u\in N$ and $l\in \{1, 2, ..., d\}$, ${}^MG$ is a diffusion generator on $M$; and $M$, $N$, and $P$ are smooth manifolds of dimension $n$, $m$, and $q$ respectively.

The usual way of ensuring that the constraint is satisfied is by ensuring that the differential evolution of the constraint is zero. Hence, if we consider the Stratonovich differential of the constraint equation \ref{equation:SDAE on manifold - constraint} to be zero, then we obtain
\begin{equation}
D_1h\circ \delta X_t + D_2h \circ\delta U_t = 0.
\end{equation}
Therefore, as long as $D_2h$ is invertible, we get
\begin{equation}
\delta U_t = -(D_2h)^{-1}D_1h \circ \delta X_t.
\end{equation}
Using equation \ref{equation:conversion_shorthand_stratonovich},
\begin{equation}
\delta U_t = -\left((D_2h)^{-1}D_1h \right) \mathbf{d}X_t  + \dfrac{1}{2}(D_2h)^{-1}D_1h \sum_{l = 1}^p{}^MG_S(\sigma_l(X_t,U_t)) dt.
\end{equation}
Therefore, following the nomenclature from \cite{suthar2021explicit}, we give the following definitions.
\begin{definition}
An SDAE in form of equation \ref{equation:SDAE on manifold} will be called an \textbf{\textit{index 1 SDAE}}, if $D_2(h)(x,u)$ is invertible. 
\end{definition}
\begin{lemma}
\label{lem:SDEforSDAE}
For an index 1 SDAE in the form of equation \ref{equation:SDAE on manifold}, the local solution (local in time) is given by the following Stratonovich SDE.
\begin{equation}
\begin{aligned}
\delta \begin{bmatrix}X_t\\ U_t \end{bmatrix} = \begin{bmatrix}
V(X_t,U_t) + \dfrac{1}{2} \sum_{l = 1}^d  \nabla^{{}^MG}_S (\sigma_l(X_t,U_t))\\ -\left((D_2h)^{-1}D_1h \right) \left(V(X_t,U_t) +\dfrac{1}{2} \sum_{l = 1}^d  \nabla^{{}^MG}_S (\sigma_l(X_t,U_t))\right)\end{bmatrix} dt\\ + \sum_{l = 1}^d \begin{bmatrix}\sigma_l(X_t,U_t)\\ -\left((D_2h)^{-1}D_1h \right)\sigma_l(X_t,U_t) \end{bmatrix}\circ dW^l_t.
\end{aligned}
\end{equation}
\end{lemma}
\noindent The problem starts when $D_2h$ is not invertible. A clear case in which $D_2h$ is not invertible is the one in which $h$ is not a function of the algebraic variable. We will call such equations as \textit{high index SDAEs}. 
\begin{definition}
\label{def:high_index}
An SDAE that does not admit algebraic variable in the algebraic equation (i.e. when $h$ is not a function of the algebraic variable) will be called {\textbf{high index}} SDAE. A generic form for high index SDAEs is given by,
\begin{subequations}
\label{equation: high index SDAE}
\begin{equation}
\label{equation: high index SDAE - SDE}
\mathbf{d}X_t = \left[V(X_t,U_t) + \dfrac{1}{2} \sum_{l = 1}^d  {}^MG(\sigma_l(X_t,U_t))\right]dt + \sum_{l = 1}^d\sigma_l(X_t,U_t) dW^l_t,
\end{equation}
\begin{equation}
\label{equation: high index SDAE - AE}
p = h(X_t)
\end{equation}
\end{subequations}
where $X_t\in M$, $U_t\in N$, $p\in P$, $h:M \to P$, $V:M\times N\to TM$ such that $\tau_M(V(\cdot, u)) = Id_M$ for all $u\in N$, $\sigma_l:M\times N\to TM$ such that $\tau_M(\sigma_l(\cdot, u)) = Id_M$ for all $u\in N$ and $l\in \{1, 2, ..., d\}$, ${}^MG$ is a diffusion generator on $M$; and $M$, $N$, and $P$ are smooth manifolds of dimension $n$, $m$, and $q$ respectively.
\end{definition}

From this definition it is clear that the solution in lemma \ref{lem:SDEforSDAE} is not valid for high index SDAEs. The following theorem from \cite{suthar2021explicit}, states the necessary condition for existence of a local solution of SDAE in Euclidean spaces.
\begin{theorem}
\label{corollary:subbundle}
Consider a high index SDAE given as 
\begin{subequations}
\label{eq:SDAE2}
\begin{equation}
\label{eq:SDAE2-SDE}
dx = f(x,u)dt + \sum_{l = 1}^d\sigma_l(x,u)dW^l_t,  
\end{equation}
\begin{equation}
\label{eq:SDAE2-AE}
0 = h(x);
\end{equation}
\end{subequations}
with $x(t)\in \mathbb{R}^n$, $u(t)\in \mathbb{R}^m$, $W_t = (W^1_t, W^2_t, ..., W^d_t)$ as d-dimensional Wiener process, $f:\mathbb{R}^n\times \mathbb{R}^m\to \mathbb{R}^n$, $\sigma_l:\mathbb{R}^n\times \mathbb{R}^m\to \mathbb{R}^n$, $h:\mathbb{R}^n\to \mathbb{R}^q$. Let $h(x)$ be $\mathcal{C}^3$, all the vector fields be $\mathcal{C}^1$, and assume that there is a solution of the SDAE. Then
\[ span(\sigma_1 (x(t,\omega),u(t,\omega)), \sigma_2 (x(t,\omega),u(t,\omega)), ..., \sigma_d (x(t,\omega),u(t,\omega)))\subset Ker(D_{x(t,\omega)}h).\]
\end{theorem}
\noindent On manifolds, this theorem has a direct extension and can be stated as the following.
\begin{proposition}
\label{corollary:subbundle}
Consider a high index SDAE given by equation \ref{equation: high index SDAE}. Assume that there is a semi-martingale on $M\times N$ that satisfies the SDAE. Then
\[ span(\sigma_1 (X_t,U_t), \sigma_2 (X_t,U_t), ..., \sigma_d (X_t,U_t))\subset Ker(D_{X_t}h).\]
\end{proposition}
\begin{proof}
Since, $h(X_t) = p$, $D_{X_t}h\circ \delta X_t  = 0$. Therefore, $\delta X_t \in Ker(D_{X_t}h)$. As
\begin{equation}
\delta X_t = \left[V(X_t,U_t) + \dfrac{1}{2} \sum_{l = 1}^p \nabla^{{}^MG}_S(\sigma_l(X_t,U_t))\right]dt + \sum_{l = 1}^p\sigma_l(X_t,U_t)\circ dW^l_t,
\end{equation}
and $\delta X_t \in Ker(D_{X_t}h)$ for all sample paths,
\[ span(\sigma_1 (X_t,U_t), \sigma_2 (X_t,U_t), ..., \sigma_d (X_t,U_t))\subset Ker(D_{X_t}h).\]
Hence, proved.
\end{proof}
\noindent Clearly, if the vector fields $\sigma_l\in \mathfrak{X}(M)$ are not functions of the algebraic variable and if \[span(\sigma_1 (x), \sigma_2 (x), ..., \sigma_d (x))\not\subset Ker(D_{x}h),\] then there is no semi-martingale $(X_t,U_t)\in M\times N$ that satisfies the SDAE. Following the nomenclature from \cite{suthar2021explicit}, we will call such high index SDAEs as \textbf{\textit{ill-posed SDAEs}}.
\begin{corollary_p}
Ill-posed SDAEs have no solution.
\end{corollary_p}
\section{Approximate solution of high index SDAEs.}
\label{section:approximate solution}
We know that for any function $f\in \mathfrak{F}(M)$ and a semi-martingale $X_t$,
\[f(X_t) - f(X_0)= \int_0^t(\textbf{d}X_s[f])(X_s).\]
If $\mathbf{d}X_t = \left[V(X_t) + \dfrac{1}{2} \sum_{l = 1}^d G(\sigma_l(X_t))\right]dt + \sum_{l = 1}^d\sigma_l(X_t) dW^l_t,$
then 
 \[f(X_t) = f(X_0) + \int_0^t\left[\left(V+ \dfrac{1}{2} \sum_{l = 1}^d G(\sigma_l)\right)[f]\right](X_s) ds + \sum_{l = 1}^d \int_0^t(\sigma_l[f] )(X_s) dW^l_s.\]
Since this is an Ito SDE on a real line, if the initial condition $X_0$ is not random, then
\begin{equation}
\label{equation:expectation of fX_t}
\mathbb{E}(f(X_{t})) = f(X_0) + \mathbb{E}\left[\int_0^{t} \left( V+ \dfrac{1}{2} \sum_{l = 1}^d G(\sigma_l)\right)[f](X_s) ds\right].
\end{equation}

For stability analysis we will assume that $f\in \mathfrak{F}(M)$ is a positive semi-definite function with $f(y) = 0$, $\left[ V[f] + \dfrac{1}{2} \sum_{l = 1}^d G(\sigma_l)[f]\right](x) < 0$ $\forall$ $x\in M\setminus y$, and $\left[ V[f] + \dfrac{1}{2} \sum_{l = 1}^d G(\sigma_l)[f]\right](y) = 0$.

If $f\in \mathfrak{F}(M)$ were geodesically convex in $M$, then using Jensen's inequality we obtain $f(\mathbb{E} (X_{{t}}))\leq \mathbb{E} (f(X_{{t}}))$ $\forall$ $t>0$. Therefore,  $\mathbb{E}(X_t)\to y$ as $t\to \infty$, i.e. $X_t$ is asymptotically stable in the mean. Unfortunately, globally convex smooth functions may not exist e.g. in case of compact manifolds. Hence, global mean asymptotic stability analysis on manifolds is difficult. On the other hand if we have a locally convex function, say convex in $U\subset M$, then the probability that $X_t\not \in U$ for some $t>0$ is non-zero. Once the process $X_t$ is outside the set $U$, we lose the convexity of the function $f$ and hence, nothing can be concluded about the mean asymptotic stability of the process $X_t$. In case $M$ is a Riemannian manifold, if we considering $f$ as the distance function from point $y\in M$, i.e. $f = \Delta_y:M\to \mathbb{R}$, then the function will be geodesically convex upto the cut-locus of the point $y$. However, the function loses its differentiability on the cut-locus. In deterministic systems, the problem of avoidance of cut-locus has been considered in \cite{simha2018global}. However, in case of stochastic systems, due to randomness the process $X_t$ may hit the cut-locus of the point $y$. Since the problem of SDAE is similar to that of control problem (where algebraic variable is similar to the control input), if we could ensure that the probability of being in neighbourhood of the cut-locus is very low, then it is reasonable to consider the distance function from point $y$ as the candidate function.
\begin{definition}
Let $\mathfrak{C}_y\subset M$ denote the cut-locus of point $y$ on a Riemannian manifold $(M,g)$. If there exists a function $f:M\to \mathbb{R}$ such that $f$ is geodesically convex on $M\setminus \mathfrak{C}_y$ then we will call the function $f$ as \textit{\textbf{approximately convex}} for $y$.
\end{definition}

Consider the high index SDAE given by equation \ref{equation: high index SDAE},
\begin{subequations}
\begin{equation*}
\mathbf{d}X_t = \left[V(X_t,U_t) + \dfrac{1}{2} \sum_{l = 1}^d  {}^MG(\sigma_l(X_t,U_t))\right]dt + \sum_{l = 1}^d\sigma_l(X_t,U_t) dW^l_t,
\end{equation*}
\begin{equation*}
p = h(X_t);
\end{equation*}
\end{subequations}
where $X_t\in M$, $U_t\in N$, $p\in P$, $h:M \to P$, $V:M\times N\to TM$ such that $\tau_M(V(\cdot, u)) = Id_M$ for all $u\in N$, $\sigma_l:M\times N\to TM$ such that $\tau_M(\sigma_l(\cdot, u)) = Id_M$ for all $u\in N$ and $l\in \{1, 2, ..., d\}$, ${}^MG$ is a diffusion generator on $M$; and $M$, $N$, and $P$ are smooth manifolds of dimension $n$, $m$, and $q$ respectively. We will assume that the manifold $P$ is equipped with a Riemannian metric $g$.
Following the definitions introduced in \cite{suthar2021explicit}, for SDAEs on manifolds, we give the following definitions.
\begin{definition}
A high index SDAE will be said to have an \textbf{{m-solution}} if there exists semi-martingale $(X_t,U_t) \in M\times N$ such that $\mathbb{E}(h(X_t)) = p$.
\end{definition}

\begin{proposition}
Let the high index SDAE be given by equation \ref{equation: high index SDAE},
\begin{subequations}
\begin{equation*}
\mathbf{d}X_t = \left[V(X_t,U_t) + \dfrac{1}{2} \sum_{l = 1}^d  {}^MG(\sigma_l(X_t,U_t))\right]dt + \sum_{l = 1}^d\sigma_l(X_t,U_t) dW^l_t,
\end{equation*}
\begin{equation*}
p = h(X_t).
\end{equation*}
\end{subequations}
Let $P$ be a Riemannian manifold.  If there exists a positive semi-definite, and approximately convex function $f:P\to \mathbb{R}$ such that $f|_{P\setminus \overline{\mathfrak{C}^\epsilon_p}}\in \mathfrak{F}^3(P\setminus \overline{\mathfrak{C}^\epsilon_p})$ for $p\in P$ and a semi-martingale $(X_t,U_t)\in M\times N$ such that \[\left(D_{X_t}hV(X_t,U_t)+ \dfrac{1}{2} \sum_{l = 1}^d D_{X_t}h\nabla^{{}^MG}_S(\sigma_l)(X_t,U_t) + \dfrac{1}{2} \sum_{l = 1}^d{}^PG_S(D_{X_t}h\sigma_l)(X_t,U_t) \right)[f(h(X_t))] = 0,\]
\[f(p) = 0,\]
and
\[h(X_0) = p;\]
then the equation \ref{equation: high index SDAE} has an m-solution. $\overline{\mathfrak{C}^\epsilon_p}$ is the closure of the $\epsilon$ - neighbourhood of cut-locus of $p$, $\mathfrak{C}^\epsilon_p$.
\end{proposition}
\begin{proof}
As \[\delta X_t = \left[ V(X_t,U_t) + \dfrac{1}{2} \sum_{l = 1}^d  \nabla^{{}^MG}_S (\sigma_l(X_t,U_t))\right]dt + \sum_{l = 1}^d\sigma_l(X_t,U_t) dW^l_t,\]
\[\delta h(X_t) = \left[ D_xhV(X_t,U_t) + \dfrac{1}{2} \sum_{l = 1}^d D_xh \nabla^{{}^MG}_S (\sigma_l(X_t,U_t))\right]dt + \sum_{l = 1}^dD_xh\sigma_l(X_t,U_t) dW^l_t.\]
 \[\therefore\mathbf{d}h(X_t) = \left[ D_xhV(X_t,U_t) + \dfrac{1}{2} \sum_{l = 1}^d D_xh \nabla^{{}^MG}_S (\sigma_l(X_t,U_t)) + \dfrac{1}{2} \sum_{l = 1}^d{}^PG_S(D_{X_t}h\sigma_l) (X_t,U_t)\right]dt\]
 \[+ \sum_{l = 1}^dD_xh\sigma_l(X_t,U_t) dW^l_t.\]
From equation \ref{equation:expectation of fX_t}, we know that if $h(X_0) = p$ and $f(p) = 0$, then
 \[\mathbb{E}f(h(X_{t})) = \mathbb{E}\int_0^{t}\left[  D_xhV(X_s,U_s) + \dfrac{1}{2} \sum_{l = 1}^d D_xh \nabla^{{}^MG}_S (\sigma_l(X_s,U_s)) + \dfrac{1}{2} \sum_{l = 1}^d{}^PG_S(D_{X_s}h\sigma_l) (X_s,U_s)\right][f(h(X_s))] ds.\]
 \[\therefore \mathbb{E}f(h(X_{t})) = 0.\]
The $\mathcal{C}^3$ smoothness of $f:P\to \mathbb{R}$ in $P\setminus \overline{\mathfrak{C}^\epsilon_p}$ is required to ensure that the drift term $D_xhV(x,u) + \dfrac{1}{2} \sum_{l = 1}^d D_xh \nabla^{{}^MG}_S (\sigma_l(x,u)) + \dfrac{1}{2} \sum_{l = 1}^d{}^PG_S(D_{x}h\sigma_l) (x,u)$ is Lipschitz.
Using Jensen's inequality, $f(\mathbb{E}h(X_{t})) = 0$. Hence, $\mathbb{E}(h(X_{t})) = p$.
\end{proof}

\begin{corollary_p}
Let the high index SDAE be given by equation \ref{equation: high index SDAE}, with $P$ as a Riemannian manifold.  If there exists a semi-martingale $(X_t,U_t)\in M\times N$ such that \[\left(D_{X_t}hV(X_t,U_t)+ \dfrac{1}{2} \sum_{l = 1}^d D_{X_t}h\nabla^{{}^MG}_S(\sigma_l)(X_t,U_t) + \dfrac{1}{2} \sum_{l = 1}^d{}^PG_S(D_{X_t}h\sigma_l)(X_t,U_t) \right)[\Delta^4_p(h(X_t))] = 0,\]
and
\[h(X_0) = p;\]
then the equation \ref{equation: high index SDAE} has an m-solution. $\Delta_p(z):P\to \mathbb{R}$ is the distance of point $z\in P$ from point $p\in P$.  
\end{corollary_p}
\begin{proof}
$\Delta^4_p(h(X_t))$ is clearly a positive semi-definite function with $\Delta^4_p(p) = 0$. The function is also approximately convex  for $p\in P$. The $\mathcal{C}^3$ smoothness of $\Delta^4_p$ follows because $\Delta_p(z) = ||exp_p^{-1}(z)||$ where $||\cdot||$ is the 2-norm in $T_pP$.
\end{proof}

As discussed before, in order to ensure that the solution has low probability of hitting the cut-locus $\mathfrak{C}_p$, we must have a probability bound on the explosion of the solution.
\begin{definition}
We define an $\epsilon$ \textbf{\textit{{bounded m-solution}}}, with probability of explosion less than or equal to $\alpha$, as an m-solution $(X_t,U_t)\in M\times N$ such that for any $t>0$, $\mathbb{P}\left(\sup_{s\in [0,t]}(\Delta_p(h(X_s))>\epsilon\right)\leq \alpha$ for some $\alpha\in (0,1]$ and $\epsilon >0$, where $\Delta_p:M\to \mathbb{R}$ such that $\Delta_p(z)$ is the distance of point $z\in P$ from the point $p\in P$. We will call $X_t$ as \textbf{\textit{local bounded m-solution}} if for any $t\in (0,a]$ for $a<\infty$, $\mathbb{P}\left(\sup_{s\in [0,t]}(\Delta_p(h(X_s))>\epsilon\right)\leq \alpha$.
\end{definition}

\begin{proposition}
\label{proposition:bounded m-solution}
Let the high index SDAE be given by equation \ref{equation: high index SDAE} such that $P$ is a Riemannian manifold.  For every $\epsilon>0$, $\alpha\in (0,1]$, there exists $b>0$ such that if there exists a semi-martingale $(X_t,U_t)\in M\times N$ that ensures $h(X_0) = p$, 
\begin{multline}
\left(D_{X_t}hV(X_t,U_t)+ \dfrac{1}{2} \sum_{l = 1}^d D_{X_t}h\nabla^{{}^MG}_S(\sigma_l)(X_t,U_t) + \dfrac{1}{2} \sum_{l = 1}^d{}^PG_S(D_{X_t}h\sigma_l)(X_t,U_t) \right)[\Delta^4_p(h(X_t))] \\=- b \Delta^4_p(h(X_t)),
\end{multline}
and $\lambda =\sup_{(x,u)\in U\times V} \sum_{l = 1}^d\left( (D_{x}h)\sigma_l(x,u)[\Delta^4_p(h(x))]\right)^2<\infty$ for $U\times V = \{(X_t(\omega),U_t(\omega))\in M\times N| t\in (0,a], \omega\in \Omega\}$, where $\Omega$ is the sample space;
then for any  $t\in (0,a]$, for some $a>0$ that makes $(0,a]$ an interval of existence,
\[\mathbb{P}\left(\sup_{s\in [0,t]}(\Delta_p(h(X_s))>\epsilon\right)\leq \alpha\]
if $b > \lambda/(\alpha \sqrt{\epsilon})$. In other words, the equation \ref{equation: high index SDAE} has a local bounded m-solution for $t\in (0,a]$.
\end{proposition}
\begin{proof}
From the previous proposition we know that
\[f(h(X_{t})) = \int_0^{t}\left[ D_xhV(X_t,U_t) + \dfrac{1}{2} \sum_{l = 1}^d D_xh \nabla^{{}^MG}_S (\sigma_l(X_t,U_t)) + \dfrac{1}{2} \sum_{l = 1}^d{}^PG_S(D_{X_t}h\sigma_l) (X_t,U_t)\right][f(h(X_s))] ds\]
\[+\int_0^t\sum_{l = 1}^d(D_{X_t}h)\sigma_l(X_s,U_s)[f(h(X_s))] dW^l_s.\]
Therefore, for $f = \Delta^4_p$, if
\begin{multline*}
\left(D_{X_t}hV(X_t,U_t)+ \dfrac{1}{2} \sum_{l = 1}^d D_{X_t}h\nabla^{{}^MG}_S(\sigma_l)(X_t,U_t) + \dfrac{1}{2} \sum_{l = 1}^d{}^PG_S(D_{X_t}h\sigma_l)(X_t,U_t) \right)[\Delta^4_p(h(X_t))] \\=- b \Delta^4_p(h(X_t)),
\end{multline*}
then 
\[d\Delta^4_p(h(X_t)) = -b \Delta^4_p(h(X_t)) dt + \sum_{l = 1}^d(D_{X_t}h)\sigma_l(X_t,U_t)[\Delta^4_p(h(X_t))] dW^l_t.\]
We know that as $h(X_0) = p$,
\[\Delta^4_p(h(X_t)) = e^{-bt}\int_0^te^{bs}\sum_{l = 1}^d(D_{X_s}h)\sigma_l(X_s,U_s)[f(h(X_s))] dW^l_s.\]
As the integral $\int_0^te^{bs}\sum_{l = 1}^d(D_{X_s}h)\sigma_l(X_s,U_s)[f(h(X_s))] dW^l_s$ is a martingale, the Doob's martingale inequality gives,
\[\mathbb{P}\left(\sup_{s\in [0,t]}(e^{bs}\Delta^4_p(h(X_s))>k)\right) \leq \dfrac{\mathbb{E}(e^{2bt}\Delta^8_p(h(X_t)))}{k^2}.\]
\[\therefore \mathbb{P}\left(\sup_{s\in [0,t]}(e^{bs}\Delta^4_p(h(X_s))>k)\right) =  \mathbb{P}\left(\sup_{s\in [0,t]}(\Delta^4_p(h(X_s))>\epsilon)\right) \leq \dfrac{\mathbb{E}(e^{2bt}\Delta^8_p(h(X_t)))}{e^{2bt}\epsilon^2},\]
where $k = e^{bt}\epsilon$.
\[\therefore \mathbb{P}\left(\sup_{s\in [0,t]}(\Delta^4_p(h(X_s))>\epsilon)\right) \leq \dfrac{\mathbb{E}\left(\int_0^te^{bs}\sum_{l = 1}^d(D_{X_s}h)\sigma_l(X_s,U_s)[f(h(X_s))] dW^l_s\right)^2}{e^{2bt}\epsilon^2}.\]
Hence using Ito-isometry property,
\[\mathbb{P}\left(\sup_{s\in [0,t]}(\Delta^4_p(h(X_s))>\epsilon)\right) \leq \dfrac{\mathbb{E}\int_0^te^{2bs}\sum_{l = 1}^d\left((D_{X_s}h)\sigma_l(X_s,U_s)[\Delta^4_p(h(X_s))]\right)^2 ds}{e^{2bt}\epsilon^2}.\]
Therefore, if $\lambda =\sup_{(x,u)\in U\times V} \sum_{l = 1}^d\left( (D_{x}h)\sigma_l(x,u)[\Delta^4_p(h(x))]\right)^2<\infty$, then
\[\dfrac{\mathbb{E}\int_0^te^{2bs}\sum_{l = 1}^d\left((D_{X_s}h)\sigma_l(X_s,U_s)[\Delta^4_p(h(X_s))]\right)^2 ds}{e^{2bt}\epsilon^2}\leq \lambda(1-e^{-2bt})/2b\epsilon^2\leq \lambda /2b\epsilon^2\quad \forall\quad t\in (0,a].\]
Since, this is true for any $0<t\leq a$, if $b > \lambda/(\epsilon^2 \alpha)$, then $\alpha< \dfrac{\lambda}{2b\epsilon^2}$ and hence
\[\mathbb{P}\left(\sup_{s\in [0,t]}(\Delta^4_p(h(X_s))>\epsilon)\right)\leq \alpha.\]
In other words if $b > \lambda/(\alpha \sqrt{\epsilon})$, then for any $0<t\leq a$
\[\mathbb{P}\left(\sup_{s\in [0,t]}(\Delta_p(h(X_s))>\epsilon)\right)\leq \alpha.\]
\end{proof}
This proposition proves that for large enough $b>0$, we can ensure that a semi-martingale $(X_t,U_t)\in M\times N$ that satisfies
\begin{subequations}
\begin{equation}
\mathbf{d}X_t = \left[V(X_t,U_t) + \dfrac{1}{2} \sum_{l = 1}^d  {}^MG(\sigma_l(X_t,U_t))\right] dt + \sum_{l = 1}^d\sigma_l(X_t,U_t) dW^l_t,
\end{equation}
\begin{multline}
\left(D_{X_t}hV(X_t,U_t)+ \dfrac{1}{2} \sum_{l = 1}^d D_{X_t}h\nabla^{{}^MG}_S(\sigma_l)(X_t,U_t) + \dfrac{1}{2} \sum_{l = 1}^d{}^PG_S(D_{X_t}h\sigma_l)(X_t,U_t) \right)[\Delta^4_p(h(X_t))]\\ =- b \Delta^4_p(h(X_t)),
\end{multline}
\end{subequations}
is a local bounded m-solution for the high index SDAE \ref{equation: high index SDAE}.
\begin{corollary_p}
\label{corollary:cutlocus_issue}
Let $\mathfrak{I}_p$ be the injectivity radius of point $p\in P$ and let $\alpha\in (0,1]$. If there exists $\epsilon>0$ such that $\epsilon <\mathfrak{I}_p$, and if there exists a semi-martingale $(X_t,U_t)\in M\times N$ such that it satisfies
\begin{subequations}
\begin{equation}
\mathbf{d}X_t = \left[V(X_t,U_t) + \dfrac{1}{2} \sum_{l = 1}^d  {}^MG(\sigma_l(X_t,U_t))\right]dt + \sum_{l = 1}^d\sigma_l(X_t,U_t) dW^l_t \; and
\end{equation}
\begin{multline}
\left(D_{X_t}hV(X_t,U_t)+ \dfrac{1}{2} \sum_{l = 1}^d D_{X_t}h\nabla^{{}^MG}_S(\sigma_l)(X_t,U_t) + \dfrac{1}{2} \sum_{l = 1}^d{}^PG_S(D_{X_t}h\sigma_l)(X_t,U_t) \right)[\Delta^4_p(h(X_t))]\\ =- b \Delta^4_p(h(X_t))
\end{multline}
\end{subequations}
for large enough $b>0$, then $(X_t,U_t)\in M\times N$ is a bounded m-solution for the high index SDAE \ref{equation: high index SDAE} such that $\mathbb{P}(h(X_t)\in \mathfrak{C}_p\; for\; some\; t\in(0,a])< \alpha$; where $(0,a]$ is the interval of existence.
\end{corollary_p}
\begin{proof}
$\mathfrak{I}_p$ is the injectivity radius from point $p\in P$. As $\mathfrak{C}_p\subset P\setminus\{z| \Delta_p(z)<\mathfrak{I}_p\}$,
\[\mathbb{P}(h(X_t)\in \mathfrak{C}_p\; for\; some\; t\in(0,a])<\mathbb{P}(\sup_{t\in(0,a]}\Delta_p(h(X_t)) > \mathfrak{I}_p).\]
As $\mathbb{P}((\sup_{t\in(0,a]}\Delta_p(h(X_t)))>\mathfrak{I}_p)\leq \mathbb{P}((\sup_{t\in(0,a]}\Delta_p(h(X_t)))>\epsilon)$ for all $\epsilon < \mathfrak{I}_p$, if $\mathbb{P}((\sup_{t\in(0,a]}\Delta_p(h(X_t)))>\epsilon)\leq \alpha$ and $\epsilon<\mathfrak{I}_p$, then $\mathbb{P}(h(X_t)\in \mathfrak{C}_p\; for\; some\; t\in(0,a])< \alpha$ is true.
\end{proof}
In proposition \ref{proposition:bounded m-solution} we have shown that if there exists a semi-martingale $(X_t,U_t)\in M\times N$ that ensures $h(X_0) = p$, and 
\begin{multline}
\left(D_{X_t}hV(X_t,U_t)+ \dfrac{1}{2} \sum_{l = 1}^d D_{X_t}h\nabla^{{}^MG}_S(\sigma_l)(X_t,U_t) + \dfrac{1}{2} \sum_{l = 1}^d{}^PG_S(D_{X_t}h\sigma_l)(X_t,U_t) \right)[\Delta^4_p(h(X_t))]\\ =- b \Delta^4_p(h(X_t));
\end{multline}
then for large enough $b>0$, \[\mathbb{P}((\sup_{t\in (0,a]}\Delta_p(h(X_t)))>\epsilon)\leq \alpha.\]
The choice of $0<b\in\mathbb{R}$, which in turn depends on the interval $(0,a]$.
In real world applications, in order to avoid the problem of hitting the cut-locus, we will choose $\alpha$ to be small enough (say $\alpha = 0.98$) and $\epsilon$ to be much smaller than $\mathfrak{I}_p$. Then from corollary \ref{corollary:cutlocus_issue}, we know that the probability of hitting the cut-locus of the point $p$ will be very low.
\section{Method for numerical computation of the approximate solution}
\label{section:numerical}
In previous section, we have demonstrated that if $Y:M\times N\to \mathbb{R}$ such that
\[Y(x,u) = b \Delta^4_p(h(x)) + \left(D_{x}hV(x,u)+ \dfrac{1}{2} \sum_{l = 1}^d D_{x}h\nabla^{{}^MG}_S(\sigma_l)(x,u) + \dfrac{1}{2} \sum_{l = 1}^d{}^PG_S(D_{x}h\sigma_l)(x,u) \right)[\Delta^4_p(h(x))],\] then we are seeking a solution $(X_t,U_t)\in M\times N$ such that $Y(X_t,U_t) = 0$. Since, $dim(N) = m\geq dim(\mathbb{R}) = 1$, implicit function theorem cannot be used to find $U_t$ unless $m = 1$.  This makes the computation of the process $U_t$ difficult. If we assume that $N$ is a Riemannian manifold with $g$ as the metric then we can use gradient descent to minimize $K = Y^2(X_t,\cdot):N\to \mathbb{R}$ i.e. 
the solution of $\dfrac{du}{ds} = -g^\sharp dK$ as $s\to \infty$ at every time $t\in (0,a]$. Off-course, the gradient descent has a problem of converging to a local minimum. Therefore, if we start with $U_0\in N$ that ensures $Y(X_0,U_0) = 0$ and consider the initial condition for  $\dot{u} = -g^\sharp dK$ as the solution of the previous time step, then we may be able to avoid local minima. From a computational perspective, implementation of gradient descent at every time step is not recommended as it will unavoidably add to the computational complexity. Moreover, as this solution is not a result of application of implicit function theorem, we do not know if the resulting process $U_t$ will be a semi-martingale. In the following proposition we address this dimensionality issue in the application of implicit function theorem.
\begin{proposition}
\label{proposition:high_index_approximate_solution}
Let the high index SDAE be given by equation \ref{equation: high index SDAE},
\begin{subequations}
\begin{equation}
\mathbf{d}X_t = \left[V(X_t,U_t) + \dfrac{1}{2} \sum_{l = 1}^d  {}^MG(\sigma_l(X_t,U_t))\right]dt + \sum_{l = 1}^d\sigma_l(X_t,U_t) dW^l_t,
\tag{\ref{equation: high index SDAE - SDE}}
\end{equation}
\begin{equation}
p = h(X_t), \tag{\ref{equation: high index SDAE - AE}}
\end{equation}
\end{subequations}
where $X_t\in M$, $U_t\in N$ with $N$ as a Riemannian manifold with metric $g$, $p\in P$, $h:M \to P$, $V:M\times N\to TM$ such that $\tau_M(V(\cdot, u)) = Id_M$ for all $u\in N$, $\sigma_l:M\times N\to TM$ such that $\tau_M(\sigma_l(\cdot, u)) = Id_M$ for all $u\in N$ and $l\in \{1, 2, ..., d\}$, ${}^MG$ is a diffusion generator on $M$; and $M$, $N$, and $P$ are smooth manifolds of dimension $n$, $m$, and $q$ respectively. Let $P$ be a Riemannian manifold.  Suppose there exists a semi-martingale $(X_t,U_t)\in M\times N$ that satisfies
\begin{subequations}
\label{equation:solution_bounded-m}
\begin{equation}
\mathbf{d}X_t = \left[V(X_t,U_t) + \dfrac{1}{2} \sum_{l = 1}^d  {}^MG(\sigma_l(X_t,U_t))\right]dt + \sum_{l = 1}^d\sigma_l(X_t,U_t) dW^l_t,
\end{equation}
\begin{equation}
\delta U_t = \alpha_0(X_t,U_t) a(X_t,U_t) dt + \sum_{l = 1}^d \alpha_l(X_t,U_t)a(X_t,U_t) dW^l_t,
\end{equation}
\end{subequations}
where
\[\alpha_0(X_t,U_t) = \dfrac{-D_1Y \left[V(X_t,U_t) + \dfrac{1}{2} \sum_{l = 1}^d  \nabla^{{}^MG}_S(\sigma_l(X_t,U_t))\right]}{D_2Y a(X_t,U_t)},\]
\[\alpha_l(X_t,U_t) = \dfrac{-D_1Y\sigma_l(X_t,U_t)}{D_2Y a(X_t,U_t)},\]
\[a\not \in Ker(D_2Y),\]
and $Y:M\times N\to \mathbb{R}$ such that
\[Y(x,u) = b \Delta^4_p(h(x)) + \left(D_{x}hV(x,u)+ \dfrac{1}{2} \sum_{l = 1}^d D_{x}h\nabla^{{}^MG}_S(\sigma_l)(x,u) + \dfrac{1}{2} \sum_{l = 1}^d{}^PG_S(D_{x}h\sigma_l)(x,u) \right)[\Delta^4_p(h(x))],\]
where $b > \lambda/(\alpha \sqrt{\epsilon})$. If there exists $U_0\in N$ such that $Y(X_0,U_0) = 0$ and if $D_2Y(X_t,U_t)\not = 0$, then $(X_t,U_t)$ is a bounded m-solution for the high index SDAE \ref{equation: high index SDAE} such that $\mathbb{P}((\sup_{t\in(0,a]}\Delta_p(h(X_t)))> \epsilon)< \alpha$ and $\epsilon$ is much smaller than $\mathfrak{I}_p$.
\end{proposition}
\begin{proof}
As $Y:M\times N\to \mathbb{R}$ such that
\[Y(x,u) = b \Delta^4_p(h(x)) + \left(D_{x}hV(x,u)+ \dfrac{1}{2} \sum_{l = 1}^d D_{x}h\nabla^{{}^MG}_S(\sigma_l)(x,u) + \dfrac{1}{2} \sum_{l = 1}^d{}^PG_S(D_{x}h\sigma_l)(x,u) \right)[\Delta^4_p(h(x))],\]
from proposition \ref{proposition:bounded m-solution}, we know that if there exists a semi-martingale $(X_t,U_t)\in M\times N$ such that it satisfies $Y(X_t,U_t) = 0$, then $\mathbb{P}((\sup_{t\in(0,a]}\Delta_p(h(X_t)))> \epsilon)< \alpha$ for large enough $b>0$.
If we assume that
\[\delta U_t = \alpha_0(X_t,U_t) a(X_t,U_t) dt + \sum_{l = 1}^d \alpha_l(X_t,U_t)a(X_t,U_t) dW^l_t,\]
for some $\lambda_i(X_t,U_t) \in \mathbb{R}$ for $i\in \{0, 1, ..., d\}$, and $a(x,\cdot)\in \mathfrak{X}(N)$; then Stratonovich differentiation of $Y(X_t,U_t)$ gives 
\[D_1Y \left[V(X_t,U_t) + \dfrac{1}{2} \sum_{l = 1}^d  \nabla^{{}^MG}_S(\sigma_l(X_t,U_t))\right]dt + \sum_{l = 1}^d D_1Y\sigma_l(X_t,U_t) dW^l_t\]
\[+ \alpha_0(X_t,U_t) D_2Y a(X_t,U_t) dt +  \sum_{l = 1}^d \alpha_l(X_t,U_t)D_2Y a(X_t,U_t) dW^l_t = 0.\]
If $a\not \in Ker(D_2Y)$, then
\[\alpha_0(X_t,U_t) = \dfrac{-D_1Y \left[V(X_t,U_t) + \dfrac{1}{2} \sum_{l = 1}^d  \nabla^{{}^MG}_S(\sigma_l(X_t,U_t))\right]}{D_2Y a(X_t,U_t)},\]
and
\[\alpha_l(X_t,U_t) = \dfrac{-D_1Y\sigma_l(X_t,U_t)}{D_2Y a(X_t,U_t)}\]
will ensure that $Y(X_t,U_t) = 0$ for any $U_0$ that satisfies $Y(X_0,U_0) = 0$.
\end{proof}
\subsection{Algorithm}
The computation of the initial condition $U_0$ is still a problem. If $N$ is a Riemannian manifold then we can use gradient descent to compute $U_0$. In any case, we will assume that we know the initial condition $U_0$ that ensures $Y(X_0,U_0) = 0$. Now, in proposition \ref{proposition:high_index_approximate_solution}, we have presented a method of computing the bounded m-solution by simply solving an SDE. However, this method depends on choosing a vector field $a(x,\cdot )\in \mathfrak{X}(N)$ such that $a\not\in Ker(D_2Y)$ and choosing $b>\dfrac{\lambda}{\alpha\sqrt{\epsilon}}$. Therefore, the problem is two folds:
\begin{itemize}
\item Computing $\lambda =\sup_{(x,u)\in U\times V} \sum_{l = 1}^d\left( (D_{x}h)\sigma_l(x,u)[\Delta^4_p(h(x))]\right)^2$ for $U\times V = \{(X_t(\omega),U_t(\omega))\in M\times N| t\in (0,a], \omega\in \Omega\}$, where $\Omega$ is the sample space. This is necessary to choose $b>\dfrac{\lambda}{\alpha\sqrt{\epsilon}}$.
\item Choosing the vector field $a\in \mathfrak{X}(M)$ such that $a\not\in Ker(D_2Y)$. 
\end{itemize}
The computation of $\lambda =\sup_{(x,u)\in U\times V} \sum_{l = 1}^d\left( (D_{x}h)\sigma_l(x,u)[\Delta^4_p(h(x))]\right)^2$ can be done using Monte-Carlo approach. However, since we are only interested in choosing large enough $b$ that ensures \[\mathbb{P}\left(\sup_{s\in [0,t]}(\Delta_p(h(X_s))>\epsilon\right)\leq \alpha,\]
algorithmically, we can double the current value of $b$ whenever $\Delta_p(h(X_t))>\epsilon$. In this way, we avoid the computation of $\lambda$.
Now for choosing the vector field $a\not\in Ker(D_2Y)$. Let $a(x,u) = g^\sharp_u D_2Y(x,u)$. Therefore, as long as $D_2Y\not =0$, we get the required solution by solving the following SDE:
\begin{subequations}
\label{equation:numerical_bounded-m-sol}
\begin{equation}
\label{equation:numerical_bounded-m-sol:SDEfor X}
\mathbf{d}X_t = \left[V(X_t,U_t) + \dfrac{1}{2} \sum_{l = 1}^d  {}^MG(\sigma_l(X_t,U_t))\right]dt + \sum_{l = 1}^d\sigma_l(X_t,U_t) dW^l_t,
\end{equation}
\begin{multline}
\delta U_t = -D_1Y \left[V(X_t,U_t) + \dfrac{1}{2} \sum_{l = 1}^d  \nabla^{{}^MG}_S(\sigma_l(X_t,U_t))\right]\dfrac{g^\sharp_{U_t}D_2Y(X_t,U_t)}{g(D_2Y,D_2Y)} dt\\ - \sum_{l = 1}^d {D_1Y\sigma_l(X_t,U_t)}\dfrac{g^\sharp_{U_t}D_2Y(X_t,U_t)}{g(D_2Y,D_2Y)} dW^l_t.
\end{multline}
\end{subequations}
This method is summarized in algorithm \ref{alg:bounded-sol}.
\begin{algorithm}
\hrule \vspace{1em}
{Numerical computation of bounded m-solution when $D_2Y(X_t,U_t) \not = 0$}\\
\label{alg:bounded-sol}
\hrule
\begin{algorithmic}[1]
\STATE{If $U_0$ is not available, then use gradient descent to compute $U_0$ such that $Y(X_0,U_0) = 0$;}
\STATE{$b = 1$;}
\STATE {$t_0 = 0$;}
\STATE{Choose $\delta t> 0$;}
\FOR{$i = 1$ \TO $i = T/\delta t$}
\STATE{Solve the SDE \ref{equation:numerical_bounded-m-sol} upto time $\delta t$ with initial condition as $(X_{t_{i-1}},U_{t_{i-1}})$;}
\STATE{$t_i = t_{i-1} + \delta t$;}
\WHILE{$\Delta_p(h(X_{t_i}))>\epsilon$}
\STATE{$b = 2b$;}
\STATE{Solve the SDE \ref{equation:numerical_bounded-m-sol} upto time $\delta t$ with initial condition as $(X_{t_{i-1}},U_{t_{i-1}})$;}
\STATE{$t_i = t_{i-1} + \delta t$;}
\ENDWHILE
\ENDFOR
\end{algorithmic}
\hrule
\end{algorithm}	
Offcourse, the method fails when $D_2Y(X_t,U_t) = 0$, in which case we will return to the using gradient descent in pseudo time at every time step. However, in this case we will not know if the resulting solution is a semi-martingale or not. Nonetheless, at every time step we will have a point $U_t\in N$ that ensures $Y(X_t,U_t) = 0$. We will call such a solution as \textit{approximate bounded m-solution}. This can be summarized in form of algorithm \ref{alg:approximate-bounded-sol}.
\begin{algorithm}
\hrule \vspace{1em}
{Numerical computation of approximate bounded m-solution.}\\
\label{alg:approximate-bounded-sol}
\hrule
\begin{algorithmic}[1]
\STATE{If $U_0$ is not available, then use gradient descent to compute $U_0$ such that $Y(X_0,U_0) = 0$;}
\STATE{$b = 1$;}
\STATE {$t_0 = 0$;}
\STATE{Choose $h> 0$;}
\STATE{Choose a positive integer $N$;}
\STATE{$\delta t = N*h$;}
\FOR{$i = 1$ \TO $i = T/\delta t$}
\STATE{$(X_{t^{0}},U_{t^{0}}) = (X_{t_{i-1}},U_{t_{i-1}})$;}
\FOR{$j = 1$ \TO $j = N$}
\IF{$D_2Y(X_{t^0},U_{t^0})\not =0$}
\STATE{Numerical time stepping for SDE \ref{equation:numerical_bounded-m-sol} with initial condition as $(X_{t^{j-1}},U_{t^{j-1}})$ and time step of $h$;}
\STATE{$t^j = t^{j-1} + h$;}
\ELSE
\STATE{Numerical time stepping for SDE \ref{equation:numerical_bounded-m-sol:SDEfor X} with initial condition as $X_{t^{j-1}}$, time step of $h$, and constant value of $U_t = U_{t^{j-1}}$.}
\STATE{$t^j = t^{j-1} + h$;}
\STATE{Use gradient descent to compute $U_{t^j}$ such that $Y(X_{t^j},U_{t^j}) = 0$, with initial guess as $U_{t^{j-1}}$;}
\ENDIF
\ENDFOR
\STATE{$t_i = t_{i-1} + \delta t$;}
\STATE{$(X_{t_{i}},U_{t_{i}}) = (X_{t^{N}},U_{t^{N}})$;}
\WHILE{$\Delta_p(h(X_{t_i}))>\epsilon$}
\STATE{$b = 2b$;}
\STATE{$(X_{t^{0}},U_{t^{0}}) = (X_{t_{i-1}},U_{t_{i-1}})$;}
\FOR{$j = 1$ \TO $j = N$}
\IF{$D_2Y(X_{t^0},U_{t^0})\not =0$}
\STATE{Numerical time stepping for SDE \ref{equation:numerical_bounded-m-sol} with initial condition as $(X_{t^{j-1}},U_{t^{j-1}})$ and time step of $h$;}
\STATE{$t^j = t^{j-1} + h$;}
\ELSE
\STATE{Numerical time stepping for SDE \ref{equation:numerical_bounded-m-sol:SDEfor X} with initial condition as $X_{t^{j-1}}$, time step of $h$, and constant value of $U_t = U_{t^{j-1}}$.}
\STATE{$t^j = t^{j-1} + h$;}
\STATE{Use gradient descent to compute $U_{t^j}$ such that $Y(X_{t^j},U_{t^j}) = 0$, with initial guess as $U_{t^{j-1}}$;}
\ENDIF
\ENDFOR
\STATE{$t_i = t_{i-1} + \delta t$;}
\STATE{$(X_{t_{i}},U_{t_{i}}) = (X_{t^{N}},U_{t^{N}})$;}
\ENDWHILE
\ENDFOR
\end{algorithmic}
\hrule
\end{algorithm}


