text:
\section{Local approximation MCMC}
\label{sec:numerics}

We assume that our target distribution has a density $\pi(x) = \exp{\mathcal{L}(x)}$ on $\mathbb{R}^d$. Our goal is to construct an algorithm that exploits regularity in $\pi$ to reduce the computational cost of simulating from the target distribution. As discussed in the Introduction, replacing target density evaluations with a continually and infinitely refined surrogate model, within MCMC, can asymptotically generate samples from the exact target distribution \citep{Conradetal2016}. Given a finite number of MCMC samples, however, replacing the target density with such an approximation introduces a \textit{surrogate bias}. In this section, we present a LA-MCMC algorithm that extends \citet{Conradetal2016} in two significant ways: (i) we use the trade-off between surrogate bias and Monte Carlo variance to develop a refinement strategy with near-optimal convergence guarantees; and (ii) we introduce a parameter that allows LA-MCMC to easily characterize distributions with heavier tails given little \textit{a priori} knowledge, while enabling additional speedups when the user has substantial \textit{a priori} knowledge.

Our approach will build a \textit{local polynomial} approximation \citep{Connetal2009,Kohler2002,Stone1977} of a function $g:\mathbb{R}^{d} \to \mathbb{R}$, where $g$ is chosen such that evaluating $\pi(x)$ is trivial given $g(x)$. For instance, $g$ might be the logarithm of the target density, $\mathcal{L}$; or it might be the log-likelihood function in a Bayesian setting, if the prior density is relatively simple to evaluate. Let the ``evaluated set'' $\mathcal{S}_n = \{x_1, \hdots, x_n\}$ comprise the set of parameter values $x_i \in \mathbb{R}^d$ at which we have evaluated $g(x_i)$. For now, we take this set as given; later, we discuss how to construct it. With an appropriate construction of $\mathcal{S}$, local approximations allow refinement to focus on regions where the target distribution has greater mass. Polynomial approximations are appealing in this setting because they can be built easily, are cheap to evaluate, and have known analytic derivatives.

\subsection{Local polynomial approximations}

We construct local polynomial approximations using weighted regression \citep{Kohler2002}. Let $\mathcal{P}$ be a polynomial space over $\mathbb{R}^{d}$. The approximation at each $x \in \mathbb{R}^{d}$ is $\hat{g}(x, \mathcal{S}_n)$ such that
\begin{equation}
    \hat{g}(x, \mathcal{S}_n) = \argmin_{m \in \mathcal{P}}{\sum_{i=1}^{n} (m(x_i)-g(x_i))^2 W(x, x_i)},
    \label{eq:local-polynomial-estimate}
\end{equation}
where $W(x, x^{\prime})$ is a locally supported kernel. Typically, we choose the $k$-nearest neighbor kernel
\begin{equation}
    W(x, x^{\prime}) = \begin{cases}
    1 & \mbox{if $x^{\prime} \in \mathcal{B}_k(x)$} \\
    0 & \mbox{otherwise,}
    \end{cases}
    \label{eq:hat-kernel}
\end{equation}
where $\mathcal{B}_k(x)$ is the smallest ball centered at $x$ containing $k$ elements of $\mathcal{S}_n$. Define a basis $\Phi = \{\phi_j\}_{j=1}^{q}$ for $\mathcal{P}$, let $\boldsymbol{\phi}(x) = [\phi_1(x), \hdots, \phi_q(x)]^T$, and define $\boldsymbol{x}_k(x, \mathcal{S}_n) = \{\tilde{x} \in \mathcal{S}_n: W(x, \tilde{x})>0\}$. The elements of $\boldsymbol{x}_k(x, \mathcal{S}_n)$ are the $k$ nearest neighbors to $x$ in $\mathcal{S}_n$. The optimal polynomial kernel estimate is
\begin{equation}
    \hat{g}(x, \mathcal{S}_n) = \boldsymbol{\phi}(x)^T \mathbf{a}(x, \boldsymbol{x}_k(x, \mathcal{S}_n))
\end{equation}
such that
\begin{equation}
    \mathbf{a}(x, \boldsymbol{x}_k) = \argmin_{\boldsymbol{\alpha} \in \mathbb{R}^{q}}{\|\mathbf{V}(\boldsymbol{x}_k) \boldsymbol{\alpha} - \mathbf{g}(\boldsymbol{x}_k)\|_{\mathbf{W}(x, \boldsymbol{x}_k)}^2},
    \label{eq:local-polynomial-estimate-matrix-form}
\end{equation}
where $\|\mathbf{x}\|_{\mathbf{A}}^2 = \mathbf{x}^T \mathbf{A} \mathbf{x}$,
\begin{subequations}
\begin{equation}
    \mathbf{W}(x, \boldsymbol{x}_k) = \diag{[W(x, x_1), \hdots, W(x, x_k)]},
\end{equation}
$\mathbf{V}(\boldsymbol{x}_k)$ is the Vandermonde matrix
\begin{equation}
    \mathbf{V}(\boldsymbol{x}_k) = \left[ \begin{array}{c}
    \boldsymbol{\phi}(x_1)^T \\
    \vdots \\
    \boldsymbol{\phi}(x_k)^T
    \end{array} \right],
\end{equation}
and
\begin{equation}
    \mathbf{g}(\boldsymbol{x}_k) = \left[ \begin{array}{c}
     g(x_1) \\
     \vdots \\
     g(x_k)
     \end{array} \right].
\end{equation}
\end{subequations}
Assuming $\mathbf{V}(\boldsymbol{x}_k)$ has full column rank and $k \geq q$, the local polynomial approximation exists and is unique \citep{Kohler2002,Stone1977}. The solution to \eqref{eq:local-polynomial-estimate-matrix-form} is
\begin{equation}
    \mathbf{a}(x, \boldsymbol{x}_k(x, \mathcal{S}_n)) = (\mathbf{V}^T \mathbf{W} \mathbf{V})^{-1} \mathbf{V}^T \mathbf{g}.
\end{equation}

\subsubsection{Error analysis}
\label{sec:erroranalysis}

Now we derive local error bounds that, for a fixed evaluated set $\mathcal{S}_n$, depend on the number of nearest neighbors, the size of the ball containing them, and the local behavior of the surrogate. We assume that the kernel $W(x, x^{\prime})$ is the hat kernel defined in \eqref{eq:hat-kernel} and that $\mathcal{P}$ is the space of polynomials of degree less than or equal to $p$;

we write $q = \text{dim}( \mathcal{P})$.

We first investigate the local behavior of the polynomials in the ball $\mathcal{B}_k(x)$. Let
\begin{equation*}
\Delta(x) = \max_{x_i \in\boldsymbol{x}_k(x,\mathcal{S}_n)}{\|x-x_i\|}
\end{equation*}
be the radius of $\mathcal{B}_k(x)$---the ball that contains $k$ points in the evaluated set $\mathcal{S}_n$. We note that the radius of $\mathcal{B}_k(x)$ is not uniform in $x \in \mathbb{R}^{d}$. Analogous to \eqref{eq:local-polynomial-estimate-matrix-form}, define $k$ Lagrange polynomials
\begin{equation}
    \lambda_j(x) = \argmin_{m \in \mathcal{P}}{\sum_{i=1}^{k} (m(x_i)-\delta_{ij})^2},
\end{equation}
where $\delta_{ij}$ is the Dirac delta. Here, the summation is over the $k$ nearest neighbors to $x$, i.e., $x_i \in \boldsymbol{x}_k(x)$. Let $\boldsymbol{\lambda}(x) = [\lambda_1(x), \hdots, \lambda_k(x)]^T$. The nearest neighbors $\boldsymbol{x}_k(x)$ are called $\Lambda$-poised in $\mathcal{B}_k(x)$ if
\begin{equation}
    \max_{x^{\prime} \in \mathcal{B}_k(x)}{\|\boldsymbol{\lambda}(x^{\prime})\|_2} = \Lambda_2(x) \leq \Lambda.
    \label{eq:poisedness}
\end{equation}
This definition of $\Lambda$-poisedness is slightly different than in \citet{Connetal2009}, which uses the infinity norm. Since the number of nearest neighbors $k$ is fixed and finite, however,

these definitions are equivalent. $\Lambda$-poisedness measures how well distributed the points are within the ball. For example, consider points $i$ and $j$ in a unit ball and quadratic polynomials. The $i^{\text{th}}$ Lagrange polynomial is close to one at point $i$ and close to zero at point $j$. (If the number of nearest neighbors $k$ is exactly the number required to interpolate ($k=q$) then the Lagrange polynomial will be exactly one or zero at these points.)

The quadratic Lagrange polynomial is narrower if the points are close together and wider if they are far apart. The wider parabola has a smaller poisedness constant than the narrower one.

 Assuming $\Lambda$-poisedness and that $g(x)$ is at least $(p+1)$ times differentiable, we have the error bound

\begin{align}
    \vert \hat{g}(x, \mathcal{S}_n) - g(x) \vert \leq \frac{k}{(p+1)!} \Lambda \Delta(x)^{p+1} \sup_{x^{\prime} \in \mathcal{B}_k(x)}{g^{(p+1)}(x^{\prime})},
    \label{eq:local-error-bound}
\end{align}

where $p$ is the polynomial degree and $g^{(p+1)}(x)$ is the $(p+1)^{\text{th}}$ derivative of $g(x)$ \citep{Connetal2009}. The radius $\Delta(x)$ decreases as points are added to $\mathcal{S}_n$ and thus the error bound decreases, assuming that points are chosen in a way that maintains $\Lambda$-poisedness.

\subsubsection{Local refinements}

Given a surrogate $\hat{g}(x, \mathcal{S}_n)$, we can improve its accuracy in a neighborhood around $x$ by performing a \textit{local refinement}. A local refinement adds a new point $x^{*}$ to the evaluated set---i.e., $\mathcal{S}_{n+1} = \mathcal{S}_n \cup \{x^{*}\}$---and thus allows the local error bound \eqref{eq:local-error-bound} to decrease. Randomly choosing a nearby point tends to form clusters \citep{RoteTichy1996}, however, and therefore refining using a random point inside $\mathcal{B}_k(x)$ fails to maintain $\Lambda$-poisedness.

We instead choose the refinement location based on the poisedness constant $\Lambda_2(x)$ defined in \eqref{eq:poisedness}. Computing $\Lambda_2(x)$ by solving the optimization problem defined in \eqref{eq:poisedness} also defines the point
\begin{equation}
    x_{\lambda}(x) = \argmin_{x^{\prime} \in \mathcal{B}_k(x)}{\|\boldsymbol{\lambda}(x^{\prime})\|_2}.
    \label{eq:poisedness-refinement-location}
\end{equation}
If $\boldsymbol{x}_k(x)$ is poorly poised in $\mathcal{B}_k(x)$ ($\Lambda_2(x) \gg 1$) then the new point $x_{\lambda}(x)$ will be relatively far from any clusters in $\boldsymbol{x}_k(x)$ and will tend to improve the poisedness of the updated set. Conversely, if $\boldsymbol{x}_k(x)$ is well poised in $\mathcal{B}_k(x)$ ($\Lambda_2(x) \approx 1$), then the updated set remains well poised. We therefore refine the surrogate by setting $\mathcal{S}_{n+1} = \mathcal{S}_n \cup \{x_{\lambda}(x)\}$.

Having discussed {how} to refine, we must also discuss \textit{when}, i.e., under what conditions, to refine. Our primary criterion will be derived from the local error bound \eqref{eq:local-error-bound}; this process is described in Section~\ref{sec:bias-variance-trade-off}. In the meantime, we discuss an additional but natural secondary criterion, which is to trigger a refinement if the poisedness constant exceeds a user-prescribed threshold $\bar{\Lambda}$ (i.e., if $\Lambda_2(x) > \bar{\Lambda}$ then set $\mathcal{S}_{n+1} = \mathcal{S}_n \cup \{x_{\lambda}(x)\}$). In practice, we find that adding $x_{\lambda}(x)$ to the evaluated set whenever a refinement is triggered tends to maintain poisedness automatically, and that the secondary criterion thus rarely triggers refinements. Computing the poisedness constant by solving \eqref{eq:poisedness} can also be computationally burdensome (although still significantly cheaper than an expensive density evaluation). We therefore often ``turn off'' this secondary refinement criterion by setting $\bar{\Lambda} = \infty$, and only compute the poisedness constant when refinements are triggered by the primary criterion.

\subsection{Sampling methods using local approximations}

Each step of a typical MCMC algorithm consists of three stages: (i) propose a new state, (ii) compute the acceptance probability, and (iii) accept or reject the proposed state. LA-MCMC adds a fourth stage---possibly refine the surrogate model---and replaces the target density evaluations in stage (ii) exclusively with cheaper surrogate evaluations; see Algorithm \ref{alg:la-mcmc}. The refinement frequency ensures that the error incurred by using a surrogate model balances the Monte Carlo variance.

\subsubsection{Bias-variance trade-off} \label{sec:bias-variance-trade-off}

We now present a heuristic for balancing bias and variance in LA-MCMC,  deferring a rigorous discussion to Section \ref{sec:theory}.

Fix a function $f$ and Markov chain $\{X_{t}\}_{t \geq 0}$ with stationary measure $\pi$. For $T \in \mathbb{N}$, denote by $\hat{\pi}_{T}(f) = \frac{1}{T} \sum_{t=1}^{T} f(X_{t})$ the usual Monte Carlo estimate of $\pi(f) = \mathbb{E}[f(X)]$---the expectation with respect to $\pi$.

Under modest conditions (see, e.g.,  \citet{meyn2012markov}), the asymptotic bias exists and satisfies

\begin{equation}
\lim_{T \rightarrow \infty} T \, \left \vert \mathbb{E}[\hat{\pi}_{T}(f)] - \pi(f) \right \vert = 0

\end{equation}
while the asymptotic variance exists and satisfies
\begin{equation}
\lim_{T \rightarrow \infty} T \, \var{(\hat{\pi}_{T}(f))}  = C_v \in (0,\infty).
\end{equation}

