Convex Nonparametric Formulation for Identification of Gradient Flows
Mohammad Khosravi Roy S. Smith
math.OC cs.SY eess.SP eess.SY

Hydrogen Non-Equilibrium Ionisation Effects in Coronal Mass Ejections
    P. Pagano<ref> A. Bemporad<ref> D.H. Mackay<ref>
          
     
=====================================================================

empty
empty


In this paper, we develop a nonparametric system identification method for the nonlinear gradient-flow dynamics. In these systems, the vector field is the gradient field of a potential energy function. This fundamental fact about the dynamics of system plays the role of a structural prior knowledge as well as a constraint in the proposed identification method. While the nature of the identification problem is an estimation in the space of functions, we derive an equivalent finite dimensional formulation, which is a convex optimization in form of a quadratic program. This gives scalability of the problem and provides the opportunity for utilizing recently developed large-scale optimization solvers. The central idea in the proposed method is representing the energy function as a difference of two convex functions and estimating these convex functions jointly. Based on necessary and sufficient conditions for function convexity, the identification problem is formulated, and then, the existence, uniqueness and smoothness of the solution is addressed. We also illustrate the method numerically for a demonstrative example.	




§ INTRODUCTION

Nonlinear dynamics are ubiquitous in nature and widely used for modeling various phenomena in physics, chemistry, biology and other fields of science and engineering <cit.>. These models are either derived from first principles or by means of fitting and estimation methods. The latter employ techniques in optimization, statistical learning theory and system identification for deriving the model from the available measurement data. Meanwhile, in many cases modeling is beyond fitting a nonlinear dynamics to the observation data. We may additionally need to incorporate specific properties into the model, which are inherent in the nature of system. These properties includes stability, passivity, positivity or other possible characteristics of the system.  For the case of linear dynamics, many of these properties are already taken into account <cit.>.  

For the nonlinear systems, a similar line of research has received extensive attention in the past decade <cit.>. For example,  
identifying stabilizable non-autonomous dynamics is investigated in <cit.>, and learning the dynamics subject to sparsity of the dynamic modes of system is discussed in <cit.>.  
For the purpose of imitation learning, the dynamics modeled in <cit.> based on Gaussian mixtures and hidden Markov models. 
Meanwhile, a similar approach is employed in <cit.> with a global stability guarantee.
In <cit.>,  a convex quadratic potential energy as well as a linear dissipative field are considered with respect to each data point, and  the dynamics are modeled as a functional weighted sum of the gradient forces and the dissipative fields.
Also, in <cit.>, the dynamics are modeled as weakly nonlinear differential equations which have a linear part for capturing the baseline behavior together with more complicated coupling dynamics for considering more complex phenomena.

An interesting class of nonlinear dynamics are gradient-flows, also known as curl-free vector fields.
In physics, these vector fields are  called conservative forces with classical examples being electric and gravitational forces.
The gradient flows are defined as the 
negative of gradient of a potential energy function.
This property of gradient flow can be used as a structural prior knowledge as well as a constraint in the identification problem. 
In <cit.>, a learning method is introduced based on the notion of vector-valued kernels which might be suitable for learning gradient flow of convex energy functions. However, the arguments in <cit.> do not provide concise theoretical guarantees. 

Inspired by recent progresses in shape-constrained and convex regression <cit.>, we introduce a  nonparametric identification method for the gradient-flow dynamics.
The introduced identification problem is originally an estimation problem in the space of functions, i.e., it is a minimization of fitting or prediction error over the hypothesis space of  convex functions. 
Meanwhile, we derive an equivalent finite dimensional convex optimization problem. 
For the sake of more transparent discussion and ease of notation, in Section <ref>, first analyzes the case of convex energy functions. The results are then extended to the cases where the knowledge of the equilibrium is available, the energy functions are strongly convex as well as (strongly) concave, and subsequently the case of general energy functions in Section <ref>.
Finally, in Section <ref>, the method is numerically illustrated on a demonstrative
example.	



§ NOTATIONS AND PRELIMINARIES

