A Short Note on Meshfree PDE-constrained Optimization using Moving Least Squares

In many fields of science and engineering including oceanography and meteorology, there are a lot of situations where we would like to fit some physical model to observation data. When the model is an explicit function such as a linear combination of polynomials, obviously a well-known least squares works. But when the model is given in the form of partial differential equations (PDEs), a bit more sophisticated method, so called PDE-constrained optimization , is required for solving the problem. One of the most common approaches of the PDE-constrained optimization in case of linear PDEs, is solving the following optimazation problem: \begin{equation} \boldsymbol{f}_{opt}= \mathop{\text{arg min}}\limits_{\boldsymbol{Lf}=\boldsymbol{l}_r} J \label{PDE_optimize} \end{equation} where, \(\boldsymbol{Lf}=\boldsymbol{l}_r\) is a set of linearlinzed algebraic equations which is obtained by descretizing constraint PDEs, \(J\) is a positive-semidefinite function to be minimized which usually include observation data and boundary conditions, and \(\boldsymbol{f}_{opt}\) is a solution which fulfills the constraint PDEs as well as minimizes \(J\). In a linear system, \eqref{PDE_optimize} can be solved using QR decomposition technique. It is interesting that \eqref{PDE_optimize} can be used as a usual PDEs solver without any change other than letting \(J\) not include the observation data. So, formulation \eqref{PDE_optimize} is useful, but is still awkward because, in the first place, discretization of PDEs is not an easy work due to bothersome grid generation.

If some gridless method is introduced, \eqref{PDE_optimize} may become more versatile and powerful, and will be seen in wider applications. From this point of view, we conducted a basic study on gridless PDE-constrained optimization, which mainly includes the following four subjects:

  1. Introducing moving least squares to PDE-constrained optimzation formulation for gridless modeling (MLS-CLS).
  2. Implementation of a test program of MLS-CLS and its verification by comparing a theoretical solution of tidal motion.
  3. Consideration on a method to supress overfitting problem when the data have some error.
  4. Reconstruction of \(M_2\) tidal motion by fitting a simplified model to other data base

MLS assumes local basis functions defined in in each sub domain \(\Omega\) and let any function \(f\) be a linear combination of the basis functions (see ). This process provides \(f\) and its space derivatives \(\frac{\partial^{i+j} f}{\partial x^i \partial y^j}\) as a linear combination of basis functions and \(f\) values on nodes inside the local domain without any information about the connectivity among nodes. In this study, we built 2D model taking 3rd order polynomials: \begin{equation} \left\{1 \; x \; y \; x^2 \; xy \; y^2 \; x^3 \;x^2y \; xy^2 \; y^3\right\} \end{equation} as the basis functions. As for the governing equation, necessary number of algebraic linear equations were put into \(\boldsymbol{Lf}=\boldsymbol{l}_r\) by moving subdomain \(\Omega\) over the whole domain \(\sum\Omega\), and boundary conditions \(\boldsymbol{b_b}f=b_b^r\) in \(\Omega_{b}\) and observation data \(\boldsymbol{d}f=d_r\) in \(\Omega_{d}\) are incorporated into \(J\) as \begin{equation} J=\sum_{\Omega_{b}}{||\boldsymbol{b_b}f-b_b^r||}+\sum_{\Omega_{d}}{||\boldsymbol{b_d}f-b_d^r||} \end{equation} A test code was implemented in OOP manner with C++, calling a lapack complex variable routine "zgglse" for solving \eqref{PDE_optimize}.

Firstly, to check the basic function, the test code was applied to a simplified tide model: \begin{equation} \left\{ \begin{array}{rcl} \frac{\partial h}{\partial t} &=& -h_0\left( \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y}\right)\\ \frac{\partial u}{\partial t} &=& -g\frac{\partial h}{\partial x} + fv\\ \frac{\partial v}{\partial t} &=& -g\frac{\partial h}{\partial y} - fu \end{array} \right. \label{PDE_tide} \end{equation} ,and was compared with a theoretical solution of semi-infinite \(y>0\) tidal motion under Coriolis force: \begin{equation} \left\{ \begin{array}{rcl} \tilde{h}(x,y) &=& h_0e^{ -\frac{fy}{\sqrt{gD}} }e^{j\omega \left( \frac{x}{\sqrt{gD}}\right)} \\ \tilde{u}(x,y) &=& {\sqrt{\frac{g}{D}}} h_0e^{ -\frac{fy}{\sqrt{gD}} }e^{j\omega \left( \frac{x}{\sqrt{gD}} \right)} \\ \tilde{v}(x,y) &=& 0 \end{array} \right. \end{equation} where, \(\tilde{h}\) ,\(\tilde{u}\) and \(\tilde{v}\) mean complex amplitude of \(h\),\(u\) and \(v\). It was confirmed that MLS-CLS method provides a solution which agrees well with the theoretical value (see ) .

Secondary, we considered an optimization to noise-rich data. Generally, a PDE system contains a large number of eigen vectors from long wave to short wave. Random error accompanied with the data influences short-wave mode in particular and lead to unexpected solution(see (a)). To suppress this over-fitting, in this study, Laplacian filter is introduced. The Laplacian filter, which is often used for edge detection in image processing, extracts rapid value change or steep curvature. Taking Laplacian into the evaluation function as a penalty term like: \( J \rightarrow J+\sum w||\nabla^2 f|| \) , this additional term suppresses undesired growth of short-wave mode and provides smooth solution without jitter as (b) shows.

Finally, we tried to reconstruct tidal current \(u\) and \(v\) using MLS-CLS method by fitting \eqref{PDE_tide} to \(h\) of \(M_2\) component from NaoTide database. Node arrangement is shown in (a), where at the points "X", \(M_2\) tide constituents were given from NaoTide, and at the points "+", \eqref{PDE_tide} were given. When the Laplacian filter was not applied, unrealistic velocity field was obtained as (b) shows. To the contrary, when the Laplacian filter applied with suitable weight, reasonable velocity field reconstruction was obtained as (c) shows.

In conclusion, we found that MLS-CLS method works as expected and may be a good option to model and solve PDE-constrained problem. Regarding the overfitting problem, it was found that extra addition of a simple Laplacian filter to \(J\) is effective to supress jitter caused by the data noise.

    Schematics of MLS
    Comparison between MLS-CLS and theoretical value

    (a) Without Laplacian filter

    (b) With Laplacian filter
    Effect of Laplacian filter

    (a) node arrangement

    (b) without Laplacian filter

    (c) With Laplacian filter
    \(u\) \(v\) reconstruction test

    (click figures to see animation or enlarge image)


    word count=
    I.FUJITA PARI