A parameter estimation method using linear response statistics: Numerical scheme

Zhang, He; Li, Xiantao; Harlim, John

doi:10.1063/1.5081744

This paper presents a numerical method to implement the parameter estimation method using response statistics that was recently formulated by the authors. The proposed approach formulates the parameter estimation problem of Itô drift diffusions as a nonlinear least-squares problem. To avoid solving the model repeatedly when using an iterative scheme in solving the resulting least-squares problems, a polynomial surrogate model is employed on appropriate response statistics with smooth dependence on the parameters. The existence of minimizers of the approximate polynomial least-squares problems that converge to the solution of the true least square problem is established under appropriate regularity assumption of the essential statistics as functions of parameters. Numerical implementation of the proposed method is conducted on two prototypical examples that belong to classes of models with a wide range of applications, including the Langevin dynamics and the stochastically forced gradient flows. Several important practical issues, such as the selection of the appropriate response operator to ensure the identifiability of the parameters and the reduction of the parameter space, are discussed. From the numerical experiments, it is found that the proposed approach is superior compared to the conventional approach that uses equilibrium statistics to determine the parameters.

In many dynamical systems, some of the model parameters may not be identifiable from their equilibrium statistics. For problems where the parameters can be identified from their two-point statistics, this paper provides a method to numerically estimate these parameters. Mathematically, we formulate the inverse problem as a well-posed dynamic constrained nonlinear least-squares problem that fits a set of linear response statistics under small external perturbations. Since solving this problem directly is computationally not feasible, we consider a polynomial surrogate modeling to approximate this least-squares problem. Convergence of the solution of the approximate problem to the solution of the original nonlinear least-squares solution is established under appropriate regularity assumptions. Supporting numerical examples on two classes of models with a wide range of applications, the Langevin dynamics and the stochastically forced gradient flows are given.

I. INTRODUCTION

Parameter estimation is ubiquitous in modeling of nature. The goal in this inverse problem is to infer the parameters from observations with sufficient accuracy so that the resulting model becomes a useful predictive tool. Existing parameter estimation techniques often belong to one of the following two classes of approaches, the maximum likelihood estimation¹ and the Bayesian inference such as Markov chain Monte-Carlo (MCMC),² depending on the availability of prior information about the parameters. Since the model parameters are usually not directly measured, the success of any inference method depends crucially on the identifiability of the parameters from the given observations. When the dependence of the observations on the parameters is implicit, that is, through the model, sensitivity analysis (see, e.g., the review article³) is often a useful practical tool to determine the parameter identifiability.

In this paper, we consider the parameter estimation of ergodic stochastic differential equations driven by Brownian noise. When the corresponding invariant measure of the dynamical systems has an explicit dependence on the parameters, then these parameters can usually be inferred from appropriate equilibrium statistical moments. A popular approach that exploits this idea is the reverse Monte Carlo method⁴ developed in chemistry. In particular, it formulates an appropriate nonlinear least-squares system of integral equations by matching the equilibrium averages of some pre-selected observables. Subsequently, Newton’s iterations are used to estimate the parameters. At each iterative step, samples at the current parameter estimate are generated for Monte-Carlo estimation to approximate the least-squares integral equations and the corresponding Jacobian matrix. In practice, this method can be rather slow due to the repeated sampling procedure. A severe limitation of this approach is that it is restrictive to inference of parameters of the equilibrium density.

This limitation can be overcome by fitting appropriate two-point statistics. As shown in our previous work,⁵ we formulated the parameter estimation problem based on the linear response statistics subjected to an external forcing, which drives the system out of equilibrium. The fluctuation-dissipation theory (FDT), a hallmark in non-equilibrium statistical physics, suggested that the changes of the average of an observable under small perturbations can be approximated by appropriate two-point statistics, called the FDT response operators.⁶ The key point is that these FDT response operators can be estimated using the available samples of the equilibrium unperturbed dynamics so long as we know the exact form of the invariant measure of the dynamics, which we will assume to be the case for a large class of problems, such as the Langevin dynamics and stochastic gradient flows. We should point out that the proposed approach relies on the validity of FDT response statistics, which has been studied rigorously for a large class of stochastic system that includes the stochastic differential equation (SDE) setup in this paper.⁷ For deterministic dynamics, one can use, for example, the statistical technique introduced in Ref. 8 to verify the validity of the FDT linear response. Following the ideas from Refs. 9–11, we developed a parameter inference method using these response operators in Ref. 5. While the method in Refs. 9–11 involves minimizing an information-theoretic functional that depends on both the mean and variance response operators, our approach fits a finite number of quantities, to be referred to as the essential statistics, which would allow us to approximate the FDT response operators, of appropriate observables (beyond just the mean and variance).

In our previous work,⁵ we showed the well-posedness of the formulation on three examples, in the sense that under infinite sampling, appropriate choices of essential statistics will ensure the identifiability of the model parameters. On one of these examples, the simplified model for turbulence,¹² which is a nonlinear system of SDEs with a Gaussian invariant measure, we were able to explicitly determine the choice of observables and external forcings (which in turn determine the FDT response operators) that allow us to identify all of the model parameters uniquely. In this case, one can directly estimate the parameters by solving the system of equations involving these essential statistics. While this result suggests that the choice of observables and external forcings is problem dependent, the FDT theory provides a guideline for choosing the appropriate two-point statistics for a well-posed parameter estimation. On another example, the Langevin dynamics (which will be discussed in Sec. IV), while the parameters can be identified by the essential statistics as shown in Ref. 5, the formulation suggests that one must include essential statistics that involve higher-order derivatives of the FDT response operators. This can be problematic because these higher-order statistics are rarely available in practice unless the available data are solutions of high-order SDE solvers. Given this practical issue, we propose to fit only the essential statistics that include the zeroth- and first-order derivatives of the FDT response operator (if possible, use only the zeroth-order derivative information). The price paid by restricting to only the lower-order essential statistical information is that the dependence of the essential statistics on some of the model parameters becomes implicit through the solutions of the dynamical model and the identifiability of these parameters becomes questionable.

In this paper, we devise a concrete numerical method to estimate the parameters by solving dynamic-constrained least-squares problems of integral equations involving these low-order essential statistics. As one would imagine, the implementation of this approach will face several challenges. First, just like in the reverse Monte-Carlo method, naive implementation of an iterative method for solving this dynamically constrained least-squares problem requires solving the model repeatedly. To avoid solving the model repeatedly in the minimization steps, we employ a polynomial surrogate model approach^13,14 on the least-squares cost function that involves the essential statistics. This approach is motivated by the fact that the cost function (or effectively the essential statistics) is subjected to sampling error and that the essential statistics have smooth dependence on the parameters, assuming that the Fokker-Planck operator of the underlying dynamics has smooth dependence on the parameters. Under appropriate regularity assumption on the essential statistics, we will provide a theoretical guarantee for the existence of minimizers of the approximate polynomial least-squares problem that converge to the solution of the true least-squares problem. The second related issue, as mentioned above, is that certain (lower-order) essential statistics might not be sensitive enough to some of the parameters, which will lead to inaccurate estimations. To ensure the practical identifiability of the parameters, we employ an empirical a priori sensitivity analysis based on the training data used in constructing the polynomial surrogate models and a posteriori local sensitivity analysis to ensure the validity of the a prior sensitivity analysis.

While the proposed technique can be used to infer all of the model parameters, in practice, one should use the polynomial approximation to infer parameters that cannot be estimated directly to avoid the curse of dimensionality. This issue is a consequence of using the polynomial expansion in approximating the least-squares problem, where the required training data set for constructing the surrogate model increases exponentially as a function of the parameter dimension. In the numerical examples below, we will apply the approximate least-squares problems to estimate parameters that cannot be directly estimated from matching equilibrium and/or two-point statistics that involve the first-order derivative of the FDT response operator.

The rest of the paper is organized as follows. In Sec. II, we briefly review the concept of essential statistics and parameter inference method using the linear response statistics developed in Ref. 5. In Sec. III, we present the proposed numerical algorithm based on the polynomial based surrogate model and discuss its convergence (with detailed proof in Appendixes A– C). In Secs. IV and V, we show applications on two nonlinear examples, a Langevin model and a stochastic gradient system with a triple-well potential, respectively. In Sec. VI, we conclude the paper with a summary and discussion.

II. A REVIEW OF ESSENTIAL STATISTICS

We begin by reviewing the parameter estimation formulation introduced in our previous paper.⁵ Consider an $n$ -dimensional Itô diffusion

d X = b (X; θ) d t + σ (X; θ) d W_{t},

(1)

where the vector field $b (X; θ)$ denotes the drift and $σ (X; θ)$ is the diffusion tensor. Both coefficients are assumed to be smooth functions in $X$ and $θ$ ⁠; $W_{t}$ represents the standard Wiener process, and the variable $θ \in D$ contains model parameters, where $D \subset R^{N}$ is the parameter domain. In this paper, we only consider the case where $D$ is bounded and for simplicity, we assume $D = {[- 1, 1]}^{N}$ ⁠. Throughout this paper, we denote the underlying true parameter value as $θ^{†}$ and the corresponding estimate as $\hat{θ}$ ⁠. We assume that Eq. (1) is ergodic (see Ref. 15 for precise conditions) with equilibrium density $p_{e q} (x; θ)$ for all $θ \in D$ ⁠, and $p_{e q} (x; θ)$ smoothly depends on both $x$ and $θ$ ⁠. When $θ = θ^{†}$ ⁠, we assume that we have the access to the explicit formula of $p_{e q} (x; θ^{†})$ as a function of $x$ only, which is shortened as $p_{e q}^{†} (x)$ ⁠. As we shall see in many applications, $p_{e q} (x; θ)$ only depends on a subset of the parameters $θ$ ⁠.

The main idea of the parameter estimation formulation introduced in Ref. 5 is to infer the parameters from the linear response operator associated with the fluctuation dissipation theory (FDT). Recall that, in the FDT setting, we are interested in an order- $δ$ (⁠ $0 < δ ≪ 1$ ⁠) external perturbation of the form $f (x, t) = c (x) δ f (t)$ ⁠, such that the solutions of the perturbed dynamics

d X^{δ} = [b (X^{δ}, θ) + c (X^{δ}) δ f (t)] d t + σ (X^{δ}, θ) d W_{t},

which would otherwise remain at the equilibrium state of the unperturbed dynamics (1) and can be characterized by a perturbed density $p^{δ} (x, t; θ)$ ⁠. The density function is governed by the corresponding Fokker-Planck equation under initial condition $p^{δ} (x, 0; θ) = p_{e q} (x; θ)$ ⁠.

By a standard perturbation technique, e.g., Ref. 1, the difference between the perturbed and unperturbed statistics of any integrable function $A (x)$ can be estimated by a convolution integral, that is,

\begin{aligned} Δ E [A] (t) & := E_{p^{δ}} [A (X)] (t) - E_{p_{e q}} [A (X)] \\ = \int_{0}^{t} k_{A} (t - s) δ f (s) d s + O (δ^{2}) . \end{aligned}

(2)

In (2), the term $k_{A} (t)$ is known as the linear response operator. The FDT formulates the linear response operator as the following two-point statistics:

\begin{aligned} k_{A} (t; θ) & := E_{p_{e q}} {A [X (t)] \otimes B [X (0); θ]} with \\ B_{i} (X; θ) & := - \frac{\partial_{X_{i}} [c_{i} (X) p_{e q} (X; θ)]}{p_{e q} (X; θ)}, \end{aligned}

(3)

where $B_{i}$ and $c_{i}$ denote the $i th$ components of $B$ and $c$ ⁠, respectively. We should point out that the validity of the FDT has been studied rigorously under mild conditions.⁷ A more explicit form is given by

k_{A} (t; θ) = \int_{R^{n}} \int_{R^{n}} A (x) \otimes B (y) ρ (x, t | y, 0) p_{e q} (y; θ) d x d y,

(4)

where $ρ$ is the solution of the Fokker-Planck equation

\frac{\partial}{\partial t} ρ = L^{*} ρ, ρ (x, 0 | y, 0) = δ (x - y) .

(5)

Here, $L$ denotes the generator of the unperturbed dynamics (1) and $L^{*}$ is its adjoint operator in the $L^{2} (R^{n})$ ⁠.

Remark 2.1

The result in Ref. 16 states that if the coefficients of linear parabolic PDEs are $C^{k}$ as functions of $θ,$ then the weak solutions are also $k$ -time differentiable with respect to $θ$ ⁠. In our case $,$ if $b, σ \in C^{k} ({[- 1, 1]}^{N}),$ then the linear response operator $k_{A} (t; θ)$ (4) is also $C^{k}$ with respect to the parameters $θ$ under the mild assumption that $p_{e q} (x; θ)$ is smooth with respect to $θ$ ⁠. In Sec. III B, we will conduct a convergence analysis on the proposed parameter estimation method to determine the necessary value of $k,$ the smoothness of the linear response operators as functions of $θ$ ⁠.

Notice that $B^{†} (X) := B (X; θ^{†})$ can be determined analytically based on our assumption of knowing the explicit formula of $p_{e q}^{†}$ ⁠. Given $t$ ⁠, the value of $k_{A} (t; θ^{†})$ can be computed using a Monte-Carlo sum based on the time series of the unperturbed system, $X$ ⁠, at $p_{e q}^{†} (x)$ ⁠. Therefore, the linear response operator can be estimated without knowing the underlying unperturbed system in (1) so long as the time series of $X$ at $p_{e q}^{†} (x)$ and the explicit formula for this equilibrium density are known.

The FDT, along with the expression of the kernel function (3), has been implemented to predict the change of the expectation of the observable $A$ ⁠.^17–22 The advantage is that the response operator is defined with respect to the equilibrium distribution, which does not involve the prior knowledge of the perturbed density $p^{δ}$ ⁠. In statistical physics, the FDT is a crucial formulation to derive transport coefficients, e.g., electron conductivity, viscosity, diffusion coefficients, etc.⁶ In the current paper, however, we propose to use the response statistics to infer the underlying true parameter value $θ^{†}$ ⁠.

Since this linear response operator is in principle infinite-dimensional, it is necessary for practical purposes to introduce a finite-dimensional approximation. The following parametric form of the response operator has been motivated by the rational approximation of the Laplace transform.²³ To explain the idea, consider the Laplace transformation of $k_{A} (t; θ)$ denoted by $K (s; θ)$ ⁠, which can be approximated by a rational function in the following form:

\begin{aligned} K (s; θ) \approx {(I - s^{- 1} β_{1} - \dots - s^{- m} β_{m})}^{- 1} (s^{- 1} α_{1} + \dots + s^{- m} α_{m}), \\ α_{i}, β_{i} \in R^{n \times n} . \end{aligned}

In the time domain, $k_{A} (t; θ)$ can be written explicitly as

\begin{aligned} k_{A} (t; θ) & \approx g_{m} (t; θ) := (\begin{matrix} I & 0 & \dots & 0 \end{matrix}) e^{t G} (\begin{matrix} α_{1} \\ α_{2} \\ ⋮ \\ α_{m} \end{matrix}), with \\ G & = (\begin{matrix} β_{1} & I \\ β_{2} & 0 & ⋱ \\ ⋮ & ⋱ & I \\ β_{m} & 0 \end{matrix}), \end{aligned}

