Collocation Methods

Introduction

Many methods for solving inverse problems involving ODE's need costly ODE solvers at each iteration to compute the loss. E.g. for L2 loss:

$\mathrm{Loss}_{\mathrm{L2}} = \left\|u - \hat{u}_{\mathrm{solver}}(u_0, t,p) \right\|_2$

Where $u$ is the observed data, $\hat{u}_{\mathrm{solver}}(u_0, t,p)$ is the ODE solver prediction of $u$, given initial conditions $u_0$, timespan $t$ and parameters $p$.

In order to reduce the computational cost, a collocation method may be used ^[3], ^[1]. Collocation methods circumvents the need for an ODE solver by estimating the derivatives of the data. The loss is computed between these estimated derivatives and the system dynamics. E.g.

$\mathrm{Loss}_{\mathrm{L2}} = \left\|\hat{u}' - f(u, t,p) \right\|_2$

Where $\hat{u}'$ is the estimated derivatives of the observed data.

JuliaSimModelOptimizer implements three different collocation methods to estimate the derivative of the data: Kernel Collocation, Spline Collocation and Noise Robust Collocation.

Kernel Collocation

KernelCollocate(; maxiters=nothing, maxtime=nothing,  optimizer = Optim.LBFGS(linesearch = BackTracking()), kernel = EpanechnikovKernel())

Theory

Assume that each snapshot $y_i$ of the time series $\mathbf{y}$ can be approximated by a polynomial of degree $p$ in the following way:

\[y_i \approx \sum_{j=0}^{p} \alpha_j (t_i - t)^{j}\]

Where $\alpha_j$ are polynomial coefficients, $t_i$ is the time at snapshot $y_i$, in the neighborhood of $t$.

Following Eq. 2.9 in ^[2], we get an estimate of $y_i'$ by

\[\hat{y}_i' = \hat{\alpha_1}\]

Where $\hat{\alpha_1}$ is the first entry in $\hat{\alpha}$, which is the vector that minimizes the weighted least squares problem:

\[\min{\alpha}\left(\left(\mathbf{y} - \mathbf{T}\alpha\right)^{\top}\mathbf{W}\left(\mathbf{y} - \mathbf{T}\alpha\right)\right)\]

Where

\[\mathbf{T}=\left[\begin{array}{cccc} 1, & t_1-t_i & \ldots & \left(t_1-t_i\right)^p \\ \vdots & \vdots & & \vdots \\ 1, & t_n-t_i & \cdots & \left(t_n-t_i\right)^p \end{array}\right]\]

and $W=\operatorname{diag}\left(K_h\left(\mathbf{t}-t_i\right)\right)$ where $K_h$ is the kernel function with bandwith $h$ and $\mathbf{t}$ is the vector of time stamps.

To calculate $\hat{\alpha}$, we use the standard weighted least squares solution:

\[\hat{\alpha} = \left(\mathbf{T}^{\top}\mathbf{W}\mathbf{T}\right)^{-1} \mathbf{T}^{\top}\mathbf{W}\mathbf{y}\]

The implementation in JuliaSimModelOptimizer uses a polynomial of degree 1 to smoothen the data and a polynomial of degree 2 to estimate the derivatives of the data.

For reference see ^[3]

For more in depth reference see Chapter 2.1-2.3 and 3.1-3.8 in ^[2]

Available Kernels

The table below shows the available kernels in JuliaSimModelOptimizer

Spline Collocation

SplineCollocate(; maxiters=nothing, maxtime=nothing, optimizer = Optim.LBFGS(linesearch = BackTracking()), interp = Datainterpolations.CubicSpline)

Noise Robust Collocation

NoiseRobustCollocate(; maxiters=nothing, maxtime=nothing, diff_iters, dx, α, optimizer = Optim.LBFGS(linesearch = BackTracking()))