Stochastic quantization approach to modeling of constrained Ito processes in discrete time

TL;DR

This paper introduces a novel stochastic quantization method for constrained Ito processes in discrete time, enabling the modeling of nonlinear dynamics with potential terms that enforce constraints, applicable to engineering problems.

Contribution

It proposes a new approach to incorporate potential terms into stochastic quantization of general Ito processes without assumptions on the drift, linking it to Kalman filtering and constrained dynamics.

Findings

The method preserves the probability distribution role of the path integral.

Derived difference equations for expectation and covariance evolution.

Numerical results show effective constraint of nonlinear dynamics via potential terms.

Abstract

Stochastic quantization in physics has been considered to provide a path integral representation of a probability distribution for Ito processes. It has been indicated that the stochastic quantization can involve a potential term, if the Ito process is limited to Langevin equation. In this paper, in order to apply the stochastic quantization to engineering problems, we propose a novel method to incorporate a potential term into stochastic quantization of the general Ito process. This method indicates that weighted distribution gives rise to the potential term for the discrete-time path integral and preserves the role of the path integral as the probability distribution, without making any assumptions on the drift term. A second order approximation on the stochastic fluctuations for the path integral gives difference equations which represent the time evolution of expectation value and covariance matrix for the stochastic processes. The difference equations explicitly derive Extended Kalman Filter and models on the constrained Ito processes by the identification of the potential function with a penalty or barrier function. The numerical simulations of the constrained stochastic systems show that the potential term can constrain the nonlinear dynamics towards a minimum or a decreasing direction of the potential function.

Equations

$$1$$

$$\frac{{\mathbf{x}}_{t+\Delta t}-{\mathbf{x}}_t}{\Delta t}-{\mathbf{f}}_t({\mathbf{x}}_t)={\mathbf{E}}_t^{-1}{\mathbf{\eta }}_t$$

$$1$$

$$\displaystyle\qquad\qquad\times\delta\Biggl{(}\bm{E}_{t}\Bigl{(}\frac{\bm{x}_{t+\varDelta t}-\bm{x}_{t}}{\varDelta t}-\bm{f}_{t}(\bm{x}_{t})\Bigr{)}-\bm{\eta}_{t}\Biggr{)}$$

$$1$$

$$1={\int }_{-\infty }^{+\infty }{d}^m{x}_{t+\Delta t}p({\mathbf{x}}_{t+\Delta t},t+\Delta t)$$

$$p({\mathbf{x}}_{t+\Delta t},t+\Delta t)$$

$$\equiv {\int }_{-\infty }^{+\infty }{d}^m{x}_t\frac{\sqrt{g_t}}{(2\pi \Delta t{)}^{m/2}}{e}^{-T_{t+\Delta t}\Delta t}P({\mathbf{x}}_t,t),$$

$$\displaystyle T_{t+\varDelta t}=\frac{1}{2}\Bigl{(}\frac{\bm{x}_{t+\varDelta t}-\bm{x}_{t}}{\varDelta t}-\bm{f}_{t}(\bm{x}_{t})\Bigr{)}^{T}$$

$$\displaystyle\qquad\qquad\qquad\times\bm{g}_{t}\Bigl{(}\frac{\bm{x}_{t+\varDelta t}-\bm{x}_{t}}{\varDelta t}-\bm{f}_{t}(\bm{x}_{t})\Bigr{)}.$$

$$p({\mathbf{x}}_{t+\Delta t},t+\Delta t)=P({\mathbf{x}}_{t+\Delta t},t+\Delta t),$$

$$\exp (-U_{t+\Delta t}({\mathbf{x}}_{t+\Delta t})\Delta t).$$

$$P({\mathbf{x}}_{t+\Delta t},t+\Delta t)$$

$$\equiv e^{-U_{t+\Delta t}({\mathbf{x}}_{t+\Delta t})\Delta t}p({\mathbf{x}}_{t+\Delta t},t+\Delta t)$$

$$={\int }_{-\infty }^{+\infty }{d}^m{x}_t\frac{\sqrt{g_t}}{(2\pi \Delta t{)}^{m/2}}$$

$$\times e^{-(T_{t+\Delta t}+U_{t+\Delta t}({\mathbf{x}}_{t+\Delta t}))\Delta t}P({\mathbf{x}}_t,t).$$

$$1={\int }_{-\infty }^{+\infty }{d}^m{x}_{t+\Delta t}P({\mathbf{x}}_{t+\Delta t},t+\Delta t).$$

$$P({\mathbf{x}}_{t+\Delta t},t+\Delta t|{\mathbf{x}}_t,t)$$

$$=\frac{\sqrt{g_t}}{(2\pi \Delta t{)}^{m/2}}\times e^{-(T_{t+\Delta t}+U_{t+\Delta t}({\mathbf{x}}_{t+\Delta t}))\Delta t}$$