(6)

where $I$ denotes the $n$ -by- $n$ identity matrix. Here, $m$ stands for the order of the rational approximation. In Ref. 23, the coefficients in the rational approximations were determined based on certain interpolation conditions, and such a rational approximation has been proven to be an excellent approximation for the time correlation function. For the current problem, $g_{m}$ will be determined using finitely many essential statistics defined as follows.

Definition 2.2

The values ${k_{A}^{(j)} (t_{i}; θ^{†})}_{i = 1, 2, \dots, K}$ are called essential statistics if they are sufficient to approximate the response operator $k_{A} (t; θ^{†})$ up to order- $m$ using (6), for $j = 0, \dots, 2 m - 1$ and $t_{1} < t_{2} \dots < t_{K},$ where $j$ indicates the order of derivatives.

Remark 2.3

In Ref. 5, we have considered the derivatives of $k_{A}$ at $t = 0$ for a Langevin dynamics model $,$ and explicit formulas involving the model parameters have been found. Theoretically $,$ all the parameters can be identified from ${k_{A}^{(j)} (0^{+}; θ^{†})},$ which to some extent shows the consistency. However $,$ in practice $,$ we use only the derivative of order no greater than one $,$ since estimating higher order derivative of $k_{A} (t, θ^{†})$ requires time series $X (t_{i})$ with sufficient accuracy $,$ which is not necessarily available. By using the essential statistics correspond to lower-order derivatives at $t_{i}$ other than zero $,$ we give up the explicit expressions of the higher-order statistics with respect to the parameters. This prompted us to solve the problem in the least-squares sense. Furthermore $,$ it requires other more practical methods $,$ such as a sensitivity analysis test $,$ to determine the identifiability of the parameters from these statistics.

Figure 1 shows the performance of such a discrete representation, where the linear response operator arises from a Langevin model subjected to a constant external forcing (25), which will be discussed in Sec. IV. With such discrete representation, we reduce our problem of inferring from an infinite-dimensional $k_{A} (t; θ^{†})$ to a finite number of essential statistics. To derive a system of equations involving both $θ$ and those essential statistics, we introduce

{\hat{k}}_{A} (t; θ) := E_{p_{e q} (θ)} {A [X (t)] \otimes B^{†} [X (0)]} .

(7)

We also define for $j = 0, 1$ ⁠,

M_{j} (t_{i}) := k_{A}^{(j)} (t_{i}; θ^{†}) = \frac{d^{j}}{d t^{j}} E_{p_{e q}^{†}} {A [X (t_{i})] \otimes B^{†} [X (0)]} .

We should stress that (7) is not the FDT response, since the $B^{†} (x)$ in (7) is defined with respect to ${p_{e q}}^{†}$ ⁠. The key idea in Ref. 5 is to estimate the true parameter values $θ^{†}$ by solving

M_{j} (t_{i}) = {\hat{k}}_{A}^{(j)} (t_{i}; θ), j \in {0, 1}, i = 1, 2, \dots, K .

(8)

Here, the term on the left-hand-side is estimated from the available sample at $p_{e q}^{†}$ ⁠, and the term on the right-hand-side will be computed from the time series generated by solving (1) for a given $θ \in D$ ⁠.

FIG. 1.

View large Download slide

The response function (solid lines) and its order-4 approximation (dashed lines).

Remark 2.4

Since $p_{e q}^{†}$ usually depends on a subset of the true parameter value $θ^{†},$ the assumption of knowing the explicit form of $p_{e q}^{†}$ as a function of $x$ only is too strong. In general $,$ we cannot even compute the left-hand-side term in (8) since the term $B^{†},$ which may depend on the unknown true value $θ^{†},$ is not available to us. We will discuss how to address such issue on two examples in Secs. IV and V, where each of them belongs to a class of models with a wide range of applications.

To ensure the solvability, we assume that the total number of equations in (8) is always greater than the dimension of the unknown parameters $θ$ ⁠. For this over-determined scenario, we propose to solve (8) in the least-squares sense, addressing the practical issue raised in Ref. 5 when the explicit solutions are not available. The main difficulty in solving the nonlinear least-squares problem comes from evaluating ${\hat{k}}_{A}^{(j)} (t_{i}; θ)$ ⁠, which often requires solving the model (1) repeatedly. In Sec. III, we will propose an efficient numerical method to address this issue.

III. AN EFFICIENT ALGORITHM FOR PARAMETER ESTIMATION

In Sec. II, we formulated a parameter estimation method as a problem of solving a nonlinear system (8) subject to a dynamical constraint in (1). In a compact form, we denote the system (8) as

f_{i} (θ) := M_{0} (t_{i}) - {\hat{k}}_{A} (t_{i}; θ) = 0, i = 1, \dots, K,

(9)

where $θ \in {[- 1, 1]}^{N}$ ⁠, $K > N$ ⁠. Notice that we have neglected the derivative constraints in (8). While similar algorithm and convergence result as shown below also hold if the first-order derivative constraints are included in the least-squares procedure, we will only numerically solve the least-squares problem for the zeroth-order derivative constraints as in (9). As we shall see in Secs. IV and V, we will use the derivative constraints to directly estimate some of the parameters.

Our goal is to solve (9) in the sense of least-squares, that is,

min_{θ \in {[- 1, 1]}^{N}} \sum_{i = 1}^{K} f_{i}^{2} (θ) .

(10)

However, as we mentioned in Sec. II, we do not necessarily have the explicit expressions for $f_{i} (θ)$ and evaluating $f_{i}$ for certain values of $θ$ requires solving the true model (1) to approximate the integral in $f_{i} (θ)$ by a Monte-Carlo sum. For example, if we apply the Gauss-Newton method to solve the nonlinear least-squares problem (10), the overall computation will become very expensive since we have to solve the model (1) repeatedly. As a remedy, we apply the idea of polynomial based surrogate model, which was introduced in Refs. 13 and 14 for solving Bayesian inference problems, to alleviate this part of computation. The main difference of our approach from Refs. 13 and 14 is that they used the orthogonal polynomials to expand the dynamical model, $X (t, θ)$ ⁠, whereas here, we will use the orthogonal polynomials to expand $f_{i} (θ)$ ⁠. While the convergence of the approximate posterior density estimate of the Bayesian inference problem is established in Ref. 14, we will end this section by analyzing the convergence of the approximate least-squares problem.

A. Polynomial based surrogate model

The key idea of our method is to approximate $f_{i} (θ)$ by a polynomial function $f_{i}^{M} (θ)$ ⁠. This is motivated by the fact that the sampling error cannot be avoided in computing the value of essential statistics, that is, there is uncertainty on the left-hand-side of (8). Thus, it is more reasonable to think of the nonlinear least-squares solution of (10) as a random variable with certain distribution $μ$ over ${[- 1, 1]}^{N}$ ⁠. Furthermore, since ${f_{i} (θ)}$ are $C^{k}$ with respect to $θ$ (Remark 2.1) whenever $b$ and $σ$ in (1) are $C^{k}$ with respect to $θ$ ⁠, polynomial chaos expansion provides a natural choice of $f_{i}^{M} (θ)$ based on the orthogonal polynomials with respect to $μ$ ⁠.

For instance, if we consider $μ \sim U [- 1, 1]$ ⁠, i.e., the uniform distribution, then the corresponding orthogonal polynomials are the Legendre polynomials ${p_{n}}$ ⁠. They are given by $p_{0} = 1$ ⁠, $p_{1} = x$ ⁠, and the recursive formula

(n + 1) p_{n + 1} (x) = (2 n + 1) x p_{n} (x) - n p_{n - 1} (x),

which form a basis for $L^{2} ([- 1, 1], μ)$ ⁠. The orthogonality of ${p_{n}}$ over $L^{2} ([- 1, 1], μ)$ is described by

\int_{- 1}^{1} p_{m} (x) p_{n} (x) d μ (x) = \frac{2}{2 n + 1} δ_{m n} .

One can normalize $p_{n}$ by introducing $P_{n} := p_{n} \sqrt{n + \frac{1}{2}}$ ⁠.

For multi-dimensional cases, one can define $P_{\vec{k}} (θ) := \prod_{j = 1}^{N} P_{k_{j}} (θ_{j})$ (⁠ $N$ represents the dimension of $θ$ ⁠), where $\vec{k} \geq \vec{0}$ is a multi-index with $\vec{k} = (k_{1}, \dots, k_{N}) \in {0, 1, 2, \dots}^{N}$ and $θ = (θ_{1}, \dots, θ_{N}) \in {[- 1, 1]}^{N}$ ⁠. We can consider an order- $M$ approximation of $f_{i} (θ)$

f_{i} (θ_{1}, \dots, θ_{N}) \approx \sum_{∥ \vec{k} ∥_{\infty} \leq M} α_{\vec{k}}^{(i)} P_{\vec{k}} (θ) =: f_{i}^{M} (θ) .

(11)

In practice, there are two approaches to determine the coefficient $α_{\vec{k}}^{(i)}$ ⁠. The first approach is the Galerkin method, which requires that the errors in the approximation in (11) to be orthogonal to the finite-dimensional subspace of $L^{2} ([- 1, 1], μ)$ generated by ${P_{\vec{k}} (θ)}_{∥ \vec{k} ∥_{\infty} \leq M}$ ⁠. In our convergence analysis, we will deduce the error estimates based on the coefficients obtained through this Galerkin (or least-squares) fitting. Alternatively, one can determine the coefficient $α_{\vec{k}}^{(i)}$ by the collocation method, that is, matching the values of $f_{i}$ at certain points of $θ$ in $D$ ⁠, denoted by $Θ$ ⁠, which leads to the following linear equations:

f_{i} (θ) = \sum_{∥ \vec{k} ∥_{\infty} \leq M} α_{\vec{k}}^{(i)} P_{\vec{k}} (θ), θ \in Θ .

(12)

Since (12) is equivalent to a polynomial interpolation, a common choice of $Θ$ is the product space of order- $M_{C}$ Chebyshev nodes, where $M_{C} \geq M + 1$ [when $M_{C} > M + 1$ ⁠, (12) turns into a linear least-squares problem]. In the numerical examples in this paper, we will use the collocation method to determine the coefficients.

Replacing $f_{i} (θ)$ in (10) by $f_{i}^{M} (θ)$ ⁠, we obtain a new least-squares problem for the order- $M$ polynomial based surrogate model

min_{θ \in {[- 1, 1]}^{N}} \sum_{i = 1}^{K} {(f_{i}^{M})}^{2} (θ) .

(13)

For easier notations, let us use $f^{M} (θ)$ and $J^{M} (θ)$ to denote the vector-valued function, ${[f_{1}^{M} (θ), \dots, f_{N}^{M} (θ)]}^{⊤}$ ⁠, and its Jacobian matrix, respectively.

In summary, we have Algorithm 1 based on the Gauss-Newton method.

Algorithm 1.

Parameter estimation: Polynomial based surrogate model

Let

Θ

be a set of collocation nodes in

{[- 1, 1]}^{N}

⁠.

1. For each

θ \in Θ

⁠, solve .

2. Estimate

f_{i} (θ)

using a Monte-Carlo sum over

i = 1, \dots, K

⁠.

3. Compute the coefficients

α_{\vec{k}}^{(i)}

by solving the linear system and obtain the approximations

f_{i}^{M} (θ)

for

i = 1, \dots, K .

4. For an initial guess

θ^{0} \in {[- 1, 1]}^{N}

and a threshold

δ > 0

⁠. Compute

θ^{k} = θ^{k - 1} - {{[J^{M} (θ^{k - 1})]}^{⊤} J^{M} (θ^{k - 1})}^{- 1} {[J^{M} (θ^{k - 1})]}^{⊤} f^{M} (θ^{k - 1}) .

(14)

while The step length

∥ θ^{k} - θ^{k - 1} ∥ \geq δ

do

Repeat .

end while

return

θ^{k}

Computationally, the most expensive step is in solving (1) on the collocation nodes $Θ$ ⁠, which can be done in parallel.

B. The convergence of the approximate solutions

In this section, we will provide the convergence analysis of the solutions of the approximate least-squares problem in (13) to the solution of the true least-squares problem in (10) as $M \to \infty$ ⁠, under appropriate conditions. We denote the solution of (13) by $θ_{M}^{*}$ ⁠, that is,

θ_{M}^{*} := \arg min_{θ \in {[- 1, 1]}^{N}} \sum_{i = 1}^{K} {(f_{i}^{M})}^{2} (θ) .

Since each $f_{i}^{M}$ is a multivariate polynomial, the minimum (13) can be attained in ${[- 1, 1]}^{N}$ ⁠; but such $θ_{M}^{*}$ may not be unique. For the case of $K = N$ and sufficiently smooth $f_{i}$ ⁠, assuming that the true parameter value $θ^{†} \in {(- 1, 1)}^{N}$ is the unique solution of the true least-squares problem in (10), we can show that there exists a sequence of minimizers of the approximated problem (13) ${θ_{M}^{*}}$ such that $θ_{M}^{*} \to θ^{†}$ as $M \to + \infty$ ⁠. To be more precise,

Before we apply numerical method to solve the problem, a fundamental question will be: assuming $θ^{†} \in {(- 1, 1)}^{N}$ ⁠, the true parameter value, is the unique solution of (10) can we find a sequence of minimizers of the approximated problem (13) ${θ_{M}^{*}}$ such that $θ_{M}^{*} \to θ^{†}$ as $M \to + \infty$ ? We are able to prove the following theorem, which partially answers this question.

Theorem 3.1

Consider the

N

-dimensional nonlinear least-squares problem and its surrogate model . In particular

,

we let

K = N

and assume

$f_{i} \in C^{ℓ} ({[- 1, 1]}^{N}),$ $ℓ = ⌈ \frac{3}{2} N ⌉ + 2$ ⁠, $i = 1, 2, \dots, N,$
Let $θ^{†}$ be the solution of (10) such that $f_{i} (θ^{†}) = 0,$ $i = 1, 2, \dots, N,$
the Jacobian matrix of $f := {(f_{1}, f_{2}, \dots, f_{N})}^{⊤}$ at $θ^{†},$ denoted by $\nabla \vec{f} (θ^{†}),$ is invertible.

Then, there exists a sequence of minimizers

{θ_{M}^{*}}

such that

θ_{M}^{*} = \underset{θ \in {[- 1, 1]}^{N}}{argmin} \sum_{i = 1}^{N} {(f_{i}^{M})}^{2} (θ)

and

θ_{M}^{*} \to θ^{†}

as

M \to + \infty

⁠. Moreover

,

the residual error

f^{M} (θ_{M}^{*})

converges to

0

as

M \to + \infty

⁠.

Remark 3.2

We set $N = K$ since the existence of the minimizers rely on the technical assumption in Lemma 3.6. The second and third assumptions provide the well-posedness of the original nonlinear least-squares problem (10). Since we discuss the parameter estimation problem under the perfect model setting $,$ the second assumption naturally holds.

