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

OCEANS 2019 MTS/IEEE SEATTLE
Isamu FUJITA
Port and Airport Research Institute
Yokosuka, JAPAN.

Background of this Study

General background

Many cases we want to have some physical model fitted to obaservation data.
PDE model : eg.
\( j\omega h = -h_0\left( \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y}\right) \)
\( j\omega u = -g\frac{\partial h}{\partial x} \)
\( j\omega v = -g\frac{\partial h}{\partial y} \)

To get \( \left(h_i, u_i, v_i \right) \) by fitting the model to data \( h_i (i=1..n)\) without solving PDEs in an explicit form...
Explicit function model:
  • Curve fitting
  • Least Mean Squares
PDE model:
  • more comlicated technique required....
  • PDE-constraint optimization

My specific interest

For acculate oil drift simulation, tidal current is important
  • Each tidal components \(M_2,S_2,O_1,K_1...\) are sinusoidal. \[ \begin{array}{ccc} M_2(x,y,t)=\tilde{M_2}(x,y)e^{j\omega t} \quad & \tilde{M_2} & \in \mathbb{C} \\ \end{array} \]
  • Database available, but rather coarse or sometimes gives only \(h\).
  • Fine \( (u,v) \) data are required.


→ Interpolation using PDE model
  • \((h)_{coarse}\) → \( (h,u,v)_{fine} \)
  • a kind of data assimilation..
  • want to do it as easy as possible....

→ Meshfree PDE-constraint optimization
→ MLS-CLS method

MLS-CLS method formulation

Meshfree PDE-constraint optimization

PDE-constraint optimization:
\[ \boldsymbol{f}_{opt}= \mathop{\text{arg min}}\limits_{\boldsymbol{Lf}=\boldsymbol{l}_r} J(f) \tag{1} \]
finding \( \boldsymbol{f}_{opt} \) which minimizes \( J \) as well as fulfills \( \boldsymbol{Lf}=\boldsymbol{l}_r \).

In this study,
  • MLS : Moving Least Squares
    for making \(\boldsymbol{L}\),\(\boldsymbol{l}_r\) \(J\) in meshfree manner.
  • CLS : Constraint Least Squares
    for solving (1) with gereralized RQ fatorization.


MLS-CLS method makes us :
  • Free from tedious mesh generation task.
  • Able to treat observation data and boundary conditions in the same way.

MLS : PDE to algebraic equation without node connection information (meshfree)