$${\partial }_{t}P({\mathbf{x}}_t,t)=$$

$$-\sum _{\mu =1}^m{\partial }_{\mu }(f_t^{\mu }({\mathbf{x}}_t)P({\mathbf{x}}_t,t))-U_t({\mathbf{x}}_t)P({\mathbf{x}}_t,t)$$

$$P({\mathbf{x}}_t,t)=e^{W({\mathbf{x}}_t)}\rho ({\mathbf{x}}_t,t)$$

$${\partial }_t\rho ({\mathbf{x}}_t,t)$$

$$=\sum _{\mu ,\nu =1}^m\frac{1}{2}{g}_t^{\mu \nu }{\partial }_{\mu }{\partial }_{\nu }\rho ({\mathbf{x}}_t,t)$$

$$-\sum _{\mu =1}^m{\partial }_{\mu }\{(f_t^{\mu }({\mathbf{x}}_t)-\sum _{\nu =1}^m{g}_t^{\mu \nu }{\partial }_{\nu }W({\mathbf{x}}_t))\rho ({\mathbf{x}}_t,t)\}$$

$$+\sum _{\mu =1}^m(f_t^{\mu }({\mathbf{x}}_t)-\sum _{\nu =1}^m{g}_t^{\mu \nu }{\partial }_{\nu }W({\mathbf{x}}_t)){\partial }_{\mu }W({\mathbf{x}}_t)\rho ({\mathbf{x}}_t,t)$$

$$\displaystyle\quad-\Bigl{[}U_{t}(\bm{x}_{t})+\frac{1}{2}\sum_{\mu,\nu=1}^{m}\{g_{t}^{\mu\nu}\partial_{\mu}\partial_{\nu}W(\bm{x}_{t})$$

$$\displaystyle\qquad\qquad\qquad+g_{t}^{\mu\nu}\partial_{\mu}W(\bm{x}_{t})\partial_{\nu}W(\bm{x}_{t})\}\Bigr{]}\rho(\bm{x}_{t},t).$$

$$f_t^{\mu }({\mathbf{x}}_t)-\sum _{\nu =1}^m{g}_t^{\mu \nu }{\partial }_{\nu }W({\mathbf{x}}_t)=0,$$

$${\partial }_t\rho ({\mathbf{x}}_t,t)$$

$$=\sum _{\mu ,\nu =1}^m\frac{1}{2}{g}_t^{\mu \nu }{\partial }_{\mu }{\partial }_{\nu }\rho ({\mathbf{x}}_t,t)-U_t({\mathbf{x}}_t)\rho ({\mathbf{x}}_t,t)$$

$$-\frac{1}{2}\sum _{\mu ,\nu =1}^m(g_t^{\mu \nu }{\partial }_{\mu }{\partial }_{\nu }W+g_t^{\mu \nu }{\partial }_{\mu }W{\partial }_{\nu }W)\rho ({\mathbf{x}}_t,t).$$

$${\int }_{-\infty }^{+\infty }{d}^m{x}_{t_N}{C}_{t_N}({\mathbf{x}}_{t_N})P({\mathbf{x}}_{t_N},t_N)$$

$$={\int }_{-\infty }^{+\infty }{d}^m{x}_t{C}_t({\mathbf{x}}_t)P({\mathbf{x}}_t,t)$$

$$C_t({\mathbf{x}}_t)$$

$$\equiv {\int }_{-\infty }^{+\infty }{d}^m{x}_{t+\Delta t}{C}_{t+\Delta t}({\mathbf{x}}_{t+\Delta t})P({\mathbf{x}}_{t+\Delta t},t+\Delta t|{\mathbf{x}}_t,t)$$

$$={\int }_{-\infty }^{+\infty }{d}^m{x}_{t_N}{C}_{t_N}({\mathbf{x}}_{t_N})$$

$$\times \prod _{i=t}^{t_N-\Delta t}{d}^m{x}_{i}P({\mathbf{x}}_{i+\Delta t},i+\Delta t|{\mathbf{x}}_i,i).$$

$$P({\mathbf{x}}_t,t)$$

$$\delta {\mathbf{x}}_t$$

Full text

1 Introduction

Nonlinear stochastic processes appear in a wide variety of research fields such as physics, biology, filtering, control, machine learning, finance. A deep understanding of such stochastic phenomena depends on the formulation of the stochastic processes and then various mathematical models has been considered.

