\section{Introduction\label{sec:intro}}

Symmetric positive-definite matrices arise in spatial statistics, Gaussian-process inference, and spatio-temporal filtering, with a wealth of application areas, including geoscience \citep[e.g.,][]{Cressie1993,Banerjee2004}, machine learning \citep[e.g.,][]{Rasmussen2006}, data assimilation \citep[e.g.,][]{Nychka2010,Katzfuss2015b}, and the analysis of computer experiments \citep[e.g.,][]{Sacks1989,Kennedy2001}.
Inference in these areas typically relies on Cholesky decomposition of the positive-definite matrices. However, this operation scales cubically in the dimension of the matrix, and it is thus computationally infeasible for many modern problems and applications, which are increasingly high-dimensional.

Countless approaches have been proposed to address these computational challenges. \citet{Heaton2017} provide a recent review from a spatial-statistics perspective, and \citet{Liu2018} review approaches in machine learning. In  high-dimensional filtering, proposed solutions include low-dimensional approximations \citep[e.g.,][]{Verlaan1995,Pham1998,Wikle1999,Katzfuss2010}, spectral methods \citep[e.g.][]{Wikle1999,Sigrist2015}, and hierarchical approaches \citep[e.g.,][]{Johannesson2003,LiAmbikasaran2014,Saibaba2015,Jurek2018}. Operational data assimilation often relies on ensemble Kalman filters \citep[e.g.,][]{Evensen1994,Burgers1998,Anderson2001,Evensen2007,Katzfuss2015b,Katzfuss2017c}, which represent distributions by samples or ensembles.

Maybe the most promising approximations for spatial data and Gaussian processes implicitly or explicitly rely on sparse Cholesky factors. The assumption of ordered conditional independence in the popular Vecchia approximation \citep{Vecchia1988} and its extensions \citep[e.g.,][]{Stein2004,Datta2016,Guinness2016a,Katzfuss2017a,Katzfuss2018,Katzfuss2020,Schafer2020} implies sparsity in the Cholesky factor of the precision matrix. \citet{Schafer2017} uses an incomplete Cholesky decomposition to construct a sparse approximate Cholesky factor of the covariance matrix. However, these methods are not generally applicable to spatio-temporal filtering, because the assumed sparsity is not preserved under filtering operations.

Here, we relate the sparsity of the Cholesky factors of the covariance matrix and the precision matrix to specific assumptions regarding ordered conditional independence.
We show that these assumptions are simultaneously satisfied for a particular Gaussian-process approximation that we call hierarchical Vecchia (HV), which is a special case of the general Vecchia approximation \citep{Katzfuss2017a} based on hierarchical domain partitioning \citep[e.g.,][]{Katzfuss2015,Katzfuss2017b}.
We show that the HV approximation can be computed using a simple and fast incomplete Cholesky decomposition.

Due to its remarkable property of implying a sparse Cholesky factor whose inverse has equivalent sparsity structure, HV is well suited for extensions to spatio-temporal filtering; this is in contrast to other Vecchia approximations and other spatial approximations relying on sparsity. We provide a scalable HV-based filter for linear Gaussian spatio-temporal state-space models, which is related to the multi-resolution filter of \citet{Jurek2018}.
Further, by combining HV with a Laplace approximation \citep[cf.][]{Zilber2019}, our method can be used for the analysis of non-Gaussian data. Finally, by combining the methods with the extended Kalman filter \citep[e.g.,][Ch.~5]{Grewal1993}, we obtain fast filters for high-dimensional, non-linear, and non-Gaussian spatio-temporal models.
For a given formulation of HV, the computational cost of all of our algorithms scales linearly in the state dimension, assuming sufficiently sparse temporal evolution.