In subdomain \(\Omega \), supposing local basis functions   \( \textcolor{red}{\boldsymbol{g}=\left\{ 1 \ x \ x^2\right\} }\) ,
any function \(f\) can be written as:
\[ f(x)=\boldsymbol{g}\left\{a\right\}=\left\{ 1 \ x \ x^2\right\} \left\{ \begin{array}{l} a_0\\ a_1\\ a_2 \end{array} \right\} \tag{1} \]
Given data \( \left( x_i, f_i \right) \) are expected to be:
\[ \left\{ f_i \right\}= \left\{ \begin{array}{c} f_0\\f_1\\...\\f_{n-1} \end{array} \right\} =\left\{ \begin{array}{ccc} 1 & x_0 & x_0^2 \\ 1 & x_1 & x_1^2 \\ 1 & x_{n-1} & x_{n-1}^2 \\ \end{array} \right\} \left\{ \begin{array}{l} a_0\\ a_1\\ a_2 \end{array} \right\} =\textcolor{red}{G}\left\{a\right\} \tag{2} \]
"Least Squares" can solve it and \(\left\{a\right\}\) is given:
\[ \left\{ a\right\}=\left( G^TG\right)^{-1}G^T \left\{ f_i \right\} = \boldsymbol{G}^{+}\left\{ f_i \right\} \tag{3} \]
\((\partial f) \), the space derivatives of \(f\) can be written with \(\left\{ f_i \right\}\) like:
\[ (\partial f ) \equiv \left\{ \begin{array}{c} f\\f'\\f'' \end{array} \right\}_{x=0} =(\partial\boldsymbol{g})\boldsymbol{G}^{+}\left\{ f_i \right\} \quad \text{where} \quad (\partial\boldsymbol{g})= \left( \begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 2 \\ \end{array} \right) \tag{4} \]
Linear PDE can be written in matrix form. For examle:
\[ f''+ \omega^2 f=0 \tag{5} \]
Example:
\[ \begin{array}{ccl} f''+ \omega^2 f=0 & \stackrel{\text{step 1}}{\Longrightarrow}& \left\{ \omega^2 \ 0 \ 1 \right\} \left\{ \begin{array}{c} f\\f'\\f'' \end{array} \right\} =0 \\ & & \\ & \stackrel{\text{step 2}}{\Longrightarrow} & \left\{ \omega^2 \ 0 \ 1 \right\} (\partial\boldsymbol{g})\boldsymbol{G}^{+}\left\{ f_i \right\}=0 \\ & & \\ & \stackrel{\text{step 3}}{\Longrightarrow}& \left\{ \omega^2 \ 0 \ 1 \right\} \left( \begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 2 \\ \end{array} \right) \boldsymbol{G}^{+} \left\{ f_i \right\}=0 \\ & & \\ & \stackrel{\text{step 4}}{\Longrightarrow} & \left\{ \omega^2 \ 0 \ 2 \right\}\boldsymbol{G}^{+} \left\{ f_i \right\}=0 \end{array} \]


Fig. MLS example
\( f \) and \( \frac{df}{dx} \)

CLS formulation

In \( \Omega \),
  • A PDE converted to algebraic equation by MLS
\( \boldsymbol{l}(\partial f)=l_r \quad \stackrel{\text{MLS}}{\Longrightarrow} \quad \boldsymbol{l}(\partial\boldsymbol{g})\boldsymbol{G}^{+}\left\{ f_i \right\}= l_r \)

In \(\Sigma\Omega\),
  • CLS problem : \[ \left\{f\right\}_{opt}= \mathop{\text{arg min}}\limits_{\boldsymbol{L}\left\{f\right\}=\boldsymbol{l}_r} J \]
    which can be solved using RQ factorization (in LAPACK)
  • Constraint condition :   \( \boldsymbol{L}\left\{f\right\}=\boldsymbol{l}_r \)
  • Evaluation function to minizie :   \( J \)
\[ \boldsymbol{L}\left\{f\right\}=\boldsymbol{l}_r \] \[ \underset{\LARGE{\begin{array}{c} || \\ \boldsymbol{L} \end{array} }}{ \left\{ \begin{array}{l} \boldsymbol{l}^0(\partial\boldsymbol{g})\boldsymbol{G}^{+0} \\ \qquad \ddots \\ \qquad \qquad \boldsymbol{l}^i(\partial\boldsymbol{g})\boldsymbol{G}^{+i}\\ \qquad \qquad \qquad \ddots\\ \qquad \qquad \qquad \qquad \boldsymbol{l}^{n-1}(\partial\boldsymbol{g})\boldsymbol{G}^{+\\n-1} \end{array} \right\} \left\{f\right\} } = \underset{\LARGE{\begin{array}{c} || \\ \boldsymbol{l}_r \end{array} }}{ \left\{ \begin{array}{c} l_r^0\\ \vdots \\ l_r^i\\ \vdots \\ l_r^{n-1}\\ \end{array} \right\} } \]

Example in 2D space

\[ \begin{array}{ccc} \small{\text{PDE}} & & \small{\text{Linear Algebraic Eq.}} \\ \frac{\partial^2 f}{\partial x^2}+ \frac{\partial^2 f}{\partial y^2} = 0 & \Longrightarrow & \quad \boldsymbol{l}(\partial f)=l_r \\ & \Longrightarrow & \quad \boldsymbol{l}(\partial\boldsymbol{g})\boldsymbol{G}^{+}\left\{ f_i \right\}= l_r \end{array} \]
where,
\[ \boldsymbol{l}=\left\{ 0 \ 0 \ 0 \ 1 \ 0 \ 1 \right\} \quad , \quad l_r=0 \] \[ (\partial f)=\left\{ \begin{array}{c} f \\ \frac{\partial f}{\partial x} \\ \frac{\partial f}{\partial y} \\ \frac{\partial^2 f}{\partial x^2} \\ \frac{\partial^2 f}{\partial x \partial y} \\ \frac{\partial^2 f}{\partial y^2} \\ \end{array} \right\} \quad , \quad (\partial \boldsymbol{g})= \left\{ \begin{array}{cccccc|l} 1 & x & y & x^2 & xy & y^2 & g\\ \hline 1 & 0 & 0 & 0 & 0 & 0 & g\\ 0 & 1 & 0 & 0 & 0 & 0 & g_{,x}\\ 0 & 0 & 1 & 0 & 0 & 0 & g_{,y}\\ 0 & 0 & 0 & 2 & 0 & 0 & g_{,xx}\\ 0 & 0 & 0 & 0 & 1 & 0 & g_{,xy}\\ 0 & 0 & 0 & 0 & 0 & 2 & g_{,yy} \end{array} \right\} \]
\[ \begin{array}{l} J &= \displaystyle \sum_{\Omega_{b}}{||\boldsymbol{l_b}f-l_{rb}||}+ \displaystyle \sum_{\Omega_{d}}{||\boldsymbol{l_d}f-l_{rd}||}\\ &\\ &= \left|\left| \left\{ \begin{array}{l} \boldsymbol{l}_b^0(\partial\boldsymbol{g})\boldsymbol{G}_b^{+0} \\ \qquad \ddots \qquad \\ \qquad \qquad \boldsymbol{l}_b^p(\partial\boldsymbol{g})\boldsymbol{G}_b^{+p}\\ \hline \boldsymbol{l}_d^0(\partial\boldsymbol{g})\boldsymbol{G}_d^{+0} \\ \qquad \ddots \\ \qquad \qquad \boldsymbol{l}_d^{q-1}(\partial\boldsymbol{g})\boldsymbol{G}_d^{+\\q-1} \end{array} \right\} \left\{f\right\} - \left\{ \begin{array}{c} l_{rb}^0\\ \vdots \\ l_{rb}^p\\ \hline l_{rd}^0\\ \vdots \\ l_{rd}^{q-1}\\ \end{array} \right\} \right|\right| \end{array} \]
In 2D system, \[ \boldsymbol{g}=\left\{ 1 \ x \ y \ x^2 \ xy \ y^2 \ ... \right\} \] \[ (\partial g) = \left\{ \begin{array}{c} \boldsymbol{g}\\ \frac{\partial \boldsymbol{g}}{\partial x} \\ \frac{\partial \boldsymbol{g}}{\partial y} \\ \frac{\partial^2\boldsymbol{g}}{\partial x^2} \\ \frac{\partial^2\boldsymbol{g}}{\partial x \partial y} \\ \frac{\partial^2\boldsymbol{g}}{\partial y^2} \\ \end{array} \right\}_{x=y=0} = \left( \begin{array}{ccc} 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 2 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 2 \end{array} \right) \] \[ (\partial f)=\left\{ \begin{array}{c} f \\ \frac{\partial f}{\partial x} \\ \frac{\partial f}{\partial y} \\ \frac{\partial^2 f}{\partial x^2} \\ \frac{\partial^2 f}{\partial x \partial y} \\ \frac{\partial^2 f}{\partial y^2} \\ ... \end{array} \right\} \]
\[ \begin{array}{lcl} \Omega &:& \text{Governing equations}\\ \Omega_b &:& \text{Boundary conditions}\\ \Omega_d &:& \text{Observation data}\\ & \Downarrow & \\ \Sigma\Omega &:& \text{Whole domain}\\ \end{array} \]

Components included in \(\boldsymbol{L}\left\{f\right\}=\boldsymbol{l}_r \) and \( J \)

Implementation in C++

Making \( \boldsymbol{L} \) and \(\boldsymbol{l}_r\) \[ \begin{array}{ccccccc} \frac{\partial h}{\partial t} & + h_0( & \frac{\partial u}{\partial x} & + & \frac{\partial v}{\partial y} & )=& 0 \\ \Downarrow & & \Downarrow & & \Downarrow & & \Downarrow \\ (1) & & (2) & & (3) & & (4) \\ \end{array} \]
ρ 
//
// Equation for Mass conservation (H)
//
auto GoverningEquation_H=[&](particle &p,ColumnMajorMatrix<dcomplex> &L,Vector<dcomplex> &L_right,double w){
  for(neighbors::iterator nit=p.neighb.begin(); nit!=p.neighb.end();nit++) {
        *L(L.cnt, Hoff+VNUM*nit->neighbor_it->id )+=ci*omega*nit->f; ←(1)
	*L(L.cnt, Uoff+VNUM*nit->neighbor_it->id )+=p.depth*nit->df_dx; ←(2)
	*L(L.cnt, Voff+VNUM*nit->neighbor_it->id )+=p.depth*nit->df_dy; ←(3)
      }
      *L_right(L.cnt)= 0.0;  ←(4)
      Column_Weight(L,L_right,L.cnt,w);
      L.cnt++;
      return L.cnt;
};

Results and Discussions

Verify MLS-CLS

Compare MLS-CLS to theoretical solution.
  • Governing PDEs: \[ \begin{array}{l} \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} + f_cv\\ \frac{\partial v}{\partial t} = -g\frac{\partial h}{\partial y} - f_cu \end{array} \]
  • Kelvin wave
    in semi-infinite domain \( \left(y>0\right) \) under Coriolis effect \[ \begin{array}{l} \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} \]
  • Nodal arrangement
  • MLS-CLS calculation result
Fig. Kelvin wave (theory)
G.E. points (+) : 480
Neumann points ( □ ) : 38
Dirichlet points (○) : 3

Fitting to data with noise

MLS-CLS gives
  • correct answer when the data have no noise (left).
But,
  • If data have some noise, mostly often the case.., an over-fitting problem appears (right).
  • PDE system contains a lot of modes from long wave to short wave. (Eq.) (Fig.)
  • Short wave mode is sensitive to data noise
  • Raw MLS-CLS tries to minimize \(J\) using short wave responding to data noise too much.
  • Some appropriation (smoother) is required to get smooth solution... → ???
Eg. 1D Heat Conduction Equation:
\[ \frac{\partial \theta}{\partial t}= \lambda \frac{\partial^2 \theta}{\partial x^2} \ \ \ \longrightarrow \ \ \ \theta(t+\Delta t)=T_r \ \theta(t) \] \[ T_r= \left[ \begin{array}{ccc} a & b & & & 0 \\ b & a & b & & \\ & b & a & \ddots & \\ & & \ddots & \ddots & b \\ 0 & & & b & a \end{array} \right] \]
Eigen values and vectors   :   \( T_r \ \theta_j= \lambda_j \theta_j \)

\[ \left\{ \ \ \ \begin{align} \lambda_j&=a+2bcos\left(\frac{j\pi}{n+1}\right) \\ \theta_j&=\left\{ sin\left( \frac{j\pi}{n+1} \right) \ \ldots \ \underbrace{sin\left( \frac{ij\pi}{n+1} \right)}_{i \ th} \ \ldots \ sin\left(\frac{nj\pi}{n+1} \right) \right\}^T \end{align} \right. \]
Eigen vectors containe various wave-length
Example of over-fitting problem
Case 1
(noiseless data)
 Case 2
(noiseful data)
Double clicks to start webGL animation on figures

Laplacian filter to suppress over-fitting

  • Random noise of the data might influence short-wave mode in particular.
  • Extract short-wave components using Laplacian filter.
  • Laplacian term added to \(J\) as a penalty term: \[ J \rightarrow J+\sum w||\nabla^2 f|| \]
  • Example( Wave propagation ).
Laplacian: \( \frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2} \) , \( K=\left[\begin{array}{ccc} 0 & 1 & 0 \\ 1 & -4 & 1 \\ 0 & 1 & 0 \\ \end{array} \right] \)
  • extracts rapid value change,in other words, steep curvature.
  • gives large value for short-wave mode
Case 1
(raw MLS-CLS)
Laplacian 		Case 2
(with Laplacian filter)

Reproduction test of \(M_2\) tidal current

Tidal motion model
\[ \begin{array}{l} \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} + f_cv\\ \frac{\partial v}{\partial t} = -g\frac{\partial h}{\partial y} - f_cu \end{array} \]

Reproduction test of \(M_2\) tidal current

  • Nodal arrangement
  • Simplified physical model
  • Result 1 (raw MLS-CLS)
  • Result 2 (Laplacian filter-applied MLS-CLS)

Fig. Calculation area
G.E. nodes ( +) : 576
Dirichlet nodes (x) : 49
Tidal motion model
\[ \begin{array}{l} \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} + f_cv\\ \frac{\partial v}{\partial t} = -g\frac{\partial h}{\partial y} - f_cu \end{array} \]
Laplacian filter : \( w=1.0 \times 10^{-2} \)
Tide level can be reproduced, but \( \left(u,v\right) \) are not realistic...
Laplacian filter : \( w=1.0 \times 10^2 \)
Both of tide level and \( \left(u,v\right) \) are produced well...

Concluding remarks

  • Meshfree PDE-constraint optimization problem was discussed.
  • MLS-CLS method works as expected and may be a good option to solve PDE-constraint problem.
  • Raw MLS-CLS suffers from an over-fitting problem if the data have noise.
  • Extra addition of a simple Laplacian filter to the evaluation function is effective to suppress undesiable short-wave-error caused by the data noise.
This presentation was powered by Reveal.js.