In physics, euclidean path integral Kleinert:2004ev has been well used as a powerful tool for modeling the stochastic processes and applied not only to physics but also to various areas. For example, the applications that had been discussed so far are as follows. Ref. RiskenFrank:1996 has explained that the path integral is able to satisfy a master equation of Langevin equation. Ref. Peliti:1985 has proposed a birth-death process represented by a path integral for an interacting system of stochastic fields on a lattice. The author of Ref. Kappen:2004 ; Kappen:2005 has discovered an unexpected relation between the path integral and stochastic nonlinear control. The relation has indicated that the Hamilton-Jacobi-Bellman (HJB) equation for Ito process involving a control input RFStengel:1994 becomes a linear partial differential equation (PDE), assuming a relation between the weight matrix of the control input and the covariance matrix of the Ito process. The solution of the PDE is given by the path integral which gives the control input. In Theodorou:2010-01 ; Theodorou:2010-02 , a task on finding a policy for a controlling robot in reinforcement learning has used the path integral method based on Kappen:2004 ; Kappen:2005 . As an application of the path integral to financial engineering Linetsky:1998 ; Baaquie:2004 , a pricing of a path-dependent option has considered the path integrals as Feynman-Kac formula satisfying a backward PDE. The path integral has involved a discount factor which deforms the stochastic dynamics of the option as the potential term.

For those stochastic models, the path integral is a solution of PDE on the stochastic process, however, in general, solving PDEs is quite difficult and solvable cases are limited. Therefore, it is natural to require a simple method which can directly define the path integral on a given stochastic process, in order to expand the applicabilities. The stochastic quantization (SQ) ParisiSourlas:1979 ; ParisiSourlas:1981 ; Parisi:1982ud ; Parisi:1998 ; Damgaard:1987rr ; Namiki:1992wf ; Cooper:1994eh ; Namiki:1993 has been known as such a direct method in physics. SQ has been originally a method to investigate an equilibrium state of Langevin dynamics in quantum physics. In Namiki:1992wf ; Namiki:1993 , it has been pointed out that SQ gives a systematic and simple procedure defining the path integral as the probability distribution for a general Ito process, without solving the corresponding PDEs. Using this fact, SQ has been applied a pricing of option in financial engineering Dash01:1988 ; Dash02:1989 ; DashBook:2004 ; Nakayama:2009 . If Ito process of a stock is given, SQ can define the path integral representation of the corresponding probability distribution which does not involve a potential term, directly. Such probability distribution can give the Feynman-Kac formula for the option which is a functional on a discount factor and a stock.

It is natural from this financial applications of SQ to consider whether SQ can define a path integral formula of probability distribution satisfying the dual forward PDE with a potential term or not. The existence of such probability distribution gives the expectation that SQ can directly model the forward process affected from the external potential. If we identify the potential with a penalty or barrier function Fletcher:2000 ; Nocedal_Wright:2006 , we will be enable to discuss various models in which the potential function constrains the Ito process as the new applications of SQ on engineering. However, how to incorporate the potential term into SQ has been not clarified for the general Ito process. If the Ito process is restricted to Langevin equation which has the drift term consisting of a partial derivative of a scalar function, Lagrangian is able to involve a potential term by scaling of the path integral Parisi:1982ud ; Parisi:1998 . The Langevin equation is quite restrictive for engineering problems, because Ito process appearing in engineering does not necessarily has such the drift term. Therefore, we require to clarify how the potential term can be incorporated in SQ without such the assumption on the drift term.

Taking into account this issue, in this paper, we will propose the introduction of weighted distribution Rao:1965_01 ; Rao:1984_01 ; Patil_Rao:1987_01 as a novel method to incorporate the potential term into the discrete-time SQ for the general Ito process. We will show the weighted distribution naturally induces the potential term in the path integral and preserves the role of the path integral as a probability distribution without imposing any assumptions on the drift term, although the weighted distribution has not been noticed in physics so far. For applications of SQ to engineering problems, we evaluate the path integral by a second order approximation on the stochastic fluctuations and then derives the difference equations. Taking an appropriate potential function, the difference equations will be applied to the modeling of filtering and control and to the numerical simulations of the control models.

This paper is organized as follows. In section 2, we will introduce SQ using weighted distribution in discrete time. The existence of the dual backward formula is mentioned. In section 3, the path integral gives the difference equations as the evolution equations of the expected value and the covariance matrix on the stochastic process, applying a second order approximation on the stochastic fluctuations. In section 4, by taking the potential function as the penalty or barrier function, the difference equations give rise to Extended Kalman Filter (EKF) Simon:2006 ; Saerkkae:2013 and models on nonlinear stochastic control. We also discuss numerical simulations of the nonlinear stochastic control. The section 5 will be devoted to the conclusions.

2 SQ and Weighted distribution

