\section{Introduction}
The ensemble Kalman filtering (EnKF) method is a widely used data
assimilation method for high-dimensional state-space problems with
nonlinear dynamics. Owing to its simple implementation and efficiency,
ensemble-based filtering methods have rapidly gained popularity in geophysical sciences with applications, for example, in weather
forecasting~\cite{kalnay2003atmospheric}, atmosphere-ocean/lake
simulations~\cite{houtekamer2005atmospheric, baracchini2020data,
  hammoud2021moving}, and oil reservoir
management~\cite{aanonsen2009ensemble, ruchi2021fast}. The EnKF
method was originally proposed by Evensen~\cite{evensen1994sequential}. Subsequently, several variants were developed~\cite{houtekamer1998data, bishop2001adaptive,
  anderson2001ensemble}. EnKF approximates the filtering distribution
using the empirical measure of its ensemble members. The
$L^p-$convergence of the EnKF method with perturbed
observations~\cite{houtekamer1998data} has been studied in the literature~\cite{mandel2011convergence, le2009large}.

A challenge in numerical ensemble-based data assimilation is
the significant increase in the numerical simulation cost with an increase in
temporal resolution. This challenge can be overcome using the Monte Carlo
method (MLMC)~\cite{giles2008}, which achieves a substantial reduction in variance by simulating pairwise-coupled realizations on a
hierarchy of temporal discretization levels. MLMC is an extremely flexible
method that has been combined with other methods and successfully implemented in various fields; for instance, the quasi-Monte Carlo~\cite{giles2009multilevel,
  kuo2017, robbe2019}, sequential Monte
Carlo~\cite{beskos2017multilevel,beskos2018multilevel,moral2017multilevel,latz2018multilevel},
inverse problems and experimental design~\cite{beck2020, schaden2020,
  litvinenko2019computation, goda2020, taverniers2020}, differential equations with
randomness~\cite{jourdain2019, detommaso2019, barth2018, khodadadian2020, beck2020Hp, badwaik2019}, limit
theorems~\cite{alay2020, hoel2019central}, importance
sampling~\cite{hammouda2020,kebaier2018coupling, fang2020},
and machine learning~\cite{lye2019}.

The multilevel EnKF (MLEnKF) method was introduced by Hoel et al.~\cite{hoel2016} for
stochastic differential equation models with discrete-time
observations, and an alternative version based on a sample average
of independent pairwise-coupled EnKF estimators was subsequently developed~\cite{hoel2020multilevel}. The main difference between these two
versions of MLEnKF is that the first version~\cite{hoel2016} uses one universal Kalman
gain to update the ensemble members on all hierarchy levels, while the
second version~\cite{hoel2020multilevel} employs one Kalman gain per independent EnKF sample in
the full MLEnKF estimator. The latter approach introduces less
correlation between all particle members of the MLEnKF estimator,
which particularly simplifies convergence analysis and paves the way for extending MLEnKF to the multi-index EnKF (MIEnKF) method introduced in this work. The MLEnKF method was extended to spatiotemporal
(infinite-dimensional state space) models~\cite{chernov2020}.
Similar multilevel techniques were combined with other
ensemble-based filtering methods; examples of such combinations are particle
filtering~\cite{jasra2017multilevel, ballesio2020}, transform particle
filtering~\cite{gregory2016multilevel, gregory2017seamless},
multigrid~\cite{moldovan2020} and the
recent extension to continuous-time (Kalman--Bucy)
filtering~\cite{chada2020multilevel}.

The successful implementation of MLMC depends on a strong pairwise
coupling of realizations on neighboring hierarchy levels.  Achieving
such a coupling for stochastic differential equations with
sufficiently smooth coefficients is relatively straightforward (once the idea is introduced). However, in more realistic problems with
low-regularity features, strong coupling can be extremely
challenging, if at all possible. An
assessment of MLEnKF in reservoir history matching~\cite{fossum2020} and for applications using
sampling resolution constraints (the so-called multi-fidelity methods)~\cite{gao2020bi,popov2020multifidelity, popov2021multifidelity} are described in the literature.

Another important method is the multi-index Monte Carlo method
(MIMC)~\cite{abdo2016}, which forms the basis of this work. MIMC consists of a multi-index hierarchy of
coupled realizations on neighboring resolutions, and can be regarded as an extension of MLMC. Many concepts related to particle-wise coupling in the proposed MIEnKF
method are common in the MIMC
method for McKean-Vlasov dynamics~\cite{abdo2018}.

The main contribution of this work is to develop the MIEnKF method as
an extension of the MLEnKF by treating both the numerical discretization and
EnKF ensemble size as degrees of freedom. The
extension introduces four coupled of EnKF estimators (i.e., a
coupling in both degrees of freedom) yielding a stronger variance
reduction than that achieved by pairwise-coupling in MLEnKF. Under certain
assumptions, this approach achieves efficiency gains (Table~\ref{table:1}) for computing the weak approximations of quantities of interest (QoI).
\begin{table}[h!]
  \centering
  \begin{tabular}{|c |c |c |c|}
    \hline
    Methods & EnKF & MLEnKF & MIEnKF \\ [0.5ex]
    \hline\hline
    Mean-squared error & $\mathcal{O}(\epsilon^2)$& $\mathcal{O}(\epsilon^2)$ &  $\mathcal{O}(\epsilon^2)$\\
    \hline
    Computational cost & $\mathcal{O}(\epsilon^{-3})$& $\mathcal{O}(\epsilon^{-2}\abs{\log(\epsilon)}^3)$ & $\mathcal{O}(\epsilon^{-2})$ \\
    \hline
  \end{tabular}
  \bigskip
  \caption{Comparison of computational costs versus errors for ensemble Kalman filtering (EnKF), multilevel EnKF (MLEnKF) and multi-index EnKF (MIEnKF) methods, cf. Section~\ref{sec:complexity}.}
  \label{table:1}
\end{table}

The rest of this work is organized as follows. In
Section~\ref{sec:problem}, the setting and
notation for filtering problem are introduced and a brief overview of the EnKF, mean-field EnKF
(MFEnKF), and MLEnKF methods is presented.  In Section~\ref{sec:mienkf}, the
MIEnKF algorithm is provided. The main theorem of MIEnKF complexity is presented in
Section~\ref{sec:complexity}.  In Section~\ref{sec:numerics}, a
numerical comparison of filtering algorithms is provided. The concluding
remarks are presented in Section~\ref{sec:conclusion}.

\section{Problem setting}

In this section, the problem setting and basic notation
commonly used in data assimilation are introduced and the EnKF, MFEnKF
and MLEnKF methods are briefly reviewed.

\label{sec:problem}
Let $\prt{\Omega, \mathcal{F}, \mathbb{P}; \{\mathcal{F}_t\}_{t\geq0}}$ be a complete
probability space equipped with a filtration $\{\mathcal{F}_t\}_{t\geq0}$ of
sub-$\sigma$-algebras of $\mathcal{F}=\mathcal{F}_{\infty}$. We denote by
$L_t^p(\Omega, \mathbb{R}^k)$ the space of $\mathcal{F}_t\backslash\mathcal{B}^k$-measurable
functions\footnote{The function $u$ is
  $\mathcal{F}_t\backslash\mathcal{B}^k$-measurable iff $u^{-1}(B)\in \mathcal{F}_t$ for all
  $B\in \mathcal{B}^k,$ where $\mathcal{B}^k$ is the Borel $\sigma-$algebra on
  $\mathbb{R}^k$.}  $u:\Omega \rightarrow \mathbb{R}^k$ with