The set of natural numbers, the set of non-negative integers, the set of real numbers, n-dimensional Euclidean space and the space of n by m real matrices are denoted by , ,  , ^n, and ^n× m respectively.
The identity matrix and zero vector in the Euclidean space are denoted by  and  respectively.
The set of symmetric positive definite matrices in ^n× n is denoted by _++^n.
For any pair of symmetric matrices X,Y∈^n× n, we write X≽Y if X-Y∈_++^n.
For a set ⊆^n, the convex hull of  is denoted by (). 
The Euclidean norm and the inner product on ^n are respectively denoted by · and ··.
For a function f, ∇ f and ∇^2f are the gradient and Hessian of f respectively.
For a convex function φ:⊆^n→, the subgradient or subderivative of f at x∈ is denoted by ∂ f(x) and defined as the set of vectors ξ∈^n satisfying the inequality 
φ(y)-φ(x)≥ξy-x,  ∀y∈. 
Let  be a set and  be a subset of . The indicator function of , denoted by _, is defined as
_(y) = 0, if y∈ and  _(y) = ∞, otherwise.


§ PROBLEM STATEMENT

Let  be a simply-connected open subset of ^n and φ:→ be an unknown function.
We call φ the potential energy function or simply, energy function.
A conservative vector field corresponding to φ is induced over the space, and the corresponding dynamics are defined as 

    [ ẋ = f(x):=-∇φ(x). ]

Starting from initial condition x_0∈ at time t=0, the vector field generates a trajectory which is denoted here by x(·;x_0).
Consider a set of initial points x_0^1,…,x_0^ and corresponding trajectories.
For any i=1,…,, let trajectory x(·;x_0^i) be sampled at time instants 0≤ t^i_1<t^i_2<⋯<t^i_n_i, where n_i∈. Let x_k^i denote x(t_k^i;x_0^i) for 1≤ k≤ n_i.

The time derivative of x(t;x^i_0) at each sampling time instant can then be estimated by simply utilizing a nonlinear regression method and subsequently obtaining the derivatives numerically or analytically. Various other techniques, e.g. see <cit.> and the references therein are also available in the literature to estimate the derivative of the trajectory. Let these estimations be denoted by y_k^i, for 1≤ k≤ n_i.
Note that y_k^i is approximately equal to f(x_k^i).
Based on these samples and estimations, we can introduce a set of data, denoted by , which contains data pairs (x_k^i,y_k^i).
More precisely,  is defined as {(x_j,y_j) | 1≤ j≤}, where :=∑_1≤ i≤n_i and, for simplicity of notation, the superscripts are dropped. 

Given the set of data , estimate the unknown vector field f in (<ref>).


This problem is a nonlinear system identification where structural prior knowledge is provided in form of (<ref>).


The dynamics in (<ref>) can be extended to the case of differential inclusions <cit.>. More precisely, one may consider ẋ∈±∂φ(x) or ẋ∈-(∂φ_1(x)-∂φ_2(x)), where φ, φ_1 and φ_2 are convex functions and ∂ is the sub-derivative operator.


The dynamics in (<ref>) can be extended to the case of differential inclusions <cit.>. More precisely, one may consider ẋ∈±∂φ(x) or ẋ∈-(∂φ_1(x)-∂φ_2(x)), where φ, φ_1 and φ_2 are convex functions and ∂ is the sub-derivative operator.

	


§ CONVEX ENERGY FUNCTIONS

In this section, we consider the case where the energy function is convex. In the following, we relax the differentiablity assumption of energy function for the sake of generality.

Let Φ denote the set of convex functions defined over ^n.
Considering the data ,
we define the loss function for the estimation problem, denoted by _Φ,, as the sum of squared errors. 
In other words, _Φ,:Φ×^ n→ is a function such that for any given convex function φ∈Φ and vectors ξ_1,…,ξ_, the value of _Φ,(φ,) is defined as 

    _Φ,(φ,) := ∑_i=1^y_i+ξ_i^2 + ∑_i=1^_∂φ(x_i)(ξ_i),

where ∈^ n is a column vector defined as :=
[ ξ_1^    …  ξ_^ ]^.
Note that for any φ∈Φ, we have that ∂φ(x) ∅, for any x, and also, ∂φ(x) = {∇φ(x)}, when φ is differentiable at x.
Accordingly, the estimation problem is naturally defined as

    [        _ξ_1,…,ξ_∈^n, φ∈Φ          _i=1^y_i+ξ_i^2,;                     s.t. ξ_i∈∂φ(x_i),   ∀ i=1,…,. ]