We will introduce SQ to define a probability distribution on the Ito process by path integral. A weight function of the weighted distribution gives an additional potential term in the path integral and preserves the role of the path integral as the probability distribution, without imposing any assumptions on the drift term. To clarify the difference with the typical SQ discussed in physics, we will investigate PDE that the probability distribution involving the weighted distribution satisfies in continuous time. It is shown that the PDE recovers the Fokker-Planck equation, if the potential term vanishes. Although the Fokker-Planck equation is able to induce an additional potential term by a scaling of the solution and by assuming that the drift term of the Ito process is proportional to a partial derivative of a scalar function, we find that the SQ with the weighted distribution can introduce the additional potential term without the scaling and the assumption.

2.1 Path integral formula of SQ

We will consider a finite-time system in discrete time. $\bm{\eta}_{t}$ is a $m$-dimensional random variable which obeys the normal (Gaussian) distribution with $E[\bm{\eta}_{t}]=0$ and $E[\bm{\eta}_{t}\bm{\eta}_{t}^{T}]=\bm{1}_{m\times m}$. The probability distribution is given by

[TABLE]

Using $\bm{\eta}_{t}$, a discrete-time Ito process is defined by

[TABLE]

where $\bm{x}_{t}$ , $\bm{f}_{t}(\bm{x}_{t})$ are $m$-dimensional vectors. $\bm{E}^{-1}_{t}$ is an $m\times m$-matrix which has the inverse matrix $\bm{E}_{t}$ with $\bm{E}^{-1}_{t}\bm{E}_{t}=\bm{1}_{m\times m}$. The Ito process (2) indicates that $\bm{x}_{t+\varDelta t}$ is a random variable and $\bm{x}_{t}$ is already given at $t$.

SQ ParisiSourlas:1979 ; ParisiSourlas:1981 ; Parisi:1982ud ; Parisi:1998 ; Damgaard:1987rr ; Namiki:1992wf ; Cooper:1994eh ; Namiki:1993 has been known as a technique which defines the probability distribution of the Ito process (2) by path integral. To define the probability distribution for Eq. (2), the SQ introduces the following delta function:

[TABLE]

where $\bm{g}_{t}\equiv\bm{E}^{T}_{t}\bm{E}_{t}>0$ and $\sqrt{g_{t}}\equiv\det\bm{E}_{t}$. We assume that an initial probability density of $\bm{x}_{t}$ is given by a positive function $P(\bm{x}_{t},t)\geq 0$ satisfying the following equation;

[TABLE]

Inserting Eq. (1) and E. (3) to Eq. (4), we obtain the following equation after the integration on $\bm{\eta}_{t}$:

[TABLE]

where

[TABLE]

$T_{t+\varDelta t}$ is the kinetic term. Eq. (6) is the path integral representation of the probability distribution for the Ito process (2) and describes an evolution of the initial distribution $P(\bm{x}_{t},t)$. If we impose

[TABLE]

Eq. (6) is the path integral solution for the Fokker-Planck equation RiskenFrank:1996 . In order to constrain the stochastic dynamics of the Ito process (2) for the applications to various engineering problems, we have to introduce external forces without the condition (8) and any assumption for the drift term of Eq. (2). By facts on the path integral Kleinert:2004ev , it is expected that introducing of an additional potential term in the path integral induces such forces. In the following discussions, we introduce weighted distribution to naturally define the potential term in the path integral.

2.2 The introduction of weighted distribution

We will consider an extension of Eq. (6) to path integral with an additional potential function without assumption for the drift term of the Ito process (2). Inspired by the weighted distribution Rao:1965_01 ; Rao:1984_01 ; Patil_Rao:1987_01 which has been considered in the area of statistics, we would like to introduce a weight function,

[TABLE]

Using this weight function, we define $P(\bm{x}_{t+\varDelta t},t+\varDelta t)$ as a weighted distribution:

[TABLE]

We impose conservation of probability for the weighted distribution as follows:

[TABLE]

This condition implies that the transition probability distribution for Eq. (11) can be defined by

[TABLE]

where we have assumed that the normalization is involved in the potential function. Therefor, Eq. (11) is interpreted as the evolution of the initial probability density Eq. (4) through the transition probability distribution Eq. (13).

As shown in ParisiSourlas:1981 ; Parisi:1982ud ; Parisi:1998 , SQ without the weighted distribution can induce a potential term by introducing an additional drift term, after a scaling of the probability density which satisfies Fokker-Planck equation. To understand that our model (11) is different with the case in ParisiSourlas:1981 ; Parisi:1982ud ; Parisi:1998 , we will consider the continuous case of Eq. (11). Considering a fluctuation $\bm{\xi}_{t}=\bm{x}_{t}-\bm{x}_{t+\varDelta t}$ and Expanding Eq. (11) around $\bm{x}_{t+\varDelta t}$ by $\bm{x}_{t}=\bm{x}_{t+\varDelta t}+\bm{\xi}_{t}$, the continuous-time limit $\varDelta t\rightarrow 0$, after the integral on $\bm{\xi}_{t}$, gives the following equation:

[TABLE]

where $g_{t}^{\mu\nu}\equiv(\bm{g}_{t}^{-1})^{\mu\nu}$. If $U_{t}(\bm{x}_{t})=0$, Eq. (14) is exactly equivalent to the Fokker-Planck equation. Therefore, Eq. (11) is a simple extension of the typical SQ. For a case of $\bm{f}_{t}=\bm{0}$ in Eq. (14), we also find that $U_{t}(\bm{x}_{t})$ has the same role as potential function in physics, because Eq. (14) becomes the euclidean Schroedinger equation. In order to understand the difference, we will consider the following scaling of the probability distribution:

[TABLE]

where $W(\bm{x}_{t})$ is a scalar function. By the scaling, Eq. (14) becomes

[TABLE]

If we take additional condition on $\bm{f}_{t}$ as follows:

[TABLE]

Eq. (16) is reduced to the following equation ParisiSourlas:1981 ; Parisi:1982ud ; Parisi:1998 :

[TABLE]

We find that Eq. (18) is formally equivalent to the euclidean Schroedinger equation. For $U_{t}(\bm{x}_{t})=0$, this fact indicates that the Fokker-Plank equation can introduce the potential function under the assumptions (17). However, this potential function and the condition (17) are quite restrictive for applications on various engineering probrems, because a realistic stochastic process does not necessarily have such the scalar function $W(\bm{x}_{t})$ satisfying Eq. (17). In our model, we can consider the path integral involving the potential function $U_{t}(\bm{x}_{t})$ for a general stochastic process, without the assumption (17). Therefore, if an appropriate potential function $U_{t}(\bm{x}_{t})$ is taken for an engineering problem, it is expected that we can constrain a general stochastic process toward a desired path for solving a task on the problem. This point is the major difference from SQ in physics.

2.3 The existence of dual backward formula

Finally, we will mention whether a dual backward formula can be defined from the forward formula (11). If the dual function at $t_{N}\equiv t+N\varDelta t$ is given by a real function $C_{t_{N}}(\bm{x}_{t_{N}})$ which is not necessarily normalized, an inner product for (11) and $C_{t}(\bm{x}_{t})$ at $t_{N}$ can be represented by

[TABLE]

where $C_{t}(\bm{x}_{t})$ has been defined as follows:

[TABLE]

The function $C_{t}(\bm{x}_{t})$ clearly satisfies the backward Kolmogolov equation (20). If $C_{t_{N}}(\bm{x}_{t_{N}})$ is taken as a functional on options, time and an appropriate normalization constant, Eq. (21) is able to be interpreted as the Feynman-Kac formula involving a discount factor Linetsky:1998 ; Baaquie:2004 proportional to $-\sum_{i=t}^{t_{N}-\varDelta t}U_{i+\varDelta t}(\bm{x}_{i+\varDelta t})\varDelta t$ by (13).

We would like to also mention a difference with a backward path integral of a nonlinear control model discussed in Kappen:2004 ; Kappen:2005 . If we assume that $C_{t_{N}}(\bm{x}_{t_{N}})$ is an exponential function on a cost function like $\exp(-\phi(\bm{x}_{t_{N}})\varDelta t)$, $C_{t}(\bm{x}_{t})$ gives a similar path integral formula on a cost-to-go function in Kappen:2004 ; Kappen:2005 . However, in discrete time, Ref. Kappen:2004 ; Kappen:2005 uses the potential function at $t$ instead of the potential function at $t+\varDelta t$. Therefore, in the discrete-time level, the model proposed in this paper is directly different because a different discretization for the potential function in term of the time variable is applied in Kappen:2004 ; Kappen:2005 .

3 Nonlinear Difference equations

We have discussed SQ with the weight distribution in the previous section. It has been indicated without assumption like Eq. (17) that the weight function induces the potential function in the path integral, preserving the role as the probability distribution. By analogy with the path integral in physics, we have expected a constrained dynamics by a forces derived from the potential function. To confirm such expectations through explicit examples involving the numerical simulations, we would like to evaluate the path integral by a second order expansion on the stochastic fluctuations. We will derive the difference equations as evolution equations of covariance matrix and expectation value for the stochastic variables.

As the initial probability distribution (4), we will take gaussian distribution with the mean value $\bm{\hat{x}}_{t|t}$ and the covariance matrix $\hat{\varSigma}_{t|t}$:

[TABLE]

where the stochastic fluctuation has been defined by

[TABLE]