$\Ex{\abs{u}^p}<\infty$. Given the initial value
$u_0 \in \cap_{p\geq2} L_0^p(\Omega, \mathbb{R}^d)$, we consider the
discrete-time filtering problem for a system of stochastic dynamics
defined by a sequence of random maps
$\Psi_n: \mathbb{R}^d \times \Omega \rightarrow \mathbb{R}^d$ and observations with
additive noise:
\begin{equation*}
  \centering
  \left\{ \begin{split}
    & u_{n+1}(\omega)=\Psi_n (u_n, \omega),\\
    & y_{n+1}(\omega)= Hu_{n+1}(\omega)+\eta_{n+1},   \; \;
  \end{split}\right.
\end{equation*}
where $\omega\in \Omega$, $n\in \mathbb{N}_0$, $H\in \mathbb{R}^{m\times d}$ is an
observation operator, $\{\eta_k\}_{k\in\mathbb{N}}$ is an independent and
identically distributed (i.i.d) sequence with
$\eta_1 \sim N(0,\Gamma)$. Hereinafter, we do not indicate the dependence on $\omega$ for ease of notation, whenever there is no confusion. Let
$Y_n:=(y_1, y_2, ..., y_n)$ denote the accumulated observation data up
to time $n$ with the convention $Y_0:=\emptyset$. The main objective
of a filtering method is to track the underlying signal $u_n$ for given
$Y_n$ by computing the (filtering) distribution of the conditional
random variable $u_n|Y_n$. The Bayesian interpretation of the
filtering problem can be expressed formally using the following iterative steps
for computing the filtering density:
\[
\begin{split}
  \mbox{ \textbf{Prediction}  } \quad \rho_{u_n|Y_{n-1}}(u) &\propto \int_{\mathbb{R}^d}\rho_{u_n|u_{n-1}}(u) \rho_{u_{n-1}|Y_{n-1}}(v)dv\\
  \mbox{ \textbf{Update}  }\quad \quad \rho_{u_n|Y_n}(u) &\propto  \exp\big(-\big| \Gamma^{-1/2}(y_n-Hu)\big|^2/2 \big) \; \rho_{u_n|Y_{n-1}}(u),
  \end{split}
\]

In the linear Gaussian setting, the exact filter distribution is
provided in a closed form using the Kalman filter, however, when $\Psi$ is
nonlinear, the filter distribution must be approximated.  Hence, we
construct an efficient MIEnKF method
that converges weakly to the mean-field EnKF in the
large-ensemble limit. In other words, for a given QoI
$\varphi: \mathbb{R}^d \rightarrow \mathbb{R}$, our method approximates
\[
\mathbb{E}^{\bar \mu_n} [\varphi(u)]=\int_{\mathbb{R}^d} \varphi(u) \bar \mu_n (du),
\]
where $\bar \mu_n$ denotes the mean-field EnKF measure at time $n$,
cf.~Section \ref{ssec:mfenkf}.

\begin{notation}
  ~\\
  \vspace{-0.4cm}

  \begin{itemize}

  \item For $f,g:(0,\infty) \to [0,\infty)$ the notation $f \lesssim  g$ implies that there exists a
    $C>0$ such that
    \[
    f(x) \le  C g(x), \quad  \forall x \in (0,\infty).
    \]

  \item The notation $f \eqsim  g$ implies that $f \lesssim g$ and $g\lesssim f$.

  \item For $d \in \mathbb{N}$,
    $|x|$ denotes the Euclidean norm of a vector $x \in \mathbb{R}^d$.
    For $\mathcal{F}\backslash  \mathcal{B}^d$-measurable functions $u:\Omega \to \mathbb{R}^d$ and $p\ge 1$,
    \[
    \norm{u}_{p} := \norm{u}_{L^p(\Omega, \mathbb{R}^d)} = \prt{\int_{\Omega} |u(\omega)|^p \, \mathbb{P}(d\omega)}^{1/p}.
    \]

  \item  $\lceil x \rceil := \min \{ z \in \mathbb{Z} \mid z \ge x\}$.
  \end{itemize}
\end{notation}

We denote by $\Psi_n^N$ the numerical discretization of the dynamics
$\Psi_n$ using $N\geq 1$ uniform timesteps, and we consider that the following assumptions are satisfied:
\begin{assumption}\label{ass:Psi}
  Let $u,v \in \cap_{p\ge 2} L^p_n(\Omega, \mathbb{R}^d)$ for any $n\in
  \mathbb{N}_0$, $p\ge 2$. Given that $\varphi$ is Lipschitz, there exists a
  constant $c_p>0$ such that for all $N \geq 1$:
  \begin{itemize}
  \item[(i)] $
    \norm{\Psi^{N}_n(u)}_p \leq c_p (1+\norm{u}_p),
    $

  \item[(ii)]
    $
    \norm{\Psi^{N}_n(u)- \Psi^{N}_n(v)}_p < c_p \norm{u-v}_p .
    $

  \end{itemize}
\end{assumption}
\begin{remark}
Assumption~\ref{ass:Psi} is mainly introduced so that it ensures the
existence of the mean-field EnKF measure $\bar \mu_n$, cf.~\cite[Appendix A]{hoel2020multilevel}
and Section~\ref{ssec:mfenkf}.
\end{remark}

\subsection{EnKF}
\label{ssec:enkf}
The EnKF method is an ensemble-based extension of the Kalman filter to
nonlinear settings. For an EnKF ensemble of size $P$, let
$v_{n,i}:=v_n(\omega_i)$ and $\hat{v}_{n,i}:=\hat{v}_n(\omega_i)$,
respectively, denote the $i-$th particle of the prediction and updated
ensemble at time $n$. Then, the EnKF algorithm with perturbed
observations and numerical dynamics $\Psi^N$ comprises the following
steps:
\begin{equation*}
  \centering \mbox{ \textbf{Prediction}  }
  \left\{\begin{split}
  v_{n+1,i} &= \Psi_n^N(\hat{v}_{n,i}), \quad i=1,2,...,P,\\
  m_{n+1} &= \frac{1}{P} \sum_{i=1}^P v_{n+1,i},\\
  C_{n+1} &=\frac{1}{P-1}\sum_{i=1}^P \prt{v_{n+1,i}-m_{n+1}} \prt{v_{n+1,i}-m_{n+1}}^\mathbf{T}.
  \end{split}\right.
\end{equation*}

\begin{equation*}
  \centering \mbox{ \textbf{Update}  }
  \left\{\begin{split}
  \tilde{y}_{n+1,i}&=y_{n+1}+\eta_{n+1,i}, \quad i=1,2,...,P,\\
  K_{n+1}&=C_{n+1}H^{\mathbf{T}}(HC_{n+1}H^{\mathbf{T}}+\Gamma)^{-1},\\
  \hat{v}_{n+1,i} &= (I-K_{n+1}H)v_{n+1,i}+K_{n+1}\tilde{y}_{n+1,i},\\
  \end{split}\right.
\end{equation*}
where $\eta_{n+1,i}$ are i.i.d.~draws from $N(0, \Gamma)$.

The updated EnKF empirical measure is defined by
\[
\mu_n^{N,P}(dv)=\frac{1}{P} \sum_{i=1}^{P} \delta(dv; \hat{v}_{n, i}),
\]
where $\delta$ is the Dirac measure centered at $\hat{v}_{n, i}$, and
the expectation of a QoI $\varphi: \mathbb{R}^d \rightarrow \mathbb{R}$ with respect
to the EnKF empirical measure is expressed as
\begin{equation}\label{enkf:measure}
  \mu_n^{N,P}[\varphi]=\frac{1}{P}\sum_{i=1}^{P} \varphi(\hat{v}_{n, i}).
\end{equation}
Note that $\mu_n^{N,P}[\varphi]$ is a random variable that depends on
parameters $N$ and $P$.
Under sufficient regularity, $\mu_n^{N,P}[\varphi] \rightarrow \bar \mu_n[\varphi]$
as $N,P\rightarrow \infty$, where $\bar \mu_n[\varphi]$
denotes the known as the expectation of $\varphi$ with respect to the
mean-field EnKF measure (Section~\ref{ssec:mfenkf}).

\subsection{MFEnKF}
\label{ssec:mfenkf}
The MFEnKF is the large-ensemble
and fine-discretization limit of EnKF. In the large-ensemble limit,
the Kalman gain becomes a deterministic matrix. Consequently, one may
view MFEnKF as an ensemble of i.i.d.~noninteracting particles, so
that it suffices to represent the resulting filtering distribution by
one particle. Let $\bar v_n$ and $\hat{\bar{v}}_n$ denote the prediction and
updated state of a mean-field particle at time $n$, respectively.  The
following algorithm defines the MFEnKF for fully non-Gaussian models:
\begin{equation*}
  \centering \mbox{ \textbf{Prediction}  }
  \left\{\begin{split}
  \bar v_{n+1} &= \Psi_n(\hat{\bar{v}}_{n})\\
  \bar{m}_{n+1} &= \Ex{\hat{\bar{v}}_{n+1}},\\
  \bar{C}_{n+1} &= \Ex{\prt{\hat{\bar{v}}_{n+1}-\bar{m}_{n+1}} \prt{\hat{\bar{v}}_{n+1}-\bar{m}_{n+1}}^\mathbf{T}}.
  \end{split}\right.
\end{equation*}

\begin{equation*}
  \centering \mbox{ \textbf{Update}  }
  \left\{\begin{split}
  \tilde{y}_{n+1}&=y_{n+1}+\tilde{\eta}_{n+1}, \\
  \bar{K}_{n+1}&=\bar{C}_{n+1}H^{\mathbf{T}}(H\bar{C}_{n+1}H^{\mathbf{T}}+\Gamma)^{-1},\\
  \hat{\bar{v}}_{n+1,i} &= (I-\bar{K}_{n+1}H)\hat{\bar{v}}_{n+1,i}+\bar{K}_{n+1}\tilde{y}_{n+1,i},\\
  \end{split}\right.
\end{equation*}
where $\tilde{\eta}_{n+1}$ is i.i.d.~draws from $N(0, \Gamma)$.

The expectation of a QoI $\varphi:\mathbb{R}^d \rightarrow \mathbb{R}$ with respect to the updated mean-field EnKF measure is given by
\[
\bar\mu_n[\varphi]=\int_{\mathbb{R}^d}\varphi(v) \bar\mu_n(dv).
\]

\subsection{MLEnKF} \label{ssec:mlenkf} The recently developed MLEnKF
algorithm~\cite{hoel2020multilevel} is a natural stepping stone en route
from EnKF to MIEnKF, for gradually adding layers to the EnKF notation.
MLEnKF is a filtering method based on a sample average
of independent and pairwise-coupled samples of EnKF estimators at
different resolution levels. Let $L\in \mathbb{N}$ denote the finest
resolution level of the estimator, let the sequences
\[
N_\ell = N_0\times 2^\ell
\quad \text{and} \quad P_\ell = P_0 \times 2^\ell \quad \text{with} \quad N_0,P_0 \in \mathbb{N}
\]
denote the numerical resolution and ensemble size on level $\ell$, respectively.
Then, the updated MLEnKF estimator at time $n$ assumes the following form:
\begin{equation}\label{eq:updMLEnKF}
\mu_n^{ML}[\varphi] = \sum_{\ell=0}^L \sum_{m=1}^{M_\ell} \frac{\Big( \mu_n^{N_\ell, P_\ell,m} - (\mu_n^{N_{\ell-1}, P_{\ell-1},1,m} + \mu_n^{N_{\ell-1}, P_{\ell-1},2,m})/2 \Big)[\varphi]}{M_\ell}
\end{equation}
where $\{M_\ell\}_{\ell=0}^L \subset \mathbb{N}$ is a decreasing sequence with $M_\ell$ representing the number of i.i.d.~and pairwise-coupled EnKF estimators
\[
\Big\{ \mu_n^{N_\ell, P_\ell,m}[\varphi] , (\mu_n^{N_{\ell-1}, P_{\ell-1},1,m}[\varphi], \mu_n^{N_{\ell-1}, P_{\ell-1},2,m}[\varphi]) \Big\}_{m=1}^{M_\ell}
\]
on level $\ell$, and $\mu_n^{N_{-1}, P_{-1},m}[\varphi] := 0$.

\subsubsection*{Pairwise coupling of EnKF estimators} For $\ell\ge 1$,
let
\[
\hat{v}_{n,i}^\ell:=\hat{v}_n^\ell(\omega_i)  \qquad i =1, \ldots, P_\ell
\]
denote an updated ensemble at time $n$ coupled
to the following two coarser-level updated ensembles:
\[
\hat{v}_{n,i}^{\ell-1,1}:= \hat{v}_n^{\ell-1,1}(\omega_i)  \qquad i =1, \ldots, P_{\ell-1}
\]
and
\[
\hat{v}_{n,i}^{\ell-1,2}:= \hat{v}_n^{\ell-1,2}(\omega_{P_{\ell-1}+i})  \qquad i =1, \ldots, P_{\ell-1}.
\]
The coupling is imposed under the particle-wise shared initial condition:
\[
\hat{v}_{0,i}^\ell = \begin{cases} \hat{v}_{0,i}^{\ell-1,1} & \text{if} \quad  i \in\{ 1, \ldots, P_{\ell-1}\}\\
  \hat{v}_{0,i-P_{\ell-1}}^{\ell-1,2} & \text{if} \quad i \in\{ P_{\ell-1}+ 1, \ldots, P_{\ell}\},
  \end{cases}
\]
and using the following prediction and update steps, where the driving noise
and the perturbed observations are also shared particle wise:
\begin{equation*}
  \centering \mbox{ \textbf{Prediction}  }
  \left\{\begin{split}
  v_{n+1,i}^{\ell-1,1} &= \Psi_n^{N_{\ell-1}}(\hat{v}_{n,i}^{\ell-1,1}),  \qquad i=1,\ldots,P_{\ell-1},\\
  v_{n+1,i}^{\ell-1,2} &= \Psi_n^{N_{\ell-1}}(\hat{v}_{n,i}^{\ell-1,2}),  \qquad i=1,\ldots,P_{\ell-1},\\
  v_{n+1,i}^{\ell} &= \Psi_n^{N_\ell}(\hat{v}_{n,i}^{\ell}),\qquad \qquad   i=1,\ldots ,P_\ell, \\
  C_{n+1}^{\ell-1,1} &=\overline{\mathrm{Cov}}[v_{n+1, 1:P_{\ell-1}}^{\ell-1,1}],\\
  C_{n+1}^{\ell-1,2} &=\overline{\mathrm{Cov}}[v_{n+1,1:P_{\ell-1}}^{\ell-1,2}],\\
  C_{n+1}^{\ell} &=\overline{\mathrm{Cov}}[v_{n+1,1:P_\ell}^{\ell}],
  \end{split}\right.
\end{equation*}

\begin{equation*}
  \begin{split}
    \overline{\mathrm{Cov}}[v_{n,1:P_\ell}^{\ell}]&:=\sum_{i=1}^{P_\ell}\frac{v_{n,i}^{\ell}(v_{n,i}^{\ell})^\mathbf{T}}{P_\ell}-\sum_{i=1}^{P_\ell}\frac{v_{n,i}^\ell}{P_\ell}\prt{\sum_{i=1}^{P_\ell}\frac{v_{n,i}^{\ell}}{P_\ell}}^\mathbf{T},\\
    \overline{\mathrm{Cov}}[v_{n,1:P_{\ell-1}}^{\ell-1,k}]&:=\sum_{i=1}^{P_{\ell-1}}\frac{v_{n,i}^{\ell-1,k}(v_{n,i}^{\ell-1,k})^\mathbf{T}}{P_{\ell-1}}-\sum_{i=1}^{P_{\ell-1}}\frac{v_{n,i}^{\ell-1,k}}{P_{\ell-1}}\prt{\sum_{i=1}^{P_{\ell-1}}\frac{v_{n,i}^{\ell-1,k}}{P_{\ell-1}}}^\mathbf{T}, \quad k=1,2.
  \end{split}
\end{equation*}
\begin{equation*}
  \centering \mbox{ \textbf{Update} }
  \left\{\begin{split}
  \tilde{y}_{n+1,i}^\ell&=y_{n+1}+\eta_{n+1,i}^\ell, \qquad \qquad \qquad \qquad \qquad \qquad \quad  i=1,\ldots,P_\ell, \\
  K_{n+1}^{\ell-1,1}&=C_{n+1}^{\ell-1,1}H^{\mathbf{T}}(HC_{n+1}^{\ell-1,1} H^{\mathbf{T}}+\Gamma)^{-1},\\
  K_{n+1}^{\ell-1,2}&=C_{n+1}^{\ell-1,2} H^{\mathbf{T}}(HC_{n+1}^{\ell-1,2} H^{\mathbf{T}}+\Gamma)^{-1},\\
  K_{n+1}^\ell&=C_{n+1}^\ell H^{\mathbf{T}}(HC_{n+1}^\ell H^{\mathbf{T}}+\Gamma)^{-1},\\
  \hat{v}_{n+1,i}^{\ell-1,1} &= (I-K_{n+1}^{\ell-1,1}H)v_{n+1,i}^{\ell-1,1}+K_{n+1}^{\ell-1,1}\tilde{y}_{n+1,i}^\ell,\quad  \qquad i=1,\ldots,P_{\ell-1},\\
  \hat{v}_{n+1,i}^{\ell-1,2} &= (I-K_{n+1}^{\ell-1,2}H)v_{n+1,i}^{\ell-1,2}+K_{n+1}^{\ell-1,2}\tilde{y}_{n+1,P_{\ell-1}+i}^\ell, \;\; i=1,\ldots,P_{\ell-1},\\
  \hat{v}_{n+1,i}^{\ell} &= (I-K_{n+1}^{\ell}H)v_{n+1,i}^{\ell}+K_{n+1}^{\ell}\tilde{y}_{n+1,i}^\ell,\qquad \quad \quad  i=1,\ldots,P_\ell,
  \end{split}\right.
\end{equation*}
where $\{\eta_{n+1,i}^\ell\}_{i=1}^{P_\ell}$ is a sequence of independent $N(0, \Gamma)$-distributed random variables.
In the above notation, the coupling between the fine-level EnKF estimator
\[
\mu_n^{N_\ell, P_\ell}[\varphi] := \sum_{i=1}^{P_\ell} \frac{\varphi(\hat{v}_{n,i}^{\ell})}{P_\ell},
\]
and the two coarse-level estimators
\[
\mu_n^{N_{\ell-1}, P_{\ell-1},k}[\varphi] := \sum_{i=1}^{P_{\ell-1}} \frac{\varphi(\hat{v}_{n,i}^{\ell-1,k})}{P_{\ell-1}}, \qquad k=1,2
\]
is obtained through particle-wise coupling
\[
\hat{v}_{n,i}^{\ell} \xleftrightarrow{\text{coupling}} \begin{cases} \hat{v}_{n,i}^{\ell-1,1} & \text{if} \quad  i \in\{ 1, \ldots, P_{\ell-1}\}, \\
    \hat{v}_{n,i-P_{\ell-1}}^{\ell-1,2} & \text{if} \quad i \in\{ P_{\ell-1}+ 1, \ldots, P_{\ell}\}.
  \end{cases}
\]

\section{MIEnKF}
\label{sec:mienkf}
This section describes the MIEnKF method by extending
the MLEnKF method in Section~\ref{ssec:mlenkf}.

To define a set of multi-index discretizations for the MIEnKF, we
first introduce $\ensuremath{\boldsymbol\ell}:=(\ell_1, \ell_2)\in \mathbb{N}_0^2$ and
shorthands $\ensuremath{\boldsymbol e}_1:=(1,0)$, $\ensuremath{\boldsymbol e}_2:=(0,1)$, and
$\ensuremath{\boldsymbol 1}:=(1,1)$. Similarly for MLEnKF, we associate sequences of
natural numbers $N_{\ell_1} = N_0\times2^{\ell_1}$ and $P_{\ell_2} = P_0\times
2^{\ell_2}$ to the number of timesteps and particles on the 2-index ``level''
$\ensuremath{\boldsymbol\ell}$.

We seek to approximate the mean-field measure $\bar{\mu}_n[\varphi]$
and we denote the discrete approximation corresponding to the 2-index
$\ensuremath{\boldsymbol\ell}$ by $\mu_n^{\ensuremath{\boldsymbol\ell}}[\varphi]:=\mu_n^{N_{\ell_1},
  P_{\ell_2}}[\varphi]$. In other words, $\mu_n^{\ensuremath{\boldsymbol\ell}}[\varphi]$ is the
EnKF estimator~\eqref{enkf:measure} computed using $N_{\ell_1}$
timesteps and $P_{\ell_2}$ particles. We define
first-order difference operators for numbers of timesteps
and particles as follows:
\begin{equation}\label{eq:differenceOperators}
  \begin{split}
    \Delta_1\mu_n^{\ensuremath{\boldsymbol\ell}}[\varphi]&=
    \begin{cases}
      \prt{\mu_n^{\ensuremath{\boldsymbol\ell}}-\mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_1}}[\varphi], \qquad \qquad \qquad \quad &\mbox{  if } \ell_1>0, \\
      \mu_n^{\ensuremath{\boldsymbol\ell}}[\varphi], &\mbox{  if } \ell_1=0
    \end{cases}\\
    \Delta_2\mu_n^{\ensuremath{\boldsymbol\ell}}[\varphi]&=
    \begin{cases}
      \prt{\mu_n^{\ensuremath{\boldsymbol\ell}}-\prt{\mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2,1}+\mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2,2}}/2}[\varphi], &\mbox{  if } \ell_2>0, \\
      \mu_n^{\ensuremath{\boldsymbol\ell}}[\varphi], &\mbox{  if } \ell_2=0
    \end{cases}
  \end{split}
\end{equation}
where $\mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2,1}[\varphi]$ and
$\mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2,2}[\varphi]$ are two i.i.d.~copies of
$\mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2}[\varphi]$. Note that $\mu_n^{\ensuremath{\boldsymbol\ell}}[\varphi]$
comprises $P_{\ell_2}$ ensemble members, whereas
$\mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2}[\varphi]$ consists of $P_{\ell_2-1}$, i.e., half as many ensemble members.
Therefore, the pair $(\mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2,1}[\varphi],\mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2,2}[\varphi])$
are introduced to uniquely couple each ensemble member from
$\mu_n^{\ensuremath{\boldsymbol\ell}}[\varphi]$ with one from the ``level''
$\ensuremath{\boldsymbol\ell} - \ensuremath{\boldsymbol e}_2$.