Note that optimization problem (<ref>) is over the set Φ, a cone in the space of functions which is an infinite-dimensional space. We investigate this problem and introduce a tractable approach for obtaining a solution for (<ref>).


 §.§ Towards Finite-Dimensional Formulation

For any convex function φ, the following holds <cit.> 

    φ(y)-φ(x)≥ξy-x,  ∀x,y∈^n, ∀ξ∈∂φ(x).

Motivated by this property of convex functions, we introduce the following optimization problem

    [                _ξ_1,…,ξ_∈^n
    θ_1,…,θ_∈                       _i=1^y_i+ξ_i^2,;                                  s.t. θ_j-θ_i≥ξ_ix_j-x_i,     ∀ i,j∈{1,…,}. ]

Define the vector ∈^ as
:=
[ θ_1   …  θ_ ]^.
Let  be the feasible set in (<ref>),  

    := {(,)∈^×^ n | θ_j-θ_i    ≥ξ_ix_j-x_i, 
      ∀ i,j=1,…,}.

Considering the optimization problem (<ref>), one can define a loss function _,:^×^ n→ as

    _,(,) := ∑_i=1^y_i+ξ_i^2 + _(,).
 
The next theorem presents the connection between optimization problems (<ref>) and (<ref>) as well as the corresponding loss functions. 


(i) Let (φ,) be a solution of (<ref>). Then, (<ref>) has a solution (,) such that 


    _,(,) = _Φ,(φ,).
 
(ii) Conversely, if (,) is a solution of (<ref>), then (<ref>) has a solution (φ,) such that 


    _Φ,(φ,) = _,(,),
 
and ξ_i∈∂φ(x_i), for all i=1,…,. 


Proof of (i): For any i∈{1,…,}, define θ_i=φ(x_i). From (<ref>), one can easily see that (,)∈, and subsequently, _(,) =0. Therefore, we have

    _,(,) =
    _i=1^y_i+ξ_i^2
    = 
    _Φ,(φ,).

Let (<ref>) have a feasible point (,) such that _,(,) <_,(,).
Therefore, we have

    ∑_i=1^y_i+ξ̃_i^2 
    =_,(,) 
    <_,(,)
    =∑_i=1^y_i+ξ_i^2.

Let function φ̃:^n→ be defined as

    [ φ̃(x) : = max_i=1,…,ξ̃_ix-x_i+θ̃_i,   ∀x∈^n. ]

One can easily see that φ̃ is a convex function, i.e., φ̃∈Φ.
Define set-valued map Ĩ:^n⇉{1,…,} as 

    Ĩ(x)={i∈{1,…,} | φ̃(x) = ξ̃_ix-x_i +
    θ̃_i, }.
 
For any x∈^n, we know that
<cit.> 

    ∂φ̃(x)= {ξ̃_i | i∈Ĩ(x)}.

Since (,)∈,  for any i=1,…,, we have

    θ̃_i
    ≥ξ̃_jx_i-x_j +
    θ̃_j,    ∀ j=1,…,.

Therefore, from (<ref>) and (<ref>), one can see

    φ̃(x_i)≥ξ̃_ix_i-x_i +
    θ̃_i= 
    θ̃_i ≥max_j=1,…,ξ̃_jx_i-x_j+θ̃_j = φ̃(x_i).

Subsequently, due to (<ref>), one can see ξ̃_̃ĩ∈∂φ̃(x_i). Therefore, (φ̃,) is a feasible point for (<ref>). 
Subsequently, due to (<ref>) and (<ref>), we have

    _Φ,(φ̃,)= ∑_i=1^y_i + ξ̃_i^2 
    = 
    _,(,)
    <
    _,(,) 
    = _Φ,(φ,),
 
which is a contradiction and we have _,(,) ≤_,(,).
Therefore (,) is a solution of (<ref>) and the proof of part (i) is concluded.

Proof of (ii):  
Let (,) be a solution of (<ref>). Define φ:^n→ as

    φ(x) : = _i=1,…,ξ_ix-x_i+θ_i.