To perform the integration in Eq. (6), we will expand $\bm{x}_{t}+\bm{f}_{t}(\bm{x}_{t})\varDelta t$ around $\bm{\hat{x}}_{t|t}$ as follows:

[TABLE]

where we have defined

[TABLE]

After the integration on $\delta\bm{x}_{t}$ ($d\delta\bm{x}_{t}=d\bm{x}_{t}$), Eq. (6) is represented by

[TABLE]

where $\delta\bm{x}_{t+\varDelta t|t}=\bm{x}_{t+\varDelta t}-\bm{\hat{x}}_{t+\varDelta t|t}$ and the prediction processes have been defined by

[TABLE]

In order to evaluate the path integral (11), we are going to assume that the potential function $U_{t+\varDelta t}$ is defined as follows:

[TABLE]

where $V^{eff}_{t+\varDelta t}$ is the effective potential, $\delta_{c}V_{t+\varDelta t}(\bm{x}_{t+\varDelta t})$ is a counter term which is decided in later discussion and $N_{t+\varDelta t}$ is a normalization. In quantum mechanics, it is known that path integrals are able to take in the quantum fluctuations around the classical paths by performing the integral on the fluctuations Kleinert:2004ev . To apply this method for evaluating Eq. (11), we will consider $V_{t+\varDelta t}$ up to the second order terms on $\delta\bm{x}_{t+\varDelta t|t}$ as follows:

[TABLE]

where we have defined the following equations:

[TABLE]

We also take the normalizations $N_{t+\varDelta t}$ as follows:

[TABLE]

where $\mathcal{N}_{t+\varDelta t}$ is defined by

[TABLE]

For a stability of the numerical simulations, we will assume

[TABLE]

Using Eq. (37), Eq. (38), Eq. (39), Eq. (40) and Eq. (97), it is found that the weighted distribution is given by

[TABLE]

where $\bm{\hat{x}}_{t+\varDelta t|t+\varDelta t}$ and $\hat{\varSigma}_{t+\varDelta t|t+\varDelta t}$ are defined as follows:

[TABLE]

Eq. (27), Eq. (28), Eq. (42) and Eq. (43) are difference equations that we would like to derive.

Finally, we will comment on the assumptions Eq. (39) and Eq. (40). If we does not assume Eq. (39) and Eq. (40), the inverse matrix in Eq. (98) and Eq. (100) is replaced by

[TABLE]

We are able to ignore the second and the third term of the inverse matrix in the continuous limit $\varDelta t\rightarrow 0$. However, for the discrete-time, this inverse matrix has a possibility of the divergence or of the negative value because of a case where $\varSigma_{\nu_{t+\varDelta t}}^{-1}$ and ${\varSigma^{\prime}}_{\nu_{t+\varDelta t}}^{-1}$ are singular or negative, in general. This causes a difficulty of the stable numerical simulations. Therefore, we will impose Eq. (39) and Eq. (40) through this paper.

4 Applications to filtering and control

In the previous section, we have derived the difference equations on the evolution of the nonlinear stochastic systems involving the potential function. The difference equation Eq. (42) implies that the expectation value $\bm{\hat{x}}_{t+\varDelta t|t}$ is affected by forces from the potential at $t+\varDelta t$. To understand this, more explicitly, note that by Eq. (98) in Appendix B, Eq. (42) can be represented as follows:

[TABLE]

where $(\hat{\varSigma}_{t+\varDelta t|t}^{-1}+\bm{H}_{t+\varDelta t}^{T}\varSigma_{\nu_{t+\varDelta t}}^{-1}\bm{H}_{t+\varDelta t}\varDelta t)^{-1}$ is positive-definite. Eq. (44) indicates that the potential term effectively gives rise to forces like classical mechanics and constrains $\bm{\hat{x}}_{t+\varDelta t|t}$ towards a local minimum or a decreasing direction of the potential function, if we identify the potential term with a penalty or barrier function Fletcher:2000 ; Nocedal_Wright:2006 . This fact allows us to expect that SQ using the weighted distribution can be applied to engineering problems like filtering or control, taking the appropriate potential function.

In this section, we will confirm such the expectations by the analytical or numerical examples. EKF is analytically derived by a quadratic potential function as the potential function representing the filtering process. The potential function involves the observations as an external sources, however we explain that the stochastic interpretation of the observations can be recovered like typical EKF. The numerical simulations are considered as the examples of control, considering the potential function as the cost function. It will be shown that the stochastic dynamics can be constrained towards a target direction by minimizing the potential function. Through this section, we take

[TABLE]

4.1 Extended Kalman Filter

To construct EKF, we consider the following quadratic potential as a penalty function:

[TABLE]