We define the four-coupled EnKF estimator using the first-order mixed difference as follows:
\begin{equation}\label{eq:mixDiffOper}
  \begin{split}
    \ensuremath{\boldsymbol \Delta} \mu_n^{\ensuremath{\boldsymbol\ell}}[\varphi]&:=\Delta_1(\Delta_2 \mu_n^{\ensuremath{\boldsymbol\ell}}[\varphi])=\Delta_2(\Delta_1 \mu_n^{\ensuremath{\boldsymbol\ell}}[\varphi])=\Delta_2(\mu_n^{\ensuremath{\boldsymbol\ell}}-\mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_1})[\varphi]\\
    &=\Bigg(\mu_n^{\ensuremath{\boldsymbol\ell}} - \Big(\mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2,1} + \mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2,2}\Big)/2\\
    &\qquad - \mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_1} + \Big(\mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1},1} + \mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1},2}\Big)/2 \Bigg)[\varphi],

  \end{split}
\end{equation}
where the pair $(\mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1},1}[\varphi],
\mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1},2}[\varphi])$ of i.i.d.~copies of
$\mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1}}[\varphi]$  is also introduced to uniquely couple
$\mu_n^{\ensuremath{\boldsymbol\ell}}[\varphi]$ with one from the ``level'' $\ensuremath{\boldsymbol\ell} -\ensuremath{\boldsymbol 1}$.
Under the assumption that $\Ex{\mu_n^{\ensuremath{\boldsymbol\ell}}[\varphi]} \rightarrow
\bar{\mu}_n[\varphi]$ as $|\ensuremath{\boldsymbol\ell}| \rightarrow \infty$,
the linearity of the expectation operator
implies that
\begin{equation}\label{eq:telescsum}
  \bar{\mu}_n[\varphi]=\sum_{ \ensuremath{\boldsymbol\ell}  \in \mathbb{N}_0^2} \Ex{\ensuremath{\boldsymbol \Delta} \mu_n^{\ensuremath{\boldsymbol\ell}}[\varphi]} =\sum_{ \ensuremath{\boldsymbol\ell}  \in \mathcal{I}} \Ex{\ensuremath{\boldsymbol \Delta} \mu_n^{\ensuremath{\boldsymbol\ell}}[\varphi]}+\sum_{ \ensuremath{\boldsymbol\ell}  \notin \mathcal{I}} \Ex{\ensuremath{\boldsymbol \Delta} \mu_n^{\ensuremath{\boldsymbol\ell}}[\varphi]},