The proof of Theorem 3.1 can be found in Appendix A. Here, we will present the proof for the one-dimensional case (⁠ $K = N = 1$ ⁠) to illustrate the main ideas. Notice that in the one-dimensional case, the formula in (11) is reduced to

f^{M} (θ) := \sum_{k = 0}^{M} α_{k} P_{k} (θ),

(15)

where ${P_{n}}$ are normalized Legendre polynomials. We start with two Lemmas, which pertain to the pointwise convergence of $f^{M}$ ⁠, and the proof can be found in Ref. 24.

Lemma 3.3

Let

f \in C^{ℓ} ([- 1, 1]),

then the corresponding least-squares approximation given by satisfies

lim_{k \to + \infty} α_{k} k^{ℓ} = 0.

Lemma 3.3 connects the smoothness of $f$ and the decay rate of the generalized Fourier coefficients $α_{k}$ ⁠. The following Lemma is a direct result of Lemma 3.3.

Lemma 3.4

Let

f \in C^{2} ([- 1, 1])

and

f^{M}

be the corresponding least-squares approximation given by , then for any

ϵ > 0,

there exists

N > 0

such that

\forall M > N

we have

∥ f - f^{M} ∥_{\infty} \leq \frac{ϵ}{\sqrt{M}} .

Thus

,

we have

f^{M} \to f

pointwise on

[- 1, 1]

⁠.

Here, $∥ \cdot ∥_{\infty}$ denotes the $L^{\infty}$ -norm on $C^{2} [- 1, 1]$ ⁠. In the proof of Proposition 3.7, we also need $∥ {(f^{M})}^{'} - f^{'} ∥_{\infty} \to 0$ as $M \to + \infty$ ⁠. To obtain such convergence, one needs higher regularity on $f$ ⁠. In particular, we have

Lemma 3.5

Let $f \in C^{4} ([- 1, 1])$ and $f^{M}$ be the corresponding least-squares approximation given by (15). We have $∥ f^{'} - {(f^{M})}^{'} ∥_{\infty} \to 0$ as $M \to + \infty$ ⁠.

Proof.

Following the same idea in the proof of Lemma 3.4,²⁴ we are going to show that