Note that φ is a convex function, i.e. φ∈Φ. 
Define set-valued map I:^n⇉{1,…,} similar to (<ref>). 
For any x∈^n, we have 
∂φ(x)
= 
{ξ_i | i∈ I(x)}.
Since (,)∈, based on a similar argument to the proof of part (ii), we have that ξ_i∈∂φ(x_i), for any i=1,…,. Subsequently, we have 

    _Φ,(φ,) = 
    ∑_i=1^y_i+ξ_i^2
    = _,(,).

Now, let (<ref>) have a feasible point (φ̃,) such that 
_Φ,(φ̃,) <_Φ,(φ,). 
For any i=1,…,, define θ̃_i = φ̃(x_i). Since φ̃ is a convex function, due to (<ref>), one can see that (,)∈. Therefore, we have 

    _,(,) 
    = ∑_i=1^y_i + ξ̃_i^2 
    = _Φ,(φ̃,) 
    <
    _Φ,(φ,) =  _,(,),
 
which is a contradiction. This shows that (φ,) is a solution of (<ref>).
This concludes the proof of part (ii).



Optimization problem (<ref>) admits a solution in form of 

    φ(x) : = _i=1,…,ξ_ix-x_i+θ_i,
 
where (,) is a solution of (<ref>). Moreover, we have


    _Φ,(φ,) = _,(,),
 
and ξ_i∈∂φ(x_i), for all i=1,…,. 


	Optimization problem (<ref>) is a quadratic program. Since  is a non-empty polyhedral cone, (<ref>) has a solution, denoted by (,). Therefore, due to Theorem <ref>, optimization problem (<ref>) admits a solution in form of (<ref>).
	The rest of theorem concludes directly from Theorem <ref> and the given proof.
	
Based on Theorem <ref>, one can solve (<ref>) instead of the main estimation problem (<ref>) and  
introduce an estimation of the energy function as in  (<ref>)
where (,) is a solution of (<ref>).
However, there are two issues to be addressed: optimization problem (<ref>) does not have a unique solution, and it is not smooth. In the following, we will address these issues.







 §.§ Uniqueness by Regularization

One can introduce a regularized version of the optimization problem (<ref>) as

    _(,)∈^×^ n_,(,) + λ (,),

where : ^×^ n→ is the regularization function and λ≥ 0 is the weight of regularization. Based on the next theorem, we introduce a suitable candidate for the regularization function.

There exist a unique ^*∈^ n and a closed and convex set Θ^*⊂^n such that

    :=_(,)∈^×^ n_,(,) = Θ^*×{^*}.



We know that  is a non-empty closed set. Take (_1,_1),(_2,_2) ∈⊂.
Since  is a convex polyhedral cone, one has 

    (,) := 1/2(_1,_1)+1/2(_2,_2)∈.
 
Moreover, since ():= ∑_1≤ i≤y_i+ξ_i^2 is a strongly convex function with ∇^2 =2, we have 

    1/2(_1) +1/2(_2)= () + 1/4_1-_2^2.

Accordingly, we should have _1=_2, otherwise the problem admits a solution with smaller cost. 
From this argument, we know that there exists a set Θ^*⊂^n such that =Θ^*×{^*}. Since  is a closed set, Θ^* is also closed. The convexity of Θ^* follows from the convexity of  and the fact that the cost function does not depend on .

Theorem <ref> says that the potential non-uniqueness of the solution is due to the term . Accordingly, we consider the regularized cost function _λ:^n×^ n→ defined as

    _λ(,) 
    := 
    ∑_i=1^y_i+ξ_i^2 + λ^2,
        ∀(,) ∈^n×^ n.
  
The next theorem characterizes the solution of the corresponding regularized optimization.

For any λ>0, the optimization problem 

    _(,) ∈_λ(,),

has a unique solution, denoted by (,).
Moreover, lim_λ→ 0(,) exists and is equal to (^*,^*) where ^* := _∈Θ^*^2.
Also, lim_λ→∞(,) exists and equals to (,^†) where ^† is the unique solution of  min_(,)∈∑_1≤ i≤y_i+ξ_i^2.


One can easily see that (,)∈ and ∇^2_λ≽λ.