\end{equation}
for any index set $\mathcal{I}\subset \mathbb{N}_0^2$.

For a properly selected index set $\mathcal{I}$, which we will specify later (Section~\ref{sec:complexity}),
the MIEnKF estimator is defined as the sample-average estimator of the
first term on the right-hand side of~\eqref{eq:telescsum}:
\begin{equation} \label{eq:MIestimator}
    \mu^{MI}_n[\varphi] := \sum_{ \ensuremath{\boldsymbol\ell}  \in \mathcal{I}}
    \sum_{m=1}^{M_{\ensuremath{\boldsymbol\ell}}} \frac{\ensuremath{\boldsymbol \Delta} \mu_n^{\ensuremath{\boldsymbol\ell},m}[\varphi]}{M_{\ensuremath{\boldsymbol\ell}}},
\end{equation}
where $\{\ensuremath{\boldsymbol \Delta} \mu_n^{\ensuremath{\boldsymbol\ell},m}[\varphi]\}_{m=1}^{M_{\ensuremath{\boldsymbol\ell}}}$ are i.i.d. copies of $\ensuremath{\boldsymbol \Delta} \mu_n^{\ensuremath{\boldsymbol\ell},m}[\varphi]$,
and
$\{ \ensuremath{\boldsymbol \Delta} \mu_n^{\ensuremath{\boldsymbol\ell},m}[\varphi]\}_{(\ensuremath{\boldsymbol\ell},m)}$ are mutually independent.

The primary motivation for sampling four-coupled EnKF estimators in
the MIEnKF estimator is that the coupling may lead to a substantial
variance reduction, thus improving the efficiency of the sampling method.
Similar to the setting of multilevel Monte Carlo estimators,
the tractability of~\eqref{eq:MIestimator} is optimized through
careful selection of the index set $\mathcal{I}$ and the number of samples
$M_{\ensuremath{\boldsymbol\ell}}$. This can be achieved by solving a constrained
optimization problem using the knapsack algorithm~\cite{abdo2016,abdo2018}.

\subsubsection*{Four-coupled EnKF estimators}
To describe the coupling between the EnKF estimators
$(\mu_n^{\ensuremath{\boldsymbol\ell}}, \mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2},
\mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_1},\mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1}})[\varphi]$, we
introduce the four-coupled updated-state ensembles
$\{(\hat{v}_{n}^{\ensuremath{\boldsymbol\ell}},
\hat{v}_{n}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2},\hat{v}_{n}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_1},\hat{v}_{n}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1}})_i\}_{i=1}^{P_{\ell_2}}$.
The second ensemble of particles is a union of two EnKF ensembles:
\[
\hat{v}_{n, i}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2,1} := \hat{v}_{n, i}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2}  \qquad i = 1, \ldots P_{\ell_2-1},
\]
and
\[
\hat{v}_{n, i}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2,2} := \hat{v}_{n, P_{\ell_2-1} + i}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2}  \qquad i = 1, \ldots P_{\ell_2-1},
\]
and the fourth ensemble is also a union of two EnKF ensembles:
\[
\hat{v}_{n, i}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1},1} := \hat{v}_{n, i}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1}}  \qquad i = 1, \ldots P_{\ell_2-1},
\]
and
\[
\hat{v}_{n, i}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1},2} := \hat{v}_{n, P_{\ell_2-1} + i}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1}}  \qquad i = 1, \ldots P_{\ell_2-1}.
\]
Further details regarding the nature of the particle-wise four-coupling of these ensembles is
provided in Section~\ref{ssec:mienkfalg}.

The empirical estimator $\mu_n^{\ensuremath{\boldsymbol\ell}} [\varphi]$ is induced by the
ensemble $\hat{v}_{n, 1:P_{\ell_2}}^{\ensuremath{\boldsymbol\ell}} :=
\{\hat{v}_{n,i}^{\ensuremath{\boldsymbol\ell}}\}_{i=1}^{P_{\ell_2}}$, meaning that it equals the
sample average of
$\{\varphi(\hat{v}_{n,i}^{\ensuremath{\boldsymbol\ell}})\}_{i=1}^{P_{\ell_2}}$.  Similarly,
$\mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_1}[\varphi]$ is induced by $\hat{v}_{n,
  1:P_{\ell_2}}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_1}:= \{\hat{v}_{n,
  i}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_1}\}_{i=1}^{P_{\ell_2}}$,
\[
\mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2}[\varphi]:=\frac{(\mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2,1} +
  \mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2,2})[\varphi]}{2}
\]
is induced by the union of two
ensembles
\[\hat{v}_{n, 1:P_{\ell_2}}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2}:=\{\hat{v}_{n,
  k}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2,1}\}_{k=1}^{P_{\ell_2-1}}\cup \{\hat{v}_{n, k}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2,2}\}_{k=1}^{P_{\ell_2-1}},
\]
and
\[
\mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1}}[\varphi]:=\frac{(\mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1},1} + \mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1},2})[\varphi]}{2}
\]
is induced by $\hat{v}_{n, 1:P_{\ell_2}}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1}}:=\{\hat{v}_{n,
  k}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1},1}\}_{k=1}^{P_{\ell_2-1}}\cup \{\hat{v}_{n,
  k}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1},2}\}_{k=1}^{P_{\ell_2-1}}$.
For consistency with~\eqref{eq:differenceOperators}, we impose the condition that
$\mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_1}[\varphi]=\mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1}}[\varphi]=0$ when
$\ell_1=0$, and $\mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2}[\varphi]=0$ when $\ell_2=0$. Figure~\ref{fig:mienkfDiag} shows a visual description
of all the couplings of the MIEnKF estimator.

Then, the MIEnKF estimator~\eqref{eq:MIestimator} can also be written as
\begin{equation}\label{eq:MIestimshort}
  \begin{split}
    \mu^{MI}_n[\varphi] := \sum_{ \ensuremath{\boldsymbol\ell} \in \mathcal{I}}
    \sum_{m=1}^{M_{\ensuremath{\boldsymbol\ell}}} \frac{\big(\mu_n^{\ensuremath{\boldsymbol\ell},m}-\mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_1,m}-\mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2,m}+\mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1},m}\big)[\varphi]}{M_{\ensuremath{\boldsymbol\ell}}}

  \end{split}
\end{equation}
where $\{(\mu_n^{\ensuremath{\boldsymbol\ell},m},\mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_1,m},\mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2,m},\mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1},m})[\varphi]\}_m$
are independent copies of the estimators $(\mu_n^{\ensuremath{\boldsymbol\ell}},\mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_1},\mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2},\mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1}})[\varphi]$
and
$\{(\mu_n^{\ensuremath{\boldsymbol\ell},m},\mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_1,m},\mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2,m},\mu_n^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1},m})[\varphi]\}_{(\ensuremath{\boldsymbol\ell}, m)}$ are mutually independent.

\begin{remark}
  For comparison, one can also represent the MLEnKF estimator~\eqref{eq:updMLEnKF} using the above 2-index notation
  \[
  \mu_n^{ML}[\varphi]=\sum_{\ell=0}^{L} \sum_{m=1}^{M_{\ell}} \frac{(\mu_n^{(\ell, \ell),m}-(\mu_n^{(\ell,\ell) - \ensuremath{\boldsymbol 1},m})[\varphi]}{M_{\ell}}.
  \]

\end{remark}
\subsection{Particle-wise four-coupling for MIEnKF}\label{ssec:mienkfalg}
We now describe how the four-coupling of EnKF estimators manifests itself
on the particle level. At time $n=0$, the fine-index update
ensemble $\{\hat{v}^{\ensuremath{\boldsymbol\ell} }_{0,i}\}_{i=1}^{P_{\ell_2}}$ comprises
independent $\mathbb{P}_{u_0|Y_0}$-distributed samples that are
particle-wise coupled to three other ensembles as
\[
\hat{v}^{\ensuremath{\boldsymbol\ell}}_{0,i}=\hat{v}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_1 }_{0,i}=\hat{v}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2
}_{0,i}=\hat{v}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1} }_{0,i} \qquad \text{for} \quad  i
=1,2\ldots,P_{\ell_2}.
\]
Given the four-coupled updated states
$(\hat{v}_{n, i}^{\ensuremath{\boldsymbol\ell}}, \hat{v}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_1 }_{n,i}, \hat{v}_{n,
  i}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2},\hat{v}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1}}_{n,i})$ at some $n\ge0$,
their next-time prediction states are computed using the dynamics
\begin{equation*}
  \begin{split}
    v^{\ensuremath{\boldsymbol\ell} }_{n+1,i}&= \Psi^{N_{\ell_1}}_n(\hat{v}^{\ensuremath{\boldsymbol\ell} }_{n,i}), \qquad
    v^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_1}_{n+1,i}= \Psi^{N_{\ell_1-1}}_n(\hat{v}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_1 }_{n,i}), \\
    v^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2}_{n+1,i} &= \Psi^{N_{\ell_1}}_n (\hat{v}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2}_{n,i}), \qquad
    v^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1}}_{n+1,i} = \Psi^{N_{\ell_1-1}}_n (\hat{v}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1}}_{n,i}),
  \end{split}
\end{equation*}
where the four particles share the same driving noise.
The sample covariance matrices and Kalman gains are expressed as follows:
\begin{equation*}
  \begin{split}
    C_{n+1}^{\ensuremath{\boldsymbol\ell}}&=\overline{\mathrm{Cov}} [v^{\ensuremath{\boldsymbol\ell}}_{n+1}], \qquad K_{n+1}^{\ensuremath{\boldsymbol\ell} }=C_{n+1}^{\ensuremath{\boldsymbol\ell}}H^\mathbf{T}(HC_{n+1}^{\ensuremath{\boldsymbol\ell}}H^\mathbf{T}+\Gamma)^{-1},\\
    C_{n+1}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_1 }&=\overline{\mathrm{Cov}}[v^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_1}_{n+1}], \! \qquad K_{n+1}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_1 }= C_{n+1}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_1  }H^\mathbf{T}(HC_{n+1}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_1 }H^\mathbf{T}+\Gamma)^{-1},
  \end{split}
\end{equation*}
and
\begin{equation*}
  \begin{split}
    C_{n+1}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2,1}&=\overline{\mathrm{Cov}} [v^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2,1}_{n+1}], \qquad K_{n+1}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2,1}=C_{n+1}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2,1}H^\mathbf{T}(HC_{n+1}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2,1}H^\mathbf{T}+\Gamma)^{-1},\\
    C_{n+1}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2,2}&=\overline{\mathrm{Cov}} [v^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2,2}_{n+1}], \qquad K_{n+1}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2,2}=C_{n+1}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2,2}H^\mathbf{T}(HC_{n+1}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2,2}H^\mathbf{T}+\Gamma)^{-1},\\
    C_{n+1}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1},1}&=\overline{\mathrm{Cov}}[v^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1},1}_{n+1}], \qquad \quad K_{n+1}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1},1}= C_{n+1}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1},1}H^\mathbf{T}(HC_{n+1}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1},1}H^\mathbf{T}+\Gamma)^{-1},\\
    C_{n+1}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1},2}&=\overline{\mathrm{Cov}}[v^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1},2}_{n+1}], \qquad  \quad K_{n+1}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1},2}= C_{n+1}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1},2}H^\mathbf{T}(HC_{n+1}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1},2}H^\mathbf{T}+\Gamma)^{-1}.
  \end{split}
