GNKI doc source and edits

docs/make.jl

@@ -77,6 +77,7 @@ pages = [
     "Examples" => examples,
     "List of default configurations" => "defaults.md",
     "Ensemble Kalman Inversion" => "ensemble_kalman_inversion.md",
+    "Gauss Newton Kalman Inversion" => "gauss_newton_kalman_inversion.md",
     "Ensemble Kalman Sampler" => "ensemble_kalman_sampler.md",
     "Unscented Kalman Inversion" => "unscented_kalman_inversion.md",
     "Learning rate schedulers" => "learning_rate_scheduler.md",

docs/src/defaults.md

@@ -69,6 +69,15 @@ failure_handler_method = SampleSuccGauss()
 accelerator = DefaultAccelerator()
 ```
 
+## `process <: GaussNewtonInversion` 
+Process documentation [here](@ref gnki)
+```julia
+scheduler = DataMisfitController(terminate_at = 1)
+localization_method = Localizers.SECNice()
+failure_handler_method = SampleSuccGauss()
+accelerator = NesterovAccelerator()
+```
+
 ## `process <: Sampler`
 Process documentation [here](@ref eks)
 ```julia

docs/src/gauss_newton_kalman_inversion.md

@@ -0,0 +1,57 @@
+# [Gauss Newton Kalman Inversion](@id gnki)
+
+### What Is It and What Does It Do?
+Gauss Netwon Kalman Inversion (GNKI) ([Chada et al, 2020](https://arxiv.org/pdf/2010.13299)), also known as the Iterative Ensemble Kalman Filter with Statistical Linearization, is a derivative-free ensemble optimizaton method based on the Gauss Newton optimization update and the Iterative Extended Kalman Filter (IExKF) ([Jazwinski, 2007]).  In the linear case and continuous limit, GNKI recovers the true posterior mean and covariance.  Empirically, GNKI performs well as an optimization algorithm in the nonlinear case.  
+
+### Problem Formulation
+
+The data ``y`` and parameter vector ``\theta`` are assumed to be related according to:
+```math
+\tag{1} y = \mathcal{G}(\theta) + \eta \,,
+```
+where ``\mathcal{G}:  \mathbb{R}^p \rightarrow \mathbb{R}^d`` denotes the forward map, ``y \in \mathbb{R}^d`` is the vector of observations, and ``\eta`` is the observational noise, which is assumed to be drawn from a ``d``-dimensional Gaussian with distribution ``\mathcal{N}(0, \Gamma_y)``. The objective of the inverse problem is to compute the unknown parameters ``\theta`` given the observations ``y``, the known forward map ``\mathcal{G}``, and noise characteristics ``\eta`` of the process.
+
+!!! note
+    GNKI relies on minimizing a loss function that includes regularization.  The user must specify a Gaussian prior distribution. See [Prior distributions](@ref parameter-distributions) to see how one can apply flexible constraints while maintaining Gaussian priors. 
+
+The optimal parameters ``\theta^*`` given relation (1) minimize the loss 
+
+ ```math
+\mathcal{L}(\theta, y) = \langle \mathcal{G}(\theta) - y \, , \, \Gamma_y^{-1} \left ( \mathcal{G}(\theta) - y \right ) \rangle + \langle m - \theta \, , \, \Gamma_{\theta}^{-1} \left ( m - \theta  \right ) \rangle,
+```
+
+where ``m`` is the prior mean and ``\Gamma_{\theta}`` is the prior covariance. 
+
+### Algorithm
+
+GNKI updates the ``j``-th ensemble member at the ``n``-th iteration by directly approximating the Jacobian with statistics from the ensemble.
+
+First, the ensemble covariance matrices are computed: 
+```math
+\begin{aligned}
+        &\mathcal{G}_n^{(j)}  = \mathcal{G}(\theta_n^{(j)}) \qquad 
+        \bar{\mathcal{G}}_n = \dfrac{1}{J}\sum_{k=1}^J\mathcal{G}_n^{(k)} \\
+        & C^{\theta \mathcal{G}}_n = \dfrac{1}{J - 1}\sum_{k=1}^{J}
+        (\theta_n^{(k)} - \bar{\theta}_n )(\mathcal{G}_n^{(k)} - \bar{\mathcal{G}}_n)^T \\
+        & C^{\theta \theta}_n = \dfrac{1}{J - 1} \sum_{k=1}^{J} 
+        (\theta_n^{(k)} - \bar{\theta}_n )(\theta_n^{(k)} - \bar{\theta}_n )^T.
+
+\end{aligned}
+```
+
+Using the ensemble covariance matrices, the update equation from ``n`` to ``n+1`` under GNKI is
+```math
+\begin{aligned}
+        & K_n = \Gamma_{\theta} G_n^T \left(G_n \Gamma_{\theta} G_n^T + \Gamma_{y}\right)^{-1} , \qquad G_n = \left(C^{\theta \mathcal{G}}_n\right)^T \left(C^{\theta \theta}_n\right)^{-1} \\
+
+        & \theta_{n+1}^{(j)} = \theta_n^{(j)} + \alpha \left\{ K_n\left(y_n^{(j)} - \mathcal{G}(\theta_n^{(j)})\right) + \left(I - K_n G_n\right)\left(m_n^{(j)} - \theta_n^{(j)}\right) \right\},  
+\end{aligned}
+```
+
+where ``y_n^{(j)} \sim \mathcal{N}(y, 2\alpha^{-1}\Gamma_y)`` and ``m_n^{(j)} \sim \mathcal{N}(m, 2\alpha^{-1}\Gamma_{\theta})``.
+
+## Creating the EKI Object
+
+An ensemble Kalman inversion object can be created using the `EnsembleKalmanProcess` constructor by specifying the ` GaussNewtonInversion()` process type.
+
+

docs/src/index.md

@@ -5,6 +5,7 @@
 Currently, the following methods are implemented in the library:
  - Ensemble Kalman Inversion (EKI) - The traditional optimization technique based on the (perturbed-observation-based) Ensemble Kalman Filter EnKF ([Iglesias, Law, Stuart, 2013](http://dx.doi.org/10.1088/0266-5611/29/4/045001)),
  - Ensemble Transform Kalman Inversion (ETKI) - An optimization technique based on the (square-root-based) ensemble transform Kalman filter  ([Bishop et al., 2001](http://doi.org/10.1175/1520-0493(2001)129<0420:ASWTET>2.0.CO;2), [Huang et al., 2022](http://doi.org/10.1088/1361-6420/ac99fa) )
+ - Gauss Newton Kalman Inversion (GNKI) [a.k.a. Iterative Ensemble Kalman Filter with Satistical Linearization] - (description here)
  - Ensemble Kalman Sampler (EKS) - also obtains a Gaussian Approximation of the posterior distribution, through a Monte Carlo integration ([Garbuno-Inigo, Hoffmann, Li, Stuart, 2020](https://doi.org/10.1137/19M1251655)),
  - Unscented Kalman Inversion (UKI) - also obtains a Gaussian Approximation of the posterior distribution, through a quadrature based integration approach ([Huang, Schneider, Stuart, 2022](https://doi.org/10.1016/j.jcp.2022.111262)),
 - Sparsity-inducing Ensemble Kalman Inversion (SEKI) - Additionally adds approximate ``L^0`` and ``L^1`` penalization to the EKI ([Schneider, Stuart, Wu, 2020](https://doi.org/10.48550/arXiv.2007.06175)).