Therefore, (<ref>) is an optimization problem with a strongly convex cost function and non-empty closed and convex feasible set. Therefore, (<ref>) has a unique solution. Similarly, since Θ^* is non-empty, closed and convex,
^* := _∈Θ^*^2 is well-defined and exists uniquely.
From the definition of (^*,^*) and  (,), one can easily see that for any λ>0, we have 

    [   _i=1^y_i+ξ_i^*^2 
    + 
    λ^2                          ≤    _i=1^y_i+ξ_λ,i^2 
    +
    λ^2;                                                     ≤ _i=1^y_i+ξ_i^*^2 
    +
    λ^*^2, ]

and subsequently, it holds that ≤^*.
Similarly, since (,)∈, one can  see that ^2≤ 4 ∑_i=1^y_i^2
and
^*^2≤ 4 ∑_i=1^y_i^2.
Now, define set ⊂^×^ n as

    :={(,) | ≤^*,
    ≤ 2(∑_i=1^y_i^2)^1/2},

which is a compact and convex set.
Define _0 and _λ as

    _0:=_(,)∈∩(,),

and

    _λ:=_(,)∈∩_λ(,),

respectively.
We know that ∩ is a compact set and _λ(,) is a continuous function with respect to (,,λ). Therefore, due to Maximum Theorem <cit.>, we know that the set-valued map λ↦_λ is upper hemicontinuous with non-empty and compact values. Moreover, one has _0={(^*,^*)} and _λ={(,)}. Subsequently, from the upper hemicontinuity of the map λ↦_λ, we have lim_λ→ 0(,) = (^*,^*). 
Replacing λ with λ^-1 and repeating same steps of the proof, one can show the last part of the theorem similarly. 
	
Given λ>0, we can define our estimator as following

    (x) : = _i=1,…,ξ_λ,ix-x_i+θ_λ,i,

where (,) is the unique solution of (<ref>).

In addition to inducing the uniqueness of solution, the regularization improves the numerical stability 
and the robustness  with respect to noise. 


If further regularization is required for improving the performance of the estimation, we can use a Tikhonov regularization <cit.> by defining 
(,) :=(,)^2.
In this case,  is also regularized and pushed towards the origin which might be not desirable.









 §.§ Smoothing the Estimator
 
The smooth version of (<ref>) is defined as

    (x)
    =
    τln(1/∑_1≤ i≤exp(1/τ[ξ_λ,ix-x_i
    +
    θ_λ,i]
    )),

where τ is a positive real scalar <cit.>.

Let the log-sum-exp function ℓ:^n→ be defined as 

    ℓ(x) := ln(∑_i=1^ne^x_i),  ∀x=(x_1,…,x_n)∈^n.
	
This function is an analytical and convex function.
The gradient and Hessian of ℓ are  

    ∇ℓ(x) = 1/^zz,

and

    ∇^2ℓ(x) = 
    1/^z(z) - 
    1/(^z)^2zz^ ,

where z = [e^x_1,…,e^x_n]^. Moreover, the following holds

    [ max_1≤ i≤ n x_i ≤ℓ(x) ≤max_1≤ i≤ n x_i + ln n,   ∀x∈^n. ]


This function is used to define a smooth approximation to  in (<ref>).
More precisely, let Ξ_λ be the matrix defined as 

    Ξ_λ:=
    [ ξ_λ,1 ξ_λ,2     …  ξ_λ, ]∈^n×,

and, for i=1,…,, define η_λ,i as η_λ,i:= θ_i-ξ_λ,ix_i and subsequently, let  be the vector defined as 

    := [ η_λ,1     …  η_λ, ]^.

Subsequently, one can see that

    (x) = τℓ(1/τ(Ξ_λ^x+))
    -
    τln.

The next corollary motivates the use of  as a smooth approximant to .

For any τ>0, the function , defined in (<ref>), is a convex and analytical function. Moreover, we have	

    ∇(x) = Ξ_λz/^z,

and

    ∇^2 (x) = 1/^zΞ_λ(z)Ξ_λ^ - 1/(^z)^2Ξ_λzz^Ξ_λ^,