The remainder of this document is organized as follows. In Section \ref{sec:sparsity}, we specify the relationship between ordered conditional independence and sparse (inverse) Cholesky factors. Then, we build up increasingly complex and general methods, culminating in non-linear and non-Gaussian spatio-temporal filters: in Section \ref{sec:hv}, we introduce hierarchical Vecchia for a linear Gaussian spatial field at a single time point; in Section \ref{sec:nongaussian}, we extend this to non-Gaussian data; and in Section \ref{sec:filter}, we consider the general spatio-temporal filtering case, including nonlinear evolution and parameter inference on unknown parameters in the model. Section \ref{sec:comparison} contains numerical comparisons to existing approaches. Section \ref{sec:conclusion} concludes.
Appendices \ref{app:graph}--\ref{app:proofs} contain proofs and further details.
Code implementing our methods and numerical comparisons is available at \url{https://github.com/katzfuss-group/vecchiaFilter}.

\section{Sparsity of Cholesky factors\label{sec:sparsity}}

We begin by specifying the connections between ordered conditional independence and sparsity of the Cholesky factor of the covariance and precision matrix.
\begin{claim}
Let $\mathbf{w}$ be a normal random vector with variance-covariance matrix $\mathbf{K}$.
\begin{enumerate}
    \item Let $\mathbf{L} = \chol(\mathbf{K})$ be the lower-triangular Cholesky factor of the covariance matrix $\mathbf{K}$. For $i>j$:
    \[ \mathbf{L}_{i,j}=0 \iff w_i \perp w_j \, | \, \mathbf{w}_{1:j-1} \]
    \item Let $\mathbf{U} = \rchol(\mathbf{K}^{-1}) = \mathbf{P} \chol( \mathbf{P}\mathbf{K}^{-1}\mathbf{P})\,\mathbf{P}$ be the Cholesky factor of the precision matrix under reverse ordering, where $\mathbf{P}$ is the reverse-ordering permutation matrix. Then $\mathbf{U}$ is upper-triangular, and for $i>j$:
    \[ \mathbf{U}_{j,i}=0 \iff w_i \perp w_j \, | \, \mathbf{w}_{1:j-1},\mathbf{w}_{j+1:i-1} \]
\end{enumerate}
\label{claim:cond-indep}
\end{claim}
The connection between ordered conditional independence and the Cholesky factor of the precision matrix is well known \citep[e.g.,][]{Rue2010}; Part 2 of our claim states this connection under reverse ordering \citep[e.g.,][Prop.~3.3]{Katzfuss2017a}. In Part 1, we consider the lesser-known relationship between ordered conditional independence and sparsity of the Cholesky factor of the covariance matrix, which was recently discussed in \citet[][Sect.~1.4.2]{Schafer2017}. For completeness, we provide a proof of Claim \ref{claim:cond-indep} in Appendix \ref{app:proofs}.

Claim \ref{claim:cond-indep} is crucial for our later developments and proofs. In Section \ref{sec:hv}, we specify a hierarchical Vecchia approximation of Gaussian processes that satisfies both types of conditional independence in Claim \ref{claim:cond-indep}; the resulting sparsity of the Cholesky factor and its inverse allows extensions to spatio-temporal filtering in Section \ref{sec:filter}.

\section{Hierarchical Vecchia for large Gaussian spatial data\label{sec:hv}}

Consider a Gaussian process $x(\cdot)$ and a vector $\mathbf{x} = (x_1, \ldots, x_n)^\top$ representing $x(\cdot)$ evaluated on a grid $\mathcal{S} = \{\mathbf{s}_1,\ldots,\mathbf{s}_n \}$ with $\mathbf{s}_i \in \mathcal{D} \subset \mathbb{R}^d$ and $x_i = x(\mathbf{s}_i)$ for $\mathbf{s}_i \in \mathcal{D}$, $i=1,\ldots,n$. We assume the following model:
\begin{align}
\label{eq:gobs} y_{i} \,|\, \mathbf{x} & \stackrel{ind}{\sim} \mathcal{N}(x_i,\tau_i^2),\qquad i \in \mathcal{I}, \\
\label{eq:gevol} \mathbf{x}  & \sim \mathcal{N}_n(\bm{\mu},\bm{\Sigma}),
\end{align}
where $\mathcal{I} \subset \{1,\ldots,n\}$ are the indices of grid points at which observations are available, and $\mathbf{y}$ is the data vector consisting of these observations $\{y_{i}: i \in \mathcal{I}\}$.
Note that we can equivalently express \eqref{eq:gobs} using matrix notation as $\mathbf{y} \,|\, \mathbf{x} \sim \mathcal{N}(\mathbf{H}\mathbf{x},\mathbf{R})$, where $\mathbf{H}$ is obtained by selecting only the rows with indices $i \in \mathcal{I}$ from an identity matrix, and $\mathbf{R}$ is a diagonal matrix with entries $\{\tau_i^2: i \in \mathcal{I}\}$.

Our interest is in computing the posterior distribution of $\mathbf{x}$ given $\mathbf{y}$, which requires inverting or decomposing an $n \times n$ matrix at a cost of $\mathcal{O}(n^3)$ if $|\mathcal{I}| = \mathcal{O}(n)$. This is computationally infeasible for large $n$.

\subsection{The hierarchical Vecchia approximation\label{sec:hv-intro}}

We now describe a hierarchical Vecchia approximation with unique sparsity and computational properties, which enable fast computation for spatial models as in \eqref{eq:gobs}--\eqref{eq:gevol} and also allow extensions to spatio-temporal filtering as explained later.

Assume that the elements of the vector $\mathbf{x}$ are hierarchically partitioned into a set $\mathcal{X}^{0:M} = \bigcup_{m=0}^M \mathcal{X}^m$, where $\mathcal{X}^m = \bigcup_{k=1}^m \bigcup_{j_k=1}^{J_k} \mathcal{X}_{j_1,\ldots,j_m}$, and $\mathcal{X}_{j_1,\ldots,j_m}$ is a set consisting of $|\mathcal{X}_{j_1,\ldots,j_m}|$ elements of $\mathbf{x}$, such that there is no overlap between any two sets, $\mathcal{X}_{j_1,\ldots,j_m} \cap \mathcal{X}_{i_1,\ldots,i_l} = \emptyset$ for $({j_1,\ldots,j_m}) \neq ({i_1,\ldots,i_l})$.
We assume that $\mathbf{x}$ is ordered according to $\mathcal{X}^{0:M}$, in the sense that if $i>j$, then $x_i \in \mathcal{X}^{m_1}$ and $x_j \in \mathcal{X}^{m_2}$ with $m_1 \geq m_2$.
As a toy example with $n=6$, the vector $\mathbf{x} = (x_1,\ldots,x_6)$ might be partitioned with $M=1$, $J_1 = 2$ as $\mathcal{X}^{0:1} = \mathcal{X}^0 \cup \mathcal{X}^1$, $\mathcal{X}^0 = \mathcal{X} = \{x_1,x_2\}$, and $\mathcal{X}^1 = \mathcal{X}_{1,1} \cup \mathcal{X}_{1,2}$, where $\mathcal{X}_{1,1}=\{x_3,x_4\}$, and $\mathcal{X}_{1,2} = \{x_5,x_6\}$. Another toy example is illustrated in Figure \ref{fig:toy}.

The exact distribution of $\mathbf{x} \sim \mathcal{N}_n(\bm{\mu},\bm{\Sigma})$ can be written as
\[
\textstyle p(\mathbf{x}) = \prod_{m=0}^M \prod_{j_1,\ldots,j_m} p(\mathcal{X}_{j_1,\ldots,j_m}|\mathcal{X}^{0:m-1}, \mathcal{X}_{{j_1,\ldots,j_{m-1}},1:j_m-1}),
\]
where the conditioning set of $\mathcal{X}_{j_1,\ldots,j_m}$ consists of all sets $\mathcal{X}^{0:m-1}$ at lower resolution, plus those at the same resolution that are previous in lexicographic ordering. The idea of \citet{Vecchia1988} was to remove many of these variables in the conditioning set, which for geostatistical applications often incurs only small approximation error due to the so-called screening effect \citep[e.g.,][]{Stein2002,Stein2011}.

\begin{figure}[!htbp]
    \centering
    \begin{subfigure}{.99\textwidth}
    \centering
    \includegraphics[width=0.9\textwidth]{plots/domain.pdf}
    \caption{Iterative domain partitioning for $n=35$ locations}
    \label{fig:domain-&-knots}
\end{subfigure}\\
\vspace{3mm}
    \begin{subfigure}{0.54\textwidth}
    \includegraphics[width=0.85\textwidth]{plots/graph.pdf}
    \caption{Directed acyclic graph (DAG)}
    \label{fig:dag}
    \end{subfigure}
    \begin{subfigure}{.44\textwidth}
    \includegraphics[trim=0mm 12mm 0mm 12mm, clip,width=0.95\textwidth]{plots/matrix.pdf}
    \caption{Sparsity pattern of the $\mathbf{U}$ matrix}
    \label{fig:sparsity-pattern}
    \end{subfigure}
    \caption{Toy example with $n=35$ of the hierarchical Vecchia approximation in \eqref{eq:vecchia-approx} with $M=2$ and $J_1=J_2=2$; the color for each set $\mathcal{X}_{j_1,\ldots,j_m}$ is consistent across (a)--(c).
    (a) Partitioning of the spatial domain $\mathcal{D}$ and the locations $\mathcal{S}$; for resolution $m=0,1,2$, locations of $\mathcal{X}^{0:m}$ (solid dots) and locations of points at finer resolutions ($\circ$).
    (b) DAG illustrating the conditional-dependence structure, with  bigger arrows for connections between vertices at neighboring levels of the hierarchy, to emphasize the tree structure.
    (c) Corresponding sparsity pattern of $\mathbf{U}$ (see Proposition \ref{prop:same-sparsity-factors}), with groups of columns/rows corresponding to different resolutions separated by pink lines, and groups of columns/rows corresponding to different $\mathcal{X}_{j_1,\ldots,j_m}$ at the same resolution separated by blue lines. }
    \label{fig:toy}
\end{figure}

Here we consider a hierarchical Vecchia (HV) approximation of the form
\begin{equation}
\textstyle \hat{p}(\mathbf{x}) = \prod_{m=0}^M \prod_{j_1,\ldots,j_m} p(\mathcal{X}_{j_1,\ldots,j_m}|\mathcal{A}_{j_1,\ldots,j_m}),
\label{eq:vecchia-approx}
\end{equation}
where $\mathcal{A}_{j_1,\ldots,j_m} = \mathcal{X} \cup \mathcal{X}_{j_1} \cup \ldots \cup \mathcal{X}_{{j_1,\ldots,j_{m-1}}}$ is the set of ancestors of $\mathcal{X}_{j_1,\ldots,j_m}$. For example, the set of ancestors of $\mathcal{X}_{2,1,2}$ is $\mathcal{A}_{2,1,2} = \mathcal{X} \cup \mathcal{X}_2 \cup \mathcal{X}_{2,1}$. Thus, $\mathcal{A}_{j_1,\ldots,j_m} = \mathcal{A}_{j_1,\ldots,j_{m-1}} \cup \mathcal{X}_{j_1,\ldots,j_{m-1}}$, and the ancestor sets are nested: $\mathcal{A}_{j_1,\ldots,j_{m-1}} \subset \mathcal{A}_{j_1,\ldots,j_m}$. We can equivalently write \eqref{eq:vecchia-approx} in terms of individual variables as
\begin{equation}
\textstyle \hat{p}(\mathbf{x}) = \prod_{i=1}^n p(x_i|\,\mathcal{C}_i),
\label{eq:vecchia-approx-ind}
\end{equation}
where $\mathcal{C}_i = \mathcal{A}_{j_1,\ldots,j_m} \cup \{x_k \in \mathcal{X}_{j_1,\ldots,j_m} \! : \, k<i \}$ for $x_i \in \mathcal{X}_{j_1,\ldots,j_m}$. The choice of the $\mathcal{C}_i$ involves a trade-off: generally, the larger the $\mathcal{C}_i$, the higher the computational cost (see Proposition \ref{prop:complexity-vf} below), but the smaller the approximation error; HV is exact when all $\mathcal{C}_i = \{x_1,\ldots,x_{i-1}\}$.

Vecchia approximations and their conditional-independence assumptions are closely connected to directed acyclic graphs \citep[DAGs;][]{Datta2016,Katzfuss2017a}. Summarizing briefly, as illustrated in Figure \ref{fig:dag}, we associate a vertex with each set $\mathcal{X}_{j_1,\ldots,j_m}$, and we draw an arrow from the vertex corresponding to $\mathcal{X}_{{i_1,\ldots,i_l}}$ to the vertex corresponding to $\mathcal{X}_{j_1,\ldots,j_m}$ if and only if $\mathcal{X}_{{i_1,\ldots,i_l}}$ is in the conditioning set of $\mathcal{X}_{j_1,\ldots,j_m}$ (i.e., $\mathcal{X}_{{i_1,\ldots,i_l}} \subset \mathcal{A}_{j_1,\ldots,j_m}$). DAGs corresponding to HV approximations always have a tree structure, due to the nested ancestor sets. Necessary terminology and notation from graph theory is reviewed in Appendix \ref{app:graph}.

In practice, as illustrated in Figure \ref{fig:domain-&-knots}, we partition the spatial field $\mathbf{x}$ into the hierarchical set $\mathcal{X}^{0:M}$ based on a recursive partitioning of the spatial domain $\mathcal{D}$ into $J_1$ regions $\mathcal{D}_{1},\ldots,\mathcal{D}_{J_1}$, each of which is again split into $J_2$ regions, and so forth, up to resolution $M$ \citep{Katzfuss2015}: $\mathcal{D}_{j_1,\ldots,j_{m-1}} = \bigcup_{j_m=1}^{J_m} \mathcal{D}_{j_1,\ldots,j_m}$, $m=1,\ldots,M$.
We then set each $\mathcal{X}_{j_1,\ldots,j_m}$ to be a subset of the variables in $\mathbf{x}$ whose location is in $\mathcal{D}_{j_1,\ldots,j_m}$: $\mathcal{X}_{j_1,\ldots,j_m} \subset \{x_i: \mathbf{s}_i \in \mathcal{D}_{j_1,\ldots,j_m}\}$. This implies that the ancestors $\mathcal{A}_{j_1,\ldots,j_m}$ of each set $\mathcal{X}_{j_1,\ldots,j_m}$ consist of the variables associated with regions at lower resolutions $m=0,\ldots,m-1$ that contain $\mathcal{D}_{j_1,\ldots,j_m}$.
Specifically, for all our numerical examples, we set $J_1=\ldots=J_M=2$, and we select each $\mathcal{X}_{j_1,\ldots,j_m}$ corresponding to the first $|\mathcal{X}_{j_1,\ldots,j_m}|$ locations in a maximum-distance ordering \citep{Guinness2016a,Schafer2017} of $\mathcal{S}$ that are contained in $\mathcal{D}_{j_1,\ldots,j_m}$ but are not already in $\mathcal{A}_{j_1,\ldots,j_m}$.

The HV approximation \eqref{eq:vecchia-approx} is closely related to the multi-resolution approximation \citep{Katzfuss2015,Katzfuss2017b}, as noted in \citet[][Sec.~2.5]{Katzfuss2017a}, which in turn is closely related to hierarchical off-diagonal low-rank (HODLR) matrices \citep[e.g.][]{Hackbusch2015,Ambikasaran2016,Saibaba2015,Geoga2018}, as noted in \citet{Jurek2018}.
However, the definition, exposition, and details provided here enable our later proofs, simple incomplete-Cholesky-based computation, and extensions to non-Gaussian data and to nonlinear space-time filtering.

\subsection{Sparsity of the hierarchical Vecchia approximation}

For all Vecchia approximations, the assumed conditional independence implies a sparse Cholesky factor of the precision matrix \citep[e.g.,][Prop.~3.3]{Datta2016,Katzfuss2017a}. The conditional-independence assumption made in our HV approximation also implies a sparse Cholesky factor of the covariance matrix, which is in contrast to many other formulations of the Vecchia approximation:
\begin{proposition}
For the HV approximation in \eqref{eq:vecchia-approx}, we have $\hat p(\mathbf{x}) = \mathcal{N}_n(\mathbf{x}|\bm{\mu},\hat\bm{\Sigma})$. Define $\mathbf{L} = \chol(\hat\bm{\Sigma})$ and $\mathbf{U} = \rchol(\hat\bm{\Sigma}^{-1}) = \mathbf{P} \chol( \mathbf{P}\hat\bm{\Sigma}^{-1}\mathbf{P})\,\mathbf{P}$, where $\mathbf{P}$ is the reverse-ordering permutation matrix.
\begin{enumerate}
    \item For $i\neq j$:
    \begin{enumerate}
        \item $\mathbf{L}_{i,j} = 0$ unless $x_j \in \mathcal{C}_i$
        \item $\mathbf{U}_{j,i} = 0$ unless $x_j \in \mathcal{C}_i$
    \end{enumerate}
    \item $\mathbf{U} = \mathbf{L}^{-\top}$
\end{enumerate}
\label{prop:same-sparsity-factors}
\end{proposition}
The proof relies on Claim \ref{claim:cond-indep}. All proofs can be found in Appendix \ref{app:proofs}. Proposition \ref{prop:same-sparsity-factors} says that the Cholesky factors of the covariance and precision matrix implied by a HV approximation are both sparse, and $\mathbf{U}$ has the same sparsity pattern as $\mathbf{L}^{\top}$. An example of this pattern is shown in Figure \ref{fig:sparsity-pattern}. Furthermore, because $\mathbf{L} = \mathbf{U}^{-\top}$, we can quickly compute one of these factors given the other, as described in Section \ref{sec:vecchia-complexity} below.

For other Vecchia approximations, the sparsity of the prior Cholesky factor $\mathbf{U}$ does not necessarily imply the same sparsity for the Cholesky factor of the posterior precision matrix, and in fact there can be substantial in-fill \citep{Katzfuss2017a}. However, this is not the case for the particular case of HV, for which the posterior sparsity is exactly the same as the prior sparsity:
\begin{proposition}
Assume that $\mathbf{x}$ has the distribution $\hat p(\mathbf{x})$ given by the HV approximation in \eqref{eq:vecchia-approx}. Let $\widetilde\bm{\Sigma} = \var(\mathbf{x}|\mathbf{y})$ be the posterior covariance matrix of $\mathbf{x}$ given data $y_{i} \,|\, \mathbf{x} \stackrel{ind}{\sim} \mathcal{N}(x_i,\tau_i^2)$, $i \in \mathcal{I} \subset \{1,\ldots,n\}$ as in \eqref{eq:gobs}. Then:
\begin{enumerate}
    \item $\widetilde\mathbf{U} = \rchol(\widetilde\bm{\Sigma}^{-1})$ has the same sparsity pattern as $\mathbf{U} = \rchol(\hat\bm{\Sigma}^{-1})$.
    \item $\widetilde\mathbf{L} = \chol(\widetilde\bm{\Sigma})$ has the same sparsity pattern as $\mathbf{L} = \chol(\hat\bm{\Sigma})$.
\end{enumerate}
\label{prop:posterior}
\end{proposition}

\subsection{Fast computation using incomplete Cholesky factorization\label{sec:vecchia-complexity}}

For notational and computational convenience, we assume now that each conditioning set $\mathcal{C}_i$ consists of at most $N$ elements of $\mathbf{x}$. For example, this can be achieved by setting $|\mathcal{X}_{j_1,\ldots,j_m}|\leq r$ with $r = N/(M+1)$. Then $\mathbf{U}$ can be computed using general expressions for the Vecchia approximation in $\mathcal{O}(nN^3)$ time \citep[e.g.,][]{Katzfuss2017a}. Alternatively, inference can be carried out using multi-resolution decompositions \citep{Katzfuss2015,Katzfuss2017b,Jurek2018} in $\mathcal{O}(nN^2)$, but these algorithms are fairly involved.

Instead, we show here how HV inference can be carried out in $\mathcal{O}(nN^2)$ time using standard sparse-matrix algorithms, including the incomplete Cholesky factorization, based on at most $nN$ entries of $\bm{\Sigma}$. Our algorithm, which is based on ideas in \citet{Schafer2017}, is much simpler than multi-resolution decompositions.

\begin{algorithm}[ht]
\caption{Incomplete Cholesky decomposition: $\ichol(\mathbf{A},\mathbf{S})$}
\KwInput{ positive-definite matrix $\mathbf{A} \in \mathbb{R}^{n \times n}$, sparsity matrix $\mathbf{S} \in \{0,1\}^{n \times n}$}
\KwResult{ lower-triangular $n\times n$ matrix $\mathbf{L}$ }
\begin{algorithmic}[1]
\FOR{$i=1$ to $n$}
\FOR{$j=1$ to $i-1$}
\IF{$\mathbf{S}_{i,j}=1$}
\STATE $\mathbf{L}_{i,j} = (\mathbf{A}_{i,j} - \sum_{k=1}^{j-1}\mathbf{L}_{i,k}\mathbf{L}_{j,k})/\mathbf{L}_{j,j}$ \\
\ENDIF
\ENDFOR
\STATE $\mathbf{L}_{i,i} = (\mathbf{A}_{i,i} - \sum_{k=1}^{i-1}\mathbf{L}_{k,k})^{1/2}$
\ENDFOR
\end{algorithmic}
\label{alg:ic0}
\end{algorithm}

The incomplete Cholesky factorization \citep[e.g.,][]{Golub2012}, denoted by $\ichol(\mathbf{A}, \mathbf{S})$ and given in Algorithm \ref{alg:ic0}, is identical to the standard Cholesky factorization of the matrix $\mathbf{A}$, except that we skip all operations that involve elements that are not in the sparsity pattern represented by the zero-one matrix $\mathbf{S}$. It is important to note that to compute $\mathbf{L}=\ichol(\mathbf{A}, \mathbf{S})$ for a large dense matrix $\mathbf{A}$, we do not actually need to form or access the entire $\mathbf{A}$; instead, to reduce memory usage and computational cost, we simply compute $\mathbf{L}=\ichol(\mathbf{A} \circ \mathbf{S}, \mathbf{S})$ based on the sparse matrix $\mathbf{A} \circ \mathbf{S}$, where $\circ$ denotes element-wise multiplication. Thus, while we write expressions like $\mathbf{L}=\ichol(\mathbf{A}, \mathbf{S})$ for notational simplicity below, this should always be read as $\mathbf{L}=\ichol(\mathbf{A} \circ \mathbf{S}, \mathbf{S})$.

For our HV approximation in \eqref{eq:vecchia-approx}, we henceforth set $\mathbf{S}$ to be a sparse lower-triangular matrix with $\mathbf{S}_{i,j}=1$ if $x_j \in \mathcal{C}_i$, and $0$ otherwise. Thus, the sparsity pattern of $\mathbf{S}$ is the same as that of $\mathbf{L}$, and its transpose is that of $\mathbf{U}$ shown in Figure \ref{fig:sparsity-pattern}.

\begin{proposition} \label{prop:chol-lemma}
Assuming \eqref{eq:vecchia-approx}, denote $\var(\mathbf{x}) = \hat\bm{\Sigma}$ and $\mathbf{L} = \chol(\hat\bm{\Sigma})$. Then, $\mathbf{L} = \ichol(\bm{\Sigma}, \mathbf{S})$.
\end{proposition}

Hence, the Cholesky factor of the covariance matrix $\hat\bm{\Sigma}$ implied by the HV approximation can be computed using the incomplete Cholesky algorithm based on the (at most) $nN$ entries of the exact covariance $\bm{\Sigma}$ indicated by $\mathbf{S}$. Using this result, we propose Algorithm \ref{alg:VF} for posterior inference on $\mathbf{x}$ given $\mathbf{y}$.

\begin{algorithm}[ht]
\caption{Posterior inference using hierarchical Vecchia: $\hv(\mathbf{y},\mathbf{S},\bm{\mu},\bm{\Sigma},\mathbf{H},\mathbf{R})$}
\KwInput{data $\mathbf{y}$; sparsity $\mathbf{S}$; $\bm{\mu}, \bm{\Sigma}$ s.t.\ $\mathbf{x} \sim \mathcal{N}_n(\bm{\mu}, \bm{\Sigma})$; obs.\ matrix $\mathbf{H}$; noise variances $\mathbf{R}$}
\KwResult{$\widetilde\bm{\mu}$ and $\widetilde\mathbf{L}$ such that $\hat{p}(\mathbf{x}|\mathbf{y}) = \mathcal{N}_n(\mathbf{x}|\widetilde\bm{\mu},\widetilde\mathbf{L}\widetilde\mathbf{L}^\top)$ }
\begin{algorithmic}[1]
\STATE $\mathbf{L} = \ichol(\bm{\Sigma}, \mathbf{S})$, using Algorithm \ref{alg:ic0}
\STATE $\mathbf{U} = \mathbf{L}^{-\top}$
\STATE $\bm{\Lambda} = \mathbf{U}\bU^\top + \mathbf{H}^\top\mathbf{R}^{-1}\mathbf{H}$
\STATE $\widetilde{\mathbf{U}} = \mathbf{P}\left(\chol(\mathbf{P}\bm{\Lambda}\mathbf{P})\right)\mathbf{P}$, where $\mathbf{P}$ is the order-reversing permutation matrix
\STATE $\widetilde{\mathbf{L}} = \widetilde{\mathbf{U}}^{-\top}$
\STATE $\widetilde{\bm{\mu}} = \bm{\mu} + \widetilde{\mathbf{L}}\widetilde{\mathbf{L}}^\top\mathbf{H}^{\top}\mathbf{R}^{-1}\left(\mathbf{y} - \mathbf{H} \bm{\mu}\right)$
\end{algorithmic}
\label{alg:VF}
\end{algorithm}

By combining the incomplete Cholesky factorization with the results in Propositions \ref{prop:same-sparsity-factors} and \ref{prop:posterior} (saying that all involved Cholesky factors are sparse), we can perform fast posterior inference:
\begin{proposition} \label{prop:complexity-vf}
Algorithm \ref{alg:VF} can be carried out in $\mathcal{O}(nN^2)$ time and $\mathcal{O}(nN)$ space, assuming that $|\mathcal{C}_i| \leq N$ for all $i=1,\ldots,n$.
\end{proposition}

\section{Extensions to non-Gaussian spatial data using the Laplace approximation\label{sec:nongaussian}}

Now consider the model
\begin{align}
\label{eq:ngobs} y_{i} \,|\, \mathbf{x} & \stackrel{ind}{\sim} g_i(y_{i} | x_{i}),\qquad i \in \mathcal{I}, \\
\label{eq:ngevol} \mathbf{x}  & \sim \mathcal{N}_n(\bm{\mu},\bm{\Sigma}),
\end{align}
where $g_i$ is a distribution from an exponential family. Using the HV approximation in \eqref{eq:vecchia-approx}--\eqref{eq:vecchia-approx-ind} for $\mathbf{x}$, the implied posterior can be written as:
\begin{equation}
    \label{eq:nongaussianpost}
\hat{p}(\mathbf{x}|\mathbf{y}) = \frac{p(\mathbf{y}|\mathbf{x}) \hat{p}(\mathbf{x})}{\int p(\mathbf{y}|\mathbf{x}) \hat{p}(\mathbf{x}) d\mathbf{x}} = \frac{(\prod_{i \in \mathcal{I}} g_i(y_i|x_i)) \hat{p}(\mathbf{x})}{\int (\prod_{i \in \mathcal{I}} g_i(y_i|x_i)) \hat{p}(\mathbf{x}) d\mathbf{x}}.
\end{equation}
Unlike in the Gaussian case as in \eqref{eq:gobs}, the integral in the denominator cannot generally be evaluated in closed form, and Markov Chain Monte Carlo methods are often used to numerically approximate the posterior. Instead, \citet{Zilber2019} proposed a much faster method that combines a general Vecchia approximation with the Laplace approximation \citep[e.g.][Sect.~3.4]{Tierney1986,Rasmussen2006}. The Laplace approximation is combined with a Gaussian approximation of the posterior, obtained by carrying out a second-order Taylor expansion of the posterior log-density around its mode. Although the mode cannot generally be obtained in closed form, it can be computed straightforwardly using a Newton-Raphson procedure, because $\log \hat{p}(\mathbf{x}|\mathbf{y}) = \log p(\mathbf{y}|\mathbf{x}) + \log \hat{p}(\mathbf{x}) + c$ is a sum of two concave functions and hence also concave (as a function of $\mathbf{x}$, under appropriate parametrization of the $g_i$).

While each Newton-Raphson update requires the computation and decomposition of the $n \times n$ Hessian matrix, the update can be carried out quickly by making use of the sparsity implied by the Vecchia approximation. To do so, we follow \citet{Zilber2019} in exploiting the fact that the Newton-Raphson update is equivalent to computing the conditional mean of $\mathbf{x}$ given pseudo-data. Specifically, at the $\ell$-th iteration of the algorithm, given the current state value $\mathbf{x}^{(\ell)}$, let us define
\begin{equation}
 \textstyle \mathbf{u}^{(\ell)} = \big[ u^{(\ell)}_i \big]_{i \in \mathcal{I}}, \qquad \text{where} \quad u^{(\ell)}_i = \frac{\partial}{\partial x} \log g_i(y_i|x)\big\vert_{x=x^{(\ell)}_i},
\label{eq:u}
\end{equation}
and
\begin{equation}
\textstyle \mathbf{D}^{(\ell)} = \diag\big(\{d_i^{(\ell)}: i \in \mathcal{I}\}\big), \qquad \text{ where} \quad d_i^{(\ell)} = -\big(\frac{\partial^2}{\partial x^2} \log g_i(y_i|x)\big)^{-1}\big\vert_{x=x^{(\ell)}_i}.
\label{eq:D}
\end{equation}
Then, we compute the next iteration's state value $\mathbf{x}^{(\ell+1)} = \mathbb{E}(\mathbf{x}|\mathbf{t}^{(\ell)})$ as the conditional mean of $\mathbf{x}$ given pseudo-data
$\mathbf{t}^{(\ell)} = \mathbf{x}^{(\ell)} + \mathbf{D}^{(\ell)} \mathbf{u}^{(\ell)}$
assuming Gaussian noise,
$t^{(\ell)}_i | \mathbf{x} \stackrel{ind.}{\sim} \mathcal{N}(x_i,d_i^{(\ell)})$,
$i \in \mathcal{I}$.
\citet{Zilber2019} recommend computing the conditional mean $\mathbb{E}(\mathbf{x}|\mathbf{t}^{(\ell)})$ based on a general-Vecchia-prediction approach proposed in \citet{Katzfuss2018}. Here, we instead compute the posterior mean using Algorithm \ref{alg:VF} based on the HV method described in Section \ref{sec:hv}, due to its sparsity-preserving properties. In contrast to the approach recommended in \citet{Zilber2019}, our algorithm is guaranteed to converge, because it is equivalent to Newton-Raphson optimization of the log of the posterior density in \eqref{eq:nongaussianpost}, which is concave.
Once the algorithm converges to the posterior mode $\widetilde\bm{\mu}$, we obtain a Gaussian HV-Laplace approximation of the posterior as
\[
\hat{p}_L(\mathbf{x}|\mathbf{y}) = \mathcal{N}_n(\mathbf{x}|\widetilde\bm{\mu},\widetilde\mathbf{L}\widetilde\mathbf{L}^\top),
\]
where $\widetilde\mathbf{L}$ is the Cholesky factor of the negative Hessian of the log-posterior at $\widetilde\bm{\mu}$.
Our approach is described in Algorithm \ref{alg:VL}. The main computational expense for each iteration of the \texttt{for} loop is carrying out Algorithm \ref{alg:VF}, and so each iteration requires only $\mathcal{O}(nN^2)$ time.

\begin{algorithm}[ht]
\caption{Hierarchical-Vecchia-Laplace inference: $\hvl(\mathbf{y},\mathbf{S},\bm{\mu},\bm{\Sigma},\{g_i\})$}
\KwInput{ data $\mathbf{y}$; sparsity $\mathbf{S}$; $\bm{\mu}, \bm{\Sigma}$ such that $\mathbf{x} \sim \mathcal{N}_n(\bm{\mu}, \bm{\Sigma})$; likelihoods $\{g_i: i \in \mathcal{I} \}$}
\KwResult{ $\widetilde\bm{\mu}$ and $\widetilde\mathbf{L}$ such that $\hat{p}_L(\mathbf{x}|\mathbf{y}) = \mathcal{N}_n(\mathbf{x}|\widetilde\bm{\mu},\widetilde\mathbf{L}\widetilde\mathbf{L}^\top)$ }
\begin{algorithmic}[1]
\STATE Initialize $\mathbf{x}^{(0)} = \bm{\mu}$
\STATE Set $\mathbf{H} = \mathbf{I}_{\mathcal{I},:}$ as the rows $\mathcal{I}$ of the $n \times n$ identity matrix $\mathbf{I}$
\FOR{$\ell = 0,1,2,\ldots$}
    \STATE Calculate $\mathbf{u}^{(\ell)}$ as in \eqref{eq:u}, $\mathbf{D}^{(\ell)}$ as in \eqref{eq:D}, and $\mathbf{t}^{(\ell)} = \mathbf{x}^{(\ell)} + \mathbf{D}^{(\ell)} \mathbf{u}^{(\ell)}$
    \STATE Calculate $[\mathbf{x}^{(\ell+1)},\widetilde\mathbf{L}]=\hv(\mathbf{t}^{(\ell)},\mathbf{S},\bm{\mu},\bm{\Sigma},\mathbf{H},\mathbf{D}^{(\ell)})$ using Algorithm \ref{alg:VF}
    \IF{$\norm{\mathbf{x}^{(\ell+1)} - \mathbf{x}^{(\ell)}}/\norm{\mathbf{x}^{(\ell)}} < \epsilon$}
    \STATE \algorithmicbreak
    \ENDIF
\ENDFOR
\RETURN $\widetilde{\bm{\mu}}=\mathbf{x}^{(\ell+1)}$ and $\widetilde\mathbf{L}$
\end{algorithmic}
\label{alg:VL}
\end{algorithm}

In the Gaussian case, when $g_i(y_i | x_i) = \mathcal{N}(y_i|a_i x_i,\tau_i^2)$ for some $a_i \in \mathbb{R}$, it can be shown using straightforward calculations that the pseudo-data $t_i = y_i/a_i$ and pseudo-variances $d_i = \tau_i^2$ do not depend on $\mathbf{x}$, and so Algorithm \ref{alg:VL} converges in a single iteration. If, in addition, $a_i=1$ for all $i=1,\ldots,n$, then \eqref{eq:ngobs} becomes equivalent to \eqref{eq:gobs}, and Algorithm \ref{alg:VL} simplifies to Algorithm \ref{alg:VF}.
For non-Gaussian data, our Laplace and Gaussian approximations introduce an additional source of error. While this error is difficult to quantify theoretically, empirical studies \citep[e.g.,][]{Bonat2016, Zilber2019} have shown that Laplace-type approximations can be very accurate and can strongly outperform sampling-based approaches such as Markov Chain Monte Carlo.

\section{Fast filters for spatio-temporal models} \label{sec:filter}

\subsection{Linear evolution} \label{sec:linear-filter}

We now turn to a spatio-temporal state-space model (SSM), which adds a temporal evolution model to the spatial model \eqref{eq:ngobs} considered in Section \ref{sec:nongaussian}. For now, assume that the evolution is linear.
Starting with an initial distribution $\mathbf{x}_0 \sim \mathcal{N}_{n}(\bm{\mu}_{0|0},\bm{\Sigma}_{0|0})$, we consider the following SSM for discrete time $t=1,2,\ldots$:
\begin{align}
y_{ti} \,|\, \mathbf{x}_t & \stackrel{ind}{\sim} g_{ti}(y_{ti} | x_{ti}), \qquad i \in \mathcal{I}_t \label{eq:linear-dynamics-obs}\\
\mathbf{x}_t \,|\, \mathbf{x}_{t-1} & \sim \mathcal{N}_n(\mathbf{E}_t\mathbf{x}_{t-1},\mathbf{Q}_t),
\label{eq:linear-dynamics}
\end{align}
where $\mathbf{y}_{t}$ is the data vector consisting of $n_t \leq n$ observations $\{y_{ti}: i \in \mathcal{I}_t\}$, $\mathcal{I}_t \subset \{1,\ldots,n\}$ contains the observation indices at time $t$, $g_{ti}$ is a distribution from the exponential family, $\mathbf{x}_t = (x_1,\ldots,x_n)^\top$ is the latent spatial field of interest at time $t$ observed at a spatial grid $\mathcal{S}$, and $\mathbf{E}_t$ is a sparse $n \times n$ evolution matrix.

At time $t$, our goal is to obtain or approximate the filtering distribution $p(\mathbf{x}_t|\mathbf{y}_{1:t})$ of $\mathbf{x}_t$ given data $\mathbf{y}_{1:t}$ up to the current time $t$. This task, also referred to as data assimilation or on-line inference, is commonly encountered in many fields of science whenever one is interested in quantifying the uncertainty in the current state or in obtaining forecasts into the future.
If the observation equations $g_{ti}$ are all Gaussian, the filtering distribution can be derived using the Kalman filter \citep{Kalman1960} for small to moderate $n$. At each time $t$, the Kalman filter consist of a forecast step that computes $p(\mathbf{x}_t|\mathbf{y}_{1:t-1})$, and an update step which then obtains $p(\mathbf{x}_t|\mathbf{y}_{1:t})$. For linear Gaussian SSMs, both of these distributions are multivariate normal.

Our Kalman-Vecchia-Laplace (KVL) filter extends the Kalman filter to high-dimensional SSMs (i.e., large $n$) with non-Gaussian data, as in \eqref{eq:linear-dynamics-obs}--\eqref{eq:linear-dynamics}.
Its update step is very similar to the inference problem in Section \ref{sec:nongaussian}, and hence it essentially consists of the HVL in Algorithm \ref{alg:VL}.
We complement this update with a forecast step, in which the moment estimates are propagated forward using the temporal evolution model. This forecast step is exact, and so the KVL approximation error is solely due to the HVL approximation at each update step. The KVL filter is given in Algorithm \ref{alg:KVLfilter}.

\begin{algorithm}[ht]
\caption{Kalman-Vecchia-Laplace (KVL) filter}
\KwInput{$\mathbf{S}$, $\bm{\mu}_{0|0}$, $\bm{\Sigma}_{0|0}$, $\{(\mathbf{y}_t,\mathbf{E}_t,\mathbf{Q}_t,\{g_{t,i}\}): t=1,2,\ldots \}$}
\KwResult{$\bm{\mu}_{t|t}, \mathbf{L}_{t|t}$, such that $\hat{p}(\mathbf{x}_t|\mathbf{y}_{1:t}) = \mathcal{N}_n(\mathbf{x}_t|\bm{\mu}_{t|t},\mathbf{L}_{t|t}\mathbf{L}_{t|t}^\top)$}
\begin{algorithmic}[1]
\STATE Compute $\mathbf{U}_{0|0} = \ichol(\bm{\Sigma}_{0|0}, \mathbf{S})$ and $\mathbf{L}_{0|0} = \mathbf{U}_{0|0}^{-\top}$
\FOR{$t = 1, 2, \dots$}

    \STATE Forecast: $\bm{\mu}_{t|t-1} = \mathbf{E}_t\bm{\mu}_{t-1|t-1}$ and $\mathbf{L}_{t|t-1} = \mathbf{E}_t\mathbf{L}_{t-1|t-1}$ \label{ln:evolL}
    \STATE For all $(i,j)$ with $\mathbf{S}_{i,j}=1$: $\bm{\Sigma}_{t|t-1; i,j} = \mathbf{L}_{t|t-1;i,:} \mathbf{L}_{t|t-1;j,:}^\top + \mathbf{Q}_{t;i,j}$ \label{ln:forecov}
    \STATE Update: $[\bm{\mu}_t, \mathbf{L}_{t|t}] = \hvl(\mathbf{y}_t, \mathbf{S}, \bm{\mu}_{t|t-1}, \bm{\Sigma}_{t|t-1}, \left\{g_{t,i}\right\}$) using Algorithm \ref{alg:VL}
    \RETURN $\bm{\mu}_{t|t}, \mathbf{L}_{t|t}$
\ENDFOR
\end{algorithmic}
\label{alg:KVLfilter}
\end{algorithm}

In Line \ref{ln:forecov}, $\mathbf{L}_{t|t-1;i,:}$ denotes the $i$th row of $\mathbf{L}_{t|t-1}$. The KVL filter scales well with the state dimension $n$. The evolution matrix $\mathbf{E}_t$, which is often derived using a forward-finite-difference scheme and thus has only a few nonzero elements in each row, can be quickly multiplied with $\mathbf{L}_{t-1|t-1}$ in Line \ref{ln:evolL}, as the latter is sparse (see Section \ref{sec:vecchia-complexity}). The $\mathcal{O}(nN)$ necessary entries of $\bm{\Sigma}_{t|t-1}$ in Line \ref{ln:forecov} can also be calculated quickly due to the sparsity of $\mathbf{L}_{t|t-1;i,:}$. The low computational cost of the \texttt{HVL} algorithm has already been discussed in Section \ref{sec:nongaussian}. Thus, assuming sufficiently sparse $\mathbf{E}_t$, the KVL filter scales approximately as $\mathcal{O}(nN^2)$ per iteration.
In the case of Gaussian data (i.e., all $g_{ti}$ in \eqref{eq:linear-dynamics-obs} are Gaussian), our KVL filter will produce essentially equivalent filtering distributions as the more complicated multi-resolution filter of \citet{Jurek2018}.

\subsection{An extended filter for nonlinear evolution}

Finally, we consider a nonlinear and non-Gaussian model, which extends \eqref{eq:linear-dynamics-obs}--\eqref{eq:linear-dynamics} by allowing nonlinear evolution operators, $\mathcal{E}_t: \mathbb{R}^n \rightarrow \mathbb{R}^n$. This results in the model
\begin{align}
\label{eq:obs} y_{ti} \,|\, \mathbf{x}_t & \stackrel{ind}{\sim} g_{ti}(y_{ti} | x_{ti}), \qquad i \in \mathcal{I}_t \\
\label{eq:evol} \mathbf{x}_t \,|\, \mathbf{x}_{t-1} & \sim \mathcal{N}_n(\mathcal{E}_t(\mathbf{x}_{t-1}),\mathbf{Q}_t).
\end{align}

Due to the nonlinearity of the evolution operator $\mathcal{E}_t$, the KVL filter in Algorithm \ref{alg:KVLfilter} is not directly applicable anymore. However, similar inference is still possible as long as the evolution is not too far from linear. Approximating the evolution as linear is generally reasonable if the time steps are short, or if the measurements are highly informative. In this case, we propose the extended Kalman-Vecchia-Laplace filter (EKVL) in Algorithm~\ref{alg:eKVLfilter}, which approximates the extended Kalman filter \citep[e.g.,][Ch.~5]{Grewal1993} and extends it to non-Gaussian data using the Vecchia-Laplace aproach. For the forecast step, EKVL computes the forecast mean as $\bm{\mu}_{t|t-1} = \mathcal{E}_t(\bm{\mu}_{t-1|t-1})$. The forecast covariance matrix $\bm{\Sigma}_{t|t-1}$ is obtained as before, after approximating the evolution using the Jacobian as $\mathbf{E}_t = \frac{\partial \mathcal{E}_t(\mathbf{y}_{t-1})}{\partial \mathbf{y}_{t-1}} \big|_{\mathbf{y}_{t-1} = \bm{\mu}_{t-1|t-1} }$.  Errors in the forecast covariance matrix due to this linear approximation can be captured in the innovation covariance, $\mathbf{Q}_t$. If the Jacobian matrix cannot be computed, it is sometimes possible to build a statistical emulator \citep[e.g.,][]{Kaufman2011} instead, which approximates the true evolution operator.

Once $\bm{\mu}_{t|t-1}$ and $\bm{\Sigma}_{t|t-1}$ have been obtained, the update step of the EKVL proceeds exactly as in the KVL filter by approximating the forecast distribution as Gaussian.

\begin{algorithm}[ht]
\caption{Extended Kalman-Vecchia-Laplace (EKVL) filter}
\KwInput{$\mathbf{S}$, $\bm{\mu}_{0|0}$, $\bm{\Sigma}_{0|0}$, $\{(\mathbf{y}_t,\mathcal{E}_t,\mathbf{Q}_t,\{g_{t,i}\}): t=1,2,\ldots \}$}
\KwResult{$\bm{\mu}_{t|t}, \mathbf{L}_{t|t}$, such that $\hat{p}(\mathbf{x}_t|\mathbf{y}_{1:t}) = \mathcal{N}_n(\mathbf{x}_t|\bm{\mu}_{t|t},\mathbf{L}_{t|t}\mathbf{L}_{t|t}^\top)$}
\begin{algorithmic}[1]
\STATE Compute $\mathbf{U}_{0|0} = \ichol(\bm{\Sigma}_{0|0}, \mathbf{S})$ and $\mathbf{L}_{0|0} = \mathbf{U}_{0|0}^{-\top}$
\FOR{$t = 1, 2, \dots$}
    \STATE Calculate $\mathbf{E}_t = \frac{\partial \mathcal{E}_t(\mathbf{x}_{t-1})}{\partial \mathbf{x}_{t-1}} \big|_{\mathbf{x}_{t-1} = \bm{\mu}_{t-1|t-1} }$ \label{ln:jacobian}
    \STATE Forecast: $\bm{\mu}_{t|t-1} = \mathcal{E}_t(\bm{\mu}_{t-1|t-1})$ and $\mathbf{L}_{t|t-1} = \mathbf{E}_t\mathbf{L}_{t-1|t-1}$
     \STATE For all $(i,j)$ with $\mathbf{S}_{i,j}=1$: $\bm{\Sigma}_{t|t-1; i,j} = \mathbf{L}_{t|t-1;i,:} \mathbf{L}_{t|t-1;j,:}^\top + \mathbf{Q}_{t;i,j}$
    \STATE Update: $[\bm{\mu}_t, \mathbf{L}_{t|t}] = \hvl(\mathbf{y}_t, \mathbf{S}, \bm{\mu}_{t|t-1}, \bm{\Sigma}_{t|t-1}, \left\{g_{t,i}\right\}$) using Algorithm \ref{alg:VL}
    \RETURN $\bm{\mu}_{t|t}, \mathbf{L}_{t|t}$
\ENDFOR
\end{algorithmic}
\label{alg:eKVLfilter}
\end{algorithm}

Similarly to Algorithm \ref{alg:KVLfilter}, EKVL scales very well with the dimension of $\mathbf{x}$, the only difference being the additional operation of calculating the Jacobian in Line \ref{ln:jacobian}, whose cost is problem dependent. Only those entries of $\mathbf{E}_t$ need to be calculated that are multiplied with non-zero entries of $\mathbf{L}_{t-1|t-1}$, whose sparsity structure is known ahead of time.

\subsection{A particle-EKVL filter in case of unknown parameters}

The distributions and matrices in model \eqref{eq:obs}--\eqref{eq:evol} may depend on parameters $\bm{\theta}_t$ at each time $t$, which we have implicitly assumed to be known thus far. We now discuss the case of a (small) number of unknown parameters $\bm{\theta}_t$. Specifically, $\bm{\mu}_{0|0}$ and $\bm{\Sigma}_{0|0}$ may depend on $\bm{\theta}_0$, and the quantities $\{g_{t,i}\}$, $\mathcal{E}_t$, and $\mathbf{Q}_t$ at each time $t$ may depend on $\bm{\theta}_t$.
There are two main approaches to simultaneous filtering for the state $\mathbf{x}_t$ and the parameters $\bm{\theta}_t$: state augmentation and Rao-Blackwellized filters \citep{Doucet2011}. The main idea behind the former is to include $\bm{\theta}_t$ in the state vector $\mathbf{x}_t$ and to modify the evolution and the model error matrices accordingly, but this approach is known to work poorly in certain cases \citep[e.g.,][]{delsole2010state,Katzfuss2017c}. Thus, following \citet{Jurek2018}, we now present a Rao-Blackwellized filter in which integration over $\mathbf{x}_t$ is performed based on our HVL approximation.

Writing $\bm{\theta}_{0:t} = (\bm{\theta}_0, \ldots, \bm{\theta}_t)$, the integrated likelihood at time $t$ is given by
\begin{equation}
\textstyle p(\mathbf{y}_{1:t}|\bm{\theta}_{0:t}) = p(\mathbf{y}_1 | \bm{\theta}_{0:1})\prod_{k=2}^t p(\mathbf{y}_{k}|\mathbf{y}_{1:k-1}, \bm{\theta}_{0:k}).
\label{eq:int-lik-approx}
\end{equation}
It is well known that
\begin{equation}
p(\mathbf{y}_{t}|\mathbf{y}_{1:t-1}, \bm{\theta}_{0:t}) =
    \frac{p(\mathbf{y}_t, \mathbf{x}_t | \mathbf{y}_{1:t-1},\bm{\theta}_{0:t})}{p(\mathbf{x}_t|\mathbf{y}_{1:t},\bm{\theta}_{0:t})} =
    \frac{p(\mathbf{y}_t | \mathbf{x}_t,\bm{\theta}_{t}) p(\mathbf{x}_t|\mathbf{y}_{1:t-1},\bm{\theta}_{0:t})}{ p(\mathbf{x}_t|\mathbf{y}_{1:t},\bm{\theta}_{0:t})},
    \label{eq:lik-decom}
\end{equation}
where $p(\mathbf{y}_t | \mathbf{x}_t,\bm{\theta}_{t})$ is available in closed form from \eqref{eq:obs}, and the forecast and filtering distributions can be approximated using the EKVL, to obtain
\begin{equation}
    \mathcal{L}_t(\bm{\theta}_{0:t}) \colonequals \hat p(\mathbf{y}_{t}|\mathbf{y}_{1:t-1}, \bm{\theta}_{0:t})  = \frac{p(\mathbf{y}_t | \mathbf{x}_t,\bm{\theta}_{t})\mathcal{N}(\bm{\mu}_{t|t-1}, \bm{\Sigma}_{t|t-1})}{\mathcal{N}(\bm{\mu}_{t|t}, \bm{\Sigma}_{t|t})}.
\label{eq:approximate-yAR}
\end{equation}
The normal densities can be quickly evaluated for given parameter values $\bm{\theta}_{0:t}$, because Algorithm \ref{alg:eKVLfilter} calculates sparse Cholesky factors of their precision matrices. For $t=1$, the term $\mathcal{L}_1(\bm{\theta}_{0:1}) \colonequals \hat p(\mathbf{y}_1|\bm{\theta}_{0:1})$ can be approximated in a similar way using $\bm{\mu}_{0|0}$ and $\bm{\Sigma}_{0|0}$.

The particle-EKVL filter is given by Algorithm \ref{alg:peKVLfilter}, assuming that the parameter priors are given by $f_0(\bm{\theta}_0)$ and then recursively by $f_t(\bm{\theta}_t|\bm{\theta}_{t-1})$.

\begin{algorithm}[ht]
\caption{Particle-EKVL filter}
\KwInput{$\mathbf{S}$, $\bm{\mu}_{0|0}$, $\bm{\Sigma}_{0|0}$, $\{(\mathbf{y}_t,\mathcal{E}_t,\mathbf{Q}_t,\{g_{t,i}\}): t=1,2,\ldots \}$, priors $\{f_t\}$, proposal distributions $\{q_t\}$, desired number of particles $N_p$}
\KwResult{
    $\{(\bm{\theta}_t^{(l)},w_t^{(l)},\bm{\mu}_{t|t}^{(l)},\mathbf{L}^{(l)}_{t|t}): l=1,\ldots,N_p\}$,
    such that $\hat p(\bm{\theta}_t,\mathbf{x}_t|\mathbf{y}_{1:t}) = \sum_{l=1}^{N_p} w_t^{(l)} \delta_{\bm{\theta}_t^{(l)}}(\bm{\theta}_t) \mathcal{N}_n(\mathbf{x}_t|\bm{\mu}_{t|t}^{(l)},\mathbf{L}_{t|t}^{(l)}\mathbf{L}_{t|t}^{(l)}{}^\top)$}
\begin{algorithmic}[1]
\FOR{$l = 1, 2, \dots, N_p$}
    \STATE Draw $\bm{\theta}_0^{(l)} \sim f_0(\bm{\theta}_0)$ and set weight $w_0^{(l)} = 1/N_p$
    \STATE Compute $\mathbf{L}_{0|0}(\bm{\theta}_0^{(l)}) = \ichol(\bm{\Sigma}_{0|0}(\bm{\theta}_0^{(l)}), \mathbf{S})$
    \STATE Compute $\bm{\mu}^{(l)}_{0|0}(\bm{\theta}_0^{(l)})$ and $\mathbf{U}_{0|0}(\bm{\theta}_0^{(l)}) = \mathbf{L}_{0|0}^{-\top}(\bm{\theta}_0^{(l)})$
\ENDFOR
\FOR{$t = 1, 2, \ldots$}
    \FOR{$l = 1, 2, \ldots, N_p$}
        \STATE Draw $\bm{\theta}_t^{(l)} \sim q_t(\bm{\theta}_t^{(l)}|\bm{\theta}_{t-1}^{(l)})$
        \STATE Calculate $\mathbf{E}_t^{(l)} = \frac{\partial \mathcal{E}_t(\mathbf{y}_{t-1}, \bm{\theta}_t^{(l)})}{\partial \mathbf{y}_{t-1}} \big|_{\mathbf{y}_{t-1} = \bm{\mu}_{t-1|t-1}(\bm{\theta}_{t-1}^{(l)}) }$ \label{ln:jacobian2}
        \STATE Forecast: $\bm{\mu}_{t|t-1}^{(l)} = \mathcal{E}_t(\bm{\mu}_{t-1|t-1}, \bm{\theta}_{t-1}^{(l)})$ and $\mathbf{L}_{t|t-1}^{(l)} = \mathbf{E}_t^{(l)}\mathbf{L}_{t-1|t-1}^{(l)}$
        \STATE For $(i,j)$ s.t.\ $\mathbf{S}_{i,j}=1$: $\bm{\Sigma}_{t|t-1; i,j}^{(l)} = \mathbf{L}_{t|t-1;i,:}^{(l)} (\mathbf{L}_{t|t-1;j,:}^{(l)})^\top + \mathbf{Q}_{t;i,j}(\bm{\theta}^{(l)}_t)$
        \STATE Update: $[\bm{\mu}_t^{(l)}, \mathbf{L}_{t|t}^{(l)}] = \hvl(\mathbf{y}_t, \mathbf{S}, \bm{\mu}_{t|t-1}^{(l)}, \bm{\Sigma}_{t|t-1}^{(l)}, \{g_{t,i}(\bm{\theta}^{(l)}_t)\}$)
        \STATE Calculate $\mathcal{L}_t(\bm{\theta}_{0:t}^{(l)})$ as in \eqref{eq:approximate-yAR}
        \STATE Update particle weight $w_t^{(l)} \propto w_{t-1}^{(l)} \mathcal{L}_t(\bm{\theta}_{0:t}^{(l)})f_t(\bm{\theta}_{t}^{(l)} | \bm{\theta}_{t-1}^{(l)}) / q_t(\bm{\theta}_{t}^{(l)} | \bm{\theta}_{t-1}^{(l)})$
        \RETURN $\bm{\mu}_{t|t}^{(l)}, \mathbf{L}_{t|t}^{(l)}, \bm{\theta}_t^{(l)}, w_t^{(l)}$
    \ENDFOR
    \STATE Resample $\{(\bm{\theta}_{t}^{(l)}, \bm{\mu}_{t|t}^{(l)}, \mathbf{L}_{t|t}^{(l)})\}_{l=1}^{N_p}$ with weights $\{w_t^{(l)}\}_{l=1}^{N_p}$ to obtain equally weighted particles \citep[e.g.,][]{Douc2005}
\ENDFOR
\end{algorithmic}
\label{alg:peKVLfilter}
\end{algorithm}

\section{Numerical comparison\label{sec:comparison}}

\subsection{Methods and criteria\label{sec:simul-general}}

We considered and compared the following methods:
\begin{description}
\item[Hierarchical Vecchia (HV):] Our methods as described in this paper.
\item[Low rank (LR):] A special case of HV with $M=1$, in which the diagonal and the first $N$ columns of $\mathbf{S}$ are nonzero, and all other entries are zero. This results in a matrix approximation $\hat\bm{\Sigma}$ that is of rank $N$ plus diagonal, known as the modified predictive process \citep[][]{Banerjee2008,Finley2009} in spatial statistics. LR has the same computational complexity as HV.
\item[Dense Laplace (DL):] A further special case of HV with $M=0$, in which $\mathbf{S}$ is a fully dense matrix of ones. Thus, there is no error due to the Vecchia approximation, and so in the non-Gaussian spatial-only setting, this is equivalent to a Gaussian Laplace approximation. DL will generally be more accurate than HV and low-rank, but it scales as $\mathcal{O}(n^3)$ and is thus not feasible for high dimension $n$.
\end{description}

For each scenario below, we simulated observations using \eqref{eq:obs}, taking $g_{t,i}$ to be each of four exponential-family distributions: Gaussian, $\mathcal{N}(x,\tau^2)$; logistic Bernoulli, $\mathcal{B}(1/(1+e^{-x}))$; Poisson, $\mathcal{P}(e^x)$; and gamma, $\mathcal{G}(a, ae^{-x})$, with shape parameter $a=2$. For most scenarios, we assumed a moderate state dimension $n$, so that DL remained feasible; a large $n$ was considered in Section \ref{sec:big-simul}.

The main metric to compare HV and LR was the difference in KL divergence between their posterior or filtering distributions and those generated by DL; as the exact distributions were not known here, we approximated this metric by the average difference in log scores \citep[dLS; e.g.,][]{Gneiting2014} over several simulations. We also calculated the relative root mean square prediction error (RRMSPE), defined as the root mean square prediction error of HV and LR, respectively, divided by the root mean square prediction error of DL. For both criteria, lower values are better.

\subsection{Spatial-only data}
\label{sec:purely-spatial}

In our first scenario, we considered spatial-only data according to \eqref{eq:ngobs}--\eqref{eq:ngevol} on a grid $\mathcal{S}$ of size $n=34 \times 34 = 1{,}156$ on the unit square, $\mathcal{D} = [0,1]^2$. We set $\bm{\mu} = \mathbf{0}$ and $\bm{\Sigma}_{i,j} = \exp(-\|\mathbf{s}_i - \mathbf{s}_j\|/0.15)$. For the Gaussian likelihood, we assumed variance $\tau^2 = 0.2$.

The comparison scores averaged over 100 simulations for the posteriors obtained using Algorithm \ref{alg:VL} are shown as a function of $N$ in Figure \ref{fig:spatial-scores}. HV (Algorithm \ref{alg:VF}) was much more accurate than LR for each value of $N$.

\begin{figure}[ht]
    \centering
    \begin{subfigure}{1.0\textwidth}
        \includegraphics[trim=0mm 6mm 0mm 0mm, clip,width=1.0\textwidth]{plots/spatial-RRMSPE.pdf}

    \end{subfigure}
    \medskip
    \begin{subfigure}{1.0\textwidth}
        \includegraphics[trim=0mm 0mm 0mm 6mm, clip,width=1.0\textwidth]{plots/spatial-dLS.pdf}

    \end{subfigure}
    \caption{Approximation accuracy for the posterior distribution $\mathbf{x}|\mathbf{y}$ for spatial data (see Section \ref{sec:purely-spatial})}
    \label{fig:spatial-scores}
\end{figure}

\subsection{Linear temporal evolution\label{sec:simul-linear}}

Next, we considered a linear spatio-temporal advection-diffusion process with diffusion parameter $\alpha = 4 \times 10^{-5}$ and advection parameter $\beta = 10^{-2}$ as in \citet{Jurek2018}. The spatial domain $\mathcal{D} = [0,1]^2$ was discretized on a grid of size $n=34 \times 34 = 1{,}156$ using the centered finite differences, and we considered discrete time points $t=1,\ldots,T$ with $T=20$. After this discretization, our model was of the form \eqref{eq:linear-dynamics-obs}--\eqref{eq:linear-dynamics}, where $\bm{\Sigma}_{0|0}=\mathbf{Q}_1=\ldots=\mathbf{Q}_T$ with $(i,j)$th entry $\exp(-\|\mathbf{s}_i - \mathbf{s}_j\|/0.15)$, and $\mathbf{E}_t$ was a sparse matrix with nonzero entries corresponding to interactions between neighboring grid points to the right, left, top and bottom. See the supplementary material of \citet{Jurek2018} for details.
At each time $t$, we generated $n_t = 0.1 n$ observations with indices $\mathcal{I}_t$ sampled randomly from $\{1,\ldots,n\}$. For the Gaussian case, we assumed variance $\tau^2=0.25$.
We used conditioning sets of size at most $N=41$ for both HV and LR; specifically, for HV, we used $J=2$ partitions at $M=7$ resolutions, with set sizes $|\mathcal{X}_{j_1,\ldots,j_m}|$ of $5, 5, 5, 5, 6, 6, 6$, respectively, for $m=0,1,\ldots,M-1$, and $|\mathcal{X}_{j_1,\ldots,j_M}| \leq 3$.

Figure \ref{fig:linear-scores} compares the scores for the filtering distributions $\mathbf{x}_t | \mathbf{y}_{1:t}$ obtained using Algorithm \ref{alg:KVLfilter}, averaged over 80 simulations. Again, HV was much more accurate than LR. Importantly, while the accuracy of HV was relatively stable over time, LR became less accurate over time, with the approximation error accumulating.

\begin{figure}[ht]
    \centering
    \begin{subfigure}{1.0\textwidth}
        \includegraphics[trim=0mm 6mm 0mm 0mm, clip,width=1.0\textwidth]{plots/linear-RRMSPE.pdf}

    \end{subfigure}
    \begin{subfigure}{1.0\textwidth}
        \includegraphics[trim=0mm 0mm 0mm 6mm, clip,width=1.0\textwidth]{plots/linear-dLS.pdf}

    \end{subfigure}
    \caption{Accuracy of filtering distributions $\mathbf{x}_t | \mathbf{y}_{1:t}$ for the advection-diffusion model in Section \ref{sec:simul-linear}}
    \label{fig:linear-scores}
\end{figure}

\subsection{Simulations using a very large \texorpdfstring{$n$}{n}} \label{sec:big-simul}

We repeated the advection-diffusion experiment from Section \ref{sec:simul-linear} on a high-resolution grid of size $n=300 \times 300 = 90{,}000$, with $n_t=9{,}000$ observations corresponding to 10\% of the grid points. In order to avoid numerical artifacts related to the finite differencing scheme, we reduced the advection and diffusion coefficients to $\alpha = 10^{-7}$ and $\beta=10^{-3}$, respectively. We set $N=44$, $M=14$, $J=2$, and $|\mathcal{X}_{j_1,\ldots,j_m}| = 3$ for $m=0,1,\ldots,M-1$, and $|\mathcal{X}_{j_1,\ldots,j_M}| \leq 2$. DL was too computationally expensive due to the high dimension $n$, and so we simply compared HV and LR based on the root mean square prediction error (RMSPE) between the true state and their respective filtering means, averaged over 10 simulations.

As shown in Figure \ref{fig:linear-scores-big}, HV was again much more accurate than LR. Comparing to Figure \ref{fig:linear-scores}, we see that the relative improvement of HV to LR increased even further; taking the Gaussian case as an example, the ratio of the RMSPE for HV and LR was around 1.2 in the small-$n$ setting, and greater than 2 in the large-$n$ setting.

\begin{figure}[ht]
    \centering
    \begin{subfigure}{1.0\textwidth}
        \includegraphics[width=1.0\textwidth]{plots/linear-RMSPE.pdf}
    \end{subfigure}
    \caption{Root mean square prediction error (RMSPE) for the filtering mean in the high-dimensional advection-diffusion model with $n=90{,}000$ in Section \ref{sec:big-simul}}
    \label{fig:linear-scores-big}
\end{figure}

\subsection{Nonlinear evolution with non-Gaussian data\label{sec:lorenz}}

Our final set of our simulations involves the most general model \eqref{eq:obs}-\eqref{eq:evol} with nonlinear evolution $\mathcal{E}_t$. Specifically, we considered a complicated model \citep[][Sect.~3]{Lorenz2005} that realistically replicates many features of atmospheric variables along a latitudinal band. A special case of this model \citep{Lorenz1996} is an important benchmark for data-assimilation techniques.
The model dynamics are described by
\begin{equation}
\textstyle    \frac{\partial }{\partial t} \tilde x_i = \frac{1}{K^2} \sum_{l=-K/2}^{K/2}\sum_{j=-K/2}^{K/2} -\tilde x_{i-2K-l}\tilde x_{i-K-j} + \tilde x_{i-K+j-l}\tilde x_{i+K+j} - \tilde x_i + F,
    \label{eq:Lorenz04-1}
\end{equation}
where
$\tilde x_{-i} = \tilde x_{n-i}$, and we used $K=32$ and $F=10$.
By solving \eqref{eq:Lorenz04-1} on a regular grid of size $n=960$ on a circle with unit circumference using a 4-th order Runge-Kutta scheme, with five internal steps of size $dt=0.005$, and setting $x_i = b\tilde x_i$ with $b=0.2$, we obtained the evolution operator $\mathcal{E}_t$. We also calculated an analytic expression for its derivative $\nabla\mathcal{E}_t$, which is necessary for Algorithm \ref{alg:eKVLfilter}.

To complete our state-space model \eqref{eq:obs}-\eqref{eq:evol}, we assumed $\mathbf{Q}_{t;i,j} = 0.2 \exp(-\|\mathbf{s}_i - \mathbf{s}_j\|/0.15)$, we randomly selected $n_t = n/10 = 96$ observation indices at each $t$, and we took the initial moments $\bm{\mu}_{0|0}$ and $\bm{\Sigma}_{0|0}$ to be the corresponding sample moments from a long simulation from \eqref{eq:Lorenz04-1}. We simulated 40 datasets from the state-space model, each at $T=20$ time steps, and for each of the four exponential-family likelihoods, using $\tau^2=0.2$ in the Gaussian case.

For each data set obtained in this way, we applied the three filtering methods described in Section \ref{sec:simul-general}. We used $N=39$, and for the EKVL filter (Algorithm \ref{alg:eKVLfilter}) we set $J=2$, $M=7$, and $\lvert \mathcal{X}_{j_1,\ldots,j_m} \rvert$ equal to 5, 5, 5, 5, 6, 6, 6, respectively, for $m=0, 1, \dots, M-1$, and $\lvert \mathcal{X}_{j_1,\ldots,j_M}\rvert \leq 1$. The average scores over the 40 simulations are shown in Figure \ref{fig:lorenz-scores}.
Our method (HV) compared favorably to the low-rank filter and provided excellent approximation accuracy as evidenced by very low RRMSPE and dLS scores.

\begin{figure}[ht]
    \centering
    \begin{subfigure}{1.0\textwidth}
        \includegraphics[trim=0mm 6mm 0mm 0mm, clip,width=1.0\textwidth]{plots/lorenz-RRMSPE.pdf}

    \end{subfigure}
    \begin{subfigure}{1.0\textwidth}
        \includegraphics[trim=0mm 0mm 0mm 12mm, clip,width=1.0\textwidth]{plots/lorenz-dLS.pdf}

    \end{subfigure}
    \caption{Accuracy of filtering distribution $\mathbf{x}_t | \mathbf{y}_{1:t}$ for the Lorenz model in Section \ref{sec:lorenz}}
    \label{fig:lorenz-scores}
\end{figure}

\section{Conclusions\label{sec:conclusion}}

We specified the relationship between ordered conditional independence and sparse (inverse) Cholesky factors. Next, we described a hierarchical Vecchia approximation and showed that it exhibits equivalent sparsity in the Cholesky factors of both its precision and covariance matrices. Due to this remarkable property, the approximation is suitable for high-dimensional spatio-temporal filtering.
The hierarchical Vecchia approximation can be computed using a simple and fast incomplete Cholesky decomposition. Further, by combining the approach with a Laplace approximation and the extended Kalman filter, we obtained scalable filters for non-Gaussian and non-linear spatio-temporal state-space models.

Our methods can also be directly applied to spatio-temporal point patterns modeled using log-Gaussian Cox processes, which can be viewed as Poisson data after discretization of the spatial domain, resulting in accurate Vecchia-Laplace-type approximations \citep{Zilber2019}.
We plan on investigating an extension of our methods to retrospective smoothing over a fixed time period. Another interesting extension would be to combine our methodology with the unscented Kalman filter \citep[][]{Julier1997} for strongly nonlinear evolution.
Finally, while we focused our attention on spatio-temporal data, our work can be extended to other applications, as long as a sensible hierarchical partitioning of the state vector can be obtained as in Section \ref{sec:hv-intro}.

Code implementing our methods and reproducing our numerical experiments is available at \url{https://github.com/katzfuss-group/vecchiaFilter}.

\FloatBarrier
\footnotesize