\end{equation*}
The perturbed observations are also particle-wise coupled, so that one obtains the updated particle states:
\begin{equation}
  \left.\begin{split}
    \tilde{y}_{n+1,i}^{\ensuremath{\boldsymbol\ell}}&=y_{n+1}+\eta_{n+1,i}^{\ensuremath{\boldsymbol\ell}}, \\
    \hat{v}_{n+1,i}^{\ensuremath{\boldsymbol\ell}}&=(I-K_{n+1}^{\ensuremath{\boldsymbol\ell} }H)v^{\ensuremath{\boldsymbol\ell}}_{n+1,i}+K_{n+1}^{\ensuremath{\boldsymbol\ell} }\tilde{y}_{n+1,i}^{\ensuremath{\boldsymbol\ell}},\\
    \hat{v}_{n+1,i}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_1}&=(I-K_{n+1}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_1}H)v_{n+1,i}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_1}+K_{n+1}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_1}\tilde{y}_{n+1,i}^{\ensuremath{\boldsymbol\ell}},
  \end{split}\right\} \quad i =1,\ldots, P_{\ell_2}.
  \label{ml:update1}
\end{equation}
\begin{equation}
  \left.\begin{split}
    \hat{v}_{n+1,i}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2,1}&=(I-K_{n+1}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2,1 }H)v^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2,1}_{n+1,i}+K_{n+1}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2,1 }\tilde{y}_{n+1,i}^{\ensuremath{\boldsymbol\ell}},\\
    \hat{v}_{n+1,i}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2,2}&=(I-K_{n+1}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2,2}H)v^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2,2}_{n+1,i}+K_{n+1}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol e}_2,2 }\tilde{y}_{n+1,i+P_{\ell_2-1}}^{\ensuremath{\boldsymbol\ell}},\\
    \hat{v}_{n+1,i}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1},1}&=(I-K_{n+1}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1},1}H)v_{n+1,i}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1},1}+K_{n+1}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1},1}\tilde{y}_{n+1,i}^{\ensuremath{\boldsymbol\ell}},\\
    \hat{v}_{n+1,i}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1},2}&=(I-K_{n+1}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1},2}H)v_{n+1,i}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1},2}+K_{n+1}^{\ensuremath{\boldsymbol\ell}-\ensuremath{\boldsymbol 1},2}\tilde{y}_{n+1,i+P_{\ell_2-1}}^{\ensuremath{\boldsymbol\ell}},
  \end{split}\right\} \quad i =1,\ldots, P_{\ell_2-1}.
  \label{ml:update2}
\end{equation}
where $\{\eta_{n+1,i}^{\ell_2}\}_{i=1}^{P_{\ell_2}}$ are i.i.d.~with
$\eta_{n+1,1}^{\ell_2} \sim N(0,\Gamma)$. Thus, the four-coupled
ensembles are correlated based on the shared initial condition,
driving noise, and perturbed observations.

\begin{figure}[h!]
  \centering
  \includegraphics[width=1.0\textwidth]{{figures/diagramNEW1.png}}\\
  \includegraphics[width=1.0\textwidth]{{figures/diagramNEW2.png}}
  \caption{One prediction-update iteration of the multi-index ensemble Kalman filtering (MIEnKF) estimator described in Section~\ref{ssec:mienkfalg}. The ovals represent four-coupled prediction-state particles, sharing the same driving noise $\omega^{\ensuremath{\boldsymbol\ell}}$ and coupled initial conditions. The respective squares represent updated-state particles sharing the same perturbed observations.

  }
  \label{fig:mienkfDiag}
\end{figure}

\section{MIEnKF complexity}
\label{sec:complexity}
This section presents a complexity result for MIEnKF and a comparison of
MIEnKF, MLEnKF and EnKF in terms of complexity.

\begin{assumption}\label{ass:Psi2}
  For some $N_0, P_0 \in \mathbb{N}$,
  \[
  N_{\ell_1} = N_0\times 2^{\ell_1} \quad \text{and} \quad P_{\ell_2} = P_0\times 2^{\ell_2} \qquad \forall \ensuremath{\boldsymbol\ell} \in \mathbb{N}_0^2.
  \]
  For any $\mathcal{B}^d \backslash \mathcal{B}$-measurable\footnote{ The mapping $\varphi$ is $\mathcal{B}^d\backslash \mathcal{B}-$measurable iff $\varphi^{-1}(B)\in \mathcal{B}^d$ for all $B\in \mathcal{B},$ where $\mathcal{B}$ and $\mathcal{B}^d$ are the $\sigma-$algebras on $\mathbb{R}$ and $\mathbb{R}^d$, respectively.}  mapping $\varphi: \mathbb{R}^d \to \mathbb{R}$ with integrability conditions $\bar \mu_n [\varphi], \mu_n^{\ensuremath{\boldsymbol\ell}}[\varphi] < \infty$  for any $n\ge 0$,
  the four-coupled EnKF estimator $\ensuremath{\boldsymbol \Delta}\mu_n^{\ensuremath{\boldsymbol\ell}}[\varphi]$ satisfies the following assumptions:
  \[
  \abs{\Ex{\ensuremath{\boldsymbol \Delta}\mu_n^{\ensuremath{\boldsymbol\ell}}[\varphi]}}  \lesssim N_{\ell_1}^{-1} P_{\ell_2}^{-1},\tag{\textbf{A1}} \label{ass:A1}
  \]
  \[
  \mathbb{V}[\ensuremath{\boldsymbol \Delta}\mu_n^{\ensuremath{\boldsymbol\ell}}[\varphi]] \lesssim N_{\ell_1}^{-2}P_{\ell_2}^{-2},\tag{\textbf{A2}} \label{ass:A2}
  \]
  and
  \[
  \mathrm{Cost}(\ensuremath{\boldsymbol \Delta}\mu_n^{\ensuremath{\boldsymbol\ell}}[\varphi]) \eqsim  N_{\ell_1}P_{\ell_2} \tag{\textbf{A3}}. \label{ass:A3}
  \footnote{ Note that \ref{ass:A3} is the trivial assumption followed from Section~\ref{ssec:mienkfalg} because the cost of the one prediction-update iteration of the MIEnKF algorithm is linear in $N_{\ell_1}$ and $P_{\ell_2}$.}
  \]
\end{assumption}

\begin{theorem}[MIEnKF complexity]\label{thm:complexityMIEnKF}
  Let Assumptions~\ref{ass:Psi} and~\ref{ass:Psi2} hold and for any
  $\epsilon >0$ consider the MIEnKF estimator with the index set
  $\mathcal{I}=\{\ensuremath{\boldsymbol\ell} \in \mathbb{N}_0^2 \mid \ell_1+ \ell_2 \leq L\}$
  and
  \begin{equation*}
    \begin{split}
      L&= \max\Big(\lceil \log \epsilon^{-1}+\log\log \epsilon^{-1} \rceil-L_0, \; 1\Big), \qquad L_0 \in \mathbb{N}_0, \\
      M_{\ensuremath{\boldsymbol\ell}} &\eqsim \epsilon^{-2}N_{\ell_1}^{-3/2}P_{\ell_2}^{-3/2}.
    \end{split}
  \end{equation*}
  Then it holds for any $n\ge0$ that
  \begin{equation}\label{eq:mienkfMSE}
  \Ex{\prt{\mu^{MI}_n[\varphi]-\bar \mu_n [\varphi]}^2}  \lesssim \epsilon^2,
  \end{equation}
  and the computational cost of the MIEnKF estimator is bounded by
  \[
  \mathrm{Cost(\mu_n^{MI}[\varphi])}\eqsim \epsilon^{-2}.
  \]
\end{theorem}
\begin{proof}
  Introducing $\pm \Ex{\mu^{MI}_n[\varphi]}$ in the mean-squared error, we obtain
  \[
  \Ex{\prt{\mu^{MI}_n[\varphi]-\bar \mu_n [\varphi]}^2}= \Var{\mu^{MI}_n[\varphi]}+\prt{\Ex{\mu^{MI}_n[\varphi]}-\bar \mu_n [\varphi]}^2.
  \]
  For the variance term, the independence of the random variables
  $\{\ensuremath{\boldsymbol \Delta} \mu_n^{\ensuremath{\boldsymbol\ell},m}[\varphi]\}_{(\ensuremath{\boldsymbol\ell}, m)}$ and~\ref{ass:A2} yield
  \[
  \Var{\mu^{MI}_n[\varphi]}= \sum_{ \ensuremath{\boldsymbol\ell}  \in \mathcal{I}} \sum_{m=1}^{M_{\ensuremath{\boldsymbol\ell}}} \frac{\Var{\ensuremath{\boldsymbol \Delta}\mu_n^{\ensuremath{\boldsymbol\ell}, m}[\varphi]}}{M_{\ensuremath{\boldsymbol\ell}}} \lesssim  \sum_{ \ensuremath{\boldsymbol\ell}  \in \mathcal{I}} M_{\ensuremath{\boldsymbol\ell}}^{-1} N_{\ell_1}^{-2}P_{\ell_2}^{-2} \lesssim \epsilon^2.
  \]
  For the squared bias term,~\ref{ass:A1} and the multi-index telescoping properties of the MIEnKF estimator imply that
  \[
  \prt{\Ex{\mu^{MI}_n[\varphi]}-\bar \mu_n [\varphi]}^2 \leq \prt{\sum_{ \ensuremath{\boldsymbol\ell}  \notin \mathcal{I}} \Ex{\ensuremath{\boldsymbol \Delta}\mu_n^{\ensuremath{\boldsymbol\ell}}[\varphi]}}^2 \lesssim  \prt{\sum_{ \ensuremath{\boldsymbol\ell}  \notin \mathcal{I}} N_{\ell_1}^{-1} P_{\ell_2}^{-1}}^2.
  \]
  The mean-squared error bound~\eqref{eq:mienkfMSE} follows by
  \begin{equation*}
    \begin{split}
      \sum_{ \ensuremath{\boldsymbol\ell}  \notin \mathcal{I}} N_{\ell_1}^{-1} P_{\ell_2}^{-1} \lesssim
      \sum_{\ell_1 + \ell_2 > L}^{\infty}2^{-(\ell_1+\ell_2)} = \sum_{k=L+1}^\infty (k+1)2^{-k} \eqsim 2^{-L}L \eqsim \epsilon
    \end{split}
  \end{equation*}
  and
  \[
  \mathrm{Cost}(\mu_n^{MI}[\varphi]) = \sum_{ \ensuremath{\boldsymbol\ell}  \in \mathcal{I}}
  M_{\ensuremath{\boldsymbol\ell}} \mathrm{Cost} (\ensuremath{\boldsymbol \Delta} \mu_n^{\ensuremath{\boldsymbol\ell}} [\varphi])  \eqsim \sum_{ \ensuremath{\boldsymbol\ell}  \in \mathcal{I}}
  M_{\ensuremath{\boldsymbol\ell}} N_{\ell_1}P_{\ell_2}  \eqsim \epsilon^{-2}.
  \]