where

    z = 
    [ exp(1/τ[ξ_λ,1x-x_1
    +
    θ_λ,1]);                            ⋮;    exp(1/τ[ξ_λ,x-x_
    +
    θ_λ,]) ].

Also, we have the following inequality

    (x)-
    τln≤(x) ≤(x),   ∀x∈^n.



By taking τ small enough,  in a uniform approximant of . 
More precisely, let ϵ be an arbitrary positive real scalar and let τ < ϵ/ln. Then (<ref>) shows that |(x)-
(x)|<ϵ, for any x∈^n.








 §.§ Further Extensions

We now introduce  extensions of the proposed estimation strategy to  other settings. 
We briefly present how the main estimation problem (<ref>) and its finite-dimensional version (<ref>) are adapted. The regularization and the smoothing procedures follow the same lines as before.





  §.§.§ Concave Energy Functions
 We know that ψ is a concave function if and only if φ:=-ψ is a convex function. Accordingly, the estimation problem (<ref>) is modified to the following optimization problem

    [        _ξ_1,…,ξ_∈^n, φ∈Φ          _i=1^y_i-ξ_i^2,;                     s.t. ξ_i∈∂φ(x_i),   ∀ i=1,…,. ]
 
Analogously to optimization problem (<ref>), 
this leads to the following finite dimensional problem

    [               _ξ_1,…,ξ_∈^n
    θ_1,…,θ_∈                      _i=1^y_i-ξ_i^2,;                                 s.t. θ_j-θ_i≥ξ_ix_j-x_i,   ∀  i,j∈{1,…,}. ]






  §.§.§ Strongly Convex Energy Functions
 
For μ>0, the convex function φ is said to be μ-strongly convex if for any x,y∈^n, we have

    φ(y)-φ(x)≥ξy-x+μ/2x-y^2,  ∀ξ∈∂φ(x).

Let Φ_μ denote the set of μ-strongly convex functions. If we know that the energy function φ belongs to Φ_μ, then the estimation problem 
(<ref>) is adapted to the following optimization problem

    [      _ξ_1,…,ξ_∈^n, φ∈Φ_μ          _i=1^y_i+ξ_i^2,;                     s.t. ξ_i∈∂φ(x_i),   ∀ i=1,…,. ]
 
In this case, the optimization problem (<ref>) is modified to 
    [                           _ξ_1,…,ξ_∈^n
    θ_1,…,θ_∈                                  _i=1^y_i+ξ_i^2,;                                             s.t. θ_j-θ_i≥ξ_ix_j-x_i +μ/2x_i-x_j^2,  ∀
    i,j∈{1,…,}. ]

In the case of μ-strongly concave energy functions, this can be adapted along the lines of the formulation in (<ref>) .





  §.§.§ Including Knowledge of Equilibrium Points
 
Let assume that we know x_0 = is an equilibrium of the dynamical system. Accordingly, we need to have ∈∂φ(). Therefore, in order to incorporate this knowledge, the estimation problem 
(<ref>) should be modified to the following,

    [       _ξ_1,…,ξ_∈^n, φ∈Φ         _i=1^y_i+ξ_i^2,;                    s.t. ξ_i∈∂φ(x_i),    i=1,…,,;                  ∈∂φ(). ]
 
Without loss of generality, we can assume that φ()=0. 
Accordingly, one can set θ_0=0 and ξ_0=.
Therefore, we modify the optimization problem (<ref>) as following 

    [                 _ξ_0,…,ξ_∈^n
    θ_0,…,θ_∈                        _i=1^y_i+ξ_i^2,;                                   s.t. θ_j-θ_i≥ξ_ix_j-x_i,   ∀  i,j∈{0,1,…,},;                                 θ_0=0,;                                  ξ_0=. ]

	


§ GENERAL ENERGY FUNCTIONS

In this section, we consider general energy functions.
For the sake of generality, the differentiability assumption of the energy function is relaxed initially.
The next theorem plays a key role in the formulation of the estimation problem.
 (i)
Let Ω⊆^n be a  convex set and φ:Ω→ be a C^2(Ω,) function with bounded Hessian, i.e., sup_x∈Ω∇^2φ(x)<∞. Then, there exist convex functions ,:Ω→ such that φ(x) = (x)-(x), for any x∈Ω.
(ii) Moreover, if Ω is convex and compact, then the Hessian is bounded and φ is decomposable as the difference of two convex functions. 