In Eq. (47), $\bm{y}_{t+1}$ is an external source, $\bm{h}_{t+1}(\bm{x}_{t+1})$ is a $m$-dimensional vector and $\varSigma_{\nu_{t+1}}$ is an $m$-dimensional positive-definite matrix. Then, for the Ito process (2), the expected dynamics is represented by the prediction process (27)-(28) and the update process (42)-(43) which becomes as follows:

[TABLE]

It is found that (27)-(28) and (48)-(49) are the formula of EKF. However, Eq. (47) treats $\bm{y}_{t+1}$ as an external source and does not assume the probability distribution. $\varSigma_{\nu_{t+1}}>0$ is a constant matrix which is not directly related with $\bm{y}_{t+1}$ as the covariance matrix. This is quite different point from typical EKF Simon:2006 ; Saerkkae:2013 in which $\bm{y}_{t}$ is treated as an observation given by $\bm{y}_{t+1}=\bm{h}_{t+1}(\bm{x}_{t+1})+\bm{\xi}_{t+1}$ for $\bm{\xi}_{t+1}\sim\mathcal{N}(\bm{0},\varSigma_{\nu_{t+1}})$. In our case, EKF can be interpreted as a process in which the potential constrains the expected motion of the stochastic process towards a direction satisfying $\bm{y}_{t+1}\simeq\bm{h}_{t+1}(\bm{\hat{x}}_{t+1|t})$ in which the potential function becomes the minimum.

We can also recover the stochastic interpretation of $\bm{y}_{t+1}$ as a special case. For example, by taking an approximation like

[TABLE]

for Eq. (47), the weight function can be reduced to

[TABLE]

where $\alpha_{t+1}$ is an arbitrary function. If we take $\alpha_{t+1}=1$,

[TABLE]

and then we can interpret Eq. (52) and Eq. (53) as the probability densities. Then, $\bm{y}_{t+1}$ obeys the Gaussian distribution. This is the standard interpretation on $\bm{y}_{t+1}$ in EKF. In general, $\alpha_{t+1}$ is not necessarily equal to one and then Eq. (52) and Eq. (53) do not have an interpretation as a probability density. Our model is able to relax the assumption of the probability distribution for the observation.

4.2 Constrained dynamics as stochastic nonlinear control

For the applications of SQ to a nonlinear stochastic control, we will rewrite the evolution equation (42) to the following form:

[TABLE]

where $\bm{u}_{t}$ is identified with the control input derived from the potential and has been defined by

[TABLE]

$B_{t}$ and $R_{t}$ are $m\times l$ and $l\times l$ matrix. we have assumed $R_{t}$ and $B_{t}R_{t}^{-1}B_{t}^{T}$ are regular and positive. Note that the covariance matrix evolves according to Eq. (43) as the forward process. This is different from a typical feedback control RFStengel:1994 and then it is not trivial whether the potential function can control the dynamics represented by Eq. (56) or not. To check how the control is realized, we will perform the numerical simulations by specifying the Ito process and the potential term. For simplicity, we take

[TABLE]

In this part, we would like to consider a stochastic nonlinear control of the stochastic process (2) with the following definitions:

[TABLE]

where $\bm{F}_{t}$ is given by Eq. (28) and $a$ and $b$ of $(\bm{F}_{t})^{a}_{~{}b}$ indicate row dimensions and column dimensions, respectively. Fig. 1 indicates sample paths of the stochastic process (2). In below discussion, we will consider two cases, ${\varSigma^{\prime}}_{\nu_{t+1}}^{-1}=0$ and ${\varSigma^{\prime}}_{\nu_{t+1}}^{-1}\neq 0$.

4.2.1 ${\varSigma^{\prime}}_{\nu_{t+1}}^{-1}=0$

In this case, the counter term vanishes, $\delta_{c}V_{t+1}=0$ and the equation (36) implies that we are able to take $\bm{l}_{t+1}(\bm{x}_{t+1})$ as a linear function on $\bm{x}_{t+1}$.

As a first example, to constrain the stochastic process towards a target $\bm{d}_{t+1}$, we will consider the following effective potential as the cost function defined by a quadratic penalty function like (46):

[TABLE]

This effective potential is a positive function and has a local minimum at $\bm{d}_{t+1}$. By the previous discussions, the principle of least action implies $\bm{\hat{x}}_{t+1|t+1}\rightarrow\bm{d}_{t+1}$, if we appropriately tune the parameter $\bm{g}_{t}$ and $\varSigma_{\nu_{t+1}}$. To numerically check this expectation, we will take the initial values as follows:

[TABLE]

Fig. 2 and Fig. 3 show sample paths of $\bm{x_{t}}$ and $\bm{u_{t}}$. These sample paths are given by the following steps. The initial control $\bm{u}_{0}$ is decided by substituting the conditions (73)-(80) for Eq. (57). $\bm{x}_{1}$ is sampled by the probability distribution (41) for the conditions. By repeating similar steps, we can obtain the sample paths. Fig. 2 and Fig. 3 indicate that the stochastic dynamics is controlled to the local minimum of $V^{eff}_{t+1}$.