When this happens, there is (asymptotically) no trade-off: the bias quickly becomes a negligible source of error.

In the context of LA-MCMC,  more care is required because of an additional bias from the surrogate model. We assume that $\pi$ has density $\pi(x) = \exp{ \mathcal{L}(x) }$, and view LA-MCMC as an ``approximate'' MCMC algorithm that tries to target\footnote{We note that our $\widehat{\pi}(x) = \exp \widehat{\mathcal{L}}(x)$ often fails to be a probability density. We resolve this technical problem in Section \ref{sec:theory}, and somewhat surprisingly our approach means that this problem has very little impact on the following heuristic calculations.} $\widehat{\pi}(x) = \exp \widehat{\mathcal{L}}(x)$, where $\widehat{\mathcal{L}}(x) \approx \mathcal{L}(x)$. In particular, we assume that the surrogate $\widehat{\mathcal{L}}$ adheres to the error bound in \eqref{eq:local-error-bound}, which is a pointwise condition of the form:
\begin{equation}
 \vert \widehat{\mathcal{L}}(x) - {\mathcal{L}}(x) \vert \leq C_{b} \Delta(x)^{p+1}
 \label{eq:transition-kernel-radius-bound}
\end{equation}
for $C_b>0$. We will show in Section \ref{sec:theory} that the bias incurred from using this surrogate model inherits a similar error bound,
\begin{equation}
    \left| \hat{\pi}_T(f) - \pi(f) \right| \leq C_{b}^{\prime} \bar{\Delta}^{p+1},
    \label{eq:expectation-error}
\end{equation}
where $C_b^{\prime}>0$ and, informally, we can think of $\bar{\Delta}^{p+1}$ as the maximum of $\Delta(x)^{p+1}/V(x)$, where $V(x) > 0$ is a penalty function to be described below. We will thus use the local radius $\Delta(x)^{p+1}$ as an \textit{error indicator}, and sequentially refine our approximation so that $\bar{\Delta} = \bar{\Delta}(T)$ decays with time $T$.

We would like to balance the additional error introduced by the surrogate model against the usual ``Monte Carlo'' error of the original Markov chain. As is typically the case when balancing different sources of error, we choose our parameters so that the upper bound on surrogate error decays at the same rate as the Monte Carlo error.

To ensure that the pointwise surrogate error is small enough, we partition the chain into

consecutive intervals called \emph{levels}, $\ell = 0, 1, 2, \ldots$, and prescribe a piecewise constant error threshold for each level $\ell$,
\begin{equation}
    \gamma_\ell(x) = \gamma_0 \ell^{-\gamma_1} V(x),
    \label{eq:error-threshold}
\end{equation}
with $\gamma_0, \gamma_1 > 0$.

Here, $V(x) > 0$ is a $\pi$-integrable penalty function that allows the error threshold to be larger in low-probability regions (e.g., in the tails of the distribution). We assume that $V(x)$ satisfies the Lyapunov inequality \eqref{eq:lyapunov}; this requirement guides how to choose the function. As a reasonable default, we set $V(x) = \exp{(\nu_0 \|x-\bar{x}\|^{\nu_1})}$, where $\nu_0>0$, $0 < \nu_1 \leq 1$, and $\bar{x}$ is a user-prescribed estimate of the centroid of the distribution (e.g., the mean or mode). During MCMC, we locally refine the surrogate model if the error indicator $\Delta(x)^{p+1}$ at the current state exceeds the current threshold $\gamma_\ell(x)$.

Level $\ell$ ends at step $T_{\ell}$, when the MCMC variance is of the same size as the squared surrogate bias:

\begin{equation}
    \left ( \gamma_0 \ell^{-\gamma_1} \right )^2 \sim C_{v} T^{-1}_\ell.
    \label{eq:bias-variance}
\end{equation}

The number of steps in the $\ell^{\text{th}}$ level is, therefore,
\begin{equation}
    T_{\ell} = \tau_0 \ell^{2\gamma_1},
\end{equation}
where $\tau_0 = C_v / \gamma_0^2$ is the length of the first level. In practice, we set the the level as a function of the MCMC step $t$,

\begin{equation}
    \ell(t) = \lfloor (t/\tau_0)^{1/(2\gamma_1)} \rfloor.
    \label{eq:leveldefn}
\end{equation}
The length of each level, $\tau_{\ell} = T_{\ell}-T_{\ell-1}$, will then decrease with $\ell$ if $\gamma_1 < 0.5$. This is undesirable because then, beyond some critical step $t^* < \infty$, the level length will become less than one. If this occurs, we will need to increment the level more than once per MCMC step. Rather than addressing this complication, we simply require that lengths of the levels strictly increase with $\ell$ by setting $\gamma_1 > 0.5$.

The constant $\tau_0$ is rarely known, and we treat it as an algorithmic parameter. Given the initial error threshold $\gamma_0>0$, error decay rate $\gamma_1>0.5$, initial level length $\tau_0 \geq 1$, and the Lyapunov function $V(x)$, we locally refine whenever the error indicator at an accepted state in the chain exceeds the decaying error threshold. See Algorithm \ref{alg:la-mcmc} for a complete summary.

\subsubsection{Tail correction}
\label{sec:tailcorrection}

There is a major technical difficulty in making the heuristic from Section \ref{sec:bias-variance-trade-off} precise: it is possible that our surrogate function has small pointwise error as in \eqref{eq:transition-kernel-radius-bound}, but fails to be globally integrable (and thus cannot possibly be an unnormalized probability density). This results in large practical and theoretical difficulties: the associated stochastic process may wander off to infinity, and the pointwise bound in \eqref{eq:transition-kernel-radius-bound} will not give any bound on the Monte Carlo error. In previous work (e.g., \citet{Conradetal2016}), we avoided this problem by working on compact state spaces or obtaining much stronger pointwise approximations than \eqref{eq:transition-kernel-radius-bound}. Here, we obtain much stronger results by slightly tweaking our algorithm's acceptance probability.

Our tweak uses the Lyapunov function of the MCMC algorithm with exact evaluations, $V(x)$, to change the acceptance probability. Rather than computing the typical Metropolis-Hastings acceptance probability targeting $\hat{\pi}$, we slightly \textit{increase} the chances of moves that would decrease the Lyapunov function and \textit{decrease} the chances of moves that would increase it. \citet{roberts1996geometric} give very general conditions under which we can easily define the Lyapunov function.

More precisely, given the current state $x_t$ and proposed state $x^{\prime}_t$, we tweak the usual Metropolis-Hastings acceptance probability with the estimate
\begin{equation}
    \log{\pi(x^{\prime}_t)} \approx \widetilde{\mathcal{L}}(x^{\prime}_t) \equiv \widehat{\mathcal{L}}(x^{\prime}_t) + Q_V(x_t, x^{\prime}_t),
    \label{eq:Lypunov-correction}
\end{equation}
where