We construct the loss function for the estimation problem, denoted by _Φ,, as the sum of squared errors. 
More precisely, the function _Φ,:Φ×Φ×^ n×^ n→ is defined as 

    _Φ,(,,,) 
    :=
    ∑_i=1^y_i
    +
    _i
    -
    _i^2
    + 
    ∑_i=1^_∂(x_i)(_i) 
    +
    ∑_i=1^_∂(x_i)(_i),

for any pair of convex functions ,∈Φ and vectors _1,…,_,_1,…,_∈^n, 
,∈^ n where ,∈^ n are column vectors respectively defined as :=
[_1^… _^ ]^
and
:=
[_1^… _^ ]^.
Accordingly, the estimation problem is defined as

    [ __1,…,_∈^n, ∈Φ
    _1,…,_∈^n, ∈Φ            _i=1^y_i+_i-_i^2,;                         s.t.     _i ∈∂(x_i),   ∀  i=1,…,,;                                   _i ∈∂(x_i),   ∀ i=1,…,. ]

Analogous to the previous section, we can introduce a finite-dimensional formulation as following

    [                         _,∈^n
    ,∈^ n                   _i=1^y_i+_i-_i^2,;                                s.t. _j-_i ≥_ix_j - x_i,   ∀ i,j∈{1,…,},;                                     _j-_i ≥_ix_j - x_i,   ∀ i,j∈{1,…,}, ]

where ,∈^n are defined respectively as =[_1,…,_]^ and
=[_1,…,_]^.
Considering this optimization problem, we define a loss function _,:^×^×^ n×^ n→  as

    _,(,,,) := _i=1^y_i+_i-_i^2 
    + _(,)+ _(,).
 
With lines of proof similar to those in Section <ref>, we formalize the connection between optimization problems (<ref>) and (<ref>).

(i) Let (,,,) be a solution of (<ref>). Then, (<ref>) has a solution (,,,) such that 

    _,(,,,)
    =
    _Φ,(,,,).

(ii) If (,,,) is a solution of (<ref>), then (<ref>) has a solution (,,,) such that 

    _Φ,(,,,) = _,(,,,),

and _i∈∂(x_i), _i∈∂(x_i), for all i=1,…,.


Optimization problem (<ref>) admits a solution in form of 

    (x) : =  (x) - (x) 
    =
    max_1≤ i_1≤_i_1x-x_i_1 
    +
    _i_1
    -
    max_1≤ i_2≤_i_2x-x_i_2
    +
    _i_2.

where (,,,) is a solution of (<ref>). Moreover, we have


    _Φ,(,,,) = _,(,,,),

and _i∈∂(x_i), _i∈∂(x_i), for all i=1,…,.

As in Section <ref>, regularization can be used for imposing uniqueness in the estimation.
Using the same arguments as those given in the proof of Theorem <ref>, one can obtain a similar conclusion and subsequently  show that the difference - is unique and the potential non-uniqueness of the solution is due to the other terms. Consequently, one can introduce the regularized cost function _λ:^n
×^n×^ n×^ n→ as 

    _λ(,,,) := ∑_i=1^y_i+_i-_i^2  
    +
    λ(^2+^2++^2),
  
and solve the following regularized optimization problem

    _(,)∈,(,)∈_λ(,,,),

where λ> 0 is the regularization weight. Similar to Theorem <ref>, one can show that, for any λ>0, the regularized optimization problem (<ref>) has a unique solution, denoted by  (,,,).
Consequently, we define  and  similar to (<ref>), and thus, 
the estimation is defined as

    (x) : =  (x) - (x) 
    =
    max_1≤ i_1≤_λ,i_1x-x_i_1 
    +
    _λ,i_1
    -
    max_1≤ i_2≤_λ,i_2x-x_i_2
    +
    _λ,i_2.

Similar to the previous section, we can smooth this estimator 
using log-sum-exp function. 
In this regard, let 
 and  be defined as in (<ref>), and then define the smooth estimator, denoted by , as 

    [ (x) = (x)-(x). ]

Note that Corollary <ref> and Corollary <ref> are valid for  and . 
Moreover, we have the following corollary for . 