\end{proof}

\begin{remark}
  For comparison, we briefly recall the cost-versus-accuracy results
  for EnKF and MLEnKF. For any $\epsilon>0$ and
  sufficient regularity
  \[
  \begin{split}
    &\norm{(\mu_n^{N,P}-\bar \mu_n)[\varphi]}_p \lesssim \epsilon, \qquad \mathrm{(EnKF)}\\
    &\norm{(\mu_n^{ML}-\bar \mu_n)[\varphi]}_p \lesssim \epsilon, \qquad \mathrm{(MLEnKF)}
  \end{split}
  \]
  with the computational cost bounded by
  \[
  \begin{split}
    &\mathrm{Cost}(\mu_n^{N,P}[\varphi])\eqsim \epsilon^{-3},\\
    &\mathrm{Cost}(\mu_n^{ML}[\varphi])\eqsim \epsilon^{-2} \abs{\log(\epsilon)}^3.
  \end{split}
  \]
  For more details, see~\cite{hoel2020multilevel}.
\end{remark}

\begin{remark}
  An alternative to~\ref{ass:A2} that is more aligned
  with assumption made for the existing convergence results for MLEnKF
is to assume that
  \[
  \norm{\ensuremath{\boldsymbol \Delta}\mu_n^{\ensuremath{\boldsymbol\ell}}[\varphi]}_p \lesssim N_{\ell_1}^{-1}P_{\ell_2}^{-1} \tag{\textbf{A2$^\star$}}, \label{ass:A2star}
  \]
  for $p\ge 2$. Then, using the same aforementioned index set $\mathcal{I}$ and $L$ and with a slight change in the
  sample size $M_{\ensuremath{\boldsymbol\ell}} \eqsim \epsilon^{-2} N_{\ell_1}^{-4/3}P_{\ell_2}^{-4/3}$, the MIEnKF estimator
  satisfies
  \begin{equation}\label{mienkf:LpError}
    \norm{(\mu_n^{MI}-\bar \mu_n)[\varphi]}_p \lesssim \epsilon,
  \end{equation}
  with the asymptotic MIEnKF cost bounded by $\mathcal{O}(\epsilon^{-2})$.
  This can be proved similarly as the case of Theorem~\ref{thm:complexityMIEnKF}, where the
  $L_p$-norm of the statistical error can be bounded using the Marcinkiewicz-Zygmund inequality:
  \begin{equation*}
    \begin{split}
      \norm{\prt{\mu^{MI}_n-\bar \mu_n}[\varphi]}_p &\leq \norm{\mu^{MI}_n[\varphi]-\mathbb{E}[\mu^{MI}_n][\varphi]}_p+\norm{\mathbb{E}[\mu^{MI}_n][\varphi]-\bar \mu_n[\varphi]}_p\\
      & \lesssim  \sum_{ \ensuremath{\boldsymbol\ell}  \in \mathcal{I}} M_{\ensuremath{\boldsymbol\ell}}^{-1/2} \norm{\ensuremath{\boldsymbol \Delta}\mu_n^{\ensuremath{\boldsymbol\ell}, m}[\varphi] }_p+ \sum_{ \ensuremath{\boldsymbol\ell}  \notin \mathcal{I}} \abs{\mathbb{E}[\ensuremath{\boldsymbol \Delta}\mu_n^{\ensuremath{\boldsymbol\ell},m}[\varphi]]}.
    \end{split}
  \end{equation*}
  We note that \ref{ass:A2} is a weaker assumption than \ref{ass:A2star} since
  \begin{equation}\label{ass:A2p=2}
    \Var{\ensuremath{\boldsymbol \Delta}\mu_n^{\ensuremath{\boldsymbol\ell}}[\varphi]} \leq \norm{\ensuremath{\boldsymbol \Delta}\mu_n^{\ensuremath{\boldsymbol\ell}}[\varphi]}_2^2 \le \norm{\ensuremath{\boldsymbol \Delta}\mu_n^{\ensuremath{\boldsymbol\ell}}[\varphi]}_p^2.
  \end{equation}
\end{remark}

\begin{remark}
  Under more general settings, Assumption~\ref{ass:Psi2} may be transformed into
  \[ \abs{\mathbb{E}[\ensuremath{\boldsymbol \Delta}\mu_n^{\ensuremath{\boldsymbol\ell}}[\varphi]] }  \lesssim N_{\ell_1}^{-\alpha_1} P_{\ell_2}^{-\alpha_2},\tag{$\mathcal{A}_1$} \label{ass:A1gen}\]
  \[\norm{\ensuremath{\boldsymbol \Delta}\mu_n^{\ensuremath{\boldsymbol\ell}}[\varphi]}_p \lesssim N_{\ell_1}^{-\beta_1}P_{\ell_2}^{-\beta_2},\tag{$\mathcal{A}_2$} \label{ass:A2gen}\]
  \[ \mathrm{Cost}(\ensuremath{\boldsymbol \Delta}\mu_n^{\ensuremath{\boldsymbol\ell}}[\varphi]) \eqsim  N_{\ell_1}^{\gamma_1}P_{\ell_2}^{\gamma_2} \tag{$\mathcal{A}_3$}, \label{ass:A3gen}\]
  for some $\alpha_1, \alpha_2, \beta_1, \beta_2, \gamma_1, \gamma_2
  >0$.  The construction of an efficient MIEnKF estimator may then
  lead to a differently shaped index set $\mathcal{I}$ and number of samples
  $\{M_{\ensuremath{\boldsymbol\ell}}\}_{\ensuremath{\boldsymbol\ell}}$, and it may even be sensible to
  consider different common ratios for the geometric sequences
  $\{N_{\ell_1}\}$ and $\{P_{\ell_2}\}$. The problem of optimizing the set $\mathcal{I}$ may be recast as a knapsack problem,
  which is a well-studied optimization problem with many available
  solution algorithms, cf.~\cite{abdo2016} and~\cite[equation~(21)]{abdo2018}.

\end{remark}

\section{Numerical examples}
\label{sec:numerics}
This section presents a numerical comparison of MIEnKF with the
EnKF and MLEnKF methods outlined in Section~\ref{sec:problem}. Three
problems will be considered: the Ornstein-Uhlenbeck (OU)
process, a stochastic differential equation (SDE) with a double-well (DW) potential,
and Langevin dynamics.

We consider SDE on the general form
\begin{equation}\label{sde:genform}
  du=-U'(u)dt+\sigma dW_t,
\end{equation}
with a constant diffusion coefficient $\sigma =0.5$
and two types of potential functions:

$(i)\; \:U(u)=u^2/2, \qquad \qquad  \qquad  \qquad   \mbox{  \textbf{(OU)}}$

$(ii)\: U(u)=u^2/4 +1/(4u^2+2), \qquad \mbox{\textbf{(DW)}}$.\\

The numerical discretizations of~\eqref{sde:genform} are computed using
the Milstein numerical scheme with uniform timestep $\Delta t = 1/N$
for any $N\ge 1$.  The observations of the process $u$
are equally spaced with
observation time interval $\tau=1$, observation operator $H=1$,
$\Gamma=0.1$ and the QoI $\varphi(x)=x$.

To numerically verify assumptions~\ref{ass:A1}
and~\ref{ass:A2star}, the following rates are estimated from
$S$ independent copies of $\ensuremath{\boldsymbol \Delta} \mu_{n}^{\ensuremath{\boldsymbol\ell}}[\varphi]$:
\begin{equation*}
  \begin{split}
    \abs{\mathbb{E}[\ensuremath{\boldsymbol \Delta}\mu_n^{\ensuremath{\boldsymbol\ell}} [\varphi]]} &\approx \Big| \sum_{i=1}^{S} \frac{\ensuremath{\boldsymbol \Delta} \mu_{n,i}^{\ensuremath{\boldsymbol\ell}}[\varphi]}{S} \Big|,\\
    \norm{\ensuremath{\boldsymbol \Delta}\mu_{n}^{\ensuremath{\boldsymbol\ell}}[\varphi]}_2& \approx\sqrt{\frac{1}{S}\sum_{i=1}^{S} \abs{\ensuremath{\boldsymbol \Delta} \mu_{n,i}^{\ensuremath{\boldsymbol\ell}}[\varphi]}^2 }.
  \end{split}
\end{equation*}
We analyze the convergence rates of the methods by computing the
time-averaged root-mean-squared error (RMSE).
\[
\mbox{RMSE}:=\sqrt{ \frac{1}{S(\mathcal{N}+1)}\sum_{i=1}^{S}\sum_{n=0}^{\mathcal{N}}\abs{\mu^{*}_{n,i} [\varphi] - \bar{\mu}_n [\varphi]}^2},
\]
where $\{\mu_{.,i}^{*}[\varphi]\}_{i=1}^{S}$ are independent copies of
$\mu_{.}^{*}[\varphi]$ for the specific methods (EnKF, MLEnKF, and
MIEnKF).

\subsection{Reference solutions and computer architecture}
Since dynamics $\Psi$ is linear for the OU problem, the reference solution $\bar{\mu}_n[\varphi]$ can be computed exactly
using the Kalman filter. However, the reference solution for the DW problem, which involves
nonlinear dynamics, must be approximated. This solution
is computed using the deterministic mean-field EnKF
algorithm, cf.~\cite[Appendix C]{hoel2020multilevel}. A
pseudoreference solution for the final test problem based on Langevin
dynamics is computed by the sample average of
$S=180$ independent simulations of the
MIEnKF estimator at the tolerance
$\epsilon=2^{-11}$ using the following parameters:
\begin{equation*}\label{par:refmienkf}
  \begin{split}
    L&= \lceil L_*+\log_2(L_*)\rceil -1, \qquad \text{with} \quad L_*= \lceil \log_2(\epsilon^{-1})\rceil -1, \\
    N_{\ell_1} &= 4 \times 2^{\ell_1}, \\
    P_{\ell_2} &= 30 \times 2^{\ell_2}, \\
    M_{\ensuremath{\boldsymbol\ell}} &=
    \begin{cases}
      6 \times \lceil \epsilon^{-2} N_{\ell_1}^{-3/2}P_{\ell_2}^{-3/2} \rceil  & \text{if} \quad \ell_1 =0 \mbox{ and } \ell_2 =0,\\
      90 \times \lceil \epsilon^{-2} N_{\ell_1}^{-3/2}P_{\ell_2}^{-3/2} \rceil & \text{if} \quad  1 \le \ell_1+\ell_2 \le  L.
    \end{cases}
  \end{split}