The second example is a constrained dynamics in $\bm{x}_{t+1}>\bm{0}$ and then the cost function is taken by a potential like barrier function:

[TABLE]

where $\bm{a}$ is a positive vector and

[TABLE]

This potential monotonically increases towards $+\infty$ for $\bm{x}_{t+1}\rightarrow\bm{0}$ and decreases to $-\infty$ for $\bm{x}_{t+1}\rightarrow+\infty$ . It is expected that the potential constrains the stochastic processes to $\bm{x}_{t+1}>\bm{0}$. Using Eq. (35), $\varSigma_{\nu_{t+1}}$ is given by

[TABLE]

$\bm{H}_{t+1}$ is given by $\bm{H}_{t+1}=-\bm{1}_{2\times 2}$. The numerical simulation gives Fig. 4 and Fig. 5 by the same condition and manner with the first example, applying Eq. (83) and

[TABLE]

Fig. 4 and Fig. 5 show that the stochastic dynamics is constrained to $\bm{x}_{t+1}>\bm{0}$, by the potential as barrier function.

4.2.2 ${\varSigma^{\prime}}_{\nu_{t+1}}^{-1}\neq 0$

We will consider a case in which the cost function is given by a double well potential defined by Eq. (69) with

[TABLE]

where

[TABLE]

The double well potential has the minima at $d^{m}_{t+1}$ or $-d^{m}_{t+1}$ by Eq. (85). The target $\bm{d}_{t+1}$ changes from $\begin{pmatrix}0&0\end{pmatrix}^{T}$ to $\begin{pmatrix}0.4&0.4\end{pmatrix}^{T}$ and then the minima of the double well potential dynamically deform, according to $\bm{d}_{t+1}$. Substituting Eq. (69) and Eq. (85) for Eq. (36), ${\varSigma^{\prime}}_{\nu_{t+1}}^{-1}$ is given by

[TABLE]

The numerical simulation gives Fig. 6 and Fig. 7 by the same condition and manner with the previous examples, replacing the corresponding conditions by

[TABLE]

We find that the potential controls the stochastic processes along the dynamically changing minima.

5 Conclusions

In this paper, we have proposed the discrete-time SQ using the weighted distribution as a novel and simple method to treat constrained dynamics of general Ito processes in engineering problems. The introduction of the weighted distribution in SQ has preserved the role of the path integral as the probability distribution and has naturally induced the potential term in the path integral, without imposing any assumption of the drift term. A second order approximation on the stochastic fluctuations to the path integral has given the difference equations representing the time-evolution of the expectation value and the covariance matrix for the stochastic systems. The difference equations have indicated that the potential term gives rise to forces constraining the stochastic dynamics. The difference equations has derived EKF by the potential function as the penalty function. Although the penalty function has involved observations as just external sources, we have also explained a case where such the stochastic interpretation of the observations is recovered as the typical definition of EKF. The potential function as a barrier or penalty function has been applied to model the nonlinear stochastic control. It has been numerically shown that the nonlinear stochastic processes can be constrained towards a minimum or a decreasing direction of the potential function.

For more realistic applications, how to accept constraints on control inputs induced from the potential term also remains to be considered. A hyperparameter optimization is also required since the stability of a given stochastic system depends on the hyperparameters. The computation on the covariance matrix has the high cost for a large-scale and high-dimensional system. Thus, how to reduce the computational cost is remaining problem to be solved in order to be more powerful tool for real world. It is also interesting to consider parameter estimation because EKF is able to estimate parameters of a machine learning model Singhal_Wu:1988 .

Acknowledgements.

The author would like to thank Koichi Hamada, Shinichi Takayanagi and Hirofumi Yamashita for valuable discussions and comments about early ideas on this paper at a seminar. The author is also grateful to Yu Higashiyama and Takashi Kawasaki for the opportunity and the support of the seminar at Fringe81 Co., Ltd..

Appendix A Matrix inversion lemma and product rules for determinant

$A$ and $D$ are invertible square matrices, and $B$ and $C$ matrices may or may not be square. The matrix inversion lemma and the product rules for determinant are given as follows Simon:2006 :

[TABLE]

Appendix B Identities

Applying the matrix inversion lemma in Appendix A, we can show the following equations:

[TABLE]

where we have defined the following equations:

[TABLE]

[TABLE]

Conflict of interest

The authors declare that they have no conflict of interest. This manuscript is an extended version of author’s presentation on stochastic process, when the author was working at Fringe81 Co., Ltd.. The manuscript is not related to the business.