Let ϵ be an arbitrary positive real scalar and let τ < ϵ/2ln. Then, due to  (<ref>),  we have |(x)-
(x)|<ϵ, for any x∈^n.
In other words, we can uniformly approximate  with  to an arbitrary accuracy by taking τ sufficiently small.


Let {h_k}_k=1^p be a set of given vector fields. 
Then, the proposed method can be extended to the case where the dynamics is sum of a gradient flow and a parametric part, i.e., 

    f(x)=-∇φ(x)+∑_k=1^p α_k h_k(x).
 
To estimate f, we modify optimization problem (<ref>) to 

    [               _∈^ n, ∈Φ
    ∈^ n, ∈Φ
    α_1,…,α_p∈ _i=1^y_i + _i - _i - _k=1^p α_k h_k(x_i)^2,;                                        s.t.                    _i ∈∂(x_i),   ∀  i=1,…,,;                                                                 _i ∈∂(x_i),   ∀ i=1,…,. ]

and apply the previous adaptations to make the problem finite dimensional and sufficiently smooth.
		



§ NUMERICAL EXPERIMENTS

In this section, we discuss a numerical example.
To this end, consider the non-convex energy function φ:^2→ defined as

    φ(x_1,x_2) = 
    ax_1^2+  
    bx_1x_2+  
    ax_2^2
    -cx_1^4-cx_2^4,

where a = 0.7, b= -0.5 and c = 0.15.
The graph of energy function is shown in Figure <ref>.
The gradient flow dynamics of (<ref>) are given by

    [ ẋ_1 = f_1(x) =  -∂_x_1φ(x_1,x_2) = -2ax_1-bx_2+4cx_1^3,; ẋ_2 = f_2(x) =  -∂_x_2φ(x_1,x_2) = -bx_1-2ax_2+4cx_2^3. ]

Since the gradient flows are curl-free, a single  trajectory does not explore the space.
Therefore, in order to collect sufficient data points for identifying the dynamics, it is required to take a  set of sufficiently rich initial points and sample the resulting trajectories.
We consider 18 initial points in 

    Ω:=[-1.9,1.9]×[-1.9,1.9]⊂^2.
 
From the resulting trajectories,  
we take 118 noisy samples with additive white Gaussian noise of zero mean and standard deviation σ_w=0.01.
In Figure <ref>, the initial points, the sampled points and the trajectories are shown by red bullets, black bullets and dotted lines, respectively.
The derivatives are estimated using MATLAB tools for curve-fitting.
From these point, 80% are randomly chosen as the data set for estimating f(x)=-∇φ(x). 
By solving (<ref>), a close to minimum norm estimation is obtained (see Theorem <ref> and Section <ref>). 
Then, the results are used to construct an estimate of the energy function as in (<ref>).
Further more,  using the log-sum-exp function a smooth version of φ  is derived as in (<ref>).
The parameters λ and τ are chosen based on a cross-validation procedure performed using  the remaining 20% of the data points. The results are λ=10^-8 and τ=0.16.
Following this, we estimate ∂_x_1φ(x_1,x_2) and ∂_x_2φ(x_1,x_2) due to the gradient of the smoothed function and Corollary <ref>.
Figure <ref> shows that the results closely fit the true values.
The calculated coefficient of determination for the estimations, also known as R-squared,  equals to 92.4%. 



	
 

§ CONCLUSION

We have introduced a nonparametric system identification method for nonlinear systems with gradient-flow dynamics. The corresponding vector field is the gradient field of a potential energy function. 
This fact is structural prior knowledge as so is used as a constraint in the proposed method. Initially, the identification problem is formulated as a minimization of the fitting or prediction error over the hypothesis space of convex functions. To give a tractable problem, an equivalent formulation is derived as a finite dimensional quadratic program. This formulation is based on two central ideas: representing the energy function as a difference of two convex functions, and a necessary and sufficient condition for convexity.
The final optimization problem approximates the convex functions jointly. The existence, uniqueness and smoothness of the solution is addressed. Finally, a numerical example is presented where a non-convex energy function is considered and based on sample data from multiple trajectories, the corresponding gradient flow is identified.	
	

IEEEtran