\end{equation*}

The numerical simulations were computed in parallel on 18 cores on
an Intel(R) Xeon(R) CPU E5-2680 v2 20-core processor with 128 GB RAM.
The computer code was written in the Julia programming
language~\cite{Julia-2017}, and it can be downloaded from
\url{https://github.com/GaukharSH/mienkf}.

\subsection{Ornstein-Uhlenbeck process}
\label{ssec:ou}
We consider the SDE~\eqref{sde:genform} with the (OU) potential
function and initial condition $u(0)\sim N(0, \Gamma)$. Convergence rates~\ref{ass:A1}
and~\ref{ass:A2star} shown in Figure~\ref{fig:OUconjRates}  were estimated by the Monte Carlo
method using $S=10^6$ independent samples of
$\ensuremath{\boldsymbol \Delta}\mu_n^{\ensuremath{\boldsymbol\ell}} [\varphi]$.  In the figure,  the left panel shows the weak and
$L_2$ convergence rates over $\mathcal{N}=10$ observation times with
$(\ell_1+ \ell_2) \in [0,7]$, and the right panel shows the ratio of
the rates to $N_{\ell_1}^{-1}P_{\ell_2}^{-1}$. The plane-like flatness
of the right panel for $(\ell_1+\ell_2) \in [1,7]$ validates the
said rate assumptions.

When conducting runtime-versus-accuracy convergence tests for an input
tolerance $\epsilon>0$, we set the parameters of the
respective methods as follows:
\begin{equation}\label{par:enkf}
  \textbf{EnKF:} \quad  P = \lceil15 \epsilon^{-2}\rceil \quad \text{and} \quad  N  = \lceil \epsilon^{-1} \rceil,
\end{equation}
\begin{equation}\label{par:mlenkf}
  \textbf{MLEnKF:} \quad \left\{
  \begin{split}
    L&= \lceil \log_2(\epsilon^{-1})\rceil-1, \\
    N_\ell &= 2 \times 2^{\ell}, \\
    P_\ell &= 10 \times 2^\ell, \\
    M_\ell &=
    \begin{cases}
      2 \times \lceil \epsilon^{-2} L^2 2^{-3}\rceil & \text{if} \quad \ell =0,\\
      \lceil \epsilon^{-2} L^2 2^{-2\ell-3}\rceil & \text{if} \quad  1 \le \ell \le  L,
    \end{cases}
  \end{split}
  \right.
\end{equation}
and
\begin{equation}\label{par:mienkf}
  \textbf{MIEnKF:} \quad \left\{
  \begin{split}
    L&= \lceil L_*+\log_2(L_*)\rceil -1, \quad
    \text{with} \quad L_*= \lceil \log_2(\epsilon^{-1})\rceil -1,\\
    N_{\ell_1} &= 4 \times 2^{\ell_1}, \\
    P_{\ell_2} &= 30 \times 2^{\ell_2}, \\
    M_{\ensuremath{\boldsymbol\ell}} &=
    \begin{cases}
      6 \times \lceil \epsilon^{-2} N_{\ell_1}^{-3/2}P_{\ell_2}^{-3/2} \rceil & \text{if} \quad \ell_1 =0 \mbox{ and } \ell_2 =0,\\
      120 \times \lceil \epsilon^{-2} N_{\ell_1}^{-3/2}P_{\ell_2}^{-3/2} \rceil  & \text{if} \quad  1 \le \ell_1+\ell_2 \le  L.
    \end{cases}
  \end{split}
  \right.
\end{equation}

For a sequence of predefined tolerances $\epsilon = [2^{-4}, 2^{-5},
  2^{-6}, 2^{-7}, 2^{-8}, 2^{-9}]$ for EnKF and MLEnKF, and  $\epsilon = [2^{-4}, 2^{-5},
  2^{-6}, 2^{-7}, 2^{-8}, 2^{-9}, 2^{-10}, 2^{-11}]$ for MIEnKF,
Figure~\ref{fig:OUconvRates} shows the runtime against the RMSE for
the three methods over observation times of
$\mathcal{N}=10$ and $\mathcal{N}=100$ estimated using $S=100$
independent runs. MIEnKF outperforms EnKF and MLEnKF
for sufficiently small tolerances, and the
complexity rate agrees with the theory.

\begin{figure}[h!]
  \centering
  \hspace*{-1cm}
  \includegraphics[width=1.2\textwidth]{{figures/mOUS1e6ParSeed.pdf}}
  \hspace*{-1cm}
  \includegraphics[width=1.2\textwidth]{{figures/sOUS1e6ParSeed.pdf}}
  \caption{\textbf{Ornstein-Uhlenbeck problem. Estimates based on $S=10^6$ independent
      runs (Section~\ref{ssec:ou}).} Top row: Numerical evidence of
    assumption~\ref{ass:A1} over $\mathcal{N}=20$ observation times when using
    $N_{\ell_1}=4\times 2^{\ell_1}$ and $P_{\ell_2}=20 \times
    2^{\ell_2}$.  Bottom row: Similar plots for the verification of assumption~\ref{ass:A2star}.}
  \label{fig:OUconjRates}
\end{figure}

\begin{figure}[h!]
  \centering
  \includegraphics[width=0.47\textwidth]{{figures/loglogMeanRatesOUT10_18core.pdf}}
  \includegraphics[width=0.47\textwidth]{{figures/loglogMeanOUT100parVar.pdf}}
  \caption{\textbf{Ornstein-Uhlenbeck problem. Estimates based on $S=100$ independent
      runs (Section~\ref{ssec:ou}).}  Comparison of the runtime versus root-mean-squared error(RMSE)
    for mean over observation times $\mathcal{N}=10$ (left) and $\mathcal{N}=100$ (right). The solid-crossed line represents MLEnKF and the
    dot-dashed line is a fitted
    $\mathcal{O}(\log(10+\mathrm{Runtime})^{1/3}\mathrm{Runtime}^{-1/2})$
    reference line. The solid-asterisk line represents the MIEnKF and the
    dotted line is a fitted $\mathcal{O}(\mathrm{Runtime}^{-1/2})$ reference
    line. The solid-bulleted line represents EnKF and the dashed line
    is a fitted $\mathcal{O}(\mathrm{Runtime}^{-1/3})$ reference line.}
  \label{fig:OUconvRates}
\end{figure}

\subsection{Double-well SDE}
\label{ssec:dw}
We consider the SDE~\eqref{sde:genform} with the DW potential
function and $u(0)\sim N(0, \Gamma)$.  Similar to the OU case,
Figure~\ref{fig:DWconjRates} provides numerical evidence of the
conjecture rates under assumptions~\ref{ass:A1} and
\ref{ass:A2star}. For the same predefined $\epsilon-$inputs with the
same degrees of freedom setting as in the example of OU, the performance of
the three methods were compared in terms of runtime against RMSE for
observation times $\mathcal{N}=10$ and $\mathcal{N}=100$ and
estimated over $S=100$ independent runs
(Figure~\ref{fig:DWconvRates}). We observe that MIEnKF outperforms EnKF and
MLEnKF for small RMSE.

\begin{figure}[h!]
  \centering
  \hspace*{-1cm}
  \includegraphics[width=1.2\textwidth]{{figures/mDWS1e6ParSeed.pdf}}
  \hspace*{-1cm}
  \includegraphics[width=1.2\textwidth]{{figures/sDWS1e6ParSeed.pdf}}
  \caption{\textbf{Double Well problem. Estimates based on $S=10^6$ independent runs (Section~\ref{ssec:dw}).}

    Top row: Numerical evidence of assumption~\ref{ass:A1}
    for $\mathcal{N}=10$ observation times when using
    $N_{\ell_1}=4\times 2^{\ell_1}$ and
    $P_{\ell_2}=20 \times 2^{\ell_2}$.  Bottom row: Similar plots for
    verifying assumption~\ref{ass:A2star}.  }
  \label{fig:DWconjRates}
\end{figure}

\begin{figure}[h!]
  \centering
  \includegraphics[width=0.47\textwidth]{{figures/loglogMeanRatesDWT10_18core.pdf}}
  \includegraphics[width=0.47\textwidth]{{figures/loglogMeanRatesDWT100_18core.pdf}}
  \caption{\textbf{Double Well problem. Estimates based on $S=100$ independent runs (Section~\ref{ssec:dw}).}
    Similar plots as those shown in Figure~\ref{fig:OUconvRates}.

  }
  \label{fig:DWconvRates}
\end{figure}

\subsection{Langevin SDE}\label{ssec:langevin}
In the last example, we consider the two-dimensional stochastic
Langevin dynamics
\begin{equation}\label{sde:Langevin}
  \begin{split}
    dX_t&=V_tdt, \\
    dV_t&=-U'(X_t)dt - \kappa V_t dt + (2\kappa T)^{1/2} dW_t,
  \end{split}
\end{equation}
where $X_t$ and $V_t$ denotes the particle position and velocity,
respectively, $U(X)$ is the previously introduced DW potential,
$\kappa=2^{-5}\times\pi^2$ is the viscosity and $T=1$ is the temperature.
To improve the pairwise coupling between particles, we used the first-order
symplectic Euler splitting scheme~\cite{muller2015improving}.
The initial conditions are provided by
$X_0 \sim N(0, \Gamma)$ and $P_0 \sim N(0, \Gamma)$ with
$X_0$ and $P_0$ being independent. Further, based on initial test runs,
the method parameters are set to
\begin{equation*}\label{par: enkf}
  \textbf{EnKF:} \quad  P = \lceil 10 \epsilon^{-2}\rceil \quad \text{and} \quad  N  = \lceil \epsilon^{-1} \rceil,
\end{equation*}
\begin{equation*}\label{par:mlenkf}
    \textbf{MLEnKF:} \quad \left\{ \begin{split}
        L&= \lceil \log_2(\epsilon^{-1})\rceil-1,\\
        N_\ell &= 2 \times 2^{\ell},\\
        P_\ell &= 8 \times 2^\ell, \;\\
    M_\ell &=
    \begin{cases}
      2 \times \lceil \epsilon^{-2} L^2 2^{-2}\rceil & \text{if} \quad \ell =0,\\
      \lceil \epsilon^{-2} L^2 2^{-2\ell-2}\rceil & \text{if} \quad  1 \le \ell \le  L,
    \end{cases}
  \end{split}\right.
\end{equation*}
and
\begin{equation*}\label{par:mienkf}
  \textbf{MIEnKF:} \quad \left\{
    \begin{split}
        L&= \lceil L_*+\log_2(L_*)\rceil -1,\; \quad
        \text{with} \quad L_*= \lceil \log_2(\epsilon^{-1})\rceil -1 \\
        N_{\ell_1} &= 4 \times 2^{\ell_1},\\
        P_{\ell_2} &= 20 \times 2^{\ell_2},\\
    M_{\ensuremath{\boldsymbol\ell}} &=
    \begin{cases}
      6 \times \lceil \epsilon^{-2} N_{\ell_1}^{-3/2}P_{\ell_2}^{-3/2} \rceil & \text{if} \quad \ell_1 =0 \mbox{ and } \ell_2 =0,\\
      50 \times \lceil \epsilon^{-2} N_{\ell_1}^{-3/2}P_{\ell_2}^{-3/2} \rceil  & \text{if} \quad  1 \le \ell_1+\ell_2 \le  L.
    \end{cases}
  \end{split}\right.
\end{equation*}

To shed some light on the importance of the temperature parameter,
Figure~\ref{fig:LangevinEvolutionT50Split} illustrates the
phase-portrait time evolution of the realization of Langevin dynamics
up to the final time $\mathcal{N}=50$ for different temperatures $T$. Damping
causes a rapid decay of the velocity from the initial value to zero
when $T=0$. For positive temperatures, thermal fluctuation leads to
more diffusive dynamics.
Figure~\ref{fig:LangevinMeanTruthObsT50SplitH11} shows the
signal-tracking performance of MIEnKF for the full observation operator
\[
H=\begin{bmatrix} 1 & 0\\ 0 &1
\end{bmatrix}
\]
and the partial observation operators $H=[1\; 0]$ or $H=[0\; 1]$, all
computed at the tolerance $\epsilon = 2^{-7}$. The method is tracking
the true state of the observed components well in all cases, but, as
is to be expected, it does not track the true state of unobserved
components with the same level of accuracy. The numerical verification
of assumptions~\ref{ass:A1} and~\ref{ass:A2star} with respect to
different observation operators is shown in
Figures~\ref{fig:LangevinConjRatesH10}
and~\ref{fig:LangevinConjRatesH01}, respectively. For a sequence of
predefined tolerances, $\epsilon = [2^{-4}, 2^{-5}, \ldots, 2^{-9}]$
for EnKF and MLEnKF and $\epsilon = [2^{-4}, 2^{-5}, \ldots,
2^{-10}]$ for MIEnKF, we compare the performance of the three methods
in terms of runtime versus RMSE. We consider
$\mathcal{N}=10$ and $\mathcal{N}=20$ observation times, the QoI
$\varphi(X,V)=X$ and $\varphi(X,V)=V$, and we use $S=90$
indpendent runs of each method to estimate both RMSE and runtime.
Figures~\ref{fig:loglogRatesLangevinH10sepT10} and~\ref{fig:loglogRatesLangevinH11sepT10}
show the results for the observation operators
\[
H=[1\; 0] \quad \text{ and } \quad H=\begin{bmatrix} 1 & 0\\ 0 &1
\end{bmatrix},
\]
respectively. The observed complexity
rates for MIEnKF are close to the theory, and the method is more
efficient than the alternatives for small tolerances in both cases.

\begin{figure}[h!]
  \centering
  \includegraphics[width=0.47\textwidth]{{figures/LangevinPhasePortT0.pdf}}
  \includegraphics[width=0.47\textwidth]{{figures/LangevinPhasePortT001.pdf}}\\
  \includegraphics[width=0.47\textwidth]{{figures/LangevinPhasePortT01.pdf}}
  \includegraphics[width=0.47\textwidth]{{figures/LangevinPhasePortT1diffscale.pdf}}\\
  \caption{\textbf{Time evolution of a solution to Langevin dynamics
      with different temperature $T$ values}. The symplectic Euler scheme is used
    up to final time $\mathcal{N}=50$. The red dot represents the initial
    value.}
  \label{fig:LangevinEvolutionT50Split}
\end{figure}

\begin{figure}[h!]
  \centering
  \includegraphics[width=0.47\textwidth]{{figures/LangevinDataAssH10.pdf}}
  \includegraphics[width=0.47\textwidth]{{figures/LangevinDataAssH01.pdf}}
   \includegraphics[width=0.47\textwidth]{{figures/LangevinDataAssFullObs.pdf}}
  \caption{\textbf{Data assimilation for Langevin dynamics
      given different observation operators $H$,
      Section~\ref{ssec:langevin}}. The blue solid line represents the
    truth, and the red stars are the observations. The cyan dashed lines
    represent MIEnKF mean. The final time $\mathcal{N}=50$ is used with the
    observation timestep $\tau=1.0$}
  \label{fig:LangevinMeanTruthObsT50SplitH11}
\end{figure}

\begin{figure}[h!]
  \centering
  \begin{subfigure}[t]{4.5in}
    \centering
    \hspace*{-1cm}
    \includegraphics[width=1.2\textwidth]{{figures/mWXlanH10par1e6.pdf}}
    \hspace*{-1cm}
    \includegraphics[width=1.2\textwidth]{{figures/mSXlanH10par1e6.pdf}}
    \caption{The particle position $X_t$}\label{fig:LangevinConjRates(a)}
  \end{subfigure}
  \\
  \begin{subfigure}[t]{4.5in}
    \centering
    \hspace*{-1cm}
    \includegraphics[width=1.2\textwidth]{{figures/mWVlanH10par1e6.pdf}}
    \hspace*{-1cm}
    \includegraphics[width=1.2\textwidth]{{figures/mSXlanH10par1e6.pdf}}
    \caption{The particle velocity $V_t$}\label{fig:LangevinConjRates(b)}
  \end{subfigure}
  \caption{\textbf{Langevin dynamics with partial observations, $H=[1 \; 0]$. Estimates based on $S=10^6$
      independent runs  (Section~\ref{ssec:langevin}).} Top row in each subfigure: Numerical
    evidence of assumption~\ref{ass:A1} for $\mathcal{N}=10$ observation
    times when using $N_{\ell_1}=4\times 2^{\ell_1}$ and
    $P_{\ell_2}=20 \times 2^{\ell_2}$.  Bottom row in each subfigure:
    Similar plots for verifying
    assumption~\ref{ass:A2star}.}\label{fig:LangevinConjRatesH10}
\end{figure}

\begin{figure}[h!]
  \centering
  \begin{subfigure}[t]{4.5in}
    \centering
    \hspace*{-1cm}
    \includegraphics[width=1.2\textwidth]{{figures/mWXlanH11.pdf}}
    \hspace*{-1cm}
    \includegraphics[width=1.2\textwidth]{{figures/mSXlanH11.pdf}}
    \caption{The particle position $X_t$}\label{fig:LangevinConjRates(a)}
  \end{subfigure}
  \\
  \begin{subfigure}[t]{4.5in}
    \centering
    \hspace*{-1cm}
    \includegraphics[width=1.2\textwidth]{{figures/mWVlanH11.pdf}}
    \hspace*{-1cm}
    \includegraphics[width=1.2\textwidth]{{figures/mSVlanH11.pdf}}
    \caption{The particle velocity $P_t$}\label{fig:LangevinConjRates(b)}
  \end{subfigure}
  \caption{\textbf{Langevin dynamics with full obervations, $H=[1\;0; 0\; 1]$. Estimates based on $S=10^6$
      independent runs (Section~\ref{ssec:langevin}).} Similar plots as those shown in
    Figure~\ref{fig:LangevinConjRatesH10}.}
  \label{fig:LangevinConjRatesH01}
\end{figure}

\begin{figure}[h!]
	\centering
	\includegraphics[width=0.47\textwidth]{{figures/loglogMeanXLanH10_T10.pdf}}
	\includegraphics[width=0.47\textwidth]{{figures/loglogMeanVLanH10_T10.pdf}}\\
	\includegraphics[width=0.47\textwidth]{{figures/loglogMeanXLanH10_T20.pdf}}
	\includegraphics[width=0.47\textwidth]{{figures/loglogMeanVLanH10_T20.pdf}}
	\caption{\textbf{Langevin dynamics with partial observations, $H=[1 \; 0]$. Estimates based on $S=90$
			independent runs (Section~\ref{ssec:langevin}).}
		\textit{Top row}: Comparison of the runtime versus RMSE for the
	    mean of the component $X$ (left) and the component $V$ (right) over $\mathcal{N}=10$ observation
		times. The
		solid-crossed line represents MLEnKF and the dot-dashed line is a
		fitted
		$\mathcal{O}(\log(10+\mathrm{Runtime})^{1/3}\mathrm{Runtime}^{-1/2})$
		reference line.  The solid-asterisk line represents the MIEnKF and the
		dotted line is a fitted $\mathcal{O}(\mathrm{Runtime}^{-1/2})$ reference
		line. The solid-bulleted line represents EnKF and the dashed line
		is a fitted $\mathcal{O}(\mathrm{Runtime}^{-1/3})$ reference line.
		\textit{Bottom row}: Similar plots for $\mathcal{N}=20$ observation
		times.}
	\label{fig:loglogRatesLangevinH10sepT10}
\end{figure}

\begin{figure}[h!]
	\centering
	\includegraphics[width=0.47\textwidth]{{figures/loglogMeanXlanFullObs_T10.pdf}}
	\includegraphics[width=0.47\textwidth]{{figures/loglogMeanVlanFullObs_T10.pdf}}\\
		\includegraphics[width=0.47\textwidth]{{figures/loglogRatesMeanXLanFullObs_T20.pdf}}
	\includegraphics[width=0.47\textwidth]{{figures/loglogRatesMeanVLanFullObs_T20.pdf}}
	\caption{\textbf{Langevin dynamics with full observations, $H=[1 \; 0; 0 \; 1]$. Estimates based on $S=90$
			independent runs (Section~\ref{ssec:langevin}).}
	Similar plots as those shown in
	Figure~\ref{fig:loglogRatesLangevinH10sepT10}.}
	\label{fig:loglogRatesLangevinH11sepT10}
\end{figure}

\section{Conclusion}\label{sec:conclusion}
We have developed a multi-index Monte Carlo based EnKF method, namely
MIEnKF, which is based on independent samples of four-coupled EnKF
estimators on a multi-index hierarchy of resolution levels. Under
Assumptions~\ref{ass:Psi} and~\ref{ass:Psi2}, we proved the efficiency
of the proposed method theoretically. Further, we numerically verified
that this method is more efficient than MLEnKF and EnKF. A natural
extension of this work is to develop MIEnKF for spatiotemporal models
discretized in both space and time, e.g., stochastic partial
differential equations. We believe that MIEnKF will often outperform
alternative methods prominently when more than two degrees of freedom
need to be discretized. Further, it will be interesting to extend the
MIEnKF coupling theory to the more general ``splitting
parameter''/``common ratio'' $\lambda = P_{\ell_2}/P_{\ell_2-1}$ for
ensemble sizes besides $\lambda=2$ considered in this study.

\medskip

{\bf Acknowledgments } This work was supported by the KAUST Office of
Sponsored Research (OSR) under Award No. URF/1/2584-01-01 and the
Alexander von Humboldt Foundation. G.~Shaimerdenova and R.~Tempone are
members of the KAUST SRI Center for Uncertainty Quantification in
Computational Science and Engineering.