\begin{equation}
    Q_V(x_t, x^{\prime}_t) = \begin{cases}
    \ \, \, \eta (\gamma(x^{\prime}_t) + \gamma(x_t))   & \mbox{if } V(x'_t) < V(x_t)\\
    - \eta(\gamma(x^{\prime}_t) + \gamma(x_t)) \ & \mbox{if } V(x'_t) \geq V(x_t),
    \end{cases}
    \label{EqQvDef}
\end{equation}
and $\gamma(\cdot) = \gamma_{\ell(t)}(\cdot)$ is defined by \eqref{eq:error-threshold} and \eqref{eq:leveldefn}. Here, $\widehat{\mathcal{L}}(x^{\prime}_t)$ is an approximation of the log-density using the local polynomial approximation. The factor $\gamma(x^{\prime}_t) + \gamma(x_t)$ causes the correction to become less apparent as the error threshold decays,

$\gamma_{\ell(t)}(\cdot) \to 0$,
and $\eta \geq 0$ is a user-defined parameter to control the correction.

\begin{algorithm}
 \begin{algorithmic}
    \State Set initial state $X_0$ and surrogate model $\widehat{\mathcal{L}}_0$ for the log-target density
    \For{$t \gets 1$ to $\infty$}
    \State Propose $X^{\prime} \sim q_t(\cdot \vert X_t)$
    \State
    \State Possibly refine at $X_t$: $(\gamma(X_t), \widehat{\mathcal{L}}_t) = $ \Call{CheckAndRefine}{$X_t$, $\widehat{\mathcal{L}}_{t-1}$, $t$}
    \State
    \State Compute error threshold at the proposed point $$\gamma(X^{\prime}) = \gamma_0 (\lfloor (t/\tau_0)^{1/(2\gamma_1)} \rfloor)^{-\gamma_1} V(X^{\prime})$$

    \State
    \State Compute $\widetilde{\mathcal{L}}_t(X^{\prime})$ via \eqref{eq:Lypunov-correction} and evaluate the acceptance probability
    \begin{equation*}
        \alpha(X_t, X^{\prime}) = \min{\left(1, \exp{(\widetilde{\mathcal{L}}_t(X^{\prime})-\widehat{\mathcal{L}}_t(X_t))} \frac{q_t(X_t \vert X^{\prime})}{q_t(X^{\prime} \vert X_t)} \right)}
    \end{equation*}
    \State
    \State Accept/reject step:
    \begin{equation*}
        X_{t+1} = \begin{cases}
        X^{\prime} & \mbox{with probability } \alpha(X_t, X^{\prime}) \\
        X_t & \mbox{else}
    \end{cases}
    \end{equation*}
    \EndFor
    \Procedure{CheckAndRefine}{$x$, $\widehat{\mathcal{L}}$, $t$}
  \State Compute the local error threshold $$\gamma(x) = \gamma_0 (\lfloor (t/\tau_0)^{1/(2\gamma_1)} \rfloor)^{-\gamma_1} V(x)$$ and (optionally) compute the poisedness constant $\Lambda_2(x)$
  \State
  \If{$\Delta(x)^{p+1}>\gamma(x)$ or (optionally) $\Lambda_2(x) > \bar{\Lambda}$}
    \State $\widehat{\mathcal{L}}^{*}$ = \Call{RefineSurrogate}{$x$, $\widehat{\mathcal{L}}$}
  \Else{}
    $\widehat{\mathcal{L}}^{*} = \widehat{\mathcal{L}}$
  \EndIf
  \Return The error threshold $\gamma(x)$ and $\widehat{\mathcal{L}}^{*}$
  \EndProcedure
  \State
  \Procedure{RefineSurrogate}{$x$, $\widehat{\mathcal{L}}$}
  \State Compute poisedness-based refinement location $x_{\lambda}(x)$ defined in \eqref{eq:poisedness-refinement-location}
  \If{$x_{\lambda}(x) \notin \mathcal{S}_n$}
    Set $\mathcal{S}_n \gets \mathcal{S}_n \cup \{x_{\lambda}(x)\}$ and $n \gets n+1$
  \EndIf
  \If{$x_{\lambda}(x) \in \mathcal{S}_n$}
    Randomly choose $X^{*} \in \mathcal{B}_k(x)$ and set $\mathcal{S}_n \gets \mathcal{S}_n \cup \{X^{*}\}$ and $n \gets n+1$
  \EndIf
  \Return The updated surrogate model $\widehat{\mathcal{L}}^{*}$ using $\mathcal{S}_n$
  \EndProcedure
 \end{algorithmic}
 \caption{Pseudocode for the LA-MCMC algorithm. LA-MCMC requires eight parameters: (i) the initial error threshold $\gamma_0>0$, (ii) the error threshold decay rate $\gamma_1>0.5$, (iii) the maximum poisedness constant $\bar{\Lambda}$, (iv) the length of the first level $\tau_0 \geq 1$, (v) the tail-correction parameter $\eta \geq 0$, (vi) the number of nearest neighbors $k$ used to construct the local polynomial surrogate, (vii) the degree of the local polynomial surrogate $p$, and (viii) a guessed Lyapunov function $V$.}
 \label{alg:la-mcmc}
\end{algorithm}


TITLE: Show that $ \int_{0}^{\infty} \frac{t}{e^{t}-1} \, \mathrm{d}t = \frac{\pi^{2}}{6} $
QUESTION [0 upvotes]: I would like to prove that
$$\int_{0}^{\infty} \frac{t}{e^{t}-1} \, \mathrm{d}t = \frac{\pi^{2}}{6} $$
My teacher suggested to use
$$ \int_{0}^{\infty} \frac{t}{e^{t}}\sum_{n\geq 0 }e^{nt} \, \mathrm{d}t $$
But I don't understand why both integrals would be equal.

REPLY [2 votes]: Preliminaries:
$\frac {1}{1-x} = \sum_\limits{n=0}^{\infty} x^n\\
\frac {1}{1-e^x} = \sum_\limits{n=0}^{\infty} (e^{x})^n = \sum_\limits{n=0}^{\infty} e^{nx}$
And for our problem:
$\frac {t}{{e^t}-1} = \frac {te^{-t}}{1-e^{-t}}  = te^{-t}\sum_\limits{n=0}^{\infty}  e^{-nx} = t\sum_\limits{n=1}^{\infty}  e^{-nx}$
Since integration is a "linear operator" we can integrate each term, and sum the results.
that is: $\int (\sum f(t)) dt = \sum (\int f(t) dt)$
$\int_0^\infty t e^{-nt} dt = \frac {1}{n^2}$
$\sum \frac {1}{n^2} = \frac {\pi^2}{6}$


TITLE: Eigenvalues of Sturm–Liouville operator
QUESTION [4 upvotes]: Can we calculate the eigenvalues and eigenfunctions of the following operator in $W^{1,2}(\mathbb{R})$?
$$-\left(\frac{1}{\cosh^2x}\right)y''-\frac{2}{\cosh^4x}y=\lambda y.$$

REPLY [7 votes]: For me, the first (and usually only) attempt with such Sturm-Liouville equations is to rewrite them as Schrodinger equations, using this transformation.
If I got everything right, ... EDIT: but I didn't in the first attempt, as hyh pointed out. The correct transformed equation is
$$
-u'' + V(t)u = \lambda u , \quad\quad V(t) =\frac{-t^2-6}{4(t^2+1)^2} ;
$$
the variables are related to the old ones by $t=\sinh x$, $u=\sqrt{\cosh x} y$.
This is a similar, but slightly more subtle situation than with my original incorrect potential. We now have $V=-1/(4t^2)-1/t^4+O(t^{-6})$. What is clear is that we have purely absolutely continuous spectrum of multiplicity $2$ on $\sigma_{ac}=[0,\infty)$. Below $\lambda =0$, we have purely discrete spectrum.
The question of whether we have finitely or infinitely many eigenvalues below zero is probably answerable, but seems a bit tricky. The potential $W=-1/(4t^2)$ would lead to finite discrete spectrum, but is exactly the borderline case in the sense that $W=-c/t^2$ for $c>1/4$ would give infinitely many eigenvalues.
EDIT 2: I'm reasonably confident now that there are only finitely many eigenvalues, though to show it properly would probably require some work. We could try to argue as follows: on the interval $2^n\le t\le 2^{n+1}$, we can bound $V$ from below by $-c_n/t^2$, $c_n=1/4+2^{-2n}$. For this modified potential, the solutions at $\lambda =0$ look like $t^{1/2}t^{i2^{-n}}=t^{1/2}e^{i2^{-n}\log t}$ (the equation for the modified potential is an Euler equation and can be solved explicitly). So the phase increases by $\approx n2^{-n}$ on the interval under consideration. Since this is summable, we expect (real valued) solutions to have only finitely many zeros, and then the claim follows from oscillation theory (see the comments).
Third and (let's hope) final EDIT: There are indeed only finitely many eigenvalues. A much better argument (than the previous paragraph) is given by Giorgio in the comments.


TITLE: Show that $f(z)=\sum_{k=1}^{n}|f_{k}(z)|$ cannot achieve a strictly local maximum inside a connected open set $D$.
QUESTION [1 upvotes]: Suppose that we have $f_{1},\cdots, f_{n}$ holomorphic in a connected open set $D$. For $$f(z):=|f_{1}(z)|+|f_{2}(z)|+\cdot+|f_{n}(z)|\ \ \text{for}\ \ z\in D.$$ I want to show that $f(z)$ cannot achieve a strict local maximum inside $D$.

I know that this must have something to do with the maximum modulus principle, especially we have some connected open set, this will enable us to directly conclude that a function $g(z)$ is constant all over $D$ if it is holomorphic in $D$ and $|g|$ achieves a local maximum at $z_{0}\in D$.
However, I cannot connect this theorem to $f(z)$ because even thought $f_{k}(z)$'s are holomorphic, the summation of moduli of them may not be holomorphic.
I know a way to show that

Let $D$ be a disc with radius $R$ and center $z_{0}$. Suppose $f_{k}(z)$'s are holomorphic on $D$ and continuous on $\partial D$ such that $f(z_{0})\geq z$ for all $z\in D\cup C$, then $f$ is constant in $D\cup C$.

I can prove this claim using Cauchy's Integral formula, but it is long. This lemma can certainly conclude the proof quickly, but I realized that I did not really use the "connectedness" of $D$.

Is there any way for me to drive a shorter proof using maximum modulus principle? Thank you!

Edit 1:
The post is duplicated, but I'd still like to post a relatively complete proof.
Thanks to Conrad, the following is a relatively complete proof:
Suppose $\omega$ is the local maximum of $f(z):=\sum_{k=1}^{n}|f_{k}(z)|$. Note that as $f_{k}(\omega)$ is a complex number, unless it is zero, we can represent (uniquely) as $$f_{k}(\omega)=|f_{k}(\omega)|e^{it_{k}},\ \ t_{k}\in [0,2\pi).$$ Take $\alpha_{k}:=e^{-itk}$, and if $f_{k}(\omega)=0$ for some $k$, we tae $\alpha_{k}:=1$.
Then, we have $$f_{k}(\omega)=\alpha_{k}^{-1}|f_{k}(\omega)|\ \text{with}\ |\alpha_{k}|=1.$$
Now, define $g(z):=\sum_{k=1}^{n}\alpha_{k}f_{k}(z).$ Then, we have \begin{align*}
|g(z)|\leq \sum_{k=1}^{n}|\alpha_{k}||f_{k}(z)|=\sum_{k=1}^{n}|f_{k}(z)|&=f(z)\\
&\leq \sum_{k=1}^{n}|f_{k}(\omega)|=f(\omega)\\
&=\sum_{k=1}^{n}f_{k}(\omega)\alpha_{k}=g(\omega)\leq|g(\omega)|\ \ \ (*).
\end{align*}
Hence, $\omega$ is also the point at which $|g|$ achieves local max. As $g$ is holomorphic and $D$ is connected open, by the global maximum modulus principle, we conclude that $g(z)=C$ for some constant $C$ and for all $z\in D.$
But then, $(*)$ implies that $$|C|\leq C\leq |C|\implies C=|C|\implies C\geq 0.$$
So we have $g(z)=g(\omega)\geq 0$. But then (*) implies that $$g(\omega)=g(z)=|g(z)|\leq f(z)\leq f(\omega)\leq g(\omega),$$ and thus we must have $$g(\omega)=g(z)=f(z)=f(\omega)=C\geq 0\ \ \text{for all}\ \ z\in D,$$ and thus the local maximum of $f$ cannot be strict.

I really appreciate the patient and inspiring help from Conrad.

REPLY [0 votes]: Thanks to Conrad's patient and inspiring help, I wrote the proof, as follows:
Suppose $\omega$ is the point at which local maximum of $f(z):=\sum_{k=1}^{n}|f_{k}(z)|$ is achieved. Note that as $f_{k}(\omega)$ is a complex number, unless it is zero, we can represent (uniquely) as $$f_{k}(\omega)=|f_{k}(\omega)|e^{it_{k}},\ \ t_{k}\in [0,2\pi).$$ Take $\alpha_{k}:=e^{-itk}$, and if $f_{k}(\omega)=0$ for some $k$, we tae $\alpha_{k}:=1$.
Then, we have $$f_{k}(\omega)=\alpha_{k}^{-1}|f_{k}(\omega)|\ \text{with}\ |\alpha_{k}|=1.$$
Now, define $g(z):=\sum_{k=1}^{n}\alpha_{k}f_{k}(z).$ Then, we have \begin{align*}
|g(z)|\leq \sum_{k=1}^{n}|\alpha_{k}||f_{k}(z)|=\sum_{k=1}^{n}|f_{k}(z)|&=f(z)\\
&\leq \sum_{k=1}^{n}|f_{k}(\omega)|=f(\omega)\\
&=\sum_{k=1}^{n}f_{k}(\omega)\alpha_{k}=g(\omega)\leq|g(\omega)|\ \ \ (*).
\end{align*}
Hence, $\omega$ is also the point at which $|g|$ achieves local max. As $g$ is holomorphic and $D$ is connected open, by the global maximum modulus principle, we conclude that $g(z)=C$ for some constant $C$ and for all $z\in D.$
But then, $(*)$ implies that $$|C|\leq C\leq |C|\implies C=|C|\implies C\geq 0.$$
So we have $g(z)=g(\omega)\geq 0$. But then (*) implies that $$g(\omega)=g(z)=|g(z)|\leq f(z)\leq f(\omega)\leq g(\omega),$$ and thus we must have $$g(\omega)=g(z)=f(z)=f(\omega)=C\geq 0\ \ \text{for all}\ \ z\in D,$$ and thus the local maximum of $f$ cannot be strict.


TITLE: Can a normal graduate discover something in mathematics?
QUESTION [4 upvotes]: When I was in High school, I really developed my interest in mathematics. From that time , I always had an ambition to become mathematician or to discover something in mathematics. Because of some reasons, I chose B.E. in electronics & communication and got a degree in that. Now I got a government job also. My work is some what less. In free time, I usually watch mathematics related videos about unsolved mathematical mysteries and sometimes I try to work on that. My favourite part is complex numbers. My most favourite thing is infinite number series problems like grandis series etc.,
From the beginning I always think that there must be something unified in mathematics like sometimes I feel like prime numbers , Fibonacci series , Pascal triangle , e, π, i are all related in some way. One day I watched a video about Riemann hypothesis which isn't solved yet. But solving it might give us better understanding of prime numbers.
So here I come to my question part, instead of wasting my time everyday, I am thinking to achieve  something. Mathematics is really that kind of item to work on it. I mean it just needs pen & paper. Also, I can get programmes to simulate my thoughts / ideas about numbers. Especially infinite series etc.,
Internet is there for me to get all information. so my question is which branch of mathematics is best for me to continue to work on that. Number theory?. I have decided yesterday that I should dedicate myself  to mathematics in next 20 years. Now my age is 24. I hope I'll really get something positive result in future. Even if I fail also, I don't mind. I told this to my friend & he laughed saying "till now no one discovered pattern of prime numbers & you think you can discover something in mathematics, even PhD holder isn't getting results".
 sometimes crazy ideas come to my mind related to maths, & immediately I'll start working on that. After hours and hours ,I'll realise that this isn't simple thing to do.
Can I really succeed in this?
What is the best way to approach?
Should I first start reading recent discoveries in maths & start working on that?
Now I can't goto university & do PhD .
With this internet, will I be able to do a reaserch in mathematics? Or else should I stop thinking like this non sense ?
Can a person like me be able to discover?

REPLY [2 votes]: Welcome to mathematics. You can have a lot of fun playing in this world.
I think you'll have more fun if you start with easier problems than the Riemann hypothesis and the patterns of the primes.
You should count any problem you solve for yourself as "research", even if you're not the first one to solve it. Perhaps begin with easy ones - these searches suggests many places:
https://www.google.com/search?q=elementary+problems+in+mathematics&ie=utf-8&oe=utf-8
https://www.google.com/search?q=recreational++problems+in+mathematics&ie=utf-8&oe=utf-8
You may find that when you solve a problem you have an idea about how to change it slightly to make a new problem you will be the first to solve. Find a place to post it on the internet as a puzzle for others and you have begun publishing your own contributions to mathematics.


TITLE: If $\sum_n a_n $ converges and $a_n > 0$, is it true that $\sum_n \sqrt{a_n}$ is also convergent?
QUESTION [1 upvotes]: Assume that $\sum_n a_n $ is convergent and $a_n > 0$, is it true that $\sum_n \sqrt{a_n}$ is also convergent? What can you comment on $∑_n a_n^2$?

I have thought that since summation an is convergent that $\lim_{n \to ∞} a_n = 0$ so $\lim_{n \to ∞} \sqrt{a_n}$ must also be $0$, hence it is convergent.
We can't really say about $a_n^2$.
Is this right? Can anyone suggest a different approach to this? Thanks.

REPLY [2 votes]: First hint: consider
$$
a_n=\frac{1}{n^2}
$$
Second hint: if $0<a<1$, then $a^2<a$. If $\sum_n a_n$ is convergent, then $0<a_n<1$ for $n>n_0$, for some $n_0$; thus the inequality $a_n^2<a^n$ holds for $n>n_0$. What can you say about a series with positive terms that's bounded above by a convergent series (from some point on)?

Comment: the condition that a sequence $(a_n)$ converges to zero doesn't imply the series $\sum_n a_n$ converges (which is your error in the first case).


text:
\begin{document}

\title{Combining the band-limited parameterization and Semi-Lagrangian Runge-Kutta integration for efficient PDE-constrained LDDMM}

\titlerunning{BL parameterization and SL-RK for PDE-LDDMM}

\author{Monica Hernandez}

\institute{Computer Sciences Department \\ Aragon Institute on Engineering Research \\ University of Zaragoza \\ mhg@unizar.es}

\date{Received: This work was finished in January 2019  / Accepted: date}

\maketitle

\begin{abstract}
The family of PDE-constrained LDDMM methods is emerging as a particularly interesting approach for physically meaningful diffeomorphic transformations.
The original combination of Gauss--Newton--Krylov optimization
and Runge-Kutta integration,
shows excellent numerical accuracy and fast convergence rate.
However, its most significant limitation is the huge computational complexity,
hindering its extensive use in Computational Anatomy applied studies.
This limitation has been treated independently by the problem formulation in the space of band-limited vector fields and Semi-Lagrangian integration.
The purpose of this work is to combine both in three variants of band-limited PDE-constrained LDDMM
for further increasing their computational efficiency.
The accuracy of the resulting methods is evaluated extensively.
For all the variants, the proposed combined approach shows a significant increment of the computational efficiency.
In addition, the variant based on the deformation state equation is positioned consistently as the best performing method
across all the evaluation frameworks in terms of accuracy and efficiency.
This work was finished on January 2019 and is currently under review.

\end{abstract}

\newpage

\section{Introduction}

Computational Anatomy is a powerful interdisciplinary field for the analysis of anatomical shape variability~\cite{Miller_04,Miller_09}.
This discipline is based on Sir D'Arcy Thompson's original ideas for explaining the similarity of the anatomical shape
of homologous species using the transformations existing between the anatomical structures~\cite{Thompson_17}.
In Computational Anatomy, shape similarity is measured from the diffeomorphic transformations estimated between the anatomies.
These transformations yield a generative model for the analysis of shape variability.
Diffeomorphisms are computed from the anatomical images using diffeomorphic registration methods~\cite{Younes_Book}.

There exists a vast literature on diffeomorphic registration methods with differences in the transformation characterization,
regularizers, image similarity metrics, optimization methods, and additional constraints~\cite{Sotiras_13}.
Although the differentiability and invertibility of the transformations constitute crucial features for Computational Anatomy applications,
the diffeomorphic constraint does not necessarily guarantee that a transformation computed with a given method is physically meaningful for the clinical domain of interest.
PDE-constrained Large Deformation Diffeomorphic Metric Mapping (PDE-LDDMM) has become relevant in the last decade
for the computation of transformations under plausible physical models of interest~\cite{Younes_07,Hart_09,Ashburner_11,Vialard_11,Mang_15,Mang_16,Mang_17,Zhang_18,Hernandez_19b}.

Our work focuses on the family of PDE-constrained LDDMM methods pioneered by Hart et al.~\cite{Hart_09} and leading to the relevant contributions
in~\cite{Mang_15,Mang_16,Polzin_16,Hernandez_19c,Mang_19}.
In this family of methods, the registration problem is approached from an optimal control perspective, where the different physical models are imposed
directly using the physical PDEs that are attached to the LDDMM variational
problem using hard constraints.
The numerical optimization is approached using gradient-descent~\cite{Hart_09,Mang_15,Polzin_16} or second-order optimization in the form
of inexact reduced Newton-Krylov methods~\cite{Mang_15,Mang_16,Hernandez_19c,Mang_19}.
In particular, the combination of Gauss--Newton--Krylov for optimization, with sophisticated multi-level preconditioners,
spectral methods for differentiation, and Runge-Kutta schemes for PDE integration,
shows excellent numerical accuracy and an extraordinarily fast convergence rate.
However, the most significant limitation of Gauss--Newton--Krylov PDE-constrained LDDMM is the huge computational complexity, which hinders
the extensive use in Computational Anatomy applied studies.
This computational complexity is due to:
\begin{enumerate}
\item The formulation of the problem in the spatial domain.
\item The large time-sampling needed for the stability of Runge-Kutta integration.
\end{enumerate}

\noindent Both issues have been treated independently in the literature yielding to PDE-constrained LDDMM methods with increased
efficiency and an assumable cost in accuracy loss.

\subsection{Computational complexity due to problem formulation}

The computational complexity due to the formulation of the problem in the spatial domain has been successfully reduced using the band-limited vector field
parameterization proposed in~\cite{Zhang_15,Zhang_18}.
LDDMM methods, and in particular PDE-constrained LDDMM, involve the action of low-pass filters in the optimization update equations of the velocities.
Therefore, the computation of the high-frequency components of high-resolution velocity fields can be omitted since these computations result equal or nearly
equal to zero by the action of the low-pass filters.
The band-limited vector field parameterization allows a reduction of the dimensionality of the problem that circumvents the high-frequency computations.

The works in~\cite{Hernandez_19,Hernandez_19c} formulate three different variants of PDE-constrained LDDMM in the space of band-limited vector fields
and perform the computations in the GPU.
Some configurations of these variants have been really successful, greatly outperforming the state of the art methods in terms of computational complexity
while keeping a competitive accuracy.

\subsection{Computational complexity due to PDE integration}

Runge-Kutta methods are explicit techniques. Hence they are only conditionally stable.
This means that the time sampling should be selected enough to preserve the Courant-Friedrichs-Lewy (CFL) condition.
For PDE-constrained LDDMM, the time sampling values that guarantee stability are usually large.
As a result, the time and memory requirements of the problem are considerably increased.
In particular, the memory requirements of PDE-constrained LDDMM is increased to limits that hinder the execution on limited memory
devices such as the GPU.
In addition, one can experience that the time-sampling needed for the non-stationary parameterization is much higher than for the
stationary parameterization, increasing the complexity of an already not particularly memory efficient configuration.
On the other side, when stability is satisfied, the accuracy of PDE-constrained LDDMM is high~\cite{Mang_15,Mang_16,Hernandez_19,Hernandez_19c}.

Semi-Lagrangian methods are semi-implicit techniques that are unconditionally stable.
Therefore, the time sampling can be selected according to accuracy rather than stability considerations.
Semi-Lagrangian methods were originally proposed in the 90's in the context of modelling weather predictions~\cite{Staniforth_91}.
In the context of diffeomorphic registration, the original LDDMM method proposed in~\cite{Beg_05} already used Semi-Lagrangian integration
for the solution of the transport equation.
The combination of Semi-Lagrangian integration with Runge-Kutta has been recently proposed for solving some time-dependent PDEs.
Runge-Kutta has shown to increase the accuracy of first-order schemes in Semi-Lagrangian integration~\cite{Guo_13}.

The computational complexity in~\cite{Mang_15,Mang_16} due to the use of Runge-Kutta schemes for PDE integration has been successfully reduced
using Semi-Lagrangian Runge-Kutta integration~\cite{Mang_17b,Mang_19} for the stationary parameterization of diffeomorphisms.
For PDE-constrained LDDMM, the selected time sampling is usually much smaller than the time sampling typically selected with explicit
schemes, yielding to a considerable reduction of the computational complexity of the problem.
On the other hand, the expected accuracy of PDE-constrained LDDMM is lower than with explicit schemes.

Beyond the computational complexity improvement through numerical schemes, Mang et al. proposed an efficient implementation of PDE-LDDMM that
exploits massive CPU based parallel computing architectures~\cite{Mang_16_3D}. The source code has been recently released with~\cite{Mang_19}.

\subsection{Our contribution}

The purpose of this work is to further increase the computational efficiency of BL PDE-constrained LDDMM by combining the two independent methodological
approaches of circumventing the huge computational complexity of PDE-constrained LDDMM and to extensively analyze the accuracy of the resulting methods.
We have implemented the band-limited methods in~\cite{Hernandez_19,Hernandez_19c} with the Semi-Lagrangian Runge-Kutta integration scheme originally
proposed in~\cite{Mang_17b} for the stationary and the non-stationary parameterization of diffeomorphisms.
The resulting methods have been evaluated in five different datasets following the evaluation frameworks in~\cite{Hernandez_19c,Klein_09,Rohlfing_12}.
To our knowledge, this is the first time that Semi-Lagrangian Runge-Kutta integration is implemented in the space of band-limited vector fields.
It is also the first time that Semi-Lagrangian Runge-Kutta integration is used in PDE-LDDMM with the non-stationary parameterization.
Moreover, our work first provides the position achieved by benchmark PDE-constrained LDDMM methods in the ranking of Klein et al. evaluation.
The best performing method of our work coincides with the best performing variant in~\cite{Hernandez_19c}, PDE-constrained LDDMM based on the deformation state equation.
The Semi-Lagrangian Runge-Kutta scheme proposed in this work has shown to outperform the Runge-Kutta scheme in~\cite{Hernandez_19c} in terms of computational efficiency
and accuracy.
Indeed, the best performing PDE-LDDMM variant in this work has recently reached the highest sensitivity (97 \% vs a baseline of 88 \%) in the classification
of stable vs progressive mild cognitive impaired conversors in the Alzheimer's Disease Neuroimaging Initiative (ADNI) database using convolutional neural networks~\cite{Ramon_20}.

\subsection{Manuscript organization}

In the following, Section~\ref{sec:PDE-LDDMM} reviews the foundations of PDE-constrained LDDMM,
with particular emphasis on the band-limited vector field parameterization.
Section~\ref{sec:SL} presents the proposed Semi-Lagrangian Runge-Kutta integration method.
Next, Section~\ref{sec:ExpSetup} details the experimental setup.
Section~\ref{sec:Results} shows the results and Section~\ref{sec:Discussion} discusses the most important highlights.
Finally, Section~\ref{sec:Conclusions} gathers the most remarkable conclusions of our work.

\section{PDE-constrained LDDMM methods}
\label{sec:PDE-LDDMM}

\subsection{Parameterization in the spatial domain}

Let $\Omega \subseteq \mathbb{R}^d$ be the image domain.
Let $Diff(\Omega)$ be the LDDMM Riemannian manifold of diffeomorphisms and $V$ the tangent space at the identity element.
$Diff(\Omega)$ is a Lie group, and $V$ is the corresponding Lie algebra.
The Riemannian metric of $Diff(\Omega)$ is defined from the scalar product in $V$
\begin{equation}
 \langle v, w \rangle_V = \langle \Cauchy v, w \rangle_{L^2} = \int_\Omega \langle \Cauchy v(x), w(x) \rangle d\Omega,
\end{equation}
where $\Cauchy = (Id - \alpha \Delta)^s, \alpha >0, s \in \mathbb{N}$
is the invertible self-adjoint differential operator associated with the differential structure of
$Diff(\Omega)$~\cite{Beg_05}.

Let $I_0$ and $I_1$ be the source and the target images.
PDE-constrained LDDMM is formulated from the minimization of the PDE-constrained variational problem

\begin{equation}
\label{eq:LS-LDDMM}
E(v) = \frac{1}{2} \int_0^1 \langle \Cauchy v_t, v_t \rangle_{L^2} dt + \frac{1}{\sigma^2} \Vert m(1) - I_1\Vert_{L^2}^2,
\end{equation}
subject to the state equation with state variable $m(t)$
\begin{equation}
\label{eq:StateEquation}
\partial_t m(t) + \nabla m(t) \cdot v_t = 0 \textnormal{ in } \Omega \times (0,1],
\end{equation}
\noindent with initial condition $m(0) = I_0$~\cite{Hart_09,Mang_15}.

The variational problem is posed in the space of time-varying smooth flows of velocity fields, ${v} \in L^2([0,1],V)$.
Given the smooth flow ${v}:[0,1] \rightarrow V$, $v_t:\Omega \rightarrow \mathbb{R}^{d}$,
the solution at time $t=1$ to the evolution equation
\begin{equation}
\label{eq:TransportEquation}
\partial_t \phi_t^{v} = v_t \circ \phi_t^{v}
\end{equation}
with initial condition $\phi_0^{v} = id$ is a diffeomorphism, $\phi_1^{v} \in Diff(\Omega)$.
The transformation $(\phi^{v}_{1})^{-1}$, computed from the minimum of $E({v})$ constrained to the state equation,
is the diffeomorphism that solves the PDE-constrained LDDMM registration problem between $I_0$ and $I_1$.

The state equation constraint in PDE-constrained LDDMM can be imposed in two more different manners, yielding three
different variants of PDE-constrained LDDMM~\cite{Hernandez_19c}. The second variant is formulated from the
minimization of Equation~\ref{eq:LS-LDDMM}, where
\begin{equation}
 m(t) = I_0 \circ \phi(t)
\end{equation}
\noindent and $\phi$ is computed from the deformation state equation
\begin{equation}
\label{eq:DefStateEquation}
\partial_t \phi(t) + D \phi(t) \cdot v_t = 0 \textnormal{ in } \Omega \times (0,1].
\end{equation}
\noindent The third variant is formulated from the minimization of Equation~\ref{eq:LS-LDDMM} subject to the
deformation state equation (Equation~\ref{eq:DefStateEquation}).

The best optimization method for PDE-constrained LDDMM is Gauss--Newton--Krylov~\cite{Mang_15,Hernandez_19c}.
The expressions of the gradient and the Hessian-vector product are derived from the augmented Lagrangian of the
energy functional and the state equation.
The update equation has the form

\begin{equation}
\label{eq:Newton}
 v^{n+1} = v^n + \epsilon \delta v^n,
\end{equation}
where $\epsilon$ is the update length and $\delta v^n$ is computed from preconditioned conjugate gradient (PCG) on the system
\begin{equation}
\label{eq:KKT}
 H_v E_{aug}( v^n ) \delta  v^n = - \nabla_v E_{aug}(v^n),
\end{equation}
with preconditioner $\Cauchy^{-1}$.

\subsection{Parameterization in the space of band-limited vector fields}

Let $\widetilde{\Omega}$ be the discrete Fourier domain truncated with frequency bounds $K_1,$ $\dots,$ $K_d$.
We denote with $\widetilde{V}$ the space of discretized band-limited vector fields on $\Omega$ with these frequency bounds.
The elements in $\widetilde{V}$ are represented in the Fourier domain as $\tilde{v}: \widetilde{\Omega} \rightarrow \mathbb{C}^d$,
$\tilde{v}(k_1, \dots, k_d)$.

\noindent The application $\iota:\widetilde{V} \rightarrow V$ denotes the natural inclusion mapping of $\widetilde{V}$ in $V$.
The application $\pi: V \rightarrow \widetilde{V}$ denotes the projection of $V$ onto $\widetilde{V}$~\cite{Zhang_15,Zhang_18}.

The space $\widetilde{V}$ of band-limited vector fields has a finite-dimensional Lie algebra structure using the truncated
convolution $\star$ in the definition of the Lie bracket~\cite{Zhang_18}.
We denote with $Diff(\widetilde{\Omega})$ to the finite-dimensional Riemannian manifold of diffeomorphisms on $\widetilde{\Omega}$
with corresponding Lie algebra $\widetilde{V}$.
The Riemannian metric in $Diff(\widetilde{\Omega})$ is defined from the scalar product
\begin{equation}
 \langle \tilde{v}, \tilde{w} \rangle_{\tilde{V}} = \langle \tilde{\Cauchy} \tilde{v}, \tilde{w} \rangle_{l^2},
\end{equation}
\noindent where $\tilde{\Cauchy}$ is the projection of operator $\Cauchy$ in the truncated Fourier domain.
Similarly, we will denote with $\tilde{\ast}$ to the projection in the truncated Fourier domain of the differential operators $\ast$
involved in the differential equations.

The band-limited PDE-constrained variational problem is given by the minimization of
\begin{equation}
\label{eq:BLEnergy}
 E(\tilde{v}) = \frac{1}{2} \int_0^1 \langle \tilde{L}\tilde{v}_t, \tilde{v}_t \rangle_{l^2} dt + \frac{1}{\sigma^2} \Vert m(1) - I_1 \Vert_{L^2}^2.
\end{equation}
The first variant is formulated from the minimization of Equation~\ref{eq:BLEnergy} subject to
\begin{equation}
\partial_t m(t) + \nabla m(t) \cdot \iota(\tilde{v}_t) = 0 \textnormal{ in } \Omega \times (0,1], \\
\end{equation}
\noindent with initial condition $m(0) = I_0$.
For the second variant, the diffeomorphism is computed from $\phi(t) = id - \iota(\tilde{u})(t)$ where
$\tilde{u}(t)$ is computed from the deformation state equation formulated in displacement field form
\begin{equation}
\label{eq:BLDefStateEquation}
\partial_t \tilde{u}(t) + \widetilde{D} \tilde{u}(t) \star \tilde{v}_t = \tilde{u}(t) \textnormal{ in } \widetilde{\Omega} \times (0,1].
\end{equation}
The third variant is formulated analogously to the spatial case from the minimization of Equation~\ref{eq:BLEnergy} subject to the
deformation state equation (Equation~\ref{eq:BLDefStateEquation})~\cite{Hernandez_19,Hernandez_19c}.

The optimization is approached using Gauss--Newton--Krylov methods in $\widetilde{V}$ with preconditioner $\tilde{\Cauchy}^{-1}$.
The update equation has the form
\begin{equation}
\label{eq:BLNewton}
 \tilde{v}^{n+1} = \tilde{v}^n + \epsilon \delta \tilde{v}^n,
\end{equation}
\noindent where $\delta \tilde{v}^n$ is computed from
\begin{equation}
\label{eq:BLKKT}
 \widetilde{(H_{\tilde{v}} E_{aug}(  \tilde{v}^n))} \delta \tilde{v}^n = - \widetilde{(\nabla_{\tilde{v}} E_{aug}(\tilde{v}^n))}.
\end{equation}

\subsection{BL PDE-constrained LDDMM equations}
\label{sec:OPDE-LDDMM}

\subsubsection{Original BL PDE-constrained LDDMM}

Originally proposed BL PDE-constrained LDDMM uses the state equation in the augmented Lagrangian for the derivation
of the state and adjoint equations and their incremental counterparts~\cite{Hernandez_19}.
The gradient and the Gauss-Newton approximation of the Hessian-vector product are given by the equations
\begin{eqnarray}
\widetilde{(\nabla_{\tilde{v}} E_{aug}(\tilde{v}))_t} = \tilde{\Momentum} \tilde{v}_t + \tilde{\lambda}(t) \star \widetilde{\nabla} \tilde{m}(t) \\
\widetilde{ (H_{\tilde{v}} E_{aug}(\tilde{v}))_t } \delta \tilde{v} (t) = \tilde{\Momentum} \delta \tilde{v}(t) + \delta \tilde{\lambda}(t) \star \widetilde{\nabla} \tilde{m}(t),
\end{eqnarray}
\noindent where the projected state variable $\tilde{m}$ and the projected adjoint variable $\tilde{\lambda}$ are computed from $\pi(m)$ and $\pi(\lambda)$, and $m$ and $\lambda$ are computed from
\begin{eqnarray}
\label{eq:BLState}
\partial_t m(t) + \nabla m(t) \cdot \iota(\tilde{v}_t) = 0 \\
\label{eq:BLAdjoint}
-\partial_t \lambda(t) - \nabla \cdot ( \lambda(t) \cdot \iota(\tilde{v}_t) ) = 0.
\end{eqnarray}
\noindent The incremental counterparts $\delta \tilde{m}$ and $\delta \tilde{\lambda}$ are the solutions of
\begin{eqnarray}
\label{eq:BLIncState}
\partial_t \delta \tilde{m} (t) + \widetilde{\nabla} \delta \tilde{m}(t) \star \tilde{v}_t + \widetilde{\nabla} \tilde{m}(t) \star \delta \tilde{v}(t) = 0\\
\label{eq:BLIncAdjoint}
-\partial_t \delta \tilde{\lambda}(t) - \widetilde{\nabla} \cdot ( \delta \tilde{\lambda}(t) \star \tilde{v}_t ) = 0
\end{eqnarray}
\noindent in the BL domain. The initial conditions are, respectively, $m(0) = I_0$, $\lambda(1) = -\frac{2}{\sigma^2}(m(1)-I_1)$, $\delta \tilde{m}(0) = 0$,
$\delta \tilde{\lambda}(1) = -\frac{2}{\sigma^2} \delta \tilde{m}(1)$.

\subsubsection{BL PDE-constrained LDDMM based on the state equation}

This method departs from the original BL PDE-constrained LDDMM by using $m(t) = I_0 \circ \phi(t)$ and $\lambda(t) = J(t) \lambda(1) \circ \psi(t)$,
where $\psi$ is the inverse of $\phi$, and $J$ is the Jacobian determinant of $\psi$~\cite{Hart_09,Hernandez_19c}.
The transformations $\phi$ and $\psi$ and the scalar field $J$ are computed from the inclusion of the truncated displacement fields ($\tilde{u}(t)$ and $\tilde{\nu}(t)$)
and the corresponding Jacobian
\begin{eqnarray}
 \phi(t) = id - \iota(\tilde{u}(t)), \\
 \psi(t) = id - \iota(\tilde{\nu}(t)), \\
 J(t) = 1 - \iota(\tilde{U}(t)),
\end{eqnarray}
\noindent where
\begin{eqnarray}
\label{eq:BLDefState}
\partial_t \tilde{u}(t) + \tilde{D} \tilde{u}(t) \star \tilde{v}_t = \tilde{v}_t \\
-\partial_t \tilde{\nu}(t) - \tilde{D} \tilde{\nu}(t) \star \tilde{v}_t = -\tilde{\nu}_t \\
-\partial_t \tilde{U}(t) - \tilde{v}_t \star \widetilde{\nabla} \tilde{U}(t) = -\widetilde{\nabla \cdot} \tilde{v} + \tilde{U}(t) \star \widetilde{\nabla \cdot} \tilde{v}(t)
\end{eqnarray}

\noindent with initial conditions $\tilde{u}(0)= 0$, $\tilde{\nu}(1) = 0$, and $\tilde{U}(1) = 0$.

The incremental state and adjoint variables are computed from
\begin{eqnarray}
\delta m(t) = \nabla I_0 \circ \phi(t) \cdot \iota(\delta \tilde{u}(t))\\
\delta \lambda(t) = J(t) \nabla \lambda(1) \circ \psi(t) \cdot \iota(\delta \tilde{\nu}(t))
\end{eqnarray}
where $\delta \tilde{u}$ and $\delta \tilde{\nu}$ are computed from
\begin{eqnarray}
\label{eq:BLIncDefState}
\partial_t \delta \tilde{u}(t) + \tilde{D} \delta \tilde{u}(t) \star \tilde{v}_t + \tilde{D} \tilde{u}(t) \star \delta \tilde{v}_t = \delta \tilde{v}_t \\
-\partial_t \delta \tilde{\nu}(t) - \tilde{D} \delta \tilde{\nu}(t) \star \tilde{v}_t - \tilde{D} \tilde{\nu}(t) \star \delta \tilde{v}_t = -\delta \tilde{\nu}_t
\end{eqnarray}

\noindent with initial conditions $\delta \tilde{u}(0) = 0$ and $\delta \tilde{\nu}(1) = 0$.

\subsubsection{BL PDE-constrained LDDMM based on the deformation state equation}

The third method is formulated from the minimization of Equation~\ref{eq:BLEnergy}
subject to the truncated displacement state equation~\cite{Hernandez_19c}
\begin{equation}
\partial_t \tilde{u}(t) + \tilde{D} \tilde{u}(t) \star \tilde{v}_t = \tilde{v}_t.
\end{equation}

The gradient and the Hessian-vector product are given by the equations
\begin{eqnarray}
\widetilde{(\nabla_{\tilde{v}} E_{aug}(\tilde{v}))_t} = \tilde{\Momentum} \tilde{v}_t + \tilde{\rho}(t) - \widetilde{D} \tilde{u}(t) \star \tilde{\rho}(t) \\
\widetilde{ (H_{\tilde{v}} E_{aug}(\tilde{v}))_t } \delta \tilde{v} (t) = \tilde{\Momentum} \delta \tilde{v}(t) + \delta \tilde{\rho}(t) - \widetilde{D} \delta \tilde{u}(t) \star \tilde{\rho}(t),
\end{eqnarray}

\noindent where the displacement state variable $\tilde{u}$, the adjoint variable $\tilde{\rho}$, and their incremental counterparts $\delta \tilde{u}$ and $\delta \tilde{\rho}$ are computed from
\begin{eqnarray}
\partial_t \tilde{u}(t) + \widetilde{D} \tilde{u}(t) \star \tilde{v}_t = \tilde{v}_t \\
\label{eq:BLDefAdjoint}
-\partial_t \tilde{\rho}(t) - \widetilde{\nabla \cdot}  (\tilde{\rho}(t) \star \tilde{v}_t ) = 0  \\
\partial_t \delta \tilde{u}(t) + \widetilde{D} \delta \tilde{u}(t) \star \tilde{v}_t + \widetilde{D} \tilde{u}(t) \star \delta \tilde{v}(t) = \delta \tilde{v}(t) \\
\label{eq:BLIncDefAdjoint}
-\partial_t \delta \tilde{\rho}(t) - \widetilde{\nabla \cdot} ( \delta \tilde{\rho}(t) \star \tilde{v}_t ) = 0
\end{eqnarray}
\noindent with initial conditions $\tilde{u}(0) = 0$, $\tilde{\rho}(1)= \pi( -\frac{2}{\sigma^2} (m(1)-I_1) \nabla m(1))$, $\delta \tilde{u}(0) = 0$,
and $\delta \tilde{\rho}(1) = \pi( -\frac{2}{\sigma^2} \delta m(1) \nabla m(1))$.

\section{Semi-Lagrangian Runge-Kutta integration}
\label{sec:SL}

\subsection{Semi-Lagrangian integration in a spatial domain}

Semi-Lagrangian (SL) integration methods~\cite{Staniforth_91} allow solving transport equations of the general form

\begin{equation}
\label{eq:transport}
 D_t u = f(u,v),
\end{equation}

\noindent where $u: \Omega^d \times [0,1] \rightarrow \IR$ is a scalar or a vector function varying in time,
and $$D_t u = \partial_t u + D u \cdot v.$$
SL methods combine the most interesting properties of Eulerian and Lagrangian schemes.
On the one hand, SL methods involve following the characteristic lines of the differential equation,
similarly to Lagrangian approaches.
On the other hand, the equation is solved on the regular grid, similarly to Eulerian approaches.
As a result, SL methods are unconditionally stable as Lagrangian schemes.
This means that the time sampling can be selected according to accuracy considerations rather than
stability considerations.
SL methods allow selecting a time sampling usually much smaller than Eulerian methods yielding a sensible
reduction of the computational complexity.

SL schemes involve two steps.
First, the departure points are computed solving the characteristic equation

\begin{equation}
 D_t X(t) = v(t, X(t)),
\end{equation}

\noindent with initial condition $X(0) = x$.
The direction of the time integration can be forward or backward,
depending on the direction of the time integration of the transport equation.
From the several methods proposed in the literature for solving the characteristic equation,
we use the approach given by Mang et al. in~\cite{Mang_17b}
\begin{eqnarray}
 X_* = x - \delta t \cdot v \\
 v_* = v \circ X_* \\
 X_* = x - 0.5 \textnormal{ } \delta t \cdot (v_* + v).
\end{eqnarray}

Second, the transport equation (Equation~\ref{eq:transport}) is solved in the Eulerian grid
\begin{equation}
\label{eq:TransportOmega}
 D_t u(X(t), t) = f(u(X(t),t), v(t, X(t)))
\end{equation}

\noindent along the characteristic line $X$.
The use of Runge-Kutta (RK) integration has been recently proposed in this step, yielding a higher-order accurate SL-RK method~\cite{Guo_13}.
The velocity field needs to be estimated at points that do not belong to the Eulerian grid.
Therefore, an interpolator is needed. Cubic interpolation is the method of choice
for SL schemes~\cite{Riishojgaard_98}.

\subsection{Semi-Lagrangian integration in a band-limited domain}

In $\tilde{\Omega}$, the transport equations are of the general form
\begin{equation}
\label{eq:BLtransport}
 \tilde{D}_t \tilde{u} = \tilde{f}(\tilde{u},\tilde{v}),
\end{equation}
\noindent where
\begin{equation}
\tilde{D}_t \tilde{u} = \partial_t \tilde{u} + \tilde{D} \tilde{u} \star \tilde{v}.
\end{equation}

\noindent The characteristics are computed from
\begin{equation}
 D_t X(t) = \iota(\tilde{v})(t, X(t)),
\end{equation}
\noindent and the transport equation is solved from
\begin{equation}
\label{eq:TransportOmegaTilde}
 \tilde{D}_t \tilde{u}(\tilde{X}(t), t) = \tilde{f}(\tilde{u}(\tilde{X}(t),t), \tilde{v}(t, \tilde{X}(t))).
\end{equation}

\subsection{Semi-Lagrangian Runge-Kutta integration in PDE-LDDMM}

In this work, SL-RK integration has been implemented in $\Omega$ and $\tilde{\Omega}$ for the spatial and band-limited versions of the three
PDE-LDDMM variants. To be able to apply SL integration, the differential equations need to be written in the shape of equations
~\ref{eq:TransportOmega} or~\ref{eq:TransportOmegaTilde}. We focus on the derivation for the BL domain $\widetilde{\Omega}$.
The derivation for the spatial domain can be performed analogously.

The state equations, the deformation state equations, and their incremental counterparts (Equations \ref{eq:BLState}, \ref{eq:BLDefState},
\ref{eq:BLIncState}, \ref{eq:BLIncDefState}) are already in the shape of equation~\ref{eq:TransportOmegaTilde} by just moving to
the right-hand-side of the equation a remaining term. For the adjoint and the incremental adjoint equations (Equations \ref{eq:BLAdjoint},
\ref{eq:BLDefAdjoint}, \ref{eq:BLIncAdjoint}, \ref{eq:BLIncDefAdjoint}) we use the identity
\begin{eqnarray}

\widetilde{\nabla \cdot} (\tilde{u} \star \tilde{v}) = \tilde{u} \widetilde{\nabla \cdot} \tilde{v} + \tilde{v} \star \widetilde{\nabla} \tilde{u} \textnormal{ in } \tilde{\Omega}
\end{eqnarray}
\noindent and move the divergence term to the right-hand-side of the transformed equation.
Table~\ref{table:DtEquations} gathers the expressions of the resulting differential equations, needed for the implementation of BL PDE-constrained LDDMM
methods in SL form.
For SL-RK, the right-hand side expressions can be directly plugged into an RK differential solver.

\setlength{\tabcolsep}{10pt}
\renewcommand{\arraystretch}{1.5}

\begin{table*}[!t]
\scriptsize
\centering
\begin{tabular}{|c|c|}
\hline
Equation in original form & Equation in SL form\\
\hline
 $\partial_t m(t) + \nabla m(t) \cdot \iota(\tilde{v}_t) = 0$ (eq.~\ref{eq:BLState}) & $D_t m(t) = 0$ \\
 $-\partial_t \lambda(t) - \nabla \cdot ( \lambda(t) \cdot \iota(\tilde{v}_t) ) = 0$ (eq.~\ref{eq:BLAdjoint}) & $-D_t \lambda(t) = \lambda(t) \nabla \cdot \iota(\tilde{v_t})$ \\
 $\partial_t \tilde{u}(t) + \widetilde{D} \tilde{u}(t) \star \tilde{v}_t = \tilde{v}_t$ (eq.~\ref{eq:BLDefState}) & $\tilde{D}_t \tilde{u}(t) = \tilde{v_t}$ \\
 $ -\partial_t \tilde{\rho}(t) - \widetilde{\nabla \cdot}  (\tilde{\rho}(t) \star \tilde{v}_t ) = 0 $ (eq.~\ref{eq:BLDefAdjoint}) & $-\tilde{D}_t \tilde{\rho}(t) = \tilde{\rho}(t) \star \widetilde{\nabla \cdot} \tilde{v_t}$ \\
 $\partial_t \delta \tilde{m} (t) + \widetilde{\nabla} \delta \tilde{m}(t) \star \tilde{v}_t + \widetilde{\nabla} \tilde{m}(t) \star \delta \tilde{v}(t) = 0$ (eq.~\ref{eq:BLIncState}) & $\tilde{D}_t \delta \tilde{m}(t) = -\widetilde{\nabla} \tilde{m}(t) \star \delta \tilde{v_t}$ \\
 $-\partial_t \delta \tilde{\lambda}(t) - \widetilde{\nabla} \cdot ( \delta \tilde{\lambda}(t) \star \tilde{v}_t ) = 0$ (eq.~\ref{eq:BLIncAdjoint}) & $-D_t \delta \tilde{\lambda}(t) = \tilde{\lambda}(t) \star \widetilde{\nabla \cdot} \tilde{v_t}$ \\
 $\partial_t \delta \tilde{u}(t) + \widetilde{D} \delta \tilde{u}(t) \star \tilde{v}_t + \widetilde{D} \tilde{u}(t) \star \delta \tilde{v}(t) = \delta \tilde{v}(t)$ (eq.~\ref{eq:BLIncDefState}) &  $\tilde{D}_t \delta \tilde{u}(t) = \delta \tilde{v}_t - \tilde{D}\tilde{u}(t) \star \delta \tilde{v}(t)$ \\
 $-\partial_t \delta \tilde{\rho}(t) - \widetilde{\nabla \cdot} ( \delta \tilde{\rho}(t) \star \tilde{v}_t ) = 0$ (eq.~\ref{eq:BLIncDefAdjoint}) & $-\tilde{D}_t \delta \tilde{\rho}(t) = \delta \tilde{\rho}(t) \star \widetilde{\nabla \cdot} \tilde{v_t}$ \\
\hline
\end{tabular}
\caption{Original PDEs involved in BL PDE-constrained LDDMM and corresponding PDEs written in SL form.}
\label{table:DtEquations}
\end{table*}

\section{Experimental setup}
\label{sec:ExpSetup}

In this work, we evaluate the performance of SL-RK integration in all the variants of PDE-constrained LDDMM (see Table~\ref{table:Methods}).
The evaluation has been performed consistently with our previous work~\cite{Hernandez_19,Hernandez_19c}, in order to show the improvement of the
proposed integration method over RK integration.
In addition, we have performed an extensive evaluation of the most memory efficient stationary methods
in the frameworks of Klein et al.~\cite{Klein_09} and Rohlfing et al.~\cite{Rohlfing_12}
in order to establish the position achieved by PDE-constrained LDDMM methods in these evaluation rankings.

\begin{table*}[!t]
\begin{center}
\scriptsize
\begin{tabular}{|c|c|c|c|}
\hline
Abbreviation & Method & Integration & Publication \\
\hline
PDE-LDDMM RK & PDE-constrained LDDMM & RK & Mang et al. \cite{Mang_15} \\
\hline
PDE-LDDMM st. eq. RK & \begin{tabular}{c} PDE-constrained LDDMM \\ based on the state equation \end{tabular} & RK & Hernandez \cite{Hernandez_19c} \\
\hline
PDE-LDDMM def. eq. RK & \begin{tabular}{c} PDE-constrained LDDMM \\ based on the deformation state equation \end{tabular} & RK & Hernandez \cite{Hernandez_19c} \\
\hline
\hline
PDE-LDDMM SL & PDE-constrained LDDMM & SL-RK & Mang et al. \cite{Mang_17b} \\
\hline
PDE-LDDMM st. eq. SL & \begin{tabular}{c} PDE-constrained LDDMM \\ based on the state equation \end{tabular} & SL-RK & this work \\
\hline
PDE-LDDMM def. eq. SL & \begin{tabular}{c} PDE-constrained LDDMM \\ based on the deformation state equation \end{tabular} & SL-RK & this work \\
\hline
\hline
BL PDE-LDDMM RK & band-limited PDE-constrained LDDMM & RK & Hernandez \cite{Hernandez_19} \\
\hline
BL PDE-LDDMM st. eq. RK & \begin{tabular}{c}band-limited PDE-constrained LDDMM \\ based on the state equation \end{tabular} & RK & Hernandez \cite{Hernandez_19c} \\
\hline
BL PDE-LDDMM def. eq. RK & \begin{tabular}{c} band-limited PDE-constrained LDDMM \\ based on the deformation state equation \end{tabular} & RK & Hernandez \cite{Hernandez_19c} \\
\hline
\hline
BL PDE-LDDMM SL & band-limited PDE-constrained LDDMM & SL-RK & this work \\
\hline
BL PDE-LDDMM st. eq. SL & \begin{tabular}{c} band-limited PDE-constrained LDDMM \\ based on the state equation \end{tabular} & SL-RK &  this work \\
\hline
BL PDE-LDDMM def. eq. SL &  \begin{tabular}{c} band-limited PDE-constrained LDDMM \\ based on the deformation state equation \end{tabular} & SL-RK & this work \\
\hline
\end{tabular}
\caption{List of the PDE-constrained LDDMM methods compared in this work. The table gathers the abbreviations used in this experimental section,
the name of the method given in Section 2, the integration scheme, and the publication where the method first appeared. }
\label{table:Methods}
\end{center}
\end{table*}

\subsection{Datasets}

\noindent We have used five different databases in our evaluation:

{\bf NIREP16} contains 16 skull-stripped brain images with the segmentation of 32 gray matter structures.
The dimension of the images is 256 $\times$ 300 $\times$ 256 with a voxel size of $ 0.7 \times 0.7 \times 0.7$ mm.
The acquisition and post-processing details can be found at the web page~\url{http://www.nirep.org}.
The most remarkable features of this dataset are the excellent image quality and the ventricle sizes that are usually small.
The geometry of the segmentations provides a specially challenging framework for deformable registration evaluation.

{\bf LPBA40} contains 40 skull-stripped brain images without the cerebellum and the brain stem.
LPBA40 is provided with the segmentation of 50 gray matter structures together with the caudate, putamen, and hippocampus.

LPBA40 protocols can be found at

\noindent \url{http://www.loni.ucla.edu/Protocols/LPBA40}.
The image quality in LPBA40 is, overall, acceptable.
The variability of the ventricle sizes is high.

{\bf IBSR18} contains 18 brain images with the segmentation of 96 cerebral structures.
The masks for skull-stripping are available with the dataset.
In addition, the release {\bf IBSR\_V2.0} skull-stripped NIFTI~\cite{Rohlfing_12} contains 18
skull-stripped brain images with the segmentation of 62 cerebral structures.

This dataset provides the segmentation of brain structures of interest for the
evaluation of image registration methods.
The image quality is low. For example, most of the images show motion artifacts.
The variability of the ventricle sizes is high.

{\bf CUMC12} contains 12 full brain images with the segmentation of 130 cerebral structures.
The masks for skull-stripping are available with the dataset.

Overall, the image quality is acceptable, although some of the images are noisy.
The contrast of the images is low.
The variability of the ventricle sizes is high.

{\bf MGH10} contains 10 full brain images with the segmentation of 106 cerebral structures.
The masks for skull-stripping are available with the dataset.

Overall, the image quality is acceptable, although some of the images are noisy.
The contrast of the images is low.
Ventricle sizes are usually all big.

\subsection{Image registration pipeline}

The evaluation consistent with our previous work was performed in a subsampled NIREP16 database.
The registrations were carried out from the first subject to every other subject in the database,
yielding to a total of 15 registrations per method.
The subsampled NIREP16 database was obtained from the resampling of the images into volumes of size
$180 \times 210 \times 180$ with a voxel size of $1.0 \times 1.0 \times 1.0$ mm after the alignment
to a common coordinate system using affine transformations.
The images were scaled between $0$ and $1$.
The affine alignment and subsampling were performed using the Insight Toolkit (ITK).
The PDE-constrained registration methods were executed directly on this dataset.

For benchmarking, we run single- and multi-resolution versions of the SyN version of ANTS diffeomorphic registration~\cite{Avants_08}
with $L^2$ image similarity (ANTS-SSD).

The evaluation in the framework of Klein et al. was performed in NIREP16, LPBA40, IBSR18, CUMC12, and MGH10 databases.
The IBSR18, CUMC12 and MGH10 images normalized with respect to the MNI152 space were used as input data.
The registrations were carried out from every subject to every other subject in each database yielding to a total of 2154 registrations per method.
The evaluation in the framework of Rohlfing et al. was performed in IBSR18 database, with a total of 306 registrations per method.

The NIREP16, LPBA40, IBSR18, CUMC12, and MGH10 images were preprocessed similarly to~\cite{Klein_09}.
In the first place, N4 bias field correction and histogram matching were applied to all the images.
To perform these preprocessing steps we used the algorithms available in ITK.
The images were scaled between $0$ and $1$.
Next, we performed an affine registration between all the image pairs.
Instead of using the affine registered images as input of our non-rigid registration methods,
we used the affine transformation as input, and it was included in the parameterization of the
diffeomorphic transformations.

Subsampled NIREP experiments were run on a cluster equipped with one NVidia Titan RTX with 24 GBS of video memory
and an Intel Core i7 with 64 GBS of DDR3 RAM.
NIREP16, LPBA40, IBSR18, CUMC12, and MGH10 experiments were run on a cluster equipped with four NVidia GeForce GTX 1080 ti with 11 GBS of video
memory and an Intel Core i7 with 64 GBS of DDR3 RAM.
The codes were developed in the GPU with Matlab 2017a and Cuda $8.0$.
Since Matlab lacks a 3D GPU cubic interpolator, we implemented in a Cuda MEX file the GPU cubic interpolator
with prefiltering proposed in~\cite{Ruijters_12}.

\subsection{Parameter configuration}

Regularization parameters were selected from a search of the optimal parameters in the registration experiments
performed in our previous work~\cite{Hernandez_19c}.
We selected the parameters $\sigma^2 = 1.0$, $\alpha = 0.0025$, and $s = 2$.

The optimization was run a maximum of 10 iterations with the stopping conditions used in~\cite{Mang_15}.
The maximum number of PCG iterations was selected equal to 5.
These parameters were selected as optimal in our previous work since the methods achieved state of the art
accuracy at a reasonable amount of time~\cite{Hernandez_19}.

The experiments were performed with band sizes of $32 \times 32 \times 32$ for BL PDE-constrained LDDMM based on
the state and on the deformation state equations, and band sizes of $40 \times 40 \times 40$ for original BL PDE-constrained LDDMM.
This selection was found as optimal for each method in our previous work~\cite{Hernandez_19,Hernandez_19c}.

For SL-RK integration, $n_t$ was selected equal to $5$ for all the methods.
For RK integration, $n_t$ was selected equal to $25$ for the BL PDE-constrained LDDMM based on the state and on the
deformation state equations, and $50$ for the spatial methods due to stability issues.
In the evaluation with LPBA40, IBSR18, CUMC12, and MGH10 datasets, $n_t=25$ showed stability issues in a considerable
number of experiments and it was raised to $50$.

ANTS-SSD was run with the following parameters for the single-resolution experiments

\noindent \texttt{\$synconvergence="[50,1e-6,10]"},

\noindent \texttt{\$synshrinkfactors="1"},

\noindent and \texttt{\$synsmoothingsigmas="3vox"}.

For the multi-resolution experiments the parameters were set to

\noindent \texttt{\$synconvergence="[50x50x50,1e-6,10]"},

\noindent \texttt{\$synshrinkfactors="4x2x1"},

\noindent and \texttt{\$synsmoothingsigmas="3x2x1vox"}.

\noindent The selection of the number of iterations was in agreement with the number of outer $\times$ inner iterations used in Gauss--Newton--Krylov
optimization.

\section{Results}
\label{sec:Results}

\subsection{Subsampled NIREP16 evaluation results}

\subsubsection{Convergence analysis}

\begin{table*}[!t]
\begin{center}
\scriptsize

\begin{tabular}{|c|c|c|c|c|c|}
\hline
Method & $n_t$ & $MSE_{rel}$ & $\Vert g\Vert_{\infty,{rel}}$ & $\min(J(\phi_{1}^v)^{-1})$ & $\max(J(\phi_{1}^v)^{-1})$ \\
\hline
St. PDE-LDDMM RK & 50 & 18.29 $\pm$ 2.83 & 0.07 $\pm$ 0.05 & 0.16 $\pm$ 0.05 & 3.70 $\pm$ 0.51 \\
St. PDE-LDDMM st. eq. RK & 30 & 18.42 $\pm$ 2.71 & 0.18 $\pm$ 0.19 & 0.07 $\pm$ 0.06 & 2.95 $\pm$ 0.41 \\
St. PDE-LDDMM def. st. eq. RK & 30 & 16.83 $\pm$ 1.53 & 0.16 $\pm$ 0.08 & 0.12 $\pm$ 0.05 & 4.12 $\pm$ 1.41 \\
\hline
St. PDE-LDDMM SL & 5 & 19.55 $\pm$ 1.76 & 0.14 $\pm$ 0.13 & 0.03 $\pm$ 0.04 & 3.25 $\pm$ 0.29 \\
St. PDE-LDDMM st. eq. SL & 5 & 18.34 $\pm$ 1.68 & 0.23 $\pm$ 0.21 & 0.03 $\pm$ 0.04 & 3.46 $\pm$ 0.69 \\
St. PDE-LDDMM def. st. eq. SL & 5 & 17.10 $\pm$ 1.50 & 0.12 $\pm$ 0.05 & 0.14 $\pm$ 0.05 & 5.02 $\pm$ 1.01 \\
\hline
NSt. PDE-LDDMM RK & 50 & 22.68 $\pm$ 5.49 & 0.26 $\pm$ 0.12 & 0.16 $\pm$ 0.05 & 3.38 $\pm$ 0.90 \\
NSt. PDE-LDDMM, st. eq. RK & 25$^*$ & 29.92 $\pm$ 2.41 & 0.32 $\pm$ 0.10 & 0.11 $\pm$ 0.07 & 2.93 $\pm$ 0.78 \\
NSt. PDE-LDDMM, def. st. eq. RK & 30 & 16.10 $\pm$ 1.70 & 0.20 $\pm$ 0.08 & 0.09 $\pm$ 0.03 & 4.55 $\pm$ 0.69 \\
\hline
NSt. PDE-LDDMM SL & 5 & 22.44 $\pm$ 6.26 & 0.51 $\pm$ 0.18 & 0.02 $\pm$ 0.02 & 3.96 $\pm$ 0.65 \\
NSt. PDE-LDDMM st. eq. SL & 5 & 21.82 $\pm$ 2.51 & 0.49 $\pm$ 0.19 & 0.01 $\pm$ 0.02 & 4.20 $\pm$ 0.94 \\
NSt. PDE-LDDMM def. st. eq. SL & 5 & 16.86 $\pm$ 1.91 & 0.13 $\pm$ 0.04 & 0.09 $\pm$ 0.02 & 7.35 $\pm$ 2.63 \\
\hline
BL St. PDE-LDDMM 40x RK & 25 & 20.99 $\pm$ 2.59 & 0.04 $\pm$ 0.03 & 0.21 $\pm$ 0.04 & 3.32 $\pm$ 0.41 \\
BL St. PDE-LDDMM st. eq. 32x RK & 25 & 18.53 $\pm$ 1.71 & 0.02 $\pm$ 0.01 & 0.11 $\pm$ 0.07 & 3.10 $\pm$ 0.35 \\
BL St. PDE-LDDMM def. st. eq. 32x RK & 25 & 17.32 $\pm$ 1.68 & 0.03 $\pm$ 0.01 & 0.13 $\pm$ 0.05 & 3.76 $\pm$ 0.49 \\
\hline
BL St. PDE-LDDMM 40x SL & 5 & 20.41 $\pm$ 1.89 & 0.02 $\pm$ 0.00 & 0.32 $\pm$ 0.02 & 9.65 $\pm$ 4.77  \\
BL St. PDE-LDDMM st. eq. 32x SL & 5 & 19.89 $\pm$ 1.76 & 0.01 $\pm$ 0.00 & 0.29 $\pm$ 0.03 & 7.94 $\pm$ 2.81 \\
BL St. PDE-LDDMM def. st. eq. 32x SL & 5 & 17.77 $\pm$ 1.66 & 0.04 $\pm$ 0.01 & 0.21 $\pm$ 0.03 & 8.23 $\pm$ 3.24 \\
\hline
BL NSt. PDE-LDDMM 40x RK & 25 & 29.30 $\pm$ 3.50 & 0.21 $\pm$ 0.08 & 0.24 $\pm$ 0.05 & 2.25 $\pm$ 0.51 \\
BL NSt. PDE-LDDMM, st. eq. 32x RK & 25 & 29.21 $\pm$ 3.96 & 0.24 $\pm$ 0.11 & 0.15 $\pm$ 0.06 & 2.71 $\pm$ 0.48 \\
BL NSt. PDE-LDDMM, def. st. eq. 32x RK & 25 & 15.68 $\pm$ 1.52 & 0.04 $\pm$ 0.01 & 0.09 $\pm$ 0.04 & 4.77 $\pm$ 0.78 \\
\hline
BL NSt. PDE-LDDMM 40x SL & 5 & 19.14 $\pm$ 2.06 & 0.11 $\pm$ 0.05 & 0.01 $\pm$ 0.07 & 5.25 $\pm$ 0.80 \\
BL NSt. PDE-LDDMM st. eq. 32x SL & 5 & 20.49 $\pm$ 1.75 & 0.17 $\pm$ 0.15 & 0.07 $\pm$ 0.03 & 4.51 $\pm$ 0.61 \\
BL NSt. PDE-LDDMM def. st. eq. 32x SL & 5 & 16.80 $\pm$ 1.57 & 0.07 $\pm$ 0.03 & 0.09 $\pm$ 0.03 & 7.37 $\pm$ 1.43 \\
\hline
\end{tabular}

\caption{ Subsampled NIREP16. Convergence results.
Mean and standard deviation of the relative image similarity error expressed in \% ($MSE_{rel}$),
the relative gradient magnitude ($\Vert g\Vert_{\infty,{rel}}$),
and the Jacobian determinant extrema associated with the transformation $(\phi^{v}_{1})^{-1}$.
St. stands for the stationary and NSt. for the non-stationary parameterization of diffeomorphisms.
(*) The $n_t$ selected for NSt. PDE-LDDMM, st. eq. with RK integration was due to a memory load greater than
the GPU maximum capacity (24 GBS). }
\label{table:QuantitativeResults}
\end{center}
\end{table*}

Table~\ref{table:QuantitativeResults} shows, averaged by the number of experiments,
the mean and standard deviation of the relative image similarity error after registration,
$$MSE_{rel} = \frac{\Vert m(1) - I_1 \Vert_{L^2}^2}{\Vert I_0 - I_1 \Vert_{L^2}^2},$$
\noindent the relative gradient magnitude,
$$\Vert g \Vert_{\infty,rel} = \frac{\Vert \nabla_{v} E({v}^n) \Vert_\infty}{\Vert \nabla_{v} E({v}^0) \Vert_\infty},$$
and the extrema of the Jacobian determinant obtained with PDE-constrained LDDMM in the subsampled NIREP16 dataset.

We have grouped the analysis of the $MSE_{rel}$ values by integration scheme, variant, image domain, and diffeomorphism
parameterization:

\begin{itemize}
 \item RK vs SL-RK. The $MSE_{rel}$ value achieved by SL-RK methods is close or improves RK methods.
 It drives our attention the bad performance of PDE-LDDMM st. eq. for the non-stationary parameterization and RK integration
 which is indeed improved by SL-RK integration.
 \item Variants. The best performing variant is PDE-LDDMM def. st. eq. This result is persistent for the different
 integration schemes, the spatial or BL parameterization, and the diffeomorphic parameterization.
 \item SP vs BL. The $MSE_{rel}$ values obtained with the spatial methods are slightly degraded by the BL methods as expected.
 The degradation is only shown in some cases.
 \item St. vs NSt. The non-stationary methods do not outperform the stationary methods in a consistent manner.
 The (out)performance depends on the variant.
\end{itemize}

For the bad-performing methods, the high values of the relative gradient indicate the stagnation of the convergence.
For the BL methods, the relative gradient was reduced to average values ranging from $0.02$ to $0.04$, which means that the optimization
was stopped in acceptable energy values in all cases.
All the Jacobians remained above zero.

\subsubsection{Evaluation}

The evaluation is based on the accuracy of the registration results for template-based segmentation.

The Dice Similarity Coefficient (DSC) is selected as the evaluation metric.
Given two segmentations $S$ and $T$, the DSC is defined as
\begin{equation}
DSC(S,T) = \frac{2 \textnormal{Vol} (S \cap T)}{\textnormal{Vol}( S ) + \textnormal{Vol}( T )}.
\end{equation}
\noindent This metric provides the value of 1 if $S$ and $T$ exactly overlap and gradually decreases towards 0
depending on the overlap of the two volumes.

Figure~\ref{fig:NIREPSubDSC} shows, in the shape of box and whisker plots, the statistical distribution of the
DSC values obtained after the registration across the 32 segmented structures.
For the single-resolution experiments, the performance of the benchmark method ANTS-SSD was under $50 \%$.
For the multi-resolution experiments, the average DSC value achieved by ANTS-SSD equals to $55.59 \%$.
We have selected this value as a baseline of good registration accuracy for methods with $L^2$-based image
similarity.

A great number of PDE-LDDMM methods showed similar values or even improved ANTS-SSD performance.
The best performing variant was our PDE-LDDMM based on the deformation state equation (boxes in pink tones).
This variant showed similar results for RK and SL-RK integration regardless the image domain and
diffeomorphism parameterization.

For the variant associated with the PDE-constrained benchmark methods~\cite{Mang_15,Mang_17b} (boxes in blue tones),
RK integration slightly outperformed SL-RK integration for the stationary parameterization.
For the non-stationary parameterization the median DSC value for RK integration was under the value
for SL-RK integration.
Similarly, RK slightly outperformed SL-RK integration for BL PDE-LDDMM and BL PDE-LDDMM based on the state
equation for the stationary parameterization.
On the contrary, SL-RK integration greatly outperformed RK integration for the non-stationary parameterization.

The variant based on the state equation (boxes in green tones) performed similarly to benchmark PDE-LDDMM variant for the
stationary parameterization for the same integration scheme.
However, it is remarkable the low performance achieved by this variant for the non-stationary parameterization
and RK integration in both image domains.

\begin{figure*} [!t]
\centering
\scriptsize
\begin{tabular}{c}
\includegraphics[angle = 0, width = 1.25 \textwidth]{./Figures/NIREPSubsampled_DSC_SIIMS2020SL.eps}
\end{tabular}
\caption{Subsampled NIREP16.
Volume overlap obtained by the registration methods measured in terms of the DSC between
the warped and the corresponding manual target segmentations.
Box and whisker plots show the distribution of the DSC values averaged over the 32 NIREP manual segmentations.
The whiskers indicate the minimum and maximum of the DSC values.
The horizontal lines in the plot indicate the first, second, and third quartiles of multi-resolution ANTS-SSD.
}
\label{fig:NIREPSubDSC}
\end{figure*}

\subsubsection{Computational complexity}

\begin{table*}[!tpb]

\begin{center}
\scriptsize

\begin{tabular}{|c|c|c|c|}
\hline
Method & Integration & VRAM (MBS) & total time (s) \\
\hline
St. PDE-LDDMM & RK & 10311 & 2281.40 $\pm$ 378.56 \\
St. PDE-LDDMM st. eq. & RK & 17229 & 877.85 $\pm$ 176.03 \\
St. PDE-LDDMM def. st. eq. & RK & 15321 & 1298.50 $\pm$ 151.19 \\
\hline
St. PDE-LDDMM & SL-RK & 3901 & 118.29 $\pm$ 4.77 \\
St. PDE-LDDMM st. eq. & SL-RK & 6267 & 139.91 $\pm$ 22.06 \\
St. PDE-LDDMM def. st. eq. & SL-RK & 5935 & 221.01 $\pm$ 2.40 \\
\hline
NSt. PDE-LDDMM & RK & 14913 & 3443.17 $\pm$ 932.45 \\
NSt. PDE-LDDMM, st. eq. & RK & 20635 & 816.79 $\pm$ 480.61 \\
NSt. PDE-LDDMM, def. st. eq. & RK & 19065 & 2222.04 $\pm$ 369.96 \\
\hline
NSt. PDE-LDDMM & SL-RK & 11155 & 191.96 $\pm$ 45.53 \\
NSt. PDE-LDDMM st. eq. & SL-RK & 12429 & 214.09 $\pm$ 48.54 \\
NSt. PDE-LDDMM def. st. eq. & SL-RK & 12065 & 355.47 $\pm$ 4.38 \\
\hline
\hline
BL St. PDE-LDDMM 40x & RK & 5709 & 312.45 $\pm$ 5.22 \\
BL St. PDE-LDDMM st. eq. 32x & RK & 4743 & 315.00 $\pm$ 2.73 \\
BL St. PDE-LDDMM def. st. eq. 32x & RK & 1819 & 377.56 $\pm$ 4.26 \\
\hline
BL St. PDE-LDDMM 40x & SL-RK & 2685 & 68.66 $\pm$ 1.32 \\
BL St. PDE-LDDMM st. eq. 32x & SL-RK & 2877 & 80.76 $\pm$ 0.35 \\
BL St. PDE-LDDMM def. st. eq. 32x & SL-RK & 2365 & 116.99 $\pm$ 1.04 \\
\hline
BL NSt. PDE-LDDMM 40x & RK & 6657 & 825.15 $\pm$ 19.93 \\
BL NSt. PDE-LDDMM, st. eq. 32x & RK & 4789 & 535.52 $\pm$ 279.65 \\
BL NSt. PDE-LDDMM, def. st. eq. 32x & RK & 1863 & 905.48 $\pm$ 25.13 \\
\hline
BL NSt. PDE-LDDMM 40x & SL-RK & 6131 & 158.79 $\pm$ 15.29 \\
BL NSt. PDE-LDDMM st. eq. 32x & SL-RK & 6171 & 159.49 $\pm$ 20.59 \\
BL NSt. PDE-LDDMM def. st. eq. 32x & SL-RK & 5577 & 249.60 $\pm$ 10.85 \\
\hline

\end{tabular}
\\
\vspace{0.1 cm}
\caption{Computational complexity. GPU peak memory usage and mean and standard deviation of the total computation
time. Experiments run in an NVidia Titan RTX with 24 GBS of video memory.

}
\label{table:Complexity}

\end{center}

\end{table*}

Table~\ref{table:Complexity} shows the VRAM peak memory reached through the computations,
and the average and standard deviation of the total computation time in the subsampled NIREP16 experiments.
For the spatial methods, SL-RK integration achieved a substantial time and memory reduction over RK integration.
The time and memory reduction achieved by SL-RK over RK integration was also considerable for the BL parameterized methods.
For the stationary parameterization, the BL parameterization further decreased the complexity of spatial SL-RK integration methods,
as expected.
However, SL-RK integration did not reduced the memory VRAM usage for the non-stationary parameterization while the total computation
time was effectively reduced.

\subsubsection{Qualitative assessment}

For a qualitative assessment of the proposed registration methods, we show the registration results obtained by
Mang et al. benchmark methods~\cite{Mang_15,Mang_17}, and
PDE-LDDMM based on the deformation state equation in a selected experiment representative of a difficult deformable registration problem.
For the non-stationary parameterization the images for a qualitative assessment were similar to the stationary parameterization.
Figure~\ref{fig:QualitativeBenchmark} shows the warped images, the difference between the warped and the target images after registration,
and the velocity field.

\subsection{NIREP16, LPBA40, IBSR18, CUMC12, and MGH10 evaluation results}

From the evaluation measurements used in Klein et al. framework, we focus on the accuracy of the registration results for
template-based segmentation.
Since~\cite{Klein_09}, this has been adopted as a widely extended criterion for non-rigid registration evaluation.
From the metrics proposed in~\cite{Klein_09}, we select the Dice Similarity Coefficient (DSC) as evaluation metric.

Figures~\ref{fig:NIREPDSC},~\ref{fig:KleinDSC1} and~\ref{fig:KleinDSC2} show the statistical distribution of the DSC values obtained after the registration
across the manually segmented structures for the five databases.
For the NIREP16 dataset, we show the results obtained with ANTS-SSD.
For the remaining databases, we include the results reported in~\cite{Klein_09} for affine registration (FLIRT),
and three diffeomorphic registration methods: Diffeomorphic Demons, SyN, and Dartel.

On the one hand, the results from NIREP16 show how PDE-LDDMM based on the deformation state equation outperformed the other variants
of PDE-constrained LDDMM.
The distribution of the method parameterized in the BL domain resulted almost identical to the distribution of the method parameterized
in the spatial domain.
As happened with subsampled NIREP16 evaluation, ANTS-SSD was among the worst performing methods.
BL PDE-LDDMM and BL PDE-LDDMM based on the state equation with RK integration exceeded the maximum VRAM capacity of our GPUs for this dataset.

\begin{figure*} [!t]
\centering
\scriptsize
\begin{tabular}{c}
\includegraphics[angle = 0, width = 0.60 \textwidth]{./Figures/NIREP16_DSC.eps}
\\
\end{tabular}
\caption{ NIREP16. Distribution of the DSC values averaged over the 32 NIREP manual segmentations in the 240 experiments.
The whiskers indicate the minimum and maximum of the DSC values.
The methods running out of 11 GBS VRAM are indicated in the plots.
The horizontal red lines indicate the first and third quartiles of BL PDE-LDDMM based on the deformation state equation.}
\label{fig:NIREPDSC}
\end{figure*}

On the other hand, the results from LPBA40, IBSR18, CUMC12, and MGH10 show that, from the PDE-constrained LDDMM methods,
the best performing method was BL PDE-LDDMM based on the deformation state equation and SL-RK integration.
The performance of the method with RK integration was slightly lower.
The performance of original PDE-LDDMM, PDE-LDDMM based on the state equation and their BL versions was significantly lower in IBSR18,
CUMC12, and MGH10 databases.
For these methods, RK integration performed slightly better than SL-RK integration.
These results are consistent with NIREP16 evaluation results.

In IBSR18, CUMC12, and MGH10 databases, our PDE-constrained LDDMM methods were not able to reach SyN or Dartel performance.
This is probably because the image similarity metrics used in these methods (Cross-Correlation and multinomial model, respectively)
favour the accuracy in template-based segmentation.
In contrast, PDE-constrained LDDMM uses SSD, which is known to restrict the performance in template-based segmentation.
However, in LPBA40 databases, our best performing PDE-constrained LDDMM methods overpass Dartel and almost reached SyN
performance, while showing a significantly reduced number of outliers.

Our methods significantly outperformed FLIRT and Diffeomorphic Demons, where the third quartile in the distribution of our
best performing method was close to the median of Demons for the four databases.
It should be noticed that Diffeomorphic Demons also uses SSD as image similarity metric.

\begin{figure*}[t]
\centering
\scriptsize
\begin{tabular}{c}
 \includegraphics[angle = 0, width=0.75\textwidth]{./Figures/LPBA40_DSC.eps} \\
 \includegraphics[angle = 0, width=0.75\textwidth]{./Figures/IBSR18_DSC_Klein_June2020.eps} \\
\end{tabular}
\caption{ LPBA40 and IBSR18.
Distribution of the DSC values averaged over the manual segmentations in the registration experiments.
The whiskers indicate the minimum and maximum of the DSC values.
The horizontal red lines indicate the first and third quartiles of BL PDE-LDDMM based on the deformation state equation.
}
\label{fig:KleinDSC1}
\end{figure*}

\begin{figure*}[t]
\centering
\scriptsize
\begin{tabular}{c}
 \includegraphics[angle = 0, width=0.75\textwidth]{./Figures/CUMC12_DSC.eps} \\
 \includegraphics[angle = 0, width=0.75\textwidth]{./Figures/MGH10_DSC.eps} \\
\end{tabular}
\caption{ CUMC12 and MGH10.
Distribution of the DSC values averaged over the manual segmentations in the registration experiments.
The whiskers indicate the minimum and maximum of the DSC values.
The horizontal red lines indicate the first and third quartiles of BL PDE-LDDMM based on the deformation state equation.
}
\label{fig:KleinDSC2}
\end{figure*}

\subsection{IBSR18 V2.0 evaluation results}

\begin{figure*}[h]
\centering
\scriptsize
\begin{tabular}{c}
 \includegraphics[angle = 0, width=0.75\textwidth]{./Figures/IBSR18_DSC_Rohlfing_Left.eps} \\
 \includegraphics[angle = 0, width=0.75\textwidth]{./Figures/IBSR18_DSC_Rohlfing_Right.eps} \\
\end{tabular}
\caption{ IBSR18 V2.0. Volume overlap obtained by the proposed registration methods.
The box and whisker plots show the distribution of the DSC values averaged over the manual segmentations for each region.
The whiskers indicate the minimum and maximum of the DSC values.
For each group of plots, the first three correspond with the spatial domain parameterization and the three last
with the band-limited domain parameterization.
WM and GM stand for white matter and grey matter, respectively.
}
\label{fig:IBSRDSC}
\end{figure*}

Figure~\ref{fig:IBSRDSC} shows the statistical distribution of the DSC values obtained by our proposed registration
methods in the regions of interest of Rohlfing et al. evaluation framework~\cite{Rohlfing_12}.
Consistently with the rest of the evaluation results, the best performing method was BL PDE-LDDMM based on the deformation
state equation, which significantly outperformed the others in the great majority of regions.

\section{Discussion}
\label{sec:Discussion}

The increase of the computational efficiency achieved by the combined over the split BL and SL-RK approaches was significant in terms of
computation time and memory.
The reduction of the memory requirements allowed us to perform the evaluation of the SL-RK PDE-constrained LDDMM methods extensively,
even in the highest-resolution level of NIREP16.

In all the evaluation frameworks, BL PDE-constrained LDDMM based on the deformation state equation with SL-RK integration resulted our best performing method.
This method achieved an identical DSC distribution compared with RK integration in the NIREP16 database.
In NIREP16 database, the method greatly outperformed ANTS-SSD.
In LPBA40, IBSR18, CUMC12, and MGH10 databases, the method outperformed Diffeomorphic Demons.
In addition, the evaluation results in the regions of interest of Rohlfing et al. corroborated its excellent performance.

For Mang et al. benchmark PDE-constrained LDDMM methods~\cite{Mang_15,Mang_17b}, our evaluation results reported a significative loss of accuracy between
RK and SL-RK integration. Interestingly, this loss of accuracy was not observed for our best performing method.

In IBSR18, CUMC12, and MGH10 databases, our PDE-constrained LDDMM methods were not able to reach SyN or Dartel performance.
This is because SSD image similarity metric restricts the performance of the methods in template-based segmentation.
This problem will be tackled in future work by formulating the PDE-constrained problem with other image similarity metrics such as
Normalized Cross Correlation, local Normalized Cross Correlation, or Mutual Information~\cite{Modersitzki_09_book}.
We expect that this change in the formulation of the problem will increase the performance results of PDE-constrained LDDMM.
In addition, it will allow us to apply these methods to other clinical applications involving multi-modal registration.

Simultaneously to the development of this work, Mang et al. released Claire software~\cite{Mang_19}.
The software is intended to exploit massive CPU based parallel computing architectures to accelerate the computation time of PDE-LDDMM.
The codes implement original PDE-constrained LDDMM with a variational extension to nearly incompressible fluids and include $H^1$ and $H^2$ regularization terms.
The software is restricted to the stationary parameterization of diffeomorphisms.
The PDE integration scheme is SL-RK.
The software includes a sophisticated multi-level preconditioner that shows to improve the convergence of PCG with respect to the original proposal in~\cite{Mang_15}.
The massive computation allows increasing the number of inner and outer iterations and use the norm of the gradient as stopping
condition for achieving an extraordinary accuracy at convergence in a simulated experiment.

In contrast to Claire, our methods of choice are intended to run completely in the VRAM of commodity GPUs ($<$ 4GBS).
We have limited the variational formulation to the one proposed in~\cite{Mang_15}, although it is straightforwardly extendible to the
nearly incompressible fluid problem.
We have limited our study to the traditional LDDMM regularizer.
Our software works for the stationary and the non-stationary parameterization of diffeomorphisms.
We have limited the preconditioner to the one proposed in~\cite{Mang_15} since we are interested on the comparison of the three different
variational variants.
We used the stopping conditions suggested in~\cite{Modersitzki_09_book} and used for PDE-LDDMM in~\cite{Mang_15,Mang_17,Hernandez_19,Hernandez_19b,Hernandez_19c}.
The variety of methods, the extensive evaluation conducted in this work, and our modest hardware capacity hindered us
the use of the inner and outer iteration values needed for achieving the stopping conditions based on the norm of the gradient
suggested in~\cite{Mang_19}.
In fact, we observed in a selected NIREP experiment that increasing the number of inner iterations in PCG resulted into a faster initial convergence that finally
stagnated in greater $MSE_{rel}$ values and lower $DSC$ scores than our considered stopping conditions.
This stagnation was also reported in~\cite{Mang_19} for the simulated experiment.
Instead, our selected inner and outer values consumed a reasonable amount of time while obtaining state of the art results for the
evaluation metrics.

Comparing Claire and our results, we believe that it would be of interest to implement our best performing variant as a part of Claire's software.
In the other direction, it would be very interesting to adopt the multi-level preconditioners in our methods.

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

In this work, we have proposed to combine the two different methodological approaches used to circumventing the huge computational complexity of
Gauss--Newton--Krylov PDE-constrained LDDMM.
In particular, we have included Semi-Lagrangian Runge-Kutta integration~\cite{Mang_17b} in the variants of band-limited PDE-constrained LDDMM proposed
in~\cite{Hernandez_19,Hernandez_19c}
for further increasing the computational efficiency of these methods.
The resulting methods have been extensively evaluated in five different datasets following three different evaluation frameworks.
To our knowledge, this is the first time that SL-RK integration is implemented in the framework of PDE-LDDMM for the non-stationary parameterization and
in the space of band-limited vector fields.
Moreover, our work first provides the position achieved by PDE-constrained LDDMM methods in the ranking of Klein et al. evaluation.

This study positions the formulation of BL PDE-constrained LDDMM based on the deformation state equation and SL-RK integration as the best performing
among all PDE-constrained LDDMM methods in terms of accuracy and efficiency.
The proposed method has reached the highest sensitivity in the classification of stable vs progressive mild cognitive impaired conversors in the Alzheimer's
Disease Neuroimaging Initiative (ADNI) database using convolutional neural networks. This result has been recently published in~\cite{Ramon_20}.

In future work, we will extend this formulation to other relevant physically meaningful LDDMM approaches such as the nearly incompressible method in~\cite{Mang_16},
and the geodesic shooting approach in~\cite{Hernandez_19b}.
We will explore the advantages of using the multi-level preconditioner in~\cite{Mang_19}.
We will adapt our methods for the use of alternative image similarity metrics that usually outperform SSD in registration evaluation rankings.
Finally, we will work in the understanding of which of the features of PDE-LDDMM allows the exceptional classification rates related with Alzheimer's disease
conversion shown in~\cite{Ramon_20}.

\section*{Acknowledgements}

The author would like to acknowledge the anonymous reviewers for their revision of the manuscript.

The author would like to give special thanks to Wen Mei Hwu from the University of Illinois
for interesting ideas in the GPU implementation of the methods,
and Nacho Navarro and Rosa Badia from the Barcelona Supercomputing Center (BSC) for his help.

This work was partially supported by the national research grant TIN2016-80347-R.
In addition, this work was supported by NVIDIA through the Polytechnical University of Catalonia /
Barcelona Supercomputing Center (UPC/BSC) GPU Center of Excellence.

\bibliographystyle{splncs}

\bibliography{abbsmall.bib,Diffeo.bib,OpticalFlow.bib,Miccai2015.bib}

config: arxiv