{{(f^{M})}^{'}}

is a Cauchy sequence under

L^{\infty}

-norm. Notice for

m > n

⁠,

∥ {(f^{m} - f^{n})}^{'} ∥_{\infty} = {∥ \sum_{k = n + 1}^{m} α_{k} P_{k}^{'} ∥}_{\infty} \leq \sum_{k = n + 1}^{m} | α_{k} | ∥ P_{k}^{'} ∥_{\infty},

which suggests that we need to find a

L^{\infty}

bound for

P_{k}^{'}

⁠. Using the recursion relation

p_{n + 1}^{'} - p_{n - 1} = (2 n + 1) p_{n}

and

∥ p_{n} ∥_{\infty} = 1

⁠, we have

∥ P_{n}^{'} ∥_{\infty} = \sqrt{n + \frac{1}{2}} ∥ p_{n}^{'} ∥_{\infty} \leq \frac{1}{2} n (n + 1) \sqrt{n + \frac{1}{2}} .

(16)

As a result,

∥ {(f^{m} - f^{n})}^{'} ∥_{\infty} \leq \frac{1}{2} \sum_{k = n + 1}^{m} | α_{k} | k (k + 1) \sqrt{k + \frac{1}{2}} .

From Lemma 3.3, we know that

f \in C^{4} ([- 1, 1])

is enough to ensure

{{(f^{M})}^{'}}

is Cauchy. Let

{(f^{M})}^{'} \to φ^{'} \in C^{0} ([- 1, 1])

in

L^{\infty}

sense. The remaining issue is whether

φ^{'} = f^{'}

⁠. This is straightforward, since by Lemma 3.4 we know

f^{M} \to f

pointwisely. Moreover,

\begin{aligned} f^{M} (θ) - f^{M} (- 1) & = \int_{- 1}^{θ} {(f^{M})}^{'} (t) d t, \\ lim_{M \to + \infty} \int_{- 1}^{θ} {(f^{M})}^{'} (t) d t & = \int_{- 1}^{θ} φ^{'} (t) d t, \forall θ \in [- 1, 1], \end{aligned}

and the convergence is uniform with respect to

x

⁠. Thus,

\forall x \in [- 1, 1]

⁠,

lim_{M \to + \infty} f^{M} (θ) - f^{M} (- 1) = f (θ) - f (- 1) = \int_{- 1}^{θ} φ^{'} (t) d t,

that is,

f^{'} = φ^{'}

⁠.

So far, we have shown that with enough regularity of $f$ ⁠, we have the uniform convergence of both $f^{M} \to f$ and ${(f^{M})}^{'} \to f^{'}$ on $[- 1, 1]$ ⁠. To show the existence of minimizer near $θ^{†}$ ⁠, the following lemma is crucial.

Lemma 3.6

Let $F : D \subset R^{n} \to R$ be continuously differentiable in the domain $D$ and suppose that there is an open ball $B (x^{0}, r) \subset D$ and a positive $γ$ such that $∥ \nabla F {(x)}^{- 1} ∥ \leq γ \forall x \in B (x^{0}, r)$ and $r > γ ∥ F (x^{0}) ∥$ ⁠. Then, $F (x) = 0$ has a solution in $B (x^{0}, r)$ ⁠.

The proof can be found in Ref. 25. In the one-dimensional case, Lemma 3.6 basically says that if a $C^{1}$ function $f$ changes with a minimum speed (⁠ $| f^{'} | > γ^{- 1}$ ⁠) on the interval $(x^{0} - r, x^{0} + r)$ and the function value at $x^{0}$ is small enough [ $f (x^{0}) \leq \frac{r}{γ}$ ], then the function must reach zero in this interval. With Lemmas 3.4–3.6, we are able to prove the following result, which is the one-dimensional analogy of Theorem 3.1.

Proposition 3.7

Consider the following one-dimensional nonlinear least-squares problem:

min_{θ \in [- 1, 1]} f^{2} (θ),

with solution

θ^{†}

⁠. If

f

is approximated by

f^{M}

given by and we assume that

$f \in C^{4} ([- 1, 1]),$
$f (θ^{†}) = 0,$
$| f^{'} (θ^{†}) | \neq 0,$

then there exists a sequence of minimizer

{θ_{M}^{*}}

such that

{(f^{M})}^{2} (θ_{M}^{*}) = min_{θ \in [- 1, 1]} {(f^{M})}^{2} (θ) =: e_{M}^{2}

and

lim_{M \to + \infty} θ_{M}^{*} = θ^{†}, lim_{M \to + \infty} | f^{M} (θ_{M}^{*}) | = 0.

Proof.

From Lemmas 3.4 and 3.5, we know that the first assumption provides

lim_{n \to + \infty} ∥ f - f^{M} ∥_{\infty} = 0, lim_{n \to + \infty} ∥ f^{'} - {(f^{M})}^{'} ∥_{\infty} = 0.

Let

F^{M} (θ) := f^{M} (θ) - e_{M}

⁠, and we are going to apply Lemma 3.6 to

F^{M}

at

θ^{†}

⁠. Since

f^{'} (θ^{†}) \neq 0

⁠, by continuity, we know that there exist positive constants

γ

and

r

such that

| f^{'} (θ) |^{- 1} \leq \frac{γ}{2}, \forall θ \in (θ^{†} - r, θ^{†} + r) \subset (- 1, 1) .

(17)

Further notice that

{(F^{M})}^{'} = {(f^{M})}^{'}

and

∥ f^{'} - {(f^{M})}^{'} ∥_{\infty} \to 0

⁠. This implies that there exists a positive constant

L_{1}

such that

\forall M > L_{1}

| {(F^{M})}^{'} {(θ)}^{- 1} | \leq γ, \forall θ \in (θ^{†} - r, θ^{†} + r) \subset (- 1, 1) .

(18)

At

θ = θ^{†}

⁠, we have

| F^{M} (θ^{†}) | = | f^{M} (θ^{†}) - e_{M} | \leq | f^{M} (θ^{†}) | + e_{M} \leq 2 | f^{M} (θ^{†}) |,

where we used the fact that

e_{M} \leq | f^{M} (θ^{†}) |

⁠, since

e_{M}^{2}

is the minimum of

{(f^{M})}^{2}

on

[- 1, 1]

⁠. As a result,

lim_{M \to + \infty} | F^{M} (θ^{†}) | \leq lim_{M \to + \infty} 2 | f^{M} (θ^{†}) | = 2 | f (θ^{†}) | = 0.

Thus, we are able to select

L_{2} > 0

such that

\forall M > L_{2}

| F^{M} (θ^{†}) | < \frac{r}{γ} .

(19)

Since we have both and

\forall M > max {L_{1}, L_{2}} =: L

⁠, applying Lemma 3.6 to

F^{M}

at

θ^{†}

⁠, we conclude that

F^{M} (θ) = 0

has a solution in

(θ^{†} - r, θ^{†} + r)

⁠. Denote such solution as

θ_{M}^{*}

⁠, which is a minimizer of the approximated least-squares problem, that is,

{(f^{M})}^{2} (θ_{M}^{*}) = e_{M}^{2} = min_{θ \in [- 1, 1]} {(f^{M})}^{2} (θ), \forall M > L .

Apply the mean value theorem to

f (θ_{M}^{*}) - f (θ^{†})

for

M > L

⁠, we have

| θ^{†} - θ_{M}^{*} | \leq | f^{'} (ξ_{M}) |^{- 1} | f (θ^{†}) - f (θ_{M}^{*}) | = | f^{'} (ξ_{M}) |^{- 1} | f (θ_{M}^{*}) |,

for some

ξ_{M} \in (θ^{†} - r, θ^{†} + r)

⁠. Thus, the stability condition in leads to

\begin{aligned} | θ^{†} - θ_{M}^{*} | & \leq \frac{γ}{2} | f (θ_{M}^{*}) | \leq \frac{γ}{2} [| f (θ_{M}^{*}) - f^{M} (θ_{M}^{*}) | + | f^{M} (θ_{M}^{*}) |] \\ \leq \frac{γ}{2} [∥ f - f^{M} ∥_{\infty} + | f^{M} (θ^{†}) |] \leq \frac{γ ϵ}{\sqrt{M}} \to 0, \end{aligned}

as

M \to \infty

and for any

ϵ > 0

⁠. Here, we have used Lemma 3.4 to bound the first term and the residual error

| f^{M} (θ_{M}^{*}) | \leq | f^{M} (θ^{†}) | = | f^{M} (θ^{†}) - f (θ^{†}) | \leq \frac{ϵ}{\sqrt{M}} .

This result guarantees the existence of minimizers of the approximate polynomial least-squares problem that converges to the solution of the true least-squares problem under a fairly strong condition, that is, $f$ is smooth enough. In practice (see Sec. IV), we shall see reasonably accurate estimates even when the essential statistics are not as regular as required by this theoretical estimate.

C. Convergence of the algorithm

In Algorithm 1, we are solving a least-squares problem (13). For the Gauss-Newton method, a standard local convergence result under a full rank assumption of the Jacobian matrix is well known (e.g., see Chapter 2 of Ref. 26 for details). To achieve such a local convergence in our numerical algorithm, we will now verify the full rank condition in the polynomial chaos setting.

For our problem, the Jacobian matrix $J^{M} (θ) \in R^{K \times N}$ is given by

{(J^{M})}_{i j} = \frac{\partial f_{i}^{M}}{\partial θ_{j}}, i = 1, \dots, K, j = 1, \dots N,

where $f_{i}^{M}$ is an element of the finite-dimensional polynomial space $Γ_{N}^{M}$ defined by

Γ_{N}^{M} := span {θ_{1}^{k_{1}} θ_{2}^{k_{2}} \dots θ_{N}^{k_{N}} | k_{i} = 0, 1, \dots, M} .

Using linear independence over the function space $Γ_{N}^{M}$ and some fundamental results in algebra, we can prove the following theorem.

Theorem 3.8

If

F := {f_{i}}_{i = 1}^{K} \subset Γ_{N}^{M}

is linearly independent and

K > max {N, {(M + 1)}^{N - 1}},

then the column space of the corresponding Jacobian matrix

J (θ)

is full rank in

{(Γ_{N}^{M})}^{K}

⁠. Furthermore

,

the set

N (J)

is defined by

N (J) := {θ \in R^{N} | rank [J (θ)] \leq N - 1}

is nowhere dense over

R^{N}

⁠.

Proof.

See Appendix B.

Remark 3.9

Here, $N (J)$ is the set of $θ$ for which $J (θ)$ is not full rank in the matrix sense. In our application $,$ $N$ is the dimension of the parameter $θ,$ $K$ is the total number of equations in (9), and $M$ is the order of approximation used in (11). In practice, since the orthogonal polynomials involved in (11) form a basis for $Γ_{N}^{M}$ ⁠, it is enough to check the rank of the coefficient matrix $A^{M}$ ⁠, where the $i th$ column consists of the coefficients $α_{\vec{k}}^{(i)}$ in (11) for all possible $\vec{k}$ ⁠. The linear independence hypothesis is indeed sensible in practice since we want to avoid solving underdetermined least-squares problems.

IV. EXAMPLE I: THE LANGEVIN MODEL

For the first example, we consider a classical model in statistical mechanics: the dynamics of a particle driven by a conservative force, a damping force, and a stochastic force. In particular, we choose the conservative force based on the Morse potential

U (x) = U_{0} [a (x - x_{0})], U_{0} (x) = ϵ (e^{- 2 x} - 2 e^{- x} + 0.01 x^{2}),

(20)

where the last quadratic term in $U_{0}$ acts as a retaining potential (also known as a confining potential), preventing the particle from moving to infinity. For this one-dimensional model, we rescale the mass to unity $m = 1$ and write the dynamics as follows:

{\begin{cases} \dot{x} = v, \\ \dot{v} = - U^{'} (x) - γ v + \sqrt{2 γ k_{B} T} \dot{W}, \end{cases}

(21)

where $\dot{W}$ is a white noise. The generator of the system (21), denoted by $L$ ⁠, is given by

L = v \frac{\partial}{\partial x} + [- U^{'} (x) - γ v] \frac{\partial}{\partial v} + γ k_{B} T \frac{\partial^{2}}{\partial v^{2}} .

(22)

The smooth retaining potential $U (x)$ guarantees the ergodicity of the Langevin system (21) (see Appendix C for details). Namely, there is an equilibrium distribution (Gibbs measure) $p_{e q} (x, v)$ ⁠, given by

p_{e q} \propto \exp {- \frac{1}{k_{B} T} [U (x) + \frac{1}{2} v^{2}]},

(23)

which is independent of $γ$ ⁠. In particular, we have $v \sim N (0, k_{B} T)$ at equilibrium. For this Langevin model, there are a total of five unknown parameters $θ := (γ, T, ϵ, a, x_{0})$ with true values $θ^{†} := (γ^{†}, T^{†}, ϵ^{†}, a^{†}, x_{0}^{†})$ ⁠.

For this specific problem, notice that all parameters except $γ$ appear in the equilibrium density function. Intuitively, this suggests that one can estimate four parameters, ${T, ϵ, a, x_{0}}$ ⁠, from equilibrium statistics and $γ$ from two-point statistics, respectively. We will present a conventional estimation method based on this guideline for diagnostic purpose. In particular, the parameter $γ$ will be directly estimated from a two-point statistic as shown in Sec. IV A. Subsequently, the remaining four parameters are obtained via one-point statistics, as shown in (27) and Sec. IV B 1.

As for the proposed approach, we will consider estimating all the parameters using the essential statistics in the most efficient manner. First, two of the parameters, ${γ, T}$ ⁠, will be estimated using the essential statistics in (27). Then, via a sensitivity analysis discussed in Sec. IV B 2, we will show that parameter $x_{0}$ is independent of the essential statistics. As a consequence, we formulate a least-square problem to estimate $θ = (ϵ, a)$ and use equilibrium statistics to estimate $x_{0}$ ⁠.

A. Reduction of the parameter space

In our previous work,⁵ we have suggested an approach to estimate $T^{†}$ and $γ^{†}$ directly from the essential statistics. To construct such essential statistics, we consider a constant external forcing $δ f$ with $δ ≪ 1$ ⁠. The corresponding perturbed system is given by

{\begin{cases} \dot{x} = v, \\ \dot{v} = - U^{'} (x) - γ v + δ f + \sqrt{2 γ k_{B} T} \dot{W}, \end{cases}

(24)

that is, $c (x) = (0, 1)$ in the FDT formula (3). By selecting the observable $A = {(0, v)}^{⊤}$ ⁠, we work with the $(2, 2)$ entry of the response operator given by

\begin{aligned} k_{A} (t; θ^{†}) & = \frac{1}{k_{B} T^{†}} E_{{p_{e q}}^{†}} [v (t) v (0)] \\ = \frac{1}{k_{B} T^{†}} \int_{R^{2}} \int_{R^{2}} v v_{0} ρ (x, v, t | x_{0}, v_{0}, 0) \\ \times p_{e q}^{†} (x_{0}, v_{0}) d x d v d x_{0} d v_{0}, \end{aligned}

(25)

where $ρ$ is the solution of corresponding Fokker-Planck equation (5) associated with the Langevin dynamics. Taking the time derivative of $E_{{p_{e q}}^{†}} [v (t) v (0)]$ for $t > 0$ ⁠, we obtain

\begin{aligned} \frac{\partial}{\partial t} E_{{p_{e q}}^{†}} [v (t) v (0)] \\ = \int_{R^{2}} \int_{R^{2}} v v_{0} L^{*} ρ (x, v, t | x_{0}, v_{0}, 0) p_{e q}^{†} (x_{0}, v_{0}) d x d v d x_{0} d v_{0} \\ = \int_{R^{2}} \int_{R^{2}} L v v_{0} ρ (x, v, t | x_{0}, v_{0}, 0) p_{e q}^{†} (x_{0}, v_{0}) d x d v d x_{0} d v_{0} \\ = \int_{R^{2}} \int_{R^{2}} [- U^{'} (x) - γ v] v_{0} ρ (x, v, t | x_{0}, v_{0}, 0) \\ \times p_{e q}^{†} (x_{0}, v_{0}) d x d v d x_{0} d v_{0} \\ = E_{{p_{e q}}^{†}} ({- U^{'} [x (t)] - γ^{†} v (t)} v (0)) . \end{aligned}

(26)

Let $t \to 0^{+}$ in (26) and recall that $v \sim N (0, k_{B} T)$ at equilibrium, we have

\begin{aligned} \begin{aligned} E_{p_{e q}^{†}} [v (t) v (0)] |_{t = 0} & = k_{B} T^{†}, \\ \frac{\partial}{\partial t} E_{p_{e q}^{†}} [v (t) v (0)] |_{t = 0^{+}} & = - γ^{†} k_{B} T^{†}, \end{aligned} \end{aligned}

(27)

where both the left-hand-sides can be computed from the sample. Thus, (27) provides direct estimates for $T^{†}$ and $γ^{†}$ ⁠. As a result, the original parameter estimation problem can be reduced into estimating $(ϵ^{†}, a^{†}, x_{0}^{†})$ in the potential function $U$ ⁠, which is a non-trivial task since the dependence of $p_{e q}$ on these parameters is quite complicated. With these estimates, $B^{†}$ becomes available since it only depends on $k_{B} T^{†}$ ⁠, which allows one to compute $M_{j} (t_{i})$ in (8) and avoid the issue pointed out in Remark 2.4.

B. Parameter estimation approaches

In this subsection, we focus on the three-dimensional parameter estimation problem with $θ := (ϵ, a, x_{0})$ ⁠. First, we review the conventional approach, which matches the equilibrium statistics of $x$ ⁠. Then, we discuss the proposed new approach using the essential statistics. In particular, we perform a sensitivity analysis to determine the identifiability of the parameters from the essential statistics. Finally, we present the numerical results including a comparison between the conventional and the new approaches.

1. Conventional method

We first look at what we can learn from the equilibrium statistics. Since the marginal distribution of $x$ at equilibrium state is proportional to $\exp [- U (x; θ) / k_{B} T]$ ⁠, a natural idea is to match the moments of $x$ with respect to this equilibrium density. In particular, one can introduce the following three equations:

E_{p_{e q}^{†}} [x^{i}] = E_{p_{e q}} [x^{i}; ϵ, a, x_{0}], i = 1, 2, 3,

(28)

where the left-hand-sides are computed from the given data, while the right-hand-sides are treated as functions of $(ϵ, a, x_{0})$ ⁠, obtained from solving the model (24). The following pre-computation simplifies Eq. (28) into a one-dimensional problem.

To begin with, we define the probability density functions

\begin{aligned} p_{e q}^{a, x_{0}} (x) & := \frac{1}{N} \exp {- \frac{1}{k_{B} T} U_{0} [a (x - x_{0})]}, \\ p_{e q}^{1, 0} (x) & := \frac{1}{N_{0}} \exp [- \frac{1}{k_{B} T} U_{0} (x)], \end{aligned}

in terms of $U$ and $U_{0}$ ⁠, respectively. With a change of variables $y := a (x - x_{0})$ ⁠, the normalizing constants $N$ and $N_{0}$ satisfy

\begin{aligned} N & = \int_{R} \exp {- \frac{1}{k_{B} T} U_{0} [a (x - x_{0})]} d x \\ = \frac{1}{a} \int_{R} \exp [- \frac{1}{k_{B} T} U_{0} (y)] d y = \frac{1}{a} N_{0} . \end{aligned}

As a result, the first two equations in (28) can be written into

\begin{aligned} E_{p_{e q}^{†}} [x] & = E_{p_{e q}}^{a, x_{0}} [x] = \frac{1}{a} E_{p_{e q}}^{1, 0} [x] + x_{0}, \\ E_{p_{e q}^{†}} [x^{2}] & = E_{p_{e q}}^{a, x_{0}} [x^{2}] = \frac{1}{a^{2}} E_{p_{e q}}^{1, 0} [x^{2}] + \frac{2 x_{0}}{a} E_{p_{e q}}^{1, 0} [x] + x_{0}^{2}, \end{aligned}

with a unique solution

a = \sqrt{\frac{{Var}_{p_{e q}}^{1, 0} [x]}{{Var}_{{p_{e q}}^{†}} [x]}}, x_{0} = E_{{p_{e q}}^{†}} [x] - E_{p_{e q}}^{1, 0} [x] \sqrt{\frac{{Var}_{{p_{e q}}^{†}} [x]}{{Var}_{p_{e q}}^{1, 0} [x]}},

(29)

where ${Var}_{p_{e q}}^{1, 0} [x]$ stands for the variance of $x$ with respect to the rescaled equilibrium density $p_{e q}^{1, 0}$ ⁠. In practice, both $E_{{p_{e q}}^{†}} [x]$ and ${Var}_{{p_{e q}}^{†}} [x]$ can be empirically estimated directly from the available data, while ${Var}_{p_{e q}}^{1, 0} [x]$ is a function of $ϵ$ only. Thus, (29) can be denoted as

a = a (ϵ); x_{0} = x_{0} (ϵ),

and (28) is reduced to a one-dimension problem for $ϵ$ ⁠,

\begin{aligned} E_{p_{e q}^{†}} [x^{3}] & = E_{p_{e q}} [x^{3}; ϵ, a, x_{0}] [= \int_{R} x^{3} p_{e q} (x; ϵ, a, x_{0}) d x], \\ s.t. a & = a (ϵ), x_{0} = x_{0} (ϵ) . \end{aligned}

(30)

Note that the left-hand-side of (30) will be estimated by Monte-Carlo averaging from the unperturbed data, whereas the right-hand-side integral will be approximated using a numerical quadrature.

2. Sensitivity analysis of the essential statistics

To ensure a successful parameter estimation, it is useful to gather some a priori knowledge of the parameter identifiability from the proposed essential statistics. This is important especially since there are non-unique choices of essential statistics to be fit and we want to ensure that we can identify the parameters from appropriate essential statistics: Recall that the essential statistics are defined based on the choices of external forcing and observable.

While local sensitivity analysis such as the pathwise derivative method described in Ref. 27 is desirable, it requires knowledge of the true parameters, which are not available to us. Essentially, the pathwise derivative method is to compute ${\bar{y}}_{θ} := E [Y (t, ω)]$ ⁠, where $Y (t, ω) := D_{θ} X (t, θ)$ and the expectation is defined with respect to the density of the Itô diffusion (written in its integral form)

X (t, θ) = X_{0} (θ) + \int_{0}^{t} b [X (s, θ), θ] d s + \int_{0}^{t} σ [X (s, θ), θ] d W_{s} .

(31)

Since we have no explicit solutions of the density $ρ (x, t; θ)$ ⁠, one way to approximate ${\bar{y}}_{θ}$ is using an ensemble average of the solutions of the following equation:

\begin{aligned} Y (t, θ) & = Y_{0} (θ) + \int_{0}^{t} {b_{X} [X (s, θ), θ] Y (s, θ) + b_{θ} [X (s, θ), θ]} d s \\ + \int_{0}^{t} {σ_{X} [X (s, θ), θ] Y (s, θ) + σ_{θ} [X (s, θ), θ]} d W_{s} . \end{aligned}

(32)

Notice that this equation depends on (31) and the unknown parameters, $θ$ ⁠. The dependence on $θ$ implies that one cannot use this method for a priori sensitivity analysis. However, it can be used for a posteriori sensitivity analysis, evaluated at the estimates, to verify the difficulties of the particular parameter regimes.

As an empirical method for a priori sensitivity analysis, we simply check the deviation of $\hat{k} (t, θ)$ as a function of the training collocation nodes, $θ \in Θ$ ⁠, which are available from Step 2 of Algorithm 1. Of course, a more elaborate global sensitivity analysis technique such as the Sobol index²⁸ can be performed as well, but we will not pursue this here. From this empirical a priori sensitivity analysis, we found that the corresponding essential statistics, $E_{p_{e q}} [v (t) v (0)]$ is not sensitive to $x_{0}$ ⁠. However, it is sensitive to $a$ and strongly sensitive to $ϵ$ ⁠. Also, the sensitivity of the low damping case (⁠ $γ = 0.5$ ⁠) is stronger than the high damping case (⁠ $γ = 0$ ⁠). To verify the validity of this empirical method, we compute the local sensitivity analysis indices ${\bar{y}}_{θ} := E [Y (t, ω)]$ ⁠, where $Y (t, ω)$ solves (32) with the true parameters (again this is not feasible in practice since the true parameters are not available).

For our Langevin model, we have the corresponding drift and diffusion coefficients

b = {[v, - U^{'} (x) - γ v]}^{⊤}, σ = {(0, \sqrt{2 γ k_{B} T})}^{⊤} .

Suppose we are interested in the dependence on parameter $x_{0}$ ⁠. Introducing the notations $x_{x_{0}} := \frac{\partial}{\partial x_{0}} x (t, x_{0})$ and $v_{x_{0}} := \frac{\partial}{\partial x_{0}} v (t, x_{0})$ ⁠, the joint system of (31) and (32) (written in differential form) is given as follows:

{\begin{cases} \dot{x} = v, \\ \dot{v} = - U^{'} (x; x_{0}) - γ v + \sqrt{2 γ k_{B} T} \dot{W}, \\ {\dot{x}}_{x_{0}} = v_{x_{0}}, \\ {\dot{v}}_{x_{0}} = - U^{″} (x; x_{0}) x_{x_{0}} - γ v_{x_{0}} - \frac{\partial}{\partial x_{0}} U^{'} (x; x_{0}) . \end{cases}

(33)

There is no stochastic term in the last equation, since the noise in the Langevin model is additive with coefficient independent of $x_{0}$ ⁠. Notice that $(x_{x_{0}}, v_{x_{0}}) \equiv (1, 0)$ is a fixed point for the last two equations in (33) because of

\begin{aligned} \frac{\partial}{\partial x_{0}} U^{'} (x; x_{0}) & = a \frac{\partial}{\partial x_{0}} U_{0}^{'} [a (x - x_{0})] \\ = - a^{2} U_{0}^{″} [a (x - x_{0})] = - \frac{\partial}{\partial x} U^{'} (x; x_{0}) . \end{aligned}

As a result, by choosing the initial condition to be $[x (0), v (0)] = (x_{0}, c)$ for an arbitrary constant $c$ ⁠, one can claim that the solution $v (t)$ of the Langevin model (21) is independent of $x_{0}$ ⁠. Thus, under this circumstance, $E_{p_{e q}^{†}} [v (t) v (0)]$ is independent of $x_{0}^{†}$ ⁠, confirming the conclusion from our empirical a priori sensitivity analysis.

While this parameter insensitivity may seem discouraging, we can use it to our advantage to simplify the problem. That is, we can assign an arbitrary value to $x_{0}$ in applying Algorithm 1 and consider estimating only $θ = (a, ϵ)$ ⁠. Once these two parameters are estimated, we can use the formula in (29) to estimate $x_{0}$ ⁠.

We also use this local sensitivity analysis to verify the validity of the a priori sensitivity analysis with respect to $ϵ$ and $a$ ⁠. Since it is difficult to analyze explicit behaviors as above, for a given realization of the Langevin model, we solve the remaining ODEs by the fourth order Runge-Kutta method. Figure 2 shows the numerical results of ${\bar{v}}_{a}$ and ${\bar{v}}_{ϵ}$ ⁠, where the average is estimated over 3000 realizations of the Langevin model. Based on the scales of ${\bar{v}}_{a}$ and ${\bar{v}}_{ϵ}$ in Fig. 2, we confirm the validity of the sensitivity of the two-point statistics $E_{p_{e q}^{†}} [v (t) v (0)]$ that was found empirically. Namely, the identifiability of $ϵ$ is stronger than $a$ and the identifiability of the low damping case (⁠ $γ = 0.5$ ⁠) is stronger than the high damping case (⁠ $γ = 5.0$ ⁠).

FIG. 2.

View large Download slide

The absolute value of ${\bar{v}}_{a}$ and ${\bar{v}}_{ϵ}$ for both $γ = 0.5$ and $γ = 5.0$ ⁠.

C. Numerical results

We now present several numerical results. Based on a time series of size $10^{7}$ (with temporal lag $h = 2 \times 10^{- 3}$ ⁠), we compare the conventional method and the proposed scheme that fits the essential statistics under $γ^{†} = 0.5$ and $5.0$ ⁠, representing a low damping case and a high damping case, respectively. The true values of the remaining four parameters are set to $(k_{B} T^{†}, ϵ^{†}, a^{†}, x_{0}^{†}) = (1, 0.2, 10, 0)$ ⁠. The data, in the form of a time series, are generated from the model (24) using an operator-splitting method,²⁹ which is a strong second-order method.

1. Low damping case

We start with the conventional method. To examine the sensitivity of the estimation for $(ϵ^{†}, a^{†}, x_{0}^{†})$ to the estimates $\hat{T}$ ⁠, as a comparison, we also include the estimation in which we use the true value $k_{B} T^{†} = 1$ in solving (30) (we name the results of this scenario as partial estimates). The results have been listed in Table I. We clearly observe that the results of the conventional method are sensitive to the value of $\hat{T}$ ⁠. To explain the sensitivity, recall that in solving (30), $a$ and $x_{0}$ are viewed as functions of $ϵ$ ⁠. This suggests that the high sensitivity would occur if

\begin{aligned} \frac{\partial ϵ}{\partial k_{B} T} & = {(\frac{\partial E_{p_{e q}} [x^{3}; T, ϵ, a (ϵ), x_{0} (ϵ)]}{\partial ϵ})}^{- 1} \\ \cdot \frac{\partial E_{p_{e q}} [x^{3}; T, ϵ, a (ϵ), x_{0} (ϵ)]}{\partial k_{B} T} ≫ 1, \end{aligned}

where $a (ϵ)$ and $x_{0} (ϵ)$ are given by (29) and the implicit function theorem has been used. By direct computation, we get

\begin{aligned} \frac{\partial}{\partial k_{B} T} E_{p_{e q}} [x^{n}] = \frac{1}{{(k_{B} T)}^{2}} {E_{p_{e q}} [x^{n} U (x)] - E_{p_{e q}} [x^{n}] E_{p_{e q}} [U (x)]}, \\ n = 1, 2, \dots . \end{aligned}

Since the explicit formula for the partial derivative with respect to $ϵ$ is not easy to compute, its value will be computed using a finite difference formula. Evaluating it at $(k_{B} T, ϵ, a, x_{0}) = (1, 0.2, 10, 0)$ ⁠, we have

\begin{aligned} \frac{\partial}{\partial k_{B} T} E_{p_{e q}} [x^{3}; T, ϵ, a (ϵ), x_{0} (ϵ)] = 9.034, \\ \frac{\partial}{\partial ϵ} E_{p_{e q}} [x^{3}; T, ϵ, a (ϵ), x_{0} (ϵ)] = 0.06396, \end{aligned}

which supports our observation.

TABLE I.

Full and partial estimates of the conventional method using the equilibrium statistics (above) and the estimates using essential statistics (below): the low damping case.

	$k_{B} T$	$γ$	$ϵ$	$a$	$x_{0}$
True	1.0000	0.50000	0.20000	10.000	0.0000
Eq. stat. full estimates	0.99903	0.50003	0.16870	10.888	0.001120
Eq. stat. partial estimates	$\dots$	0.50003	0.20039	10.005	0.011263
Essential stat estimates	0.99903	0.50003	0.20004	9.9741	0.008265

	$k_{B} T$	$γ$	$ϵ$	$a$	$x_{0}$
True	1.0000	0.50000	0.20000	10.000	0.0000
Eq. stat. full estimates	0.99903	0.50003	0.16870	10.888	0.001120
Eq. stat. partial estimates	$\dots$	0.50003	0.20039	10.005	0.011263
Essential stat estimates	0.99903	0.50003	0.20004	9.9741	0.008265

Next, we report the estimations obtained from the new approach using essential statistics, and the results are listed in Table I. The improvement in accuracy for $ϵ^{†}$ and $a^{†}$ is noticeable compared with the full estimates of the conventional method. Here, we used a total of $20$ essential statistics given by $k_{A} (t_{i}; θ^{†})$ for $t_{i} = 0.1 i$ ⁠, $i = 1, 2, \dots, 20$ ⁠. We applied Algorithm 1 to solve an ensemble of two-dimensional (⁠ $a$ and $x_{0}$ ⁠) nonlinear least-squares problems based on $300$ uniformly generated random initial conditions. Figure 3 (left) shows the contour plot of the cost function together with all of the estimates. Notice that except for few outliers, most of the estimates lie along the low value of the contour of the cost function. The estimate $\hat{ϵ}$ reported is the average of all the estimates (excluding the outliers). To satisfy the equilibrium constraint, $a (\hat{ϵ})$ given by (29) is used as the estimates of $a^{†}$ ⁠.

FIG. 3.

View large Download slide

The contour plot of the cost function together with the estimates: low damping case (left) and high damping case (right).

The results from this test indicate that for the conventional method, it is difficult to obtain very accurate estimates for $ϵ^{†}$ ⁠, unless $T^{†}$ can be estimated accurately, e.g., by using longer time series. In contrast, the proposed approach of using essential statistics is much less sensitive to the error in $\hat{T}$ ⁠. This approach, however, requires a pre-computing step as noted in Algorithm 1. In our numerical test, we solved the model (21) on 64 collocation nodes to evaluate the essential statistics over different values of $(ϵ, a)$ ⁠, that is, an order $M_{C} = 8$ Chebyshev nodes were used to construct the $Θ$ used in (12). However, as we have previously alluded to, this can be done in parallel, and it will not become a serious issue as long as the dimension of the parameter is relatively small. As for the value of $M$ in (12), we picked $M = 6$ ⁠. (The same scheme was applied to the high damping case.)

Another interesting issue arises when the damping parameter is large. Since the conventional method fully relies on one-point statistics with respect to equilibrium density, it is important to have high-quality independent samples. In Fig. 4(left), we show the time correlation of $x$ for both the low damping regimes (⁠ $γ = 0.5$ ⁠) and high damping regime (⁠ $γ = 5.0$ ⁠). We observe that in the latter case, the auto-correlation function of $x$ decays much slower, indicating a strong correlation among the samples with small lags. In this case, the estimates from the conventional method will deteriorate due to the difficulty in obtaining high quality independent sampling.

FIG. 4.

View large Download slide

The time auto-correlation of $x$ (left) and $v$ (right).

2. High damping case

In this section, we focus on the high damping regime $γ^{†} = 5.0$ ⁠. We will also make a connection between Fig. 2 and the estimation using the essential statistics. As shown in Fig. 4 (left), we no longer have a fast decay of the time auto-correlation of $x$ ⁠. As a result, without changing the sample size, significant error will occur in estimating the moment of $x$ used in the conventional method. This can be clearly seen from the results listed in Table II. In particular, besides suffering from the sensitivity to $\hat{T}$ ⁠, the error in estimating $E_{p_{e q}^{†}} [x^{n}]$ for $n = 1, 2, 3$ also leads to inaccurate estimates for $ϵ^{†}$ ⁠, $a^{†}$ ⁠, and $x_{0}^{†}$ ⁠. We should point out an interesting fact, that is, although the marginal distribution of $x$ at the equilibrium state is independent of $γ$ ⁠, a large value of $γ$ causes difficulties in estimating the moments of $x$ in practice.

TABLE II.

Full and partial estimates of the conventional method using the equilibrium statistics (above) and the estimates using essential statistics (below): the high damping case.

	$k_{B} T$	$γ$	$ϵ$	$a$	$x_{0}$
True	1.0000	5.0000	0.20000	10.000	0.0000
Eq. stat. full estimates	0.99949	4.9993	0.69579	5.4717	0.12990
Eq. stat. partial estimates	$\dots$	4.9993	0.25931	8.9295	−0.0034460
Essential stat estimates	0.99949	4.9993	0.21160	9.8663	−0.01822

	$k_{B} T$	$γ$	$ϵ$	$a$	$x_{0}$
True	1.0000	5.0000	0.20000	10.000	0.0000
Eq. stat. full estimates	0.99949	4.9993	0.69579	5.4717	0.12990
Eq. stat. partial estimates	$\dots$	4.9993	0.25931	8.9295	−0.0034460
Essential stat estimates	0.99949	4.9993	0.21160	9.8663	−0.01822

Table II and Fig. 3 (right) show the numerical results using the essential statistics. Although the relative error is not as small as in the low damping case, it is still lower than the estimates from the conventional method. Here is how our method is implemented: Similar to the low damping case, (29) is used in estimating the value of $a^{†}$ based on the estimates $\hat{ϵ}$ ⁠. The 20 essential statistics here are given by $k_{A} (t_{i}; θ^{†})$ for $t_{i} = 0.04 i + 0.2$ ⁠, $i = 1, 2, \dots, 20$ [suggested by the time auto-correlation of $v$ shown in Fig. 4 (right)], which are in a much shorter time interval. Compared to the low damping case, the loss of accuracy can be verified based on the sensitivity analysis (Fig. 2), which indicated that the parameters $(a, x_{0})$ in this regime are less identifiable when $γ$ is large.

V. EXAMPLE II: A GRADIENT SYSTEM WITH A TRIPLE-WELL POTENTIAL

Our next example is a gradient system driven by white noise. We consider a two-dimensional stochastic system as follows:

d x (t) = - C \nabla V (x) d t + \sqrt{2 k_{B} T} d W_{t},

(34)

where $W_{t}$ is a two-dimensional Wiener process, $V$ is a potential energy, and $C$ is a matrix defined by

C = (\begin{matrix} 1 & - d \\ d & 1 \end{matrix}), d \in (- 1, 1) .

The generator $L$ of the system (34) is given by

L f = - (V_{x_{1}} - d V_{x_{2}}) \frac{\partial f}{\partial x_{1}} - (d V_{x_{1}} + V_{x_{2}}) \frac{\partial f}{\partial x_{2}} + k_{B} T Δ f .

(35)

We choose a triple-well potential function $V$ similar to the model in Ref. 30,

\begin{aligned} V (x_{1}, x_{2}) & = - v (x_{1}^{2} + x_{2}^{2}) - (1 - γ) v [{(x_{1} - 2 a)}^{2} + x_{2}^{2}] \\ - (1 + γ) v [{(x_{1} - a)}^{2} + {(x_{2} - a \sqrt{3})}^{2}] \\ + 0.2 [{(x_{1} - a)}^{2} + {(x_{2} - a / \sqrt{3})}^{2}], \end{aligned}

(36)

with

v (z) = 10 \exp (\frac{1}{z^{2} - a^{2}}) \cdot χ_{(- a, a)} (z), z \in R,

where $χ_{(- a, a)} (z)$ denotes the characteristic function over the interval $(- a, a)$ ⁠. Notice that the matrix $C$ is positive definite. The additional quadratic term $0.2 [{(x_{1} - a)}^{2} + {(x_{2} - a / \sqrt{3})}^{2}]$ in the triple-well potential (36) is, again, a smooth retaining potential. It is well known¹ that the triple-well model (34) yields an equilibrium distribution given by

p_{e q} (x) \propto \exp (- \frac{V (x)}{k_{B} T}) x \in R^{2},

(37)

which is independent of parameter $d$ ⁠, the off-diagonal element of $C$ ⁠.

In the numerical tests, we set $(d^{†}, a^{†}, k_{B} T^{†}, γ^{†}) = (0.5, 1, 1.5, 0.25)$ as the true values of the parameters. Figure 5 shows the contour plot of the potential and the scatter plot of the time series under this set of parameters. To generate the data from (34), we applied the weak trapezoidal method introduced in Ref. 31, which is a weak second-order method. In the remainder of this section, we will choose appropriate essential statistics from which ${d, T}$ can be directly estimated and the remaining parameters, ${a, γ}$ ⁠, will be estimated using an appropriate least-square problem.

FIG. 5.

View large Download slide

The contour plot of the potential (left) and the scatter plot of the solution with $10^{3}$ sample points (right).

A. The essential statistics and reduction of the parameter space

Similar to the Langevin model, we consider an external forcing that is constant in $x$ ⁠. Subsequently, the perturbed dynamics is given by

d x (t) = [- C \nabla V (x) + f (t) δ] d t + \sqrt{2 k_{B} T} d W_{t},

(38)

where $| δ | ≪ 1$ ⁠. If we select $A (x) := x$ as the observable, the corresponding linear response operator reads

k_{A} (t; θ) = E_{p_{e q}} {A [x (t)] \otimes B [x (0)]},

where $B = {(B_{1}, B_{2})}^{⊤}$ with

B_{i} (x_{1}, x_{2}) = \frac{1}{k_{B} T} \frac{\partial}{\partial x_{i}} V (x_{1}, x_{2}), i = 1, 2,

given by (3). As a result, the entries of the linear response operator $k_{A} (t; θ)$ satisfy

k_{(i, j)} (t; θ) := \frac{1}{k_{B} T} E_{p_{e q}} {x_{i} (t) V_{x_{j}} [x_{1} (0), x_{2} (0)]}, i, j = 1, 2,

(39)

where $V_{x_{i}}$ denotes the partial derivative of $V$ with respect to $x_{i}$ ⁠. Using integration by parts, one can show that

\begin{aligned} \int_{R} x_{i} V_{x_{i}} (x_{1}, x_{2}) \exp [- V (x_{1}, x_{2}) / k_{B} T] d x_{j} = δ_{i j} k_{B} T \\ \int_{R} \exp [- V (x_{1}, x_{2}) / k_{B} T] d x_{j}, i, j = 1, 2, \end{aligned}

which leads to

k_{i, j} (0; θ) = δ_{i j} .

(40)

This is known as the equipartition of the energy in statistical mechanics.

Notice that the response operators $k_{(i, j)} (t; θ^{†})$ are not accessible due to the fact that the function $V_{x_{i}} (x, θ^{†})$ depends on unknown true parameter values $(k_{B} T^{†}, a^{†}, γ^{†})$ ⁠. This is precisely one of the issues raised in Remark 2.4. For gradient flow systems, this problem can be overcome by introducing a linear transformation to another set of two-point statistics. We define

m_{i, j} (t; θ) := E_{p_{e q}} [x_{i} (t) x_{j} (0)], i, j = 1, 2

(41)

and consider their time derivatives. Following the same calculations that led to (26), we obtain the following identities for $t > 0$ ⁠:

\begin{aligned} \frac{d}{d t} m_{i, j} (t; θ) = k_{B} T [- k_{j, i} (- t; θ) + {(- 1)}^{i - 1} d \cdot k_{j, 3 - i} (- t; θ)], \\ i, j = 1, 2. \end{aligned}

(42)

It is more helpful to rewrite (42) into a linear system for $t > 0$ ⁠,

k_{B} T (\begin{matrix} - 1 & d & 0 & 0 \\ 0 & 0 & - 1 & d \\ - d & - 1 & 0 & 0 \\ 0 & 0 & - d & - 1 \end{matrix}) (\begin{matrix} k_{1, 1} (- t) \\ k_{1, 2} (- t) \\ k_{2, 1} (- t) \\ k_{2, 2} (- t) \end{matrix}) = (\begin{matrix} m_{1, 1}^{'} (t) \\ m_{1, 2}^{'} (t) \\ m_{2, 1}^{'} (t) \\ m_{2, 2}^{'} (t) \end{matrix}),

(43)

where the coefficient matrix is non-singular since $d \in (- 1, 1)$ ⁠. This linear relationship suggests that one can consider fitting $m_{i, j}$ in place of $k_{i, j}$ ⁠, since the former is numerically accessible.

Let $t \to 0^{+}$ and apply (40) to the first and the third equations of (43) under $θ = θ^{†}$ ⁠, we obtain the following estimates for $d^{†}$ and $T^{†}$ ⁠:

k_{B} T^{†} = - m_{1, 1}^{'} (0^{+}; θ^{†}); d^{†} = - m_{2, 1}^{'} (0^{+}; θ^{†}) / k_{B} T^{†},

(44)

where $m_{i, j}^{'} (0^{+}; θ^{†})$ is computable from the available sample. Equation (44) reduces the problem into a two-dimensional problem of estimating $(a^{†}, γ^{†})$ only.

B. Numerical results

As we have pointed out earlier, the linear response operator, $k_{(i, j)} (t; θ^{†})$ (39), is not directly accessible. Motivated by the linear relation in (43), we consider to infer the values $(a^{†}, γ^{†})$ from the two-point statistics $m_{i, j} (t; θ^{†})$ (41). In particular, we choose $m_{1, 1} (t_{i}; θ^{†})$ with $t_{i} = 0.1 i$ ⁠, $i = 1, \dots, 20$ ⁠, as the essential statistics in our numerical test.

In light of the fact that $m_{1, 1} (t; θ)$ is the two-point statistics of $x_{1}$ ⁠, we can apply the empirical a priori sensitivity test by checking the dependence of $\hat{k} (t, θ)$ over the training collocation nodes $θ = (a, γ) \in Θ$ ⁠, which are available to us from Step 2 of Algorithm 1. In this numerical experiment, the training used 25 collocation nodes, $θ \in Θ$ ⁠. To validate this empirical sensitivity test, we also perform the local sensitivity analysis as described in Sec. IV B 2 to parameters $a$ and $γ$ (which, again, is not possible in practice since the true parameters are unknown). The results are shown in Fig. 6 (left). From the scales of ${\bar{x}}_{a}$ and ${\bar{x}}_{γ}$ (average over $2000$ realizations of $x$ ⁠), we can see that this two-point statistic, $m_{1, 1} (t; θ^{†})$ ⁠, is quite sensitive to the parameters $γ$ and $a$ ⁠, with stronger dependence on the parameter $a$ ⁠.

FIG. 6.

View large Download slide

Local sensitivity test of $x_{1}$ with respect to parameter $a$ and $γ$ (left) and the contour plot of the cost function together with the estimates (right).

Based on a time series of size $4 \times 10^{6}$ (with time lag $h = 1 \times 10^{- 3}$ ⁠), we implemented Algorithm 1 with $M = 4$ and collocation nodes, $Θ$ ⁠, constructed from order $M_{C} = 5$ Chebyshev nodes in (12).

The estimates are shown in Table III. We also show in Fig. 6 (right) the contour plot of the cost function together with the estimates based on $300$ uniformly generated initial guesses. Notice that the true parameter value does not lie on the lowest contour value and these discrepancies are due to the surrogate modeling and quality of the samples. However, the estimates $\hat{a}$ and $\hat{γ}$ are still reasonably accurate, as reported in Table III. We should point out that the estimates $\hat{a}$ and $\hat{γ}$ reported in this table are the average of all the estimates (excluding the outliers), while the estimates $k_{B} \hat{T}$ and $\hat{d}$ are obtained by solving (44). The contour plot also confirms the sensitivity analysis, which suggested that estimating $a^{†}$ is easier than $γ^{†}$ ⁠.

TABLE III.

The estimation results: the Triple-Well model.

	$a$	$γ$	$k_{B} T$	$d$
True	1.0000	0.2500	1.5000	0.5000
Estimates	1.0249	0.2463	1.4946	0.5077

One interesting question is whether the approximate values of the parameters can reproduce the essential statistics, which is useful in predicting non-equilibrium averages in the presence of external forces. Using the estimates reported in Table III, we compare the resulting two-point statistics used in estimating the parameters, $m_{1, 1} (t; \hat{θ})$ ⁠, with the corresponding true statistics, $m_{1, 1} (t; θ^{†})$ (see Fig. 7, left). In Fig. 7 (right), we also compare the estimated response operator, $k_{1, 1} (t, \hat{θ})$ ⁠, and the true response operator, $k_{1, 1} (t, θ^{†})$ ⁠. Excellent agreement is found.

FIG. 7.

View large Download slide

The recovery of the response statistics $E_{p_{e q}} [x_{1} (t) x_{1} (0)]$ (left) and the two-point statistics $E_{p_{e q}} {x_{1} (t) V_{x_{1}} [x_{1} (0), x_{2} (0)]}$ (right). All the trajectories have been normalized such that they start at point $(0, 1)$ ⁠.

VI. SUMMARY AND FURTHER DISCUSSION

This paper presented a parameter estimation method for stochastic models in the form of Itô drift diffusions. To infer parameters that do not appear in the equilibrium density function, we proposed to fit two-point essential statistics formulated using the fluctuation dissipation theory. Building upon the framework established in our previous work,⁵ we formulated the problem as a nonlinear least-squares problem subjected to a dynamical constraint when the parameters cannot be estimated directly by solving the corresponding statistical constraints. To avoid expensive computational cost in evaluating the essential statistics at each iteration when the Gauss-Newton method is used, we proposed to solve an approximate least-squares problem based on the polynomial surrogate modeling approach. This approach is motivated by the fact that the sampling error cannot be avoided in computing the value of essential statistics and the essential statistics are smooth functions of the parameters. We guaranteed the existence of minimizers of the approximate least-squares problem that converge to the solution of the true least-squares problem that involves the essential statistics under the assumption that these statistics are smooth functions of the parameters. We also showed that the polynomial approximate least-squares problem has a Jacobian that is full rank almost everywhere, which implies the local convergence of Gauss-Newton solutions. We tested the proposed methods on two examples that belong to two large classes of stochastic models—a Langevin dynamics model and a stochastic gradient system. In general, we expect that the parameter estimation procedure should be carried out as follows:

Reduce the parameter space by direct estimation using appropriate statistics that are easily computable, whenever this is possible.
Based on the observable of interest, the functional form of the equilibrium density, the external forcing term, and the available data, identify appropriate essential statistics for the remaining parameters using the a priori sensitivity analysis test. In our implementation, we compare the essential statistics computed on the training parameter values (e.g., collocation nodes). A more elaborate global sensitivity analysis technique such as the Sobol index²⁸ can be used to determine the parameter identifiability.
For the remaining parameters, formulate a nonlinear least-squares problem using the appropriate essential statistics at some training parameter values.
Apply Algorithm 1 to obtain the estimates $\hat{θ}$ of the true parameter values $θ^{†}$ ⁠.
Apply a sensitivity analysis, such as the local sensitivity analysis as discussed in Sec. IV B 2, as a posteriori confidence check for the estimates. If the parameters are found to be insensitive to the selected response statistics, go back to Step 3 and choose a different set of response statistics.

One of the restrictions with our formulation is to be able to compute $B^{†}$ in (7), which may require knowledge of (some of) the true values $θ^{†}$ and this is not feasible in general. As a result, the non-linear system (8) cannot be evaluated without knowing the true parameters. For the Langevin dynamics, we found that the parameter $k_{B} T^{†}$ in $B^{†}$ can be estimated with the equilibrium statistics, variance of $v$ ⁠. For the gradient flow problem, there exists an invertible linear transformation between the linear response statistics that are not computable and other two-point statistics that can be estimated without knowing the true parameter values. Since numerical Algorithm 1 does not rely on the FDT formulation, one can apply it on any available two-point statistics, e.g., time auto-correlations of the solution. We are not aware of a general routine that bypasses this difficulty. However, as we have demonstrated, at least for general Langevin dynamics models and stochastic gradient systems, this idea will work.

While the proposed method requires the knowledge of the equilibrium density function of the unperturbed system, one can relax this condition and estimate it from the data. For low-dimensional problems, one can use the non-parametric kernel density estimation method.^32,33 If the invariant distribution of the dynamical system can be characterized by a density function defined on a smooth manifold (embedded in a high-dimensional phase space) that the data lie on or are close to, then one can estimate the density function by applying the theory of the kernel embedding of distribution³⁴ on a Hilbert space with basis functions defined on the data manifold. Analogous to the conditional density estimation proposed in Refs. 37 and 38, these data-driven basis functions can be obtained via the diffusion maps algorithm.^35,36 Another alternative to nonparametric methods is to consider a parametric density function estimation using the moment constrained maximum entropy principle.³⁹ For densities with moderate, say five to seven, dimensions, one can use a recently developed Equation-By-Equation (EBE) method⁴⁰ for solving the moment constrained maximum entropy problems. The source code of the EBE scheme is available at Ref. 41.

One of our ultimate goals is to use this method for parameter estimation of molecular modeling, in which a Langevin type of model is usually used. The potential energy typically involves a large set of parameters, e.g., the stiffness constants associated with bond stretching and deformation of bond angles, as well as the coefficients in the van der Waals and screened electro-static interactions. The selection of damping coefficients is also non-trivial; see the review paper⁴² for various perspectives. Compared to the existing methods, the novelty of our approach is that we use the response statistics to formulate the parameter estimation problem, which can reveal parameters that do not appear in the equilibrium density. In addition, the polynomial surrogate model provides an efficient mean to solve the nonlinear least-squares problem.

To achieve this goal, however, there are some remaining challenges. For large parameter space, we suspect that one needs to combine the existing methods to compute the parameters associated with the equilibrium density as one way to reduce the problem (as suggested in Step 1 above) with the proposed method to estimate the remaining parameters. Second, the proposed method requires high quality and possibly large amount of samples for accurate evaluation of the essential statistics. However, since the essential statistics depend on the choice of the observable and external forcing, approximating the high-dimensional integral can be avoided so long as the dimensions of the ranges of these two functions are small. Finally, the underlying model might be subject to modeling error. For example, the model may be derived from a multiscale expansion or just empirically postulated. The formulation of the response statistics in the presence of modeling error will be investigated in separate works.

ACKNOWLEDGMENTS

This research was partially supported by the National Science Foundation (NSF) (Grant No. DMS-1619661). X.L. is also supported by the NSF (Grant No. DMS-1522617). J.H. is also supported by the Office of Naval Research (ONR) (Grant No. N00014-16-1-2888).

APPENDIX A: PROOF OF THEOREM 3.1

In Sec. III B, we presented the convergence result Theorem 3.1 and included the proof for the one-dimensional case. Here, we show the proof for the general $N$ -dimensional case. We will first generalize Lemmas 3.3–3.5 to the $N$ -dimensional case.

Lemma A.1

Let

f \in C^{ℓ} ({[- 1, 1]}^{N})

⁠. For the coefficient of the corresponding least-squares approximation

α_{\vec{k}},

we have

lim_{λ \to + \infty} α_{\vec{k} + λ {\vec{e}}_{i}} {(k_{i} + λ)}^{ℓ} = 0,

where

k_{i}

denotes the

i th

component of

\vec{k}

and

{\vec{e}}_{i}

is the

i th

unit vector.

Proof.

To apply Lemma 3.3, we introduce the single-variable function

g (θ_{i}) := \int_{{[- 1, 1]}^{N - 1}} f (θ) \frac{P_{\vec{k}} (θ)}{P_{k_{i}} (θ_{i})} d μ (θ_{1}) \dots d μ (θ_{i - 1}) d μ (θ_{i + 1}) \dots d μ (θ_{N}) .

It is easy to check that

g (θ_{i}) \in C^{ℓ} ([- 1, 1])

and

\int_{- 1}^{1} g (θ_{i}) P_{k_{i} + λ} (θ) d μ (θ_{i}) = α_{\vec{k} + λ {\vec{e}}_{i}}^{\vec{j}} .

We consider the least-squares approximation of

g (θ_{i})

and apply Lemma 3.3 to the corresponding coefficients. We are able to draw the conclusion.

Lemma A.2

Let $f \in C^{ℓ} ({[- 1, 1]}^{N})$ with $ℓ > \frac{3}{2} N$ and consider its least-squares approximation $f^{M}$ given by (11). Then, we have $∥ f - f^{M} ∥_{\infty} \to 0$ as $M \to + \infty$ ⁠.

Proof.

From the proof of Lemma 3.4, it is enough to show that

{f^{M}}

is a Cauchy sequence with respect to

L^{\infty}

-norm. For

m > n

⁠, we have

\begin{aligned} ∥ f^{m} - f^{n} ∥_{\infty} & \leq \sum_{λ = n + 1}^{m} \sum_{{\vec{k}; ∥ \vec{k} ∥_{\infty} = λ}} | α_{\vec{k}} | ∥ P_{\vec{k}} ∥_{\infty} \\ \leq \sum_{λ = n + 1}^{m} {(λ + \frac{1}{2})}^{\frac{N}{2}} \sum_{{\vec{k}; ∥ \vec{k} ∥_{\infty} = λ}} | α_{\vec{k}} |, \end{aligned}

where we have used the fact that for

∥ \vec{k} ∥_{\infty} = λ

⁠,

∥ P_{\vec{k}} ∥_{\infty} = {∥ \prod_{i = 1}^{N} P_{k_{i}} ∥}_{\infty} = \prod_{i = 1}^{N} {∥ P_{k_{i}} ∥}_{\infty} \leq \prod_{i = 1}^{N} \sqrt{k_{i} + \frac{1}{2}} \leq {(λ + \frac{1}{2})}^{\frac{N}{2}} .

We further notice the multi-index set

{\vec{k}; ∥ \vec{k} ∥_{\infty} = λ}

has total

{(λ + 1)}^{N} - λ^{N} = O (λ^{N - 1})

elements, and by Lemma A.1, we know

| α_{\vec{k}} | = o (λ^{- ℓ}), λ = ∥ \vec{k} ∥_{\infty},

for

λ

large enough. Thus, for large enough

m > n

⁠, we have

∥ f^{m} - f^{n} ∥_{\infty} \leq C \sum_{λ = n + 1}^{m} {(λ + \frac{1}{2})}^{\frac{N}{2}} λ^{N - 1} λ^{- ℓ},

where

C

is a constant that only depends on

N

and

ℓ

⁠. Finally, with

ℓ > \frac{3}{2} N

⁠, we know

{f^{M}}

is indeed a Cauchy sequence w.r.t. to

L^{\infty}

-norm, which completes the proof.

To generalize Lemma 3.5, it is enough to consider the partial derivative case, that is, finding the condition for $∥ f_{x_{i}}^{M} - f_{x_{i}} ∥_{\infty} \to 0$ as $n \to + \infty$ ⁠.

Lemma A.3

For

f \in C^{ℓ} ({[- 1, 1]}^{N})

with

ℓ > \frac{3}{2} N + 2,

we consider its least-squares approximation

f^{M}

⁠, which is given by . Then, we have

lim_{M \to + \infty} {∥ \frac{\partial}{\partial x_{i}} (f^{M} - f) ∥}_{\infty} = 0, i = 1, 2, \dots, N .

Proof.

Similar to the proof of Lemma A.2, we first show that

{f_{x_{i}}^{M}}

is a Cauchy sequence. For

m > n

⁠, we have

\begin{aligned} ∥ f_{x_{i}}^{m} - f_{x_{i}}^{n} ∥_{\infty} & \leq \sum_{λ = n + 1}^{m} \sum_{{\vec{k}; ∥ \vec{k} ∥_{\infty} = λ}} | α_{\vec{k}} | {∥ \frac{\partial}{\partial_{x_{i}}} P_{\vec{k}} ∥}_{\infty} \\ \leq \sum_{λ = n + 1}^{m} \frac{1}{2} λ (λ + 1) {(λ + \frac{1}{2})}^{\frac{N}{2}} \sum_{{\vec{k}; ∥ \vec{k} ∥_{\infty} = λ}} | α_{\vec{k}} |, \end{aligned}

where we have used the inequality . By Lemma A.1, similar to the proof of Lemma A.2, we know

f \in C^{ℓ} ([- 1, 1])

provides

∥ f_{x_{i}}^{m} - f_{x_{i}}^{n} ∥_{\infty} \leq C \sum_{λ = n + 1}^{m} \frac{1}{2} λ (λ + 1) {(λ + \frac{1}{2})}^{\frac{N}{2}} λ^{N - 1} λ^{- ℓ},

where

C

is a constant that only depends on

N

and

ℓ

⁠. Finally, with

ℓ > \frac{3}{2} N + 2

⁠, we are able to conclude that

{f_{x_{i}}^{M}}

is a Cauchy sequence with respect to

L^{\infty}

-norm. With the same trick in the proof of Lemma 3.5, one can show that

f_{x_{i}}^{M}

indeed converges to

f_{x_{i}}

uniformly on

[- 1, 1]

⁠.

Now, we are ready to prove Theorem 3.1.

Proof.

First, by Lemmas A.2 and A.3, the regularity assumption simply provides

lim_{M \to + \infty} ∥ f^{M} - f ∥_{\infty} = 0, lim_{M \to + \infty} ∥ \nabla f^{M} - \nabla f ∥_{\infty} = 0.

To construct a similar function

F^{M}

used in the proof of Proposition 3.7, we introduce the notation

min_{θ \in {[- 1, 1]}^{N}} \sum_{i = 1}^{N} {(f_{i}^{M})}^{2} (θ) =: \sum_{i = 1}^{N} {(e_{i}^{M})}^{2}, e^{M} := {(e_{1}^{M}, e_{2}^{M}, \dots, e_{N}^{M})}^{⊤}

and define

F^{M} := f^{M} - e^{M}

⁠. First, we show the existence of a sequence of solutions,

{θ_{M}^{*}}

⁠, that are close to

θ^{†}

by checking the hypothesis of Lemma 3.6 to

F^{M}

at

θ^{†}

⁠. Since

\nabla f (θ^{†})

is invertible, by continuity, for some

r > 0

small enough, there exists positive constants

γ

such that

∥ \nabla f {(θ)}^{- 1} ∥ \leq \frac{γ}{2}, \forall θ \in B (θ^{†}, r) .

Later, we will specify

r

for the convergence of

θ_{M}^{*} \to θ^{†}

as

M \to \infty

⁠.

We further notice that

\nabla F^{M} = \nabla f^{M}

and

\nabla f^{M}

converges to

\nabla f

uniformly. This implies that there exists

L_{1}

such that

\forall M > L_{1}

∥ \nabla F^{M} {(θ)}^{- 1} ∥ \leq γ, \forall θ \in B (θ^{†}, r) .

(A1)

By the definition of

e^{M}

⁠, we know at

θ = θ^{†}

∥ F^{M} (θ^{†}) ∥ = ∥ f^{M} (θ^{†}) - e^{M} ∥ \leq 2 ∥ f^{M} (θ^{†}) ∥,

which leads to

lim_{M \to + \infty} ∥ F^{M} (θ^{†}) ∥ \leq lim_{M \to + \infty} 2 ∥ f^{M} (θ^{†}) ∥ = 2 ∥ f (θ^{†}) ∥ = 0.

Thus, there exists

L_{2}

such that

\forall M > L_{2}

∥ F^{M} (θ) ∥ < \frac{r}{γ} .

(A2)

Let

M > max {L_{1}, L_{2}} =: L

⁠, and combine conditions and , by Lemma 3.6 to

F^{M}

at

θ^{†}

⁠, we conclude that

F^{M} (θ) = 0

has a solution in

B (θ^{†}, r)

⁠. We denote such solution as

θ_{M}^{*}

\forall M > L

⁠, which is a minimizer of the approximated least-squares problem .

Our next task is to show that

θ_{M}^{*}

converges to

θ^{†}

⁠. Notice that in the multi-dimensional case, we do not have the mean-value theorem. As a remedy, we introduce the following matrix-valued function:

J (θ; θ^{†}) := \int_{0}^{1} \nabla f [θ^{†} + t (θ - θ^{†})] d t .

From the definition, we know that

J (θ^{†}; θ^{†}) = \nabla f (θ^{†})

⁠, which is invertible and

\begin{aligned} f (θ) - f (θ^{†}) & = \int_{0}^{1} \nabla f [θ^{†} + t (θ - θ^{†})] (θ - θ^{†}) d t \\ = J (θ; θ^{†}) (θ - θ^{†}) . \end{aligned}

(A3)

Since

J {(θ; θ^{†})}^{- 1}

exists at

θ = θ^{†}

and it is a continuous function of

θ

close to

θ^{†}

⁠, we are able to select

r > 0

[the same

r

in and ] such that

J {(θ; θ^{†})}^{- 1}

exists and

∥ J {(θ; θ^{†})}^{- 1} ∥ < Δ, \forall θ \in O (θ^{†}, r),

for some positive constant

Δ

⁠. Let

θ = θ_{M}^{*}

in , and we obtain

\begin{aligned} ∥ θ_{M}^{*} - θ^{†} ∥ & \leq ∥ J {(θ_{M}^{*}; θ^{†})}^{- 1} ∥ ∥ f (θ_{M}^{*}) - f (θ^{†}) ∥ < Δ ∥ f (θ_{M}^{*}) ∥ \\ \leq Δ [∥ f (θ_{M}^{*}) - f^{M} (θ_{M}^{*}) ∥ + ∥ f^{M} (θ_{M}^{*}) ∥] \\ \leq Δ (∥ f - f^{M} ∥_{\infty} + ∥ f^{M} (θ^{†}) ∥) . \end{aligned}

Thus,

lim_{M \to + \infty} ∥ θ_{M}^{*} - θ^{†} ∥ \leq lim_{M \to + \infty} Δ [∥ f - f^{M} ∥_{\infty} + ∥ f^{M} (θ^{†}) ∥] = 0.

As for the residual error

f^{M} (θ_{M}^{*})

⁠, we know

lim_{M \to + \infty} ∥ f^{M} (θ_{M}^{*}) ∥ \leq lim_{M \to + \infty} ∥ f^{M} (θ^{†}) ∥ = 0.

This concludes the proof.

APPENDIX B: PROOF OF THEOREM 3.8

The full rank condition of the Jacobian matrix is usually required in the local convergence theorem of Newton-like methods. We are going to discuss this issue over a finite dimensional function space $Γ_{N}^{M}$ defined by

Γ_{N}^{M} := span {θ_{1}^{k_{1}} θ_{2}^{k_{2}} \dots θ_{N}^{k_{N}} | k_{i} = 0, 1, \dots, M} .

Since the zero element in $Γ_{N}^{M}$ is the function that is identically zero, we can introduce the following definition of linear independence.

Definition B.1

{f_{1}, \dots, f_{n}} \subset Γ_{N}^{M}

is said to be a linearly independent set if

\sum_{i = 1}^{n} c_{i} f_{i} \equiv 0 \Rightarrow c_{i} = 0, \forall i = 1, \dots, n .

We define the differential operator $D_{θ_{i}} : Γ_{N}^{M} ⟶ Γ_{N}^{M}$ ⁠, which is as follows:

D_{θ_{i}} f := \frac{\partial}{\partial θ_{i}} f, i = 1, \dots N .

Now, let us deduce some relevant properties of such an operator.

Lemma B.2

The operators ${D_{θ_{i}}}_{i = 1}^{N}$ satisfy

$rank (D_{θ_{i}}) = M {(M + 1)}^{N - 1};$
$\forall i \neq j,$ there exists a permutation $P_{i j}$ satisfying $D_{θ_{j}} = P_{i j}^{⊤} D_{θ_{i}} P_{i j}$ ⁠.

Proof.

We first prove

rank (D_{θ_{1}}) = M {(M + 1)}^{N - 1}

⁠. It is enough to show that the dimension of the null space of

D_{θ_{1}}

is

{(M + 1)}^{N - 1}

⁠, that is,

\dim [\ker (D_{θ_{1}})] = {(M + 1)}^{N - 1}

⁠. With the help of multi-indices, one can see that

{θ^{\vec{k}} | ∥ \vec{k} ∥_{\infty} \leq M}

provides a natural basis for

Γ_{N}^{M}

⁠, and

D_{θ_{1}} (θ^{\vec{k}}) \equiv 0 \Leftrightarrow k_{1} = 0,

which yields

\dim [\ker (D_{θ_{1}})] = {(M + 1)}^{N - 1}

⁠. We further notice that any 2-cycle

(i, j)

in the symmetric group

S_{N}

induces a permutation

P_{i j}

over the basis

P_{i j} : \dots θ_{i}^{k_{i}} \dots θ_{j}^{k_{j}} \dots \to \dots θ_{j}^{k_{j}} \dots θ_{i}^{k_{i}} \dots .

This permutation leads to the identity

D_{θ_{j}} = P_{i j}^{⊤} D_{θ_{i}} P_{i j}

⁠, which implies that

rank (D_{θ_{i}}) = rank (D_{θ_{1}}) = M {(M + 1)}^{N - 1}

for

i = 1, \dots, N

⁠.

The following lemma will also be useful to prove our main result.

Lemma B.3

Let

\sum_{i = 1}^{N} λ_{i} D_{θ_{i}}

be a non-trivial linear combination of

{D_{θ_{i}}}_{i = 1}^{N},

then

rank (\sum_{i = 1}^{N} λ_{i} D_{θ_{i}}) \geq M {(M + 1)}^{N - 1} .

Proof.

Let

Λ = \sum_{i = 1}^{N} λ_{i} D_{θ_{i}}

be a non-trivial linear combination of

{D_{θ_{i}}}_{i = 1}^{N}

⁠, i.e.,

\exists λ_{j} \neq 0

⁠. Since

Λ

is also a linear operator over

Γ_{N}^{M}

⁠, we have

\begin{aligned} rank (Λ) & = \dim [\ker {(Λ^{*})}^{⊥}] = {(M + 1)}^{N} - \dim (\ker Λ^{*}) \\ = {(M + 1)}^{N} - \dim [\ker (Λ)] . \end{aligned}

We notice that the image of

θ^{\vec{k}}

can be decomposed into

\begin{aligned} Λ (θ^{\vec{k}}) & = λ_{j} D_{θ_{j}} (θ_{1}^{k_{1}} \dots θ_{j}^{k_{j}} \dots θ_{N}^{k_{N}}) + \sum_{i \neq j} λ_{i} D_{θ_{i}} (θ_{1}^{k_{1}} \dots θ_{j}^{k_{j}} \dots θ_{N}^{k_{N}}) \\ = λ_{j} k_{j} θ_{1}^{k_{1}} \dots θ_{j}^{k_{j} - 1} \dots θ_{N}^{k_{N}} + θ_{j}^{k_{j}} R (θ_{1}, \dots, θ_{j - 1}, θ_{j + 1}, \dots, θ_{N}), \end{aligned}

where

R := \sum_{i \neq j} λ_{i} D_{θ_{i}} (θ_{1}^{k_{1}} \dots θ_{j - 1}^{k_{j - 1}} θ_{j + 1}^{k_{j + 1}} \dots θ_{N}^{k_{N}})

is a polynomial independent with

θ_{j}

⁠. It is easy to see that if

k_{j} \neq 0

⁠, the image cannot be identically zero. This implies that

\ker (Λ) \subset span {θ^{\vec{k}}, k_{j} = 0}

⁠; thus,

rank (Λ) \geq M {(M + 1)}^{N - 1}

⁠.

To gain intuitions, one can consider the following example with $N = 2$ ⁠. One can select the following ordered basis for $Γ_{2}^{M}$

\begin{aligned} B & := {1, θ_{2}, θ_{2}^{2}, \dots, θ_{2}^{M}, θ_{1}, θ_{1} θ_{2}, \dots,; θ_{1} θ_{2}^{M}, \\ θ_{1}^{2}, \dots, θ_{1}^{M} θ_{2}^{M}} . \end{aligned}

We can derive the matrix representation of $D_{θ_{1}}$ and $D_{θ_{2}}$

D_{θ_{1}} = (\begin{matrix} 0 \\ I & 0 \\ 2 I & ⋱ \\ ⋱ & ⋱ \\ M I & 0 \end{matrix}), D_{θ_{2}} = (\begin{matrix} E \\ E \\ ⋱ \\ ⋱ \\ E \end{matrix}),

where $I \in R^{M + 1 \times M + 1}$ is the identity matrix and $E \in R^{M + 1 \times M + 1}$ is defined as follows:

E := (\begin{matrix} 0 \\ 1 & 0 \\ 2 & ⋱ \\ ⋱ & ⋱ \\ M & 0 \end{matrix}) .

The linear combination $D_{θ_{1}} + λ D_{θ_{2}}$ satisfies

D_{θ_{1}} + λ D_{θ_{2}} = (\begin{matrix} λ E \\ I & λ E \\ 2 I & ⋱ \\ ⋱ & ⋱ \\ M I & λ E \end{matrix})

and rank $(D_{θ_{1}} + λ D_{θ_{2}}) \geq M (M + 1)$ for $\forall λ \in R$ ⁠.

We now turn to our discussion of the full rank condition of the Jacobian matrix,

J (θ) := {(\frac{\partial f_{i}}{\partial θ_{j}})}_{i = 1, \dots, K, j = 1, \dots, N},

(B1)

where $f_{i} \in Γ_{N}^{M}$ and $K \geq N$ ⁠. Let $f$ be the column vector ${[f_{1} (θ), \dots, f_{K} (θ)]}^{⊤}$ ⁠, then $f \in {(Γ_{N}^{M})}^{K}$ ⁠, which is the $K$ -time product space of $Γ_{N}^{M}$ ⁠. We can rewrite the Jacobian matrix in the following way:

J (θ) = (D_{θ_{1}} f, \dots, D_{θ_{N}} f),

where $D_{θ_{i}} f$ is defined by operating $D_{θ_{i}}$ on each component of $f$ ⁠. With the above observation, we are ready to prove Theorem 3.8.

Proof of Theorem 3.8.

By the definition, we know that the column space of

J (θ)

is generated by

{D_{θ_{1}} f, \dots, D_{θ_{N}} f} \subset {(Γ_{N}^{M})}^{K} .

Considering a non-trivial linear combination of those vectors denoted by

\sum_{i = 1}^{N} c_{i} D_{θ_{i}} f

⁠, then we have

\begin{aligned} \sum_{i = 1}^{N} c_{i} D_{θ_{i}} f \equiv 0 & \Leftrightarrow (\sum_{i = 1}^{N} c_{i} D_{θ_{i}}) f \equiv 0 \\ \Leftrightarrow {f_{i}}_{i = 1}^{K} \subset \ker (\sum_{i = 1}^{N} c_{i} D_{θ_{i}}) . \end{aligned}

By Lemma B.3, we know that the dimension of

\ker (\sum_{i = 1}^{N} c_{i} D_{θ_{i}})

is less than or equal to

{(M + 1)}^{N - 1}

⁠, while

\dim [span (F)] > {(M + 1)}^{N - 1}

⁠. Therefore, the column space of

J (θ)

has full rank in

{(Γ_{N}^{M})}^{K}

⁠. Furthermore, this implies that

p (θ) := det [J^{⊤} (θ) J (θ)]

is not a zero polynomial. Since

N (J)

corresponds with the root of

p (θ)

⁠, by standard results in algebra, we know that

N (J)

is a finite set for

N = 1

and a nowhere dense set for

N \geq 2

⁠.

APPENDIX C: PROOF OF THE ERGODICITY OF THE LANGEVIN DYNAMICS (21)

With a change of variables [ $q = a (x - x_{0})$ ⁠, $p = a v$ ], we reduce (21) into

\begin{aligned} {\begin{cases} \dot{q} = p, \\ \dot{p} = - U_{0}^{'} (q) - γ p + \sqrt{2 γ k_{B} T} \dot{W}, \end{cases} U_{0} (q) = ϵ (e^{- 2 q} - 2 e^{- q} + 0.01 q^{2}) . \end{aligned}

(C1)

To verify the ergodicity of (C1), we apply Theorem 3.2 in Ref. 15, which requires the potential $U_{0}$ satisfying the following two conditions:

$U_{0} (q) \geq 0$ for all $q \in R$ ⁠.
There exists an $α > 0$ and $β \in (0, 1)$ such that
$\frac{1}{2} U_{0}^{'} (q) q \geq β U_{0} (q) + γ^{2} \frac{β (2 - β)}{8 (1 - β)} q^{2} - α .$
(C2)

For the first condition, we simply notice that adding a constant term in $U_{0}$ does not affect the dynamic or the equilibrium distribution. Replacing $U_{0}$ in (C1) by

U_{0} (q) = ϵ (e^{- 2 q} - 2 e^{- q} + 0.01 q^{2} + 1),

(C3)

we have

ϵ (e^{- 2 q} - 2 e^{- q} + 0.01 q^{2} + 1) = ϵ [{(e^{- q} - 1)}^{2} + 0.01 q^{2}] \geq 0.

For the second condition, substituting $U_{0}$ in (C2) by (C3), the inequality becomes $\forall q \in R$

\begin{aligned} - q e^{- 2 q} + q e^{- q} - β e^{- 2 q} + 2 β e^{- q} \\ + [0.01 (1 - β) - \frac{γ^{2} β (2 - β)}{8 ϵ (1 - β)}] q^{2} \geq - α + β . \end{aligned}

It is enough to show that

the function $f (q) = - q e^{- 2 q} + q e^{- q} - β e^{- 2 q} + 2 β e^{- q}$ has a lower bound, and
for any given $γ, ϵ > 0$ ⁠, $\exists β \in (0, 1)$ such that
$0.01 (1 - β) - \frac{γ^{2} β (2 - β)}{8 ϵ (1 - β)} > 0.$
(C4)

Since $f (q)$ is continuous, by checking its limits on both sides

\begin{aligned} lim_{q \to + \infty} f (q) & = 0, \\ lim_{q \to - \infty} f (q) & = lim_{q \to - \infty} - (q + β) [{(e^{- q} - \frac{1}{2})}^{2} - \frac{1}{4}] \\ + β e^{- q} = + \infty, \end{aligned}

one can see that $f (q)$ must have a lower bound. Given $ϵ, γ > 0$ ⁠, (C4) can be satisfied by taking $β \in (0, 1 - \sqrt{1 - \frac{0.08}{0.08 + γ^{2} / ϵ}})$ ⁠. Thus, one can conclude that both two conditions can be satisfied and the Langevin dynamics (21) is indeed ergodic under all possible values of the parameters.

REFERENCES

1.

G. A.

Pavliotis

,

Stochastic Processes and Applications

(

Springer

,

2016

).

Google Scholar

2.

D.

Gamerman

and

H.

Lopes

, Markov Chain Monte Carlo: Stochastic Simulation for Bayesian Inference, 2nd ed., Chapman & Hall/CRC Texts in Statistical Science (Taylor & Francis, 2006).

3.

E.

Borgonovo

and

E.

Plischke

,

Eur. J. Oper. Res.

248

,

869

(

2016

).

https://doi.org/10.1016/j.ejor.2015.06.032

Crossref

4.

A. P.

Lyubartsev

and

A.

Laaksonen

,

Phys. Rev. E

52

,

3730

(

1995

).

https://doi.org/10.1103/PhysRevE.52.3730

Crossref

5.

J.

Harlim

,

X.

Li

, and

H.

Zhang

,

J. Stat. Phys.

168

,

146

(

2017

).

https://doi.org/10.1007/s10955-017-1788-9

Crossref

6.

M.

Toda

,

R.

Kubo

, and

N.

Hashitsume

,

Statistical Physics II. Nonequilibrium Statistical Mechanics

(

Springer

,

1983

).

Google Scholar

7.

M.

Hairer

and

A. J.

Majda

,

Nonlinearity

23

,

909

(

2010

).

https://doi.org/10.1088/0951-7715/23/4/008

Crossref

8.

G.

Gottwald

,

J.

Wormell

, and

J.

Wouters

,

Physica D Nonlinear Phenom.

331

,

89

(

2016

).

https://doi.org/10.1016/j.physd.2016.05.010

Crossref

9.

D.

Qi

and

A. J.

Majda

,

J. Atmos. Sci.

73

,

4609

(

2016

).

https://doi.org/10.1175/JAS-D-16-0192.1

Crossref

10.

A.

Majda

and

D.

Qi

,

J. Nonlinear Sci.

26

,

233

(

2016

).

https://doi.org/10.1007/s00332-015-9274-5

Crossref

11.

A.

Majda

and

D.

Qi

,

SIAM Rev.

60

,

491

(

2018

).

https://doi.org/10.1137/16M1104664

Crossref

12.

A.

Majda

,

Introduction to Turbulent Dynamical Systems in Complex Systems

(

Springer

,

2016

).

Google Scholar

Crossref

13.

Y. M.

Marzouk

,

H. N.

Najm

, and

L. A.

Rahn

,

J. Comput. Phys.

224

,

560

(

2007

).

https://doi.org/10.1016/j.jcp.2006.10.010

Crossref

14.

Y.

Marzouk

and

D.

Xiu

,

Commun. Comput. Phys.

6

,

826

(

2009

).

https://doi.org/10.4208/cicp

Crossref

15.

J. C.

Mattingly

,

A. M.

Stuart

, and

D. J.

Higham

,

Stoch. Process. Their Appl.

101

,

185

(

2002

).

https://doi.org/10.1016/S0304-4149(02)00150-3

Crossref

16.

J. R.

Singler

,

Math. Comput. Model.

47

,

422

(

2008

).

https://doi.org/10.1016/j.mcm.2007.02.022

Crossref

17.

C. E.

Leith

,

J. Atmos. Sci.

32

,

2022

(

1975

).

https://doi.org/10.1175/1520-0469(1975)032¡2022:CRAFD¿2.0.CO;2

18.

A.

Gritsun

,

G.

Branstator

, and

A.

Majda

,

J. Atmos. Sci.

65

,

2824

(

2008

).

https://doi.org/10.1175/2007JAS2496.1

Crossref

19.

R.

Abramov

and

A.

Majda

,

J. Nonlinear Sci.

18

,

303

(

2008

).

https://doi.org/10.1007/s00332-007-9011-9

Crossref

20.

R.

Abramov

and

A.

Majda

,

J. Atmos. Sci.

66

,

286

(

2009

).

https://doi.org/10.1175/2008JAS2813.1

Crossref

21.

A.

Majda

and

X.

Wang

,

Commun. Math. Sci.

8

,

187

(

2010

).

https://doi.org/10.4310/CMS.2010.v8.n1.a10

Crossref

22.

R.

Abramov

and

A.

Majda

,

J. Phys. Oceanogr.

42

,

243

(

2011

).

https://doi.org/10.1175/JPO-D-11-052.1

Crossref

23.

L.

Ma

,

X.

Li

, and

C.

Liu

,

J. Chem. Phys.

145

,

204117

(

2016

).

https://doi.org/10.1063/1.4967936

Crossref

PubMed

24.

E.

Isaacson

and

H. B.

Keller

,

Analysis of Numerical Methods

(

Courier Corporation

,

2012

).

Google Scholar

25.

J. M.

Ortega

and

W. C.

Rheinboldt

,

Iterative Solution of Nonlinear Equations in Several Variables

(

SIAM

,

1970

), Vol. 30.

Google Scholar

26.

C. T.

Kelley

,

Iterative Methods for Optimization

(

SIAM

,

1999

).

Google Scholar

Crossref

27.

P.

Sheppard

,

M.

Rathinam

, and

M.

Khammash

,

J. Chem. Phys.

136

,

034115

(

2012

).

https://doi.org/10.1063/1.3677230

Crossref

PubMed

28.

I. M.

Sobol

,

Math. Model. Comput. Exp.

1

,

407

(

1993

).

29.

A.

Telatovich

and

X.

Li

, e-print arXiv:1706.04237 (2017).

30.

A.

Hannachi

and

A.

O’Neill

,

Q. J. R. Meteorol. Soc.

127

,

939

(

2001

).

https://doi.org/10.1002/(ISSN)1477-870X

Crossref

31.

D. F.

Anderson

and

J. C.

Mattingly

,

Commun. Math. Sci.

9

(

1

),

301

–

318

(

2009

).

https://doi.org/10.4310/CMS.2011.v9.n1.a15

Crossref

32.

M.

Rosenblatt

,

Ann. Math. Stat.

27

,

832

(

1956

).

https://doi.org/10.1214/aoms/1177728190

Crossref

33.

E.

Parzen

,

Ann. Math. Stat.

33

,

1065

(

1962

).

https://doi.org/10.1214/aoms/1177704472

Crossref

34.

A.

Smola

,

A.

Gretton

,

L.

Song

, and

B.

Schölkopf

, in International Conference on Algorithmic Learning Theory (Springer, 2007), pp. 13–31.

35.

R. R.

Coifman

and

S.

Lafon

,

Appl. Comput. Harmon. Anal.

21

,

5

(

2006

).

https://doi.org/10.1016/j.acha.2006.04.006

Crossref

36.

T.

Berry

and

J.

Harlim

,

Appl. Comput. Harmon. Anal.

40

,

68

(

2016

).

https://doi.org/10.1016/j.acha.2015.01.001

Crossref

37.

T.

Berry

and

J.

Harlim

,

Mon. Wea. Rev.

145

,

2833

(

2017

).

https://doi.org/10.1175/MWR-D-16-0428.1

Crossref

38.

S.

Jiang

and

J.

Harlim

, e-print arXiv:1804.03272.

39.

E.

Jaynes

,

Phys. Rev.

106

,

620

(

1957

).

https://doi.org/10.1103/PhysRev.106.620

Crossref

40.

W.

Hao

and

J.

Harlim

,

Commun. Appl. Math. Comput. Sci.

13

,

189

(

2018

).

https://doi.org/10.2140/camcos

Crossref

41.

W.

Hao

and

J.

Harlim

, “Supplementary material: MATLAB software for the equation-by-equation method for solving the maximum entropy problem,” see https://github.com/whao2008/EBE (2018).

42.

W.

Noid

,

J. Chem. Phys.

139

,

090901

(

2013

).

https://doi.org/10.1063/1.4818908

Crossref

PubMed

2019

Author(s)

A parameter estimation method using linear response statistics: Numerical scheme

I. INTRODUCTION

II. A REVIEW OF ESSENTIAL STATISTICS

III. AN EFFICIENT ALGORITHM FOR PARAMETER ESTIMATION

A. Polynomial based surrogate model

B. The convergence of the approximate solutions

C. Convergence of the algorithm

IV. EXAMPLE I: THE LANGEVIN MODEL

A. Reduction of the parameter space

B. Parameter estimation approaches

1. Conventional method

2. Sensitivity analysis of the essential statistics

C. Numerical results

1. Low damping case

2. High damping case

V. EXAMPLE II: A GRADIENT SYSTEM WITH A TRIPLE-WELL POTENTIAL

A. The essential statistics and reduction of the parameter space

B. Numerical results

VI. SUMMARY AND FURTHER DISCUSSION

ACKNOWLEDGMENTS

APPENDIX A: PROOF OF THEOREM 3.1

APPENDIX B: PROOF OF THEOREM 3.8

APPENDIX C: PROOF OF THE ERGODICITY OF THE LANGEVIN DYNAMICS (21)

REFERENCES

Citing articles via

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

A parameter estimation method using linear response statistics: Numerical scheme

I. INTRODUCTION

II. A REVIEW OF ESSENTIAL STATISTICS

III. AN EFFICIENT ALGORITHM FOR PARAMETER ESTIMATION

A. Polynomial based surrogate model

B. The convergence of the approximate solutions

C. Convergence of the algorithm

IV. EXAMPLE I: THE LANGEVIN MODEL

A. Reduction of the parameter space

B. Parameter estimation approaches

1. Conventional method

2. Sensitivity analysis of the essential statistics

C. Numerical results

1. Low damping case

2. High damping case

V. EXAMPLE II: A GRADIENT SYSTEM WITH A TRIPLE-WELL POTENTIAL

A. The essential statistics and reduction of the parameter space

B. Numerical results

VI. SUMMARY AND FURTHER DISCUSSION

ACKNOWLEDGMENTS

APPENDIX A: PROOF OF THEOREM 3.1

APPENDIX B: PROOF OF THEOREM 3.8

APPENDIX C: PROOF OF THE ERGODICITY OF THE LANGEVIN DYNAMICS (21)

REFERENCES

Citing articles via

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only