TITLE: Explicit description of the "simplicial tensor product" of chain complexes
QUESTION [12 upvotes]: Recall that there is an equivalence of categories (Dold-Kan) $$N:\mathrm{s}\mathbf{Ab}\simeq \operatorname{Ch}_{\geq 0}(\mathbf{Ab}):\Gamma$$ between simplicial abelian groups and (connective) chain complexes, where $N$ sends a simplicial abelian group to its associated normalized chain complex.
Using this equivalence of categories, we can, by transport of structure, give an unorthodox tensor product on the category of chain complexes.  We may define this by the formula $X\otimes_\Delta Y=N(\Gamma(X)\otimes \Gamma(Y)),$ where the tensor product on the righthand side is the tensor product (taken pointwise) of simplicial abelian groups.
Then my question: Is there an explicit description of this tensor product in terms of the chain complexes themselves?

REPLY [10 votes]: The bad news is that in degree $n$, this tensor product has $3^n$ terms.
The functor $\Gamma$ can be roughly described as follows.  If we write $[n]$ for the ordered set $0 < 1 < \cdots < n$, then
$$
\Gamma(C)_n = \bigoplus_{k} \bigoplus_{\phi\colon [n] \twoheadrightarrow [k]} C_k.
$$
The face maps have two characters.  The map $d_i$ for $i > 0$ simply deletes the element $i$ from the ordered set $[n]$, and reindexes; if the resulting map $[n-1] \to [k]$ is no longer surjective, the corresponding factor maps to zero.  By contrast, the map $d_0$ deletes $0$ and reindexes, but if the corresponding map $\phi$ is no longer surjective its image is isomorphic to $[k-1]$, and we apply the boundary map.
When you take the tensor product of $\Gamma(C)$ and $\Gamma(D)$ levelwise, you get a direct sum indexed by pairs of surjections $[n] \twoheadrightarrow [p]$ and $[n] \twoheadrightarrow [q]$.
The functor $N$ then will take the quotient of this by the subcomplex of degenerate ones; those where the maps $[n] \twoheadrightarrow [p]$ and $[n] \twoheadrightarrow [q]$ factors through a surjection $[n] \twoheadrightarrow [m]$.  In practice, the pairs which are not degenerate are those for which the map $[n] \to [p] \times [q]$ is injective.
As a result, we have that
$$
N(\Gamma(C) \otimes \Gamma(D))_n = \bigoplus_\phi C_p \otimes D_q
$$
where the sum is indexed by injections $[n] \to [p] \times [q]$ where composing with either projection is surjective.
In practice, you can index this direct sum by $n$-tuples of strings of elements from$\{N, NE, E\}$, representing a path of length $n$; $p$ and $q$ are determined by the height and width of the path.
Unfortunately, a chain complex isn't very useful without its differential, and that's more complicated to describe.  The boundary map is the alternating sum of face maps; each face map deletes $i$ from the ordered set $0 < \cdots < n$ and reindexes.  If $i > 0$ and one of the resulting projections to $[p]$ or $[q]$ is no longer surjective, the corresponding factor maps to zero; if $i = 0$ then one or both of the maps to $[p]$ or $[q]$ misses zero, the appropriate image(s) are isomorphic to $[p-1]$ or $[q-1]$, and we apply the boundary map on those factors.
UPDATE: One way to write this is in a group homology style.  Then we can view elements of $\Gamma(C)_n$ as decorated with these $n$-tuple paths, the width and height determine which group we land in, and the boundary map takes an alternating sum of deleting commas with the understanding that $NN = EE = 0$.  So, for example, taking the boundary of this element in $C_3 \otimes D_2$:
$$
d(a \otimes b)_{(N,E,N,NE)} = (da) \otimes b_{(E,N,NE)} - a \otimes b_{(NE,N,NE)} + a \otimes b_{(N,NE,NE)} - a \otimes b_{(N,E,NNE)}.
$$
(The last term involving $NNE$ is dropped because it is zero.  The first term had a boundary map applied to the first factor because it was $N$.)


text:
\section{Multifidelity model management in statistical inference}
\label{sec:MM:InvProb}
In a Bayesian setting, inverse problems are cast in a statistical formulation where the unknown input is modeled as a random variable and is described by its posterior distribution \cite{tarantola_inverse_1982,TarantolaBook,KaipoBook,ANU:7701764}.  MCMC is a popular way to sample from the posterior distribution. Statistical inference raises several computational challenges, including the design of MCMC sampling schemes, the construction of approximate models that can reduce the costs of MCMC sampling, and the development of alternatives to MCMC sampling such as variational approaches \cite{UQDef}. Detailed discussions of these and many other important aspects of statistical inference can be found in the literature, e.g., \cite{RobertBook,KaipoBook,ANU:7701764,MCBookLiu}. We focus here on a few specific aspects of statistical inference in which multifidelity methods have been used. In particular, we survey multifidelity methods that use a two-stage formulation of MCMC, where a candidate sample has to be first accepted by a low-fidelity model before it is passed on to be either accepted or rejected by the high-fidelity model. Section~\ref{sec:MM:InvProb:Bayesian} describes the problem setup. Section~\ref{sec:MM:InvProb:TwoStageMCMC} describes a two-stage MCMC framework and Section~\ref{sec:MM:InvProb:Adaptive} discusses a framework where a low-fidelity model is adapted while it is evaluated in an MCMC algorithm. A Bayesian approach to model and correct the error of low-fidelity models is discussed in Section~\ref{sec:MM:InvProb:Error}.

\subsection{Bayesian framework for inference}
\label{sec:MM:InvProb:Bayesian}
Consider the high-fidelity model $\fhigh$ that maps inputs $\dv \in \Zcal$ onto outputs $\ov \in \Ycal$. Let $p_0$ be a prior density that describes the input $\dv$ before any measurements. Let further $\obs$ be noisy observational data with the stochastic relationship
\begin{equation}
\obs = \fhigh(\dv) + \bfeps\,,
\label{eq:MM:InvProb:Bayesian:InvProb}
\end{equation}
where $\bfeps$ is a random vector that captures the measurement error, noise, and other uncertainties of the observation $\obs$. In the following, the random vector $\bfeps$ is modeled as a zero-mean Gaussian with covariance $\bfSigma_{\bfeps} \in \mathbb{R}^{d \times d}$. Define the data-misfit function as
\[
\Phi(\dv) = \frac{1}{2}\left\|\bfSigma_{\bfeps}^{-\frac{1}{2}}\left(\fhigh(\dv) - \obs\right)\right\|^2\,,
\]
with the norm $\|\cdot\|$. The likelihood function $L: \Zcal \to \mathbb{R}$ is then proportional to
\begin{equation}
L(\obs|\dv) \propto \exp(-\Phi(\dv))\,.
\label{eq:MM:InvProb:LikelihoodHi}
\end{equation}
An evaluation of the likelihood $L$ entails an evaluation of the high-fidelity model $\fhigh$. The posterior probability density is
\[
p(\dv|\obs) \propto L(\obs|\dv)p_0(\dv)\,.
\]

We note that there are many challenging and important questions regarding the formulation of an inverse problem in the Bayesian setting. For example, the selection of the prior is a delicate and problem-dependent question, see for example \cite[Section~3.3]{KaipoBook} for a discussion of prior models. We do not address those issues here, but note that our stated goal of a multifidelity formulation is to recover a solution of the outer-loop problem (here the inverse problem) that retains the accuracy of the high-fidelity formulation. Thus, the multifidelity approaches described below will inherit the choices made in the high-fidelity formulation of the Bayesian inverse problem.

\begin{algorithm}[t]
\caption{Metropolis-Hastings}\label{alg:MCMC-MH}
\begin{algorithmic}[1]
\Procedure{MetropolisHastings}{$L, p_0, q, \nrpts$}
\State Choose a starting point $\dv_0$
\For{$i = 1, \dots, \nrpts$}
\State Draw candidate $\dv^*$ from proposal $q(\cdot | \dv_{i - 1})$
\State Compute acceptance probability
\[
\alpha(\dv_{i - 1}, \dv^*) = \min\left\{1, \frac{q(\dv_{i - 1} | \dv^*)L(\obs|\dv^*)p_0(\dv^*)}{q(\dv^*|\dv_{i - 1})L(\obs|\dv_{i - 1})p_0(\dv_{i-1})}\right\}
\]
\State Set the sample $\dv_i$ to
\[
\dv_i = \begin{cases}
\dv^*\,,\qquad &\text{with probability }\alpha(\dv_{i - 1}, \dv^*)\,,\\
\dv_{i - 1}\,,\qquad &\text{with probability }1 - \alpha(\dv_{i - 1}, \dv^*)
\end{cases}
\]
\EndFor
\State \Return $\dv_1, \dots, \dv_{\nrpts}$
\EndProcedure
\end{algorithmic}
\end{algorithm}

\paragraph{Exploring the posterior}
The solution of the inference problem is explored by drawing samples from the posterior distribution. The posterior samples can then be used to estimate the input with the maximum posterior density and to estimate expectations of functions of interest $h: \Zcal \to \mathbb{R}$ with respect to the posterior distribution. For example, one could be interested in the expected value of $h$ over the posterior distribution
\begin{equation}
\mathbb{E}[h] = \int_{\Zcal} h(\dv)p(\dv|\obs)\mathrm d\dv\,.
\label{eq:MM:InvProb:ExpH}
\end{equation}
MCMC methods are a popular way to sample from the posterior distribution, which have been studied extensively \cite{roberts1997,haario_adaptive_2001,Marzouk20091862,GirolamiHamiltonianMC,cotter2013,GhattasStochNewton,0266-5611-30-11-114015,Cui2016109}; see also the books \cite{GilksMCMC,MCBookLiu}. One example of an MCMC method is the Metropolis-Hastings algorithm, which is an iterative scheme that draws candidate samples from a proposal distribution and then accepts the candidate sample with a probability that depends on the ratio of the posterior at the current candidate sample and the posterior at the sample of the previous iteration.

\paragraph{Metropolis-Hastings algorithm}
Algorithm~\ref{alg:MCMC-MH} summarizes the Metropolis-Hastings algorithm. Inputs are the likelihood $L$, the prior density $p_0$, a proposal density $q$, and the number of samples $\nrpts$. The proposal density $q$ is used to draw the next candidate sample. A typical choice for the proposal density is a Gaussian distribution that is centered at the sample of the previous iteration. In each iteration $i = 1, \dots, \nrpts$, a candidate sample $\dv^*$ is drawn from the proposal $q(\cdot | \dv_{i - 1})$ that depends on the sample $\dv_{i - 1}$ of the previous iteration. The candidate sample $\dv^*$ is accepted $\dv_i = \dv^*$ with the probability $\alpha(\dv^*, \dv_{i - 1})$, which depends on the ratio of the likelihood at the candidate sample $\dv^*$ and the sample $\dv_{i - 1}$ of the previous iteration. If the candidate sample is rejected, then $\dv_i = \dv_{i - 1}$. Algorithm~\ref{alg:MCMC-MH} returns the samples $\dv_1, \dots, \dv_{\nrpts}$.

\paragraph{Metropolis-Hastings algorithm in practice}
Several techniques are necessary to make the Metropolis-Hastings algorithm practical. For example, the samples generated in the first iterations are usually discarded (burn-in) because they are strongly influenced by the initial starting point $\dv_0$. We refer to the literature \cite{RobertBook} for more details and further practical considerations.

\paragraph{Efficiency of MCMC sampling}
Once samples $\dv_1, \dots, \dv_{\nrpts}$ are drawn with an MCMC algorithm, they can be used to, e.g., estimate the expectation \eqref{eq:MM:InvProb:ExpH} as
\[
\bar{h} = \frac{1}{\nrpts}\sum_{i = 1}^{\nrpts} h(\dv_i)\,.
\]
MCMC generates correlated samples $\dv_1, \dots, \dv_{\nrpts}$. The efficiency of MCMC sampling can therefore be measured by the effective sample size for a given computational budget with respect to the estimator $\bar{h}$. To define the effective sample size, consider the samples $\dv_1, \dv_2, \dv_3, \dots$ drawn with an MCMC algorithm and define the integrated autocorrelation time as
\[
\tau_{\text{int}}(h) = \frac{1}{2} + \sum_{j = 1}^{\infty}\rho_j\,,
\]
where $\rho_j = \operatorname{corr}(h(\dv_1), h(\dv_{j + 1}))$ is the correlation between $h(\dv_1)$ and $h(\dv_{j + 1})$. The effective sample size $\nrpts_{\text{eff}}(h)$ is
\[
\nrpts_{\text{eff}}(h) = \frac{\nrpts}{2 \tau_{\text{int}}(h)}\,,
\]
such that
\[
\Var[\bar{h}] \approx \frac{\Var[h]}{\nrpts_{\text{eff}}(h)}\,,
\]
see \cite[p.~125f]{MCBookLiu} for a detailed derivation and further references.
This means that there are at least two ways to improve the efficiency of sampling with MCMC \cite{cui_data-driven_2014}: (1)~Increase the effective sample size for a given number of MCMC iterations with, e.g., adaptive MCMC \cite{gilks_adaptive_1994,gilks_adaptive_1998,roberts_coupling_2007,haario_adaptive_2001}. (2)~Increase the number of MCMC iterations for a given computational budget, so that more samples can be generated for a given budget with, e.g., two-stage MCMC \cite{christen_markov_2005,Fox97samplingconductivity}.

\subsection{Two-stage Markov chain Monte Carlo}
\label{sec:MM:InvProb:TwoStageMCMC}
Two-stage MCMC methods aim to increase the number of MCMC iterations for a given computational budget. In many applications, the Metropolis-Hastings algorithm, and MCMC in general, requires many iterations to produce an acceptable effective sample size. Each iteration means a likelihood evaluation, which means a high-fidelity model $\fhigh$ evaluation in the case of the likelihood $L$ as defined in \eqref{eq:MM:InvProb:LikelihoodHi}. Two-stage MCMC methods employ delayed acceptance or delayed rejection strategies that use multifidelity models to reduce the number of samples evaluated using the expensive high-fidelity model.

The work \cite{christen_markov_2005,Fox97samplingconductivity} proposes a two-stage delayed acceptance MCMC sampling. The candidate sample $\dv^*$ has to be accepted with the likelihood induced by a low-fidelity model first, before $\dv^*$ is passed on to be either accepted or rejected with the likelihood induced by the high-fidelity model.

Multi-stage MCMC methods based on delayed rejection, in contrast to delayed acceptance, have been proposed in \cite{tierney_adaptive_1999,mira_metropolis-hastings_2001,green_delayed_2001}. A candidate sample that is rejected in one stage is retried with a different proposal distribution at a subsequent stage. Ref.~\cite{tierney_adaptive_1999} suggests using an independence sampler based on a density that is thought to be a good approximation of the posterior density on the first stage. In case the density corresponding to the independence sampler is a poor approximation, a random walk is used on the second stage. Another approach discussed in \cite{tierney_adaptive_1999} is to use a proposal based on a local quadratic approximation of the log posterior density, which is restricted to smaller and smaller neighborhoods of the current point in a trust-region fashion. In \cite{green_delayed_2001}, the proposal at the first stage is a normal distribution. The proposal at the second stage is a normal distribution with the same mean but with a higher variance. In principle, the proposals in the context of delayed rejection MCMC, just as in delayed acceptance MCMC, can also be derived from models with different fidelities and costs (although we are not aware of such implementations in the literature that pre-date the work by \cite{christen_markov_2005,Fox97samplingconductivity}).

\paragraph{Two-stage delayed acceptance MCMC algorithm}
The two-stage MCMC method introduced in \cite{christen_markov_2005} is summarized in Algorithm~\ref{alg:MCMC-2Stage}. Inputs are the likelihood $L^{\hi}$ and $L^{\lo}$, corresponding to the high- and low-fidelity model respectively, the prior density $p_0$, the proposal $q$, and the number of samples $\nrpts$. Consider one of the iterations $i = 1, \dots, \nrpts$. The first stage of Algorithm~\ref{alg:MCMC-2Stage} (lines 4--6) proceeds as the Metropolis-Hastings algorithm with a candidate sample $\dv^{\lo}$ drawn from the proposal distribution, except that the likelihood of the low-fidelity model is used, instead of the likelihood of the high-fidelity model. The result of the first stage is a sample $\dv^{\hi}$, which either is $\dv^{\hi} = \dv^{\lo}$ (accept) or $\dv^{\hi} = \dv_{i - 1}$ (reject). In the second stage of the algorithm (lines 7--9), the high-fidelity model is used to either accept or reject the candidate sample $\dv^{\hi}$ of the first stage. If the sample $\dv^{\hi}$ is accepted in the second stage, then the algorithm sets $\dv_i = \dv^{\hi}$. If the sample $\dv^{\hi}$ is rejected in the second stage, then $\dv_i = \dv_{i - 1}$. Note that in case the first stage rejected the sample $\dv^{\lo}$, i.e., $\dv^{\hi} = \dv_{i - 1}$, no high-fidelity model evaluation is necessary in the second stage because the high-fidelity model output at $\dv_{i - 1}$ is available from the previous iteration.
The proposal distribution $Q$ in the second stage depends on the likelihood $L^{\lo}$ of the low-fidelity model. Note that $\delta_{\dv_{i - 1}}$ is the Dirac mass at $\dv_{i - 1}$. Algorithm~\ref{alg:MCMC-2Stage} returns the samples $\dv_1, \dots, \dv_{\nrpts}$.

\paragraph{Efficiency of two-stage delayed acceptance MCMC}
The key to the efficiency of Algorithm~\ref{alg:MCMC-2Stage} is that no high-fidelity model evaluation is necessary in the second stage if the candidate sample $\dv^{\lo}$ was rejected at the first stage. Since the rejection rate of MCMC is typically high, many high-fidelity model evaluations are saved by the two-stage MCMC. We refer to \cite{christen_markov_2005} for an asymptotic analysis and the convergence properties of the two-stage MCMC. The two-stage MCMC uses a model management based on filtering because only candidate samples accepted with the low-fidelity model are passed on to the high-fidelity model, see Section~\ref{sec:MM:strategies}.

\paragraph{Other multifidelity MCMC approaches}
In \cite{efendiev_preconditioning_2006}, two-stage MCMC is seen as a preconditioned Metropolis-Hastings algorithm. The low-fidelity model in \cite{efendiev_preconditioning_2006} is a coarse-scale model of a high-fidelity multiscale finite volume model. We also mention \cite{higdon_bayesian_2002}, where multiple MCMC chains from coarse- (low-fidelity) and fine-scale (high-fidelity) models are coupled using a product chain. The work \cite{MultiMCMC} couples multilevel Monte Carlo and MCMC to accelerate the estimation of expected values with respect to a posterior distribution. The low-fidelity models form a hierarchy of coarse-grid approximations of the high-fidelity model.

\begin{algorithm}[t]
\caption{Two-stage MCMC}\label{alg:MCMC-2Stage}
\begin{algorithmic}[1]
\Procedure{2StageMCMC}{$L^{\hi}, L^{\lo}, p_0, q, \nrpts$}
\State Choose a starting point $\dv_0$
\For{$i = 1, \dots, \nrpts$}
\State Draw candidate $\dv^{\lo}$ from proposal $q(\cdot | \dv_{i - 1})$
\State Compute acceptance probability
\[
\alpha^{\lo}(\dv_{i - 1}, \dv^{\lo}) = \min\left\{1, \frac{q(\dv_{i - 1} | \dv^{\lo})L^{\lo}(\obs|\dv^{\lo})p_0(\dv^{\lo})}{q(\dv^{\lo}|\dv_{i - 1})L^{\lo}(\obs|\dv_{i - 1})p_0(\dv_{i-1})}\right\}
\]
\State Set the candidate sample $\dv^{\hi}$ to
\[
\dv^{\hi} = \begin{cases}
\dv^{\lo}\,,\qquad &\text{with probability }\alpha^{\lo}(\dv_{i - 1}, \dv^{\lo})\,,\\
\dv_{i - 1}\,,\qquad &\text{with probability }1 - \alpha^{\lo}(\dv_{i - 1}, \dv^{\lo})
\end{cases}
\]
\State Set the distribution $Q$ to
\begin{multline*}
Q(\dv | \dv_{i - 1}) = \alpha^{\lo}(\dv_{i - 1}, \dv)q(\dv|\dv_{i - 1}) + \\\left(1 - \int_{\Zcal}\alpha^{\lo}(\dv_{i - 1}, \dv)q(\dv | \dv_{i - 1})\mathrm d\dv\right)\delta_{\dv_{i - 1}}(\dv)
\end{multline*}
\State Compute acceptance probability
\[
\alpha^{\hi}(\dv_{i - 1}, \dv^{\hi}) = \min\left\{1, \frac{Q(\dv_{i - 1} | \dv^{\hi})L^{\hi}(\obs|\dv^{\hi})p_0(\dv^{\hi})}{Q(\dv^{\hi}|\dv_{i - 1})L^{\hi}(\obs|\dv_{i - 1})p_0(\dv_{i-1})}\right\}
\]
\State Set sample $\dv_i$ to
\[
\dv_i = \begin{cases}
\dv^{\hi}\,,\qquad &\text{with probability }\alpha^{\hi}(\dv_{i - 1}, \dv^{\hi})\,,\\
\dv_{i - 1}\,,\qquad &\text{with probability }1 - \alpha^{\hi}(\dv_{i - 1}, \dv^{\hi})
\end{cases}
\]
\EndFor
\State \Return $\dv_1, \dots, \dv_{\nrpts}$
\EndProcedure
\end{algorithmic}
\end{algorithm}

\subsection{Markov chain Monte Carlo with adaptive low-fidelity models}
\label{sec:MM:InvProb:Adaptive}
In \cite{cui_data-driven_2014} an algorithm for combining high- and low-fidelity models in MCMC sampling is presented. The low-fidelity model is used in the first step of a two-stage MCMC approach to increase the acceptance rate of candidates in the second step, where the high-fidelity model is used to either accept or reject the sample. Additionally, the high-fidelity model outputs computed in the second step are used to adapt the low-fidelity model. In that sense, the approach uses a model management based on a combination of adaptation and filtering.

The low-fidelity model is a projection-based reduced model in \cite{cui_data-driven_2014}, see Section~\ref{sec:intro:models}. The low-fidelity model is constructed in an offline phase from an initial reduced basis. At each MCMC iteration, the low-fidelity model is used in a first stage to generate a certain number of samples with the Metropolis-Hastings algorithm. The goal is to generate so many samples with the low-fidelity model that the initial sample and the last sample are almost uncorrelated. Then, in the second stage of the MCMC iteration, the last sample generated with the low-fidelity model is used as a candidate sample and the acceptance probability is computed using the high-fidelity model. Similarly to the two-stage MCMC methods discussed above, this algorithm aims to increase the acceptance probability at the second stage; however, the high-fidelity model output is used to improve the low-fidelity model after a sample has been accepted, and is not discarded. In that way, the authors of \cite{cui_data-driven_2014} improve the low-fidelity model during the MCMC iterations, and consequently increase the acceptance probability in the second stage of the MCMC iterations. The authors show that under certain technical conditions, the Hellinger distance $d_{\text{Hell}}(p_{\hi}, p_{\lo})$ between the posterior distribution $p_{\hi}$ corresponding to the high-fidelity model $\fhigh$ and the approximate posterior distribution $p_{\lo}$ induced by the reduced model $\fl$ decreases as the reduced model $\fl$ is improved during the MCMC iterations. Thus, for any $\epsilon > 0$, the Hellinger distance between the high-fidelity and low-fidelity posterior distribution can be bounded
\[
d_{\text{Hell}}(p_{\hi}, p_{\lo}) \leq \epsilon
\]
after a sufficiently large number of MCMC iterations. The resulting low-fidelity model is data-driven because it uses information provided by the observation $\obs$, rather than only prior information. Using samples from this adaptive two-stage MCMC approach yields unbiased Monte Carlo estimators, see the analysis in \cite{cui_data-driven_2014}.

\subsection{Bayesian estimation of low-fidelity model error}
\label{sec:MM:InvProb:Error}
In \cite{KaipoBook,kaipio_statistical_2007}, a multifidelity approach for Bayesian inference is presented that relies on a Bayesian approximate error model of the low-fidelity model derived from the high-fidelity model. The inference is performed with the low-fidelity model but an additional term is introduced that quantifies the error between the low-fidelity and the high-fidelity model. The error term relies on the high-fidelity model and is used to correct for the error introduced by the low-fidelity model. The approach follows the Bayesian paradigm and considers the error term as a random variable, which becomes another noise term in the Bayesian formulation of the inverse problem. Within our categorization, this approach uses a model management based on model adaptation.

In \cite{KaipoBook,kaipio_statistical_2007}, the relationship $\obs = \fhigh(\dv) + \bfeps$ is formulated as
\[
\obs = \fl(\dv) + \left(\fhigh(\dv) - \fl(\dv)\right) + \bfeps\,,
\]
which can be rewritten as
\[
\obs = \fl(\dv) + \bfe(\dv) + \bfeps\,,
\]
with the error term $\bfe(\dv)$. If $\fhigh$ and $\fl$ are linear in $\dv$, the sum of the covariance of the error $\bfe(\dv)$ and noise term $\bfeps$ can derived explicitly under certain assumptions \cite{KaipoBook,kaipio_statistical_2007}. If $\fhigh$ and $\fl$ are nonlinear, a Gaussian approximation of the covariance is proposed in \cite{KaipioErrorDOT}.

The Bayesian approximate error model is used in \cite{KaipioErrorDOT,kolehmainen_approximation_2012} for a diffuse optical tomography problem, where the low-fidelity model corresponds a coarse-grid finite-element approximation and the high-fidelity model is a fine-grid approximation. In \cite{tarvainen_approximation_2010,tarvainen_corrections_2010}, diffuse optical tomography is considered with a simplified-physics low-fidelity model. The error introduced by the low-fidelity model is corrected with the Bayesian approximate error modeling approach. Correcting for errors due to truncation of the computational domain and unknown domain shape are investigated in \cite{lehikoinen_approximation_2007,nissinen_compensation_2011}. The Bayesian approximate error model approach was verified on real data in \cite{nissinen_bayesian_2008}. Extensions to time-dependent problems are presented in \cite{huttunen_approximation_2007,huttunen_approximationB_2007}.

We note that estimating the error term $\bfe$ is related to estimating the model inadequacy of the low-fidelity model with respect to the high-fidelity model. Often, Gaussian-process models are used to attempt to correct model inadequacy \cite{kennedy_bayesian_2001,kennedy_bayesian_2002,higdon_combining_2004,oakley_bayesian_2002,bayarri_framework_2007}. The corrected low-fidelity model is typically cheap to evaluate, because the Gaussian-process model is cheap, and so the corrected low-fidelity model can be used in MCMC sampling with the aim of decreasing the costs per MCMC iteration compared to using the high-fidelity model. The work \cite{KaipoBook,kaipio_statistical_2007} can be seen as dealing with model inadequacy in the specific context of inverse problems.

config: arxiv


text:
\section{Locally-implicit Lax-Wendroff discontinuous Galerkin}
\label{sec:LxW-DG}
The Lax-Wendroff method \cite{article:LxW60} is a fully discrete numerical method for
hyperbolic conservation laws (e.g., of the form \cref{eqn:conslaw}) based on
the the Cauchy-Kovalevskaya
\cite{article:Kovaleskaya1875} procedure to convert temporal derivatives into spatial derivatives.
For example, in the case of conservation law \cref{eqn:conslaw},
we begin with a Taylor series in time:
\begin{equation}
 \vec{q}(t+\Delta t, x) = \vec{q}(t,x) + \Delta t \vec{q}_{,t}(t,x) + \frac{1}{2} \Delta t^2 \vec{q}_{,t,t}(t,x) +
 \ldots,
\end{equation}
and then replace all time derivatives by spatial derivatives:
\begin{equation}
\label{eqn:LxW_time_derivs}
 \vec{q}_{,t} = - \vec{f}\left(\vec{q}\right)_{,x}, \quad
 \vec{q}_{,t,t} = - \vec{f}\left(\vec{q}\right)_{,t,x} = -\left[ \mat{f'}\left(\vec{q}\right) \, \vec{q}_{,t} \right]_{,x} =
 \left[ \mat{f'}\left(\vec{q}\right) \vec{f}\left(\vec{q}\right)_{,x} \right]_{,x}, \quad \ldots,
\end{equation}
which results in the following:
\begin{equation}
\vec{q}(t+\Delta t ,x) = \vec{q} - \Delta t \vec{f}\left(\vec{q}\right)_{,x} + \frac{1}{2} \Delta t^2 \left[ \mat{f'}\left(\vec{q}\right) \vec{f}\left(\vec{q}\right)_{,x} \right]_{,x}
+ \ldots,
\end{equation}
where on the right-hand side we have suppressed the evaluation at $(t,x)$.
The final step is to truncate the Taylor series at some finite number of terms,
and then replace all spatial derivatives by some discrete spatial derivative operators.
The above Lax-Wendroff formalism \cite{article:LxW60} has been used in conjunction with
with a variety of spatial discretizations, including finite volume  \cite{article:Le97},
weighted essentially non-oscillatory (WENO) \cite{article:TitaToro02}, and discontinuous Galerkin \cite{article:Qiu05} operators.

In this work, we are concerned with the discontinuous Galerkin version
of Lax-Wendroff \cite{article:Qiu05}; and in particular, we make use
of the reformulation of Gassner et al. \cite{article:GasDumHinMun2011} of the Lax-Wendroff discontinuous Galerkin (LxW-DG) scheme
in terms of a locally-implicit prediction step, followed by an explicit correction step.
The key advantage of this formulation is that we do not need to explicitly compute the
partial derivatives as shown in \cref{eqn:LxW_time_derivs}; and instead, the locally-implicit
solver automatically produces discrete versions of these derivatives.
The next challenge is to efficiently solve the nonlinear algebraic equations that arise from the
locally-implicit prediction step; we solve these equations by again following Gassner et al. \cite{article:GasDumHinMun2011} and making use of a Picard fixed point iteration.
One key difference in this work is that we formulate the prediction step in terms of primitive
variables, which pays dividends when we develop limiters (see \cref{sec:limiters}).
We will follow the notational conventions of Guthrey and Rossmanith \cite{article:GuRo2017}
developed for locally-implicit and regionally-implicit LxW-DG schemes.

In the remainder of this section we develop the basic locally-implicit Lax-Wendroff
discontinuous Galerkin (LxW-DG) scheme. We begin with a quick review of
the DG spatial discretization \cref{sec:dg-fem-space},
followed by a detailed description of the prediction (\cref{subsec:prediction_step}) and
correction (\cref{subsec:corr_step}) phases in one time-step of the proposed locally-implicit LxW-DG method.
The discussion in this section is on the scheme {\it without} limiters;
limiters will receive our full attention in \cref{sec:limiters}.

\subsection{\it DG-FEM spatial discretization}
\label{sec:dg-fem-space}
We discretize system \cref{eqn:conslaw} in space via the discontinuous Galerkin (DG) method, which
was first introduced by Reed and Hill \cite{article:ReedHill73} for neutron transport, and then fully developed for time-dependent hyperbolic conservation laws in a series of papers by Bernardo Cockburn, Chi-Wang Shu, and collaborators (see \cite{article:CoShu98} and references therein for details).

The computational domain
is a finite interval on the real line: $\Omega = [x_{\text{low}},x_{\text{high}}] \subset \reals$.
Let ${\mathbb P}\left(\mdeg \right)$ denote the set of polynomials from $\reals$ to $\reals$
with maximal polynomial degree $\mdeg$.
On the mesh of $\melems$ elements we define the {\it broken} finite element space:
\begin{equation}
\label{eqn:broken_space}
    \WS^{h} := \left\{ \vec{w}^h \in \left[ L^{\infty}(\Omega) \right]^{\meq}: \,
    \vec{w}^h \bigl|_{\Tm_i} \in \left[ {\mathbb P} \left(\mdeg \right) \right]^{\meq} \, \, \forall \Tm_i \right\},
\end{equation}
where $h:=\Delta x = (x_{\text{high}}-x_{\text{low}})/\melems$ is the uniform grid spacing,
$\meq$ is the number of conserved variables, and $\mdeg$ is the
maximal polynomial degree in the finite element representation.
The above expression means that $\vec{w} \in \WS^h$ has
$\meq$ components, each of which when restricted to some element $\Tm_i$
is a polynomial in ${\mathbb P}\left(\mdeg \right)$, and no continuity is assumed
across element faces.

The computational mesh is comprised of $\melems$ elements, which we denote by
\begin{equation}
{\Tm_i} = \left[ x_{i} - \frac{\Delta x}{2},  x_{i} + \frac{\Delta x}{2} \right] \quad \text{for} \quad
i=1,\ldots,\melems.
\end{equation}
For convenience, we define the canonical variable, $\xi$, which on each element is related
the physical coordinate, $x$, as follows:
\begin{equation}
\label{eqn:local_xi}
x = x_{i} + \left( \frac{\Delta x}{2} \right) \xi, \qquad \xi \in [-1,1].
\end{equation}
Using the canonical variable, we define on each element the following Legendre polynomial
basis:
\begin{equation}
\label{eqn:phi_basis_1d}
\vec{\Phi} = \left( 1, \, \sqrt{3} \xi, \, \frac{\sqrt{5}}{2} \left( 3 \xi^2 - 1 \right), \, \frac{\sqrt{7}}{2}\left(5 \xi^3 - 3 \xi \right), \,
\frac{\sqrt{9}}{8} \left(35 \xi^4 - 30 \xi^2 + 3 \right), \, \cdots \right),
\end{equation}
which also satisfies the following three-term recurrence relationship
for $k\ge3$:
\begin{equation}
\Phi_k(\xi) = \left( \frac{\sqrt{(2k-3)(2k-1)}}{(k-1)} \right) \xi \Phi_{k-1}(\xi)
       - \left( \frac{(k-2)\sqrt{2k-1}}{(k-1)\sqrt{2k-5}} \right) \Phi_{k-2}(\xi),
\end{equation}
where $\Phi_1 = 1$ and $\Phi_2 = \sqrt{3} \xi$.

The approximate solution on each element at time $t=t^n$:
\begin{equation}
\label{eqn:corrector_ansatz}
   \vec{q}^h\left(t^n, \, x_i +  \frac{\Delta x}{2} \, \xi \right) := \vec{\Phi}\left(\xi \right)^T \mat{Q_i^n},
   \quad \text{for} \quad \xi \in [-1,1],
\end{equation}
where
\begin{align}
\vec{\Phi}\left(\xi \right): [-1,1] \mapsto \reals^{\mcorr} \qquad \text{and} \qquad
\mat{Q_i^n} \in \reals^{\mcorr \times \meq}.
\end{align}
Note that we are denoting the number of basis functions by $\mcorr$ to signify that this refers to the number of basis functions in the {\it correction step}.
If the exact solution, $\vec{q}(t^n,x)$, is known {\it a priori}, then we can compute the corresponding coefficients, ${Q}$, via $L^2$-projection:
\begin{equation}
\mat{Q^n_i} = \frac{1}{2} \int_{-1}^{1} \vec{\Phi}(\xi) \left[ \vec{q}\left(t^n, x_i + \frac{\Delta x}{2} \, \xi \right) \right]^T \, d\xi
\approx \frac{1}{2} \sum_{a=1}^{\morder} \omega_a \, \vec{\Phi}(\mu_a) \left[ \vec{q}\left(t^n, x_i + \frac{\Delta x}{2} \mu_a \right) \right]^T,
\end{equation}
where $\morder = \mcorr$ is the maximum achievable order of accuracy, and $\omega_a$ and $\mu_a$ for $a=1,\ldots,\morder$ are the weights and abscissas of the $\morder$-point Gauss-Legendre quadrature rule.
In practice, the only solution that is known {\it a priori} is the initial condition: $\vec{q}(0,x)$.
For all subsequent time, the coefficients ${Q}$ must be computed by a numerical procedure, which in
in this work will be the locally-implicit Lax-Wendroff scheme described in detail below.

\subsection{\it Prediction step}
\label{subsec:prediction_step}
The numerical update of the proposed scheme is divided into two distinct parts: (1) the prediction step and (2)
the correction step. In the prediction step we will {\it not} enforce consistency of the numerical method
with the underlying conservation laws \cref{eqn:conslaw}; and therefore, we have a significant amount of freedom
in how this portion of the update can be accomplished. In particular, one freedom which we will exercise
is the choice of variables used in the prediction step (e.g., conservative, primitive, or entropy variables). For simplicity of discussion, we will simply refer to the choice of variables as the {\it primitive variables}. We denote these variables by $\alpha$ and assume (without loss of generality) that they satisfy quasilinear equation \cref{eqn:quasilinear_prim}.

The prediction step is entirely local on each element; and therefore, without loss of generality, we focus our attention on element $\Tm_i$ over time interval $[t^n, t^{n+1}]$, where $t^{n+1} = t^n + \Delta t$. On this
{\it space-time} element, we introduce the local spatial variable $\xi$ as defined by \cref{eqn:local_xi}, the local temporal variable
$\tau$ as defined by
\begin{equation}
t = t^{n}+\frac{\Delta t}{2} (1+\tau), \qquad \tau \in [-1,1],
\end{equation}
and rewrite \cref{eqn:quasilinear_prim} as follows:
\begin{equation}
\label{eqn:alpha_theta}
\vec{\alpha}_{,\tau} = \vec{\Theta}\left(\vec{\alpha} \right) :=  - \nu \mat{B}\left(\vec{\alpha} \right) \, \vec{\alpha}_{,\xi}, \qquad \nu:=\frac{\Delta t}{\Delta x}.
\end{equation}
This equation can also be written in component form as
\begin{equation}
\label{eqn:alpha_theta_comp}
{\alpha}_{m,\tau} = {\Theta}_m\left(\vec{\alpha} \right) =
- \nu \sum_{k=1}^{\meq} {B}_{mk} \left(\vec{\alpha} \right) \, {\alpha}_{k,\xi} \quad \text{for} \quad m=1,\ldots,\meq.
\end{equation}

We introduce a space-time Legendre basis on each element:
\begin{equation}
\label{eqn:spacetime_basis}
\vec{\Psi}(\tau, \xi): [-1,1]^2 \mapsto \reals^{\mpred}, \quad \mpred = \frac{\mcorr(\mcorr+1)}{2}, \quad {\Psi}_{\ell}(\tau, \xi) = {\Phi}_{\ell_\tau}(\tau) \, {\Phi}_{\ell_\xi}(\xi),
\end{equation}
which is orthonormal on $[-1,1]^2$:
\begin{equation}
\label{eqn:psi_orthonormality}
\frac{1}{4} \int_{-1}^{1} \int_{-1}^{1} \vec{\Psi} \, \vec{\Psi}^T \, d\tau \, d\xi =
  \mat{\mathbb I} \in \reals^{\mpred\times\mpred}.
\end{equation}
Note that we are denoting the number of space-time basis functions by $\mpred$ to signify that this refers to the number of basis functions in the {\it prediction step}.
We catalog, at least up to fifth-order accuracy, how the one-dimensional indices $\ell_{\tau}$ and $\ell_{\xi}$ vary with
the index $\ell$ in \cref{table:spacetime_index}.

Using these basis functions, we write the predicted solution as follows:
\begin{equation}
\label{eqn:pred_ansatz}
\vec{\alpha}^{\text{ST}}\left(t^{n}+\frac{\Delta t}{2} (1+\tau), \, x_i +\frac{\Delta x}{2} \, \xi \right)   := \vec{\Psi}\left(\tau,\xi\right)^T \mat{W^{n+\half}_i},
\qquad \mat{W_i^{n+\half}} \in \reals^{\mpred \times \meq},
\end{equation}
for $(\tau,\xi) \in [-1,1]^2$, where $\mat{W}$ represents the matrix of unknown  coefficients.

\begin{table}[!t]
\begin{center}
\begin{Large}
\begin{tabular}{|cc||cc||cc||cc||cc|}
\hline
{\normalsize $\ell$} & {\normalsize $\left(\ell_\tau, \ell_\xi \right)$} &
{\normalsize $\ell$} & {\normalsize $\left(\ell_\tau, \ell_\xi \right)$} &
{\normalsize $\ell$} & {\normalsize $\left(\ell_\tau, \ell_\xi \right)$} &
{\normalsize $\ell$} & {\normalsize $\left(\ell_\tau, \ell_\xi \right)$} &
{\normalsize $\ell$} & {\normalsize $\left(\ell_\tau, \ell_\xi \right)$} \\
\hline\hline
{\normalsize 1} & {\normalsize (1,1)} & {\normalsize 4} & {\normalsize (1,3)} & {\normalsize 7} & {\normalsize (1,4)} & {\normalsize 10} & {\normalsize (4,1)} & {\normalsize 13} & {\normalsize (3,3)} \\
{\normalsize 2} & {\normalsize (1,2)} & {\normalsize 5} & {\normalsize (2,2)} & {\normalsize 8} & {\normalsize (2,3)} & {\normalsize 11} & {\normalsize (1,5)} & {\normalsize 14} & {\normalsize (4,2)} \\
{\normalsize 3} & {\normalsize (2,1)} & {\normalsize 6} & {\normalsize (3,1)} & {\normalsize 9} & {\normalsize (3,2)} & {\normalsize 12} & {\normalsize (2,4)} & {\normalsize 15} & {\normalsize (5,1)} \\
\hline
\end{tabular}
\caption{Index conversion table for the space-time basis, $\psi$, as defined in
\cref{eqn:spacetime_basis}. Shown here are the basis element up to fifth-order of accuracy.
\label{table:spacetime_index}}
\end{Large}
\end{center}
\end{table}

Before describing how to compute the space-time coefficients, $
\mat{W}$, for the primitive variables, $\vec{\alpha}$,
we need to address one small issue: before the prediction step, the solution at time $t=t^n$ is given
only in terms of conservative variables (see equation \cref{eqn:corrector_ansatz}).
In order to convert the conservative variable coefficients  from \cref{eqn:corrector_ansatz}
to primitive variable coefficients,
\begin{equation}
\label{eqn:prim_ic}
   \vec{\alpha}^h\left(x_i +  \frac{\Delta x}{2} \, \xi \right) := \vec{\Phi}\left(\xi \right)^T \mat{A_i^n}
   \quad \text{for} \quad \xi \in [-1,1], \quad i=1,\ldots,\melems,
\end{equation}
we apply a simple $L^2[-1,1]$ projection:
\begin{equation}
\mat{A_i^n} = \frac{1}{2} \sum_{a=1}^{\morder} \omega_a \, \vec{\Phi}\left(\mu_a \right)  \left[ \vec{\alpha}\left(
\vec{\Phi}\left(\mu_a \right)^T \mat{Q_i^n} \right) \right]^T \in \reals^{\mcorr \times \meq},
\end{equation}
where $\vec{\alpha}(\vec{q}): \reals^\meq \mapsto \reals^\meq$ gives the relationship between conservative
and primitive variables.

Next, an algebraic equation for the solution of the unknown coefficients in ansatz \cref{eqn:pred_ansatz} is
obtained by multiplying \cref{eqn:alpha_theta} by $\vec{\Psi}$, integrating over
$(\tau,\xi) \in [-1,1]^2$,  integrating-parts {\it only} in $\tau$ and {\it not} in $\xi$, and
making use of ansatz \cref{eqn:pred_ansatz}:
\begin{equation}
\label{eqn:nonlinear_pred}
\begin{split}
\mat{L} \, \, \vec{W^{n+\half}_{i \, (:,m)}} &= \frac{1}{4}\int_{-1}^{1} \int_{-1}^{1}
  {\Theta}_m\left( \vec{\Psi}^T \mat{W^{n+\half}_i} \right) \vec{\Psi} \, d\tau \, d\xi \\ &+
  \left[ \frac{1}{4} \int_{-1}^{1} \, \vec{{\Psi}} \left(-1,\xi \right)  \, \vec{\Phi} \left(\xi \right)^T \, d\xi
 \right] \, \vec{A^n_{i \, (:,m)}},
 \end{split}
\end{equation}
for $m=1,\ldots,\meq$, where
\begin{gather}
\vec{W^{n+\half}_{i \, (:,m)}} \in \reals^\mpred : \, \text{$m^{\text{th}}$
column of $\mat{W^{n+\half}_i}$} \quad \text{(all coefficients for equation $m$)}, \\
\vec{A^{n}_{i \, (:,m)}} \in \reals^\mcorr : \, \text{$m^{\text{th}}$
column of $\mat{A^{n}_i}$} \quad \text{(all coefficients for equation $m$)}, \\
\mat{L} := \frac{1}{4} \int_{-1}^{1} \int_{-1}^{1} \vec{\Psi} \, \vec{\Psi}^T_{,\tau} \, d\tau \, d\xi
+ \frac{1}{4} \int_{-1}^{1} \vec{\Psi}_{|_{\tau=-1}} \vec{\Psi}_{|_{\tau=-1}}^T \, d\xi \in \reals^{\mpred \times \mpred}.
\end{gather}

We note that system \cref{eqn:nonlinear_pred} represents something akin to a single block-Jacobi iteration of a fully implicit spacetime DG approach \cite{article:KlaVegVen2006,article:Sudirham2006}, and is a set of  nonlinear algebraic equations that must be solved
independently on each space-time element. There are several techniques, including Newton's method for systems, that could be used solve these equations. However, following Gassner et al. \cite{article:GasDumHinMun2011}, we make use of an even simpler fixed-point iteration: the {\it Picard} iteration.

After replacement of the
space-time integration with Gauss-Legendre quadrature, we can write the Picard iteration as
\begin{equation}
\label{eqn:picard}
\begin{split}
\vec{W^{n+\half}_{i \, (:,m)}} \leftarrow & \frac{1}{4} \sum_{a=1}^{\morder}
\sum_{b=1}^{\morder} \, \omega_a  \, \omega_b \,
  \vec{\hat{\Psi}}\left(\mu_b, \, \mu_a \right)
   \, {\Theta}_m\left( \vec{\Psi}\left(\mu_b, \mu_a \right)^T \mat{W^{n+\half}_i} \right) \\
   + & \frac{1}{4} \sum_{b=1}^{\morder} \omega_{b} \, \vec{\hat{\Psi}} \left(-1,\xi_b \right)  \, \vec{\Phi} \left(\xi_b \right)^T \vec{A^{n}_{i \, (:,m)}},
   \end{split}
\end{equation}
for $m=1,\ldots,\meq$, where $\vec{\hat{\Psi}} = \mat{L}^{-1} \vec{\Psi}$ and
 $\omega_a$ and $\mu_a$ for $a=1,\ldots,\morder$ are the weights and abscissas of the $\morder$-point Gauss-Legendre quadrature rule.
This iteration -- like all fixed point iterations -- requires some appropriate initial guess; we explain how this
is done in \cref{subsec:full_algorithm}.

The two main advantages of the Picard iteration over Newton's method are:
(1) it is Jacobian-free (the only inverse that must be
computed is of $\mat{L}$, which is independent of the solution); and (2) the iteration converges to sufficient high-order accuracy after exactly $\morder$ iterations, obviating the need to compute residuals.
The basic principle of this approach is that each iteration improves the quality of the guess by one order of accuarcy: one iteration gives first-order accuracy, two iterations gives second-order
accuracy, etc$\ldots$,
up to the maximum possible order of accuracy\footnote{Technically, the predicted solution is not even consistent with the underlying partial differential equation. What we mean by
{\it high-order accuracy} in this context is what happens when the predicted solution is fed to the correction step. See \cref{subsec:corr_step} for more details.}: $\morder$ \cite{article:GasDumHinMun2011}.

\subsection{\it Correction step}
\label{subsec:corr_step}
The prediction step as outlined above is clearly not sufficient to produce a consistent numerical approximation
of hyperbolic conservation law \cref{eqn:conslaw}. Without having done a proper integration-by-parts in the spatial variable, the predicted solution on each space-time element is completely decoupled from all other space-time elements, which is inconsistent with the underlying partial differential equation.
Fortunately, there is a simple remedy that makes the solution not only {\it consistent} with conservation law \cref{eqn:conslaw}, but in fact  {\it high-order accurate} (under the assumption of sufficiently smooth solutions).
We refer to this remedy as the {\it correction step}, which is a single forward Euler-like step that makes use of the predicted solution.

To enact the correction step, we take hyperbolic conservation law \cref{eqn:conslaw}, multiply
by the spatial basis functions $\vec{\Phi}$ (see \cref{eqn:phi_basis_1d}),
integrate over $(\tau,\xi) \in [-1,1]^2$, make use of ansatz \cref{eqn:corrector_ansatz},
apply integration-by-parts on the spatial variable, and replace all exact integration
by Gauss-Legendre quadrature:
\begin{gather}
\label{eqn:correction_update}
\begin{split}
\mat{Q^{n+1}_i} = \mat{Q^{n}_i} &+  \frac{\nu}{2} \sum_{a=1}^{\morder}
\sum_{b=1}^{\morder} \, \omega_a  \, \omega_b \, \vec{\Phi}_{,\xi}\left(\mu_a \right) \left[
 \vec{f} \left( \vec{\Psi}\left(\mu_b,\mu_a\right)^T \mat{W^{n+\half}_{i}} \right) \right]^T \\
&- {\nu}\left( \, \vec{\Phi}(1) \left[ \vec{{\mathcal F}^{n+\half}_{i+\half}} \right]^T
-  \vec{\Phi}(-1) \left[ \vec{{\mathcal F}^{n+\half}_{i-\half}} \right]^T \, \right),
\end{split}
\end{gather}
where $\omega_a$ and $\mu_a$ for $a=1,\ldots,\morder$ are the weights and abscissas of the $\morder$-point Gauss-Legendre quadrature rule.
The time-integrated numerical fluxes are defined using the predicted solution and the Rusanov \cite{article:Ru61} time-averaged flux:
\begin{gather}
\label{eqn:correction_update_flux_1}
  \vec{{\mathcal F}^{n+\half}_{i-\half}} := \frac{1}{2} \sum_{a=1}^{\morder}
  \omega_a \, \vec{\mathcal F} \left( \mu_a \right),
\end{gather}
where the numerical flux at each temporal quadrature point is given by
\begin{gather}
  \label{eqn:correction_update_flux_2}
   \vec{\mathcal F} \left( \tau \right) := \frac{1}{2} \left( \vec{f}\left( \vec{W_{\text{R}}}(\tau) \right) +
  \vec{f} \left( \vec{W_{\text{L}}}(\tau) \right) \right)  -
       \frac{1}{2} \bigl|\lambda(\tau)\bigr| \Bigl( \vec{q} \left( \vec{W_{\text{R}}}(\tau) \right)
         - \vec{q} \left( \vec{W_{\text{L}}}(\tau) \right) \Bigr),
\end{gather}
where
\begin{gather}
         \label{eqn:correction_update_flux_3}
  \vec{W_{\text{L}}}(\tau) := \vec{\Psi}\left(\tau,1\right)^T \mat{W^{n+\half}_{i-1}}, \qquad
  \vec{W_{\text{R}}}(\tau) := \vec{\Psi}\left(\tau,-1\right)^T \mat{W^{n+\half}_{i}},
\end{gather}
and $\bigl|\lambda(\tau)\bigr|$ is a local bound on the spectral radius of $\mat{A}\left(\vec{q} \right)$ in the neighborhood of interface $x=x_{i-\half}$ and at time $\tau$.


text:
\begin{document}

\begin{abstract}
Mixing rates and decay of correlations for dynamics defined by potentials with summable variations are well understood, but little is known for non-summable variations. In this paper, we exhibit upper bounds for these quantities in the case of dynamics defined by potentials with square summable variations. We obtain these bounds as corollaries of a new block coupling inequality between pair of dynamics starting with different histories. As applications of our results, we prove a new weak invariance principle and a Chernoff-type inequality.
\end{abstract}

\maketitle

\section{Introduction}
Let $A$ be a countable set, called {\it alphabet}. Consider a measurable function $\phi: A\times A^{\{-1,-2, \ldots \}}\to \R$ such that $\sum_{a\in A}e^{\phi(a,x)}=1$ for all $x\in A^{\{-1,-2, \ldots \}}$. The function $\phi$ is called a {\it potential}, and the probability kernel $g:= e^{\phi}$ (also known as {\it $g$-function}) is a natural generalization of Markov kernels.
Let $\eta=(\eta_n)_{\Z}$ be the canonical projections on $A^{\Z}$, \emph{i.e.}, for all $x\in A^{\Z}$ and all $n\in \Z$, $\eta_n(x)=x_n$. For $y \in A^{\{0,-1, \ldots \}}$, let $\mu^y$ be the probability measure on $A^{\Z}$ such that $\mu^y[B\times A^{\N}]=\delta_{y} [B]$
for measurable $B\times A^{\N}$ and for $n\geq 0$, $\mu^y[\eta_{n+1}=a\mid (\eta_n,\eta_{n-1},\dots)=x]=e^{\phi(a,x)}$
for every $a\in A$ and $\mu^y$-a.e.\ $x$ in $A^{\{-1,-2, \ldots \}}$. Let $T$ denotes the shift operator on $A^{\{0,-1, \ldots \}}$. We indicate by $d_{\text{TV}}$ the total variation distance, that is, if $P$ and $Q$ are two probability measures on the same $\sigma$-algebra $\mathcal{F}$,
\[
d_{\text{TV}}(P,Q)=\sup_{F\in \mathcal{F}}|P[F]-Q[F]|.
\]

In this paper, we obtain upper bounds, respectively, for the relaxation rate
\begin{equation*}
L(n) :=\sup_{y,z} d_{\text{TV}}(\mu^y[\eta_n \in \cdot\;],\mu^z[\eta_n \in \cdot\;]),
\end{equation*}
the mixing rate
\begin{equation*}
M(n) := \sup_{y,z} d_{\text{TV}}(\mu^y[(\eta_j)_{j\geq n} \in \cdot\;],\mu^z[(\eta_j)_{j\geq n} \in \cdot\;]),
\end{equation*}
and the decay rate of correlations
\begin{equation*}
\rho_{f,\hat{f}}(n) := \Big| \int f\circ T^n \;{\hat f} d\tilde{\mu}-\int f d\tilde{\mu}\int {\hat f} d\tilde{\mu} \Big|
\end{equation*}
when $\tilde{\mu}$ is the unique shift-invariant measure compatible with $\phi$ (see the next section for the definition of {\it compatibility}) and $f, \hat{f}$ are suitable functions (see Theorem \ref{Maintheo2}). \cite{bressaud/fernandez/galves/1999a} and  \cite{pollicott2000rates} obtained upper bounds for $L(n), M(n)$ and $\rho_{f,\hat{f}}(n)$  for potentials of summable variation and finite alphabet $A$. Our contribution is twofold. We obtain upper bounds for $L(n), M(n),$ and $\rho_{f,\hat{f}}(n)$  when the variation rate $ \textrm{var}_k(\phi)$ decays as $\mathcal{O}(k^{-(1/2+\delta')})$ for any $ \delta' > 0$. Moreover, our results also holds for countably infinite alphabet $A$. Theorem~\ref{Maintheo} is our main result, showing a new upper bound for the coupling error between $\mu^y$ and $\mu^z$. Corollary \ref{Cor1} answers a question posed in \cite{johansson2008square}, in which the authors ask for a bound for $L(n)$ when the variation of $\phi$ is not summable. Corollary \ref{Cor2} shows a bound for $M(n)$, which  cannot be achieved by simply using the union bound and Corollary \ref{Cor1}. The result is new even for the case of summable variation. The interest in $L(n)$ stands from the fact that it is the natural generalization of mixing times for Markov chains. \cite{GGT2018} showed that $M(n)$ converges to $0$ only when the potential is square summable $\tilde{\mu}$-a.s., hence Corollary \ref{Cor2} covers the main cases of interest. Theorem \ref{Maintheo2} gives an upper bound for the speed of decay of correlations, extending Theorem 1 in \cite{bressaud/fernandez/galves/1999a}. \cite{johansson/oberg/pollicott/2012} showed, when the alphabet $A$ is finite, that there is a unique shift-invariant measure $\tilde{\mu}$ compatible with $\phi$ when $\textrm{var}_k(\phi)  \in \smallO(k^{-(1/2)})$. Moreover,  \cite{berger2018non} showed that whenever $\textrm{var}_k(\phi)  \in \mathcal{O}(k^{-(1/2-\delta)})$ for any $\delta > 0$, there exists a regular $\phi$ that exhibits multiple compatible shift-invariant measures. Therefore, Theorem \ref{Maintheo2} also covers the main variation rates of interest under uniqueness of the compatible shift-invariant measure.  We illustrate the application of our inequalities in three cases. The first application proves a novel weak invariance principle for additive functionals of dynamics with non-summable variations. The second application shows that we can obtain Chernoff-type bounds for averages of random variables when the variation of $\phi$ is not summable. The third example illustrates how we can apply our results on a Poisson autoregression model, which is popular in applied works.

The proof technique is based on a renewal equation and coupling inequalities. These ideas were developed in \cite{coelho_quas_1998, bressaud/fernandez/galves/1999a, comets_fernandez_ferrari_2002}. We improve on the coupling bounds obtained in \cite{bressaud/fernandez/galves/1999a} by using a coupling between blocks of coordinates, instead of one coordinate at a time.  A block coupling idea was used in \cite{johansson/oberg/pollicott/2012} to obtain sharp conditions for uniqueness of the equilibrium measure for $\phi$ on a finite alphabet $A$, but no rate of mixing was obtained. A difference between \cite{johansson/oberg/pollicott/2012} and our approach is that we upper bound the block coupling using different renewal processes leading to a distinct renewal equation. This new renewal equation allows us to upper bound the speed of decay of coupling inequality even when the variation is not summable (see Theorem \ref{Maintheo}).

\section{Definitions}
Let the alphabet $A$ be a countable set, $\mathcal{X}=A^{\Z}$ and $\mathcal{X}_{-} = A^{\Z_-}$ where $\Z_-=\{0,-1,-2,\dots\}$. We endow $\mathcal{X}$ and $\mathcal{X}_-$ with the product topology and its corresponding Borel $\sigma$-algebra. The topologies and $\sigma$-algebras considered on subsets of $\mathcal{X}$ and $\mathcal{X}_-$ will always be the trace topologies and $\sigma$-algebras. We denote by $x_i$ the $i$-th coordinate of $x \in \mathcal{X}$ and for $-\infty< i \leq j<\infty$ we write $x^{-i}_{-j}:=(x_{-i},\ldots, x_{-j})$, $x^{-i}_{-\infty}:=(x_{-i},x_{-i-1},\dots)$ and $x_i^{\infty}:=(\dots,x_{i+1},x_{i})$. For $x\in \mathcal{X}$ and $y \in \mathcal{X}_{-}$, a \emph{concatenation} $x^{0}_{-i}y$ is a new sequence $z\in \mathcal{X}_{-}$ with $z^{0}_{-i} = x^{-1}_{-i}$ and $z^{-i-1}_{-\infty} = y$. Let $\o$ be the neutral element of the concatenation operation, that is $\o x=x$ for all $x\in \mathcal{X}_{-}$.  Note that we are using the convention, consistent with the concatenation operation, that when we scan an element $x\in\mathcal{X}$ from the left to the right we go further into the past.

Consider a measurable function $\phi: \mathcal{X}_{-}\to \R$.  We say that $\phi$ is {\it normalized} if it satisfies
\begin{equation*}
\sum_{a\in A}e^{\phi(ax)}=1
\end{equation*}
for all $x\in \mathcal{X}_{-}$. To a normalized function $\phi$, also called {\it potential}, we can associate a {\it probability kernel} (or simply a \emph{kernel}) $g$ on the alphabet $A$ by defining $g=e^{\phi}$. The \emph{variation rate} (or \emph{continuity rate}) of order $k\geq 1$ of $\phi$ is defined by
\begin{equation*}
\textrm{var}_{k}(\phi):=\sup_{z \in \mathcal{X}}\sup_{x,y\in \mathcal{X}_-}\sum_{b\in A}\big|\phi(bz_{-k}^{-1}x)-\phi(bz_{-k}^{-1}y)\big|.
\end{equation*}
We also define, for $k\geq 1$, the {\it $\chi^2$-variation rate} of order $k$ of $\phi$ as
\begin{equation*}
\chi^2_{k}(\phi)=\sup_{z\in \mathcal{X}}\sup_{x,y\in \mathcal{X}_-}\sum_{b\in A}\frac{\big(e^{\phi(bz_{-k}^{-1}x)}-e^{\phi(bz_{-k}^{-1}y)}\big)^2}{e^{\phi(bz_{-k}^{-1}y)}}.
\end{equation*}
A potential $\phi$ is said {\it regular} if $\lim_{k\to \infty}\text{var}_k(\phi)=0$ (continuity) and if it exists $\varepsilon>0$ such that $\phi(x)\geq \varepsilon$, for all $x\in \mathcal{X}_{-}$ (positivity). Observe this last property can be satisfied only when $A$ is finite.
We always have that $ \textrm{var}^2_{k}(e^\phi)\leq \chi^2_{k}(\phi)$ for any $k\geq 1$. When $\phi$ is regular $\textrm{var}^2_k(\phi)$ and $\chi^2_k(\phi)$ are comparable, that is, there exist positive constants $K_1$ and $K_2$ such that $K_1\textrm{var}^2_k(\phi)\leq \chi^2_k(\phi)\leq  K_2\textrm{var}^2_k(\phi)$. The $\chi^2$-variation introduced in this work will be particularly useful to study asymptotic properties of positive kernels on infinite $A$ (cf. Section \ref{Poisson}).

Let $\eta=(\eta_n)_{\Z}$ be the canonical projections on $\mathcal{X}$, that is for all $x\in \mathcal{X}$, $\eta_n(x)=x_n$ for all $n\in \Z$. We say that a probability measure $\mu$ on $\mathcal{X}$ is \emph{compatible} with the potential $\phi$ if there exits a probability measure $P$ on $\mathcal{X}_-$ such that
\[
\mu[\eta_{-\infty}^{0}\in B]=P[B]
\]
for all $B\subset \mathcal{X}_-$ measurable and if for $n\geq 1$
\begin{equation*}\label{compa}
\mu[\eta_{n+1}=a\mid \eta_{-\infty}^{n}=x]=e^{\phi(ax)}
\end{equation*}
for every $a\in A$ and $\mu$-a.e.\ $x$ in $\mathcal{X}_{-}$.

A shift-invariant measure $\tilde{\mu}$ compatible with $\phi$ (not necessarily normalized)  can be characterized via variational principle, in which context is called equilibrium state for $\phi$. The measure $\tilde{\mu}$ compatible with a normalized $\phi$ is also called $g$-measure in ergodic theory \cite{keane/1972}.  In probability literature, $g$-measures are known  as chains of complete connection \cite{doeblin/fortet/1937, iosifescu/grigorescu/1990}, chains of infinite order \cite{harris/1955}, \cite{keane/1972}, random-step Markov process \cite{kalikow}, and uniform martingale \cite{kalikow}. Compatible measures that are not necessarily shift-invariant are called $g$-chains \cite{johansson/oberg/pollicott/2012} or  stochastic chains of unbounded memory \cite{GGT2018}. When there is more than one shift-invariant measure compatible with $\phi$, we say that there is a \emph{phase transition}, otherwise we say that the shift-invariant compatible measure is \emph{unique}.

\section{Results}
\label{Results}

 In this paper, we will work under the following assumption.
 \begin{assump}[A]
 \label{Assump1}
 $\phi$ is a potential on $\mathcal{X}_-$ such that for all $k\geq 1$,
 \[
 \chi^2_{k}(\phi)\leq \frac{C}{k^{1+\delta}}
 \]
for some $C>0$ and $\delta>0$.
 \end{assump}

\begin{remark}
When $\phi$ is regular, Assumption (A) implies
 \[
 \emph{var}_{k}(\phi)\leq \frac{C'}{k^{\frac{1+ \delta}{2}}},
 \]
for some $C' > 0$ and $\delta > 0$. Hence, $ \emph{var}_k(\phi)$ is not summable when $\delta \in (0, 1]$.
\end{remark}

Now, consider $\mathcal{X}\times \mathcal{X}$ with the projection maps $\hat{\eta}=(\hat{\eta}_n)_{n\in \Z}$ and $\hat{\omega}=(\hat{\omega}_n)_{n\in \Z}$ such that for $(x,y)\in \mathcal{X}\times \mathcal{X}$, $\hat{\eta}_n(x,y)=x_n$ and $\hat{\omega}_n(x,y)=y_n$ for all $n\in \Z$. Let us also denote by $\mathcal{C}(\phi)$ the set of probability measures $P$ on $\mathcal{X}\times \mathcal{X}$ such that the pushforward measures $\hat{\eta}_*P$ and $\hat{\omega}_*P$ are compatible with $\phi$.
We also introduce the process $X=(X_n)_{n\geq 1}$, such that for all $n\geq 1$,
\begin{equation*}
X_n=\ind\{\exists j\in [M_n, M_{n+1}), \hat{\eta}_j \neq \hat{\omega}_j\},
\end{equation*}
where $(M_n)_{n\geq 1}$ is a fixed strictly increasing sequence of natural numbers such that $M_1=1$. We denote by $\mathcal{P}(\phi)$ the set of probability measures that are compatible with $\phi$.
Here is our main result followed by two corollaries.
\begin{theo}
\label{Maintheo}
Consider Assumption (A) and let $M_n=\lfloor n^{\beta} \rfloor$ for $\beta\geq 1$ and $\beta>1/\delta$. For all $\mu, \nu\in \mathcal{P}(\phi)$, there exists $\P\in \mathcal{C}(\phi)$ such that $\hat{\eta}_*\P=\mu$, $\hat{\omega}_*\P=\nu$ and for  $n\geq 1$,
\begin{equation*}
\P[X_n=1]\leq \frac{C_1}{n^{\frac{\beta \delta+1}{2}}}
\end{equation*}
where $C_1$ is a positive constant depending on $C, \delta$, and $\beta$.
\end{theo}

\begin{coro}
\label{Cor1}
If $\delta>1$, we have for all $\mu, \nu\in \mathcal{P}(\phi)$ and $n\geq 1$
\begin{equation*}
d_{\emph{TV}}((\eta_n)_*\mu,(\eta_n)_*\nu)\leq \frac{C_1}{n^{\frac{1+\delta}{2}}}
\end{equation*}
where $C_1$ is the same as in Theorem \ref{Maintheo}.\\
$ $

If $\delta\in (0,1]$, we have for all $\mu, \nu\in \mathcal{P}(\phi)$, $n\geq 1$ and $\delta'<\delta$,
\begin{equation*}
d_{\emph{TV}}((\eta_n)_*\mu,(\eta_n)_*\nu)\leq \frac{C_2}{n^{\delta'}},
\end{equation*}
where $C_2$ is a positive constant that depends on $C$, $\delta$, and $\delta'$.
\end{coro}

\begin{coro}
\label{Cor2}
Let $\mu, \nu\in \mathcal{P}(\phi)$. For all $\delta'<\delta$ there exists $\P\in \mathcal{C}(\phi)$ such that $\hat{\eta}_*\P=\mu$, $\hat{\omega}_*\P=\nu$ and
\begin{equation*}
d_{\emph{TV}}((\eta_n^\infty)_*\mu,(\eta_n^\infty)_*\nu) \leq \frac{C_3}{n^{\frac{\delta'}{2}}},
\end{equation*}
where $C_3$ is a positive constant that depends on $C$, $\delta$, and $\delta'$.\\
\end{coro}

\begin{remark}
When $\delta>1$ and $A$ is finite, Corollary \ref{Cor1} recovers the rate obtained in \cite[Theorem 1]{bressaud/fernandez/galves/1999a}.
\end{remark}

\begin{remark}
When $\delta>1$ and $A$ is finite, we can use \cite[Theorem 1]{bressaud/fernandez/galves/1999a} and the union bound to obtain
\begin{equation*}
d_{\emph{TV}}((\eta_n^\infty)_*\mu,(\eta_n^\infty)_*\nu) \leq \frac{C_4}{n^{\frac{\delta-1}{2}}},
\end{equation*}
where $C_4 > 0$. Hence, the result in Corollary \ref{Cor2} gives a sharper upperbound, even when the potential is summable and the alphabet $A$ is finite.
\end{remark}

We now look at the correlations decay for the shift-invariant measure compatible with a potential $\phi$. For this, we need the following definitions. Consider the shift operator $T: \mathcal{X}_{-}\to \mathcal{X}_{-}$ such that for all $x\in  \mathcal{X}_{-}$, $Tx=Tx_{-\infty}^0=x_{-\infty}^{-1}$. For non-constant $\phi$, let us consider the seminorm
\[
\|f\|_{\phi}=\sup_{k\geq 1}\frac{\textrm{var}_k(f)}{\textrm{var}_k(e^\phi)},
\]
and the subspace of $\mathcal{C}(\mathcal{X}_-,\R)$ defined by
$$V_{\phi}=\Big\{f \in \mathcal{C}(\mathcal{X}_-,\R):\|f\|_{\phi}<\infty   \Big\}.$$
\begin{theo}
\label{Maintheo2}
Consider Assumption (A) and assume that a shift-invariant probability measure $\tilde{\mu}$ compatible with $\phi$ exists. Let $f\in L^1(\tilde{\mu})$ and $\hat{f}\in V_{\phi}$.

If $\delta>1$, we have for all $n\geq 1$
\begin{equation*}
\Big| \int f\circ T^n\; {\hat f} \,d\tilde{\mu}-\int f \,d\tilde{\mu}\int {\hat f} \,d\tilde{\mu} \Big|\leq \frac{C_5}{n^{\frac{1+\delta}{2}}}\|f\|_1 \|\hat{f}\|_{\phi}
\end{equation*}
where $C_5$ is a positive constant that depends on $C$ and $\delta$.

If $\delta\in (0,1]$, we have for all $n\geq 1$ and $\delta'<\delta$,
\begin{equation*}
\Big| \int f\circ T^n\; {\hat f}\, d\tilde{\mu}-\int f\, d\tilde{\mu}\int {\hat f} \,d\tilde{\mu} \Big|\leq \frac{C_6}{n^{\delta'}}\|f\|_1 \|\hat{f}\|_{\phi},
\end{equation*}
where $C_6$ is a positive constant that depends on $C$, $\delta$, and $\delta'$.
\end{theo}

\begin{remark}
When $A$ is finite, continuity of $\phi$ guarantees the existence of a compatible shift-invariant measure. When $A$ is not finite, the existence of a shift-invariant compatible measure is not immediate. Sufficient conditions for existence of shift-invariant compatible measures when $A$ is infinite are given in \cite{fernandez/maillard/2005, johansson/oberg/pollicot/2007}.
Whenever a shift-invariant compatible measure exists, Assumption (A) guarantees the uniqueness of $\tilde{\mu}$ in Theorem \ref{Maintheo2} \cite{johansson/oberg/pollicot/2007}. See Section \ref{Poisson} for a concrete example.
\end{remark}

\begin{remark}
\cite{bressaud/fernandez/galves/1999a} and  \cite{pollicott2000rates} also studied the correlations decay for the general case in which $\phi$ is not normalized. The idea is to reduce the problem to the study of a normalized function $\psi$ that is cohomologous to $\phi$. When the rate of variation of $\phi$ is summable, such $\psi$ exists, although with a continuity rate that in the worst case could loose a factor of $1/n$. Hence, in general, the price to pay when $\phi$ is not normalized is to require higher regularity in the continuity rate of $\phi$ to guarantee the decay of correlations. We do not know of any corresponding result for non-normalized $\phi$ with square summable continuity rate.
\end{remark}

\section{Technical lemmas}
Here we collect some results that we will use to prove Theorem \ref{Maintheo}. We first recall the definitions of the Kullback-Leibler and Pearson $\chi^2$ divergences. Let $P$ and $Q$ be two probabilities on some discrete space $\mathcal{Y}$.
\[
D_{\text{KL}}(P||Q)=\sum_{y\in \mathcal{Y}}P(y)\ln\Big(\frac{P(y)}{Q(y)}\Big)
\]
and
\[
D_{\chi^2}(P||Q)=\sum_{y\in \mathcal{Y}}\frac{(P(y)-Q(y))^2}{Q(y)}.
\]
It is well known that $D_{\text{KL}}(P||Q)\leq D_{\chi^2}(P||Q)$ (cf. \cite[eq. 5]{sason2016f}).
\begin{lemma}
\label{chideu}
Let $x, y\in \mathcal{X}_-$ and $\mu,\nu \in \mathcal{P}(\phi)$ such that $\mu[\eta_{-\infty}^{0}\in\cdot\;]=\delta_x(\cdot)$ and\\ $\nu[\eta_{-\infty}^{0}\in\cdot\;]=\delta_y(\cdot)$.  For all $n\geq 1$, $0\leq k\leq n-1$ and all $a, b, c\in \mathcal{X}$, we have that
\begin{align}
\label{Div}
D_{\emph{KL}}\Big(\mu\Big[\eta_{M_{n}}^{M_{n+1}-1}\in\; \cdot\; \Big | \eta_1^{M_n-1}= a_{M_{n-k}}^{M_n-1} b_{1}^{M_{n-k}-1}\Big]\nonumber\\
\Big|\Big| \nu \Big[\eta_{M_{n}}^{M_{n+1}-1}\in \;\cdot \;\Big | &\eta_1^{M_n-1}=a_{M_{n-k}}^{M_n-1}c_{1}^{M_{n-k}-1}\Big]\Big)\nonumber\\
\leq \sum_{j=M_n}^{M_{n+1}-1}\chi^2_{j-M_{n-k}}(\phi).
\end{align}
\end{lemma}

\noindent
{\it Proof.} Let us simply denote by $D$ the left-hand term of inequality (\ref{Div}). We have by the chain rule property of the Kullback-Leibler divergence \cite[Theorem 2.5.3.]{cover/2006}
\begin{align*}
D
=&\sum_{i=M_n}^{M_{n+1}-1}D_{\text{KL}}\Big( \mu\Big[\eta_i\in \cdot \;\Big | \eta_1^{i-1}=z_{M_n}^{i-1}a_{M_{n-k}}^{M_n-1} b_{1}^{M_{n-k}-1}\Big]  \nonumber\\
&\phantom{******************}\Big|\Big| \nu\Big[\eta_i\in \cdot \;\Big | \eta_1^{i-1}=z_{M_n}^{i-1}a_{M_{n-k}}^{M_n-1}c_{1}^{M_{n-k}-1}\Big]    \Big)\nonumber\\
=&:\sum_{i=M_n}^{M_{n+1}-1}D_i.
\end{align*}
Then, we use the well known bound,
\begin{align*}
D_i&\leq D_{\chi^2}\Big(\mu\Big[\eta_i\in \cdot \;\Big | \eta_1^{i-1}=z_{M_n}^{i-1}a_{M_{n-k}}^{M_n-1} b_{1}^{M_{n-k}-1}\Big]\nonumber\\
&\phantom{***************} \Big|\Big| \nu\Big[\eta_i\in \cdot \;\Big | \eta_1^{i-1}=z_{M_n}^{i-1}a_{M_{n-k}}^{M_n-1}c_{1}^{M_{n-k}-1}\Big]\Big)\nonumber\\
&\leq \chi^2_{i-M_n}(\phi)
\end{align*}
to conclude the proof.
\qed

\begin{lemma}
\label{Lemalg}
For $\alpha> 1$ and $0<a<b$, we have that
\begin{equation}
\label{Alg}
\frac{(b+1)^{\alpha}-a^{\alpha}}{b^{\alpha}-a^{\alpha}}\geq \frac{(b+1)^{\alpha-1}-a^{\alpha-1}}{b^{\alpha-1}-a^{\alpha-1}}.
\end{equation}
\end{lemma}
\noindent
{\it Proof.} By algebraic computations, we obtain that $(\ref{Alg})$ is equivalent to
\[
\Big(\frac{b}{a}\Big)^{\alpha-1}\geq 1+(b-a)\Big(1-\Big(\frac{b}{b+1}\Big)^{\alpha-1}\Big).
\]
This last inequality is obtained from the Bernoulli inequality $(1+x)^{r}\geq 1+rx$, for $r>0$ and $x>-1$, observing that
\[
\Big(\frac{b}{a}\Big)^{\alpha-1}=\Big(1+\frac{b-a}{a}\Big)^{\alpha-1}\geq 1+(\alpha-1)\frac{b-a}{a}
\]
and
\[
\Big(\frac{b}{b+1}\Big)^{\alpha-1}=\Big(1-\frac{1}{b+1}\Big)^{\alpha-1}\geq 1-(\alpha-1)\frac{1}{b+1}.
\]
\qed

Define for all $\delta>0$, $\beta\geq 1$, $k\geq 3$ and $n\geq k+1$
\[
\Delta^n_k:=(n^{\beta}-(n-k)^{\beta}-2)^{-\delta}-((n+1)^{\beta}-(n-k)^{\beta})^{-\delta}.
\]

\begin{lemma}
\label{HJ1}
For all $\delta>0$, $\beta\geq 1$ and $k\geq 3$, $\Delta_k^n$
is a non increasing function of $n\geq k+1$.
\end{lemma}
\noindent
{\it Proof.}
The statement of the lemma is trivial for $\beta=1$.
For $\beta>1$, consider the function $f:[4,\infty)\to \R^+$ defined by
\[
f(x)=(x^{\beta}-(x-k)^{\beta}-2)^{-\delta}-((x+1)^{\beta}-(x-k)^{\beta})^{-\delta}.
\]
 In order to prove the result, it is enough to show that the derivative of $f$ is negative.
Since
\begin{align*}
f'(x)&=-\delta \beta\Big[\Big(x^{\beta}-(x-k)^{\beta}-2\Big)^{-\delta-1}\Big(x^{\beta-1}-(x-k)^{\beta-1}\Big)\nonumber\\
&\phantom{**}-\Big((x+1)^{\beta}-(x-k)^{\beta}\Big)^{-\delta-1}\Big((x+1)^{\beta-1}-(x-k)^{\beta-1}\Big)\Big]
\end{align*}
it is enough to show that
\[
\frac{(x+1)^{\beta}-(x-k)^{\beta}}{x^{\beta}-(x-k)^{\beta}-2}\geq \frac{(x+1)^{\beta-1}-(x-k)^{\beta-1}}{x^{\beta-1}-(x-k)^{\beta-1}}.
\]
But this last inequality follows from Lemma \ref{Lemalg}.
\qed

\begin{lemma}
\label{HJ2}
For all $\delta>0$, $\beta\geq 1$ and $k\geq 3$, we have that
\begin{equation*}
\Delta_k^{k+1}\leq 6\frac{2^{\beta} 4^{\delta}\beta}{k^{\delta\beta+1}}.
\end{equation*}
\end{lemma}
\noindent
{\it Proof.} Observe that
\begin{align*}
\Delta_k^{k+1}&=\int_{(k+1)^{\beta}-3}^{(k+2)^{\beta}-1}\frac{1}{x^{1+\delta}}dx\\
&\leq \frac{(k+2)^{\beta}-(k+1)^{\beta}+2}{((k+1)^{\beta}-3)^{1+\delta}}\\
&\leq \frac{\beta (k+2)^{\beta-1}+2}{((k+1)^{\beta}-3)^{1+\delta}}\\
&\leq 6\frac{2^{\beta} 4^{\delta}\beta}{k^{\delta\beta+1}}
\end{align*}
for $k\geq 3$.\qed

Finally, we recall the following lemma in \cite{bressaud/fernandez/galves/1999a} that gives an estimate for the renewal sequence that will appear in the proof of Theorem \ref{Maintheo}. We state the lemma using a notation that is adapted to our purpose.
\begin{lemma}[ Proposition 2 item (iv) in \cite{bressaud/fernandez/galves/1999a}]
\label{renewal}
Let $(f_k)_{k \geq 1}$ be a sequence of positive real numbers such that $\sum_{k=1}^\infty f_k < 1$. Supose that $(u_k)_{k \geq 1}$ is a sequence with $u_0 = 1$ and satisfies the renewal equation
\begin{equation*}
u_n = \sum_{k=1}^n f_k u_{n-k}.
\end{equation*}
If $f_n \leq c_1/n^{1+\alpha}$, for some $\alpha >0$, then $u_n \leq c_2/n^{1+\alpha}$, where $c_1$ and $c_2$ are constants that depend on $(f_k)_{k \geq 1}$.
\end{lemma}

\section{Proof of Theorem \ref{Maintheo}}
\label{Lil}
Observe that it is enough to prove Theorem \ref{Maintheo} when $\mu, \nu$ are such that there exist $x, y\in \mathcal{X}_-$ and $\mu[\eta_{-\infty}^{0}\in\cdot\;]=\delta_x(\cdot)$, $\nu[\eta_{-\infty}^{0}\in\cdot\;]=\delta_y(\cdot)$.

We now construct the coupling of $\mu$ and $\nu$, that we call $\P^{x,y}$, as follows. We start by defining
$$\P^{x,y}[\hat{\eta}_{-\infty}^{0}\in\cdot\;, \hat{\omega}_{-\infty}^{0}\in\cdot\;]=\delta_x \otimes\delta_y.$$ Then, for all $n\geq 1$,  given the pasts $\hat{\eta}_{-\infty}^{M_n-1}$ and  $\hat{\omega}_{-\infty}^{M_n-1}$, we maximally couple $\hat{\eta}_{M_n}^{M_{n+1}-1}$ and $\hat{\omega}_{M_n}^{M_{n+1}-1}$ to complete the construction of $\P^{x,y}$.

Next, we show that $\P^{x,y}$ satisfies the inequality in Theorem \ref{Maintheo}. For all $n\geq  1$ and $0\leq k\leq n-1$, define
\begin{equation*}
q^n_k=\sup_{x,y,a,b \in\mathcal{X}}\P^{x,y}[X_{n}=1\mid X_{n-k}^{n-1}=0, \hat{\eta}^{M_{n-k}-1}_{1} = a_{1}^{M_{n-k}-1}, \hat{\omega}^{M_{n-k}-1}_{1} = b_{1}^{M_{n-k}-1}],
\end{equation*}
with the convention that if $k>l$ then elements of the form $a_k^l$ are dropped from the conditional. The shorthand notation $X_{n-k}^{n-1}=0$ means that $X_{n-k} = 0, \ldots, X_{n-1} = 0$. Observe that for all $n\geq  1$ and $0\leq k\leq n-2$ we have $q^n_{k} \geq q^n_{k+1}$.

We start by proving the following
\begin{lemma}
\label{HJ0}
Suppose that $(\chi^2_n(\phi))_{n\geq 1}\in \ell^1$. Then, there exists $\eps>0$ such that for all $k\geq 0$ and all $n\geq k+1$,
\begin{equation}
\label{Bret}
q_k^n\leq\sqrt{1-\exp\left(-\sum_{k=1}^{\infty}\chi_k^2(\phi)\right)}\leq 1-\eps.
\end{equation}
For $k\geq 1$ and $n\geq k+1$, we also have
\begin{equation}
\label{Pins}
q_k^n\leq \sqrt{\frac{1}{2}\sum_{j=M_n}^{M_{n+1}-1}\chi^2_{j-M_{n-k}}(\phi)}.
\end{equation}
\end{lemma}
\noindent
{\it Proof.}
Inequality (\ref{Bret}) is a direct consequence of the Bretagnolle-Huber inequality (cf. \cite[eq. 4]{sason2016f}) and Lemma \ref{chideu}. Inequality (\ref{Pins}) is a direct consequence of the Pinsker inequality (cf. \cite[eq. 1]{sason2016f}) and again Lemma \ref{chideu}.
\qed

Now, on some probability space $(\Omega, \mathcal{F}, P)$, consider the random process $Y=(Y_n)_{n\geq 0}$ with values in $\{0,1\}$ such that $Y_0=1$ and  for $n\geq 1$ and $0\leq k\leq n-1$,
\begin{equation*}
P[Y_n=1\mid Y^{n-1}_{n-k}=0, Y_{n-k-1}=1, Y_1^{n-k-2}]=q^n_k.
\end{equation*}
For all $m\geq 1$ and $a,b \in \{0,1\}^{m}$ we say that $a\geq b$ if $a_i \geq b_i$ for $i \in \{1,\ldots, m\}$. By construction, for all $n\geq 2$, $a,b \in \{0,1\}^{n-1}$, and $a\geq b$, we have $P[Y_n = 1 | Y^{n-1}_1 = a] \geq \P^{x,y}[X_n = 1 | X^{n-1}_1 = b]$. Therefore, by applying Strassen's Theorem on stochastic domination \cite{lindvall1999strassen} inductively on $n$, we can construct a coupling measure $Q$ such that, for  $a,b \in \{0,1\}^{n-1}$, and $a\geq b$, we have $Q[Y_n \geq X_n| Y^{n-1}_1 = a, X^{n-1}_1 = b] = 1$. Therefore, for all $n\geq 1$, we have $Q[Y_n \geq X_n] = 1$, which implies that
\begin{equation}\label{dominance}
P[Y_n = 1] \geq \P^{x,y}[X_n = 1],
\end{equation}
 for all $n \geq 1$.

Now, consider the process $Z=(Z_n)_{n\geq 0}$ with values in $\{0,1\}$ such that $Z_0=1$ and  for $n\geq 1$ and $0\leq k\leq n-1$,
\begin{equation}
\label{GHJ1}
P[Z_n=1\mid Z^{n-1}_{n-k}=0, Z_{n-k-1}=1, Z_1^{n-k-2}]=b_k:=\sup_{n\geq k+1} q_k^n.
\end{equation}
Observe that, for all $k \geq 0$, we have $b_k \geq b_{k+1}$. Using the same argument used to show \eqref{dominance}, we have that $P[Z_n = 1] \geq P[Y_n = 1]$, for all $n\geq 1$.
Also, by (\ref{Pins}), and Lemmas \ref{HJ1} and \ref{HJ2}, we have  that for $k\geq 3$ and $n\geq k+1$,
\begin{align*}
2(q_k^n)^2&\leq \sum_{j=M_n}^{M_{n+1}-1}\chi^2_{j-M_{n-k}}(\phi)\\
&\leq C\sum_{j=\lfloor n^{\beta}\rfloor}^{\lfloor(n+1)^{\beta}\rfloor-1}\frac{1}{(j-\lfloor (n-k)^{\beta}\rfloor)^{1+\delta}}\\
&=C\sum_{j=\lfloor n^{\beta}\rfloor-\lfloor (n-k)^{\beta}\rfloor}^{\lfloor(n+1)^{\beta}\rfloor-\lfloor (n-k)^{\beta}\rfloor-1}\frac{1}{j^{1+\delta}}\\
&\leq C\int_{n^{\beta}-(n-k)^{\beta}-2}^{(n+1)^{\beta}-(n-k)^{\beta}}\frac{1}{x^{1+\delta}}dx\\
&=C\frac{\Delta_k^n}{\delta}\\
&\leq C\delta^{-1}\Delta_k^{k+1}\\
&\leq 6C\frac{2^{\beta} 4^{\delta}\beta\delta^{-1}}{k^{\delta\beta+1}}.
\end{align*}
Using (\ref{Bret}), we obtain that $b_k\leq (3C\frac{2^{\beta} 4^{\delta}\beta\delta^{-1}}{k^{\delta\beta+1}})^{1/2}\wedge (1-\eps)$ for all $k\geq 3$ and
$b_k\leq 1-\eps$ for $0\leq k\leq 2$.

Next, let $f_i:= b_{i-1}\prod_{k=0}^{i-2}(1-b_k)$
for $i\geq 1$ (with the convention that $\prod_{j=0}^{-1}=1$)
and
$u_i:= P[Z_i = 1]$
for $i \geq 0$.
We have that
$$P[Z_n = 1] = \sum_{k= 1}^nP[Z_n = 1, Z^{n-1}_{n-k+1} = 0| Z_{n-k} = 1]P[Z_{n-k} = 1],$$
hence the following renewal equation holds
$$u_n = \sum_{k = 1}^nf_{k}u_{n-k}.$$
By definition, we have that $\sum_{k= 1}^\infty f_k = 1-\prod_{k = 1}^\infty(1-b_k)$. If $\beta>\delta^{-1}$, we have $\sum_{k = 0}^\infty b_k < \infty$. Hence $\sum_{k= 1}^\infty f_k < 1$. Moreover, when $\beta>\delta^{-1}$ we have that $b_k \leq c_1k^{-(\delta\beta+1)/2}$, for some positive constant $c_1$ that depends on $C, \delta$ and $\beta$. From Lemma \ref{renewal}, we have that, for all $n\geq 1$,
\begin{align*}
u_n\leq \frac{C_1}{n^{\frac{\delta\beta+1}{2}}}
\end{align*}
where $C_1$ is a positive constant that depends on $C, \delta$ and $\beta$.
Because $P[Z_n = 1] \geq \P^{x,y}[X_n = 1]$,  we obtain that for all $n\geq 1$,
\begin{equation*}
\P^{x,y}[X_n=1]\leq \frac{C_1}{n^{\frac{\delta\beta+1}{2}}}
\end{equation*}
for $\beta>\delta^{-1}$.
\qed

\section{Proofs of Corollaries \ref{Cor1} and \ref{Cor2}}
\subsection{Proof of Corollary  \ref{Cor1}}

We first observe that by definition of the total variation distance, for $k\in [M_n, M_{n+1})$,
\begin{equation*}
d_{\text{TV}}((\eta_k)_*\mu,(\eta_k)_*\nu)\leq \P[\hat{\eta}_k\neq \hat{\omega}_k]\leq \P[X_n=1].
\end{equation*}
Then, by Theorem \ref{Maintheo}, we obtain that
\begin{equation*}
\P[\hat{\eta}_k\neq \hat{\omega}_k]\leq \frac{C_1}{n^{\frac{\beta \delta+1}{2}}}.
\end{equation*}
for all $\beta\geq 1$ and $\beta>\delta^{-1}$. If $\delta>1$, just take $\beta=1$. In this case $k=n$, thus we obtain Corollary \ref{Cor1}. If $\delta\in (0,1)$, since $k\leq (n+1)^{\beta}$, we have $n\geq k^{1/\beta}-1$. This leads to
\begin{equation*}
\P[\hat{\eta}_k\neq \hat{\omega}_k]\leq \frac{C_2}{k^{\frac{\beta \delta+1}{2\beta}}}
\end{equation*}
for all $k\geq 1$, where $C_2$ is a positive constant that depends on $C, \delta$, and $\beta$. Now, observe that  for any $0< \delta'<\delta$, we can choose $\beta$ such that $\beta\geq 1$, $\beta>\delta^{-1}$ and $\frac{\beta \delta+1}{2\beta}\geq \delta'$.
\qed
\subsection{Proof of Corollary  \ref{Cor2}}

Consider $k\in [M_n, M_{n+1})$.  Let
$$\theta=\inf\{n\geq 1: \hat{\eta}_k=\hat{\omega}_k,\; \text{for all}\; k\geq n\},$$ with the convention that $\inf\emptyset=\infty$. We start by observing that
\begin{align*}
\P[\theta>k]\leq \P\Big[\bigcup_{j\geq n} \{ X_j=1\}\Big]\leq \sum_{j\geq n} \P[X_j=1].
\end{align*}
By Theorem \ref{Maintheo}, we obtain that
\begin{align*}
\P[\theta>k]\leq C_1\sum_{j\geq n} \frac{1}{n^{\frac{\beta \delta+1}{2}}}\leq \frac{C_1'}{n^{\frac{\beta\delta-1}{2}}}
\end{align*}
for $\beta\geq 1$, $\beta>\delta^{-1}$ and $C_1'$ a positive constant that depends on $C, \delta$, and $\beta$. Since $n\geq k^{1/\beta}-1$, we obtain that
\begin{align*}
\P[\theta>k]\leq \frac{C_3}{k^{\frac{\beta\delta-1}{2\beta}}}
\end{align*}
for all $k\geq 1$ and $C_3$ a positive constant that depends on $C, \delta$, and $\beta$. Finally, notice that for all $\delta'<\delta$, we can choose $\beta$ large enough such that $\frac{\beta\delta-1}{2\beta}\geq \delta'$. We finish the proof by using the coupling inequality (cf. \cite{thorisson/2000})
\begin{align*}
d_{\text{TV}}((\eta_n^\infty)_*\mu,(\eta_n^\infty)_*\nu) \leq \P[\theta>k].
\end{align*}
\qed

\section{Proof of Theorem \ref{Maintheo2}}
Consider $M_m=\lfloor m^{\beta} \rfloor$, for $m\ge 1$. For each $x, y\in \mathcal{X}_-$, we consider a probability space $(\Omega, \mathcal{F}, \P^{x,y})$ that supports the random elements $\tilde{\eta}$, $\tilde{\omega}$ and $\tilde{Z}$ defined as follows: $\tilde{\eta}_*\P^{x,y}$, $\tilde{\omega}_*\P^{x,y} \in \mathcal{P}(\phi)$ and $\tilde{\eta}_{-\infty}^{0}=x,  \tilde{\omega}_{-\infty}^{0}=y$. Also, for all $m\geq 1$ given the pasts $\tilde{\eta}_{-\infty}^{M_m-1}$ and  $\tilde{\omega}_{-\infty}^{M_m-1}$, $\tilde{\eta}_{M_m}^{M_{m+1}-1}$ and $\tilde{\omega}_{M_m}^{M_{m+1}-1}$ are maximally coupled. Under $\P^{x,y}$, the process $\tilde{Z}$ has the same law as the process $Z$ defined in Section \ref{Lil} and verifies $\tilde{Z}_m\geq \tilde{X}_m:=\ind\{\exists j\in [M_m, M_{m+1}), \tilde{\eta}_j \neq \tilde{\omega}_j\}$ for all $m\geq 1$ (this is indeed possible since $Z$ stochastically dominates $X$, see Section \ref{Lil}). We denote by $\E^{x,y}$ the expectation with respect to $\P^{x,y}$.

Fix some $n\in \N$ and let $k$ be such that $n\in [M_{k-1}, M_{k})$. We will show that
\begin{equation}
\label{RIC}
\Big| \int f\circ T^n \;{\hat f} d\tilde\mu-\int f d\tilde\mu\int {\hat f} d\tilde\mu \Big|\leq c_1 P[Z_k=1]
\end{equation}
for some positive constant $c_1$ that depends only on $C$, $\delta$ and $\beta$.
From this point, Theorem \ref{Maintheo2} is easily obtained following the proof of Corollary \ref{Cor1}.
To obtain (\ref{RIC}), we follow the argument developed in \cite{bressaud/fernandez/galves/1999a}, Section 5.
Using (3.7) in  \cite{bressaud/fernandez/galves/1999a}, we first observe that
\begin{align*}
\Big| \int f\circ T^n \;{\hat f} d\tilde\mu-\int f d\tilde\mu\int {\hat f} d\tilde\mu\Big | \leq \|f\|_1\sup_{x,y} \E^{x,y}\Big[\big|\hat{f}(\tilde{\eta}_{-\infty}^n)-\hat{f}(\tilde{\omega}_{-\infty}^{n})  \big|\Big].
\end{align*}
For $k\geq 1$, let
\[
\theta_k=\inf\{0\leq m\leq k:\tilde{Z}_{k-m}=1\}.
\]
We have that
\begin{align*}
\E^{x,y}\Big[\big|\hat{f}(\tilde{\eta}_{-\infty}^n)-\hat{f}(\tilde{\omega}_{-\infty}^{n})  \big|\Big]
&=\E^{x,y}\Big[\sum_{j=0}^k\ind\{\theta_k=j\}\big|\hat{f}(\tilde{\eta}_{-\infty}^n)-\hat{f}(\tilde{\omega}_{-\infty}^{n})  \big|\Big]\nonumber\\
&\leq \|\hat{f}\|_{\phi} \sum_{j=0}^k \textrm{var}_{n-M_{k-j}+1}(e^{\phi})\P^{x,y}[\theta_k=j].
\end{align*}
Observe that for all $0\leq j\leq k$,
\begin{align*}
\P^{x,y}[\theta_k=j]&=P[Z_k=0,\dots,Z_{k-j+1}=0,Z_{k-j}=1]\nonumber\\
&=\prod_{l=1}^j(1-b_{k-j+l})P[Z_{k-j}=1]
\end{align*}
where $b_{k-j+l}$ are from (\ref{GHJ1}).
Now, observe that for all $i\geq 1$ we have
\begin{align*}
P[Z_i=1]&=\sum_{k=1}^ib_{i-k}P[Z_{i-1}=0\mid Z_k^{i-1}=0, Z_{k-1}=1]\nonumber\\
&=\sum_{k=1}^i b_{i-k}\prod_{l=1}^{i-k}(1-b_{i-k-l})P[Z_{k-1}=1].
\end{align*}
Thus, we have that for all $i\geq 1$
\begin{align*}
P[Z_i=1]&=\sum_{k=1}^if_kP[Z_{i-k}=1]
\end{align*}
with $f_k:=b_{k-1}\prod_{l=0}^{k-2}(1-b_l)$, $k\geq 1$.
From this, we obtain that
\begin{align*}
\E^{x,y}\Big[\big|\hat{f}(\tilde{\eta}_{-\infty}^n)-\hat{f}(\tilde{\omega}_{-\infty}^{n})  \big|\Big]&\leq
\|\hat{f}\|_{\phi}\Bigg( \textrm{var}_{n-M_{k}+1}(e^\phi)\sum_{l=1}^k f_lP[Z_{k-l}=1]\nonumber\\
&\;\;\;\;\;+ \sum_{l=1}^k \textrm{var}_{n-M_{k-l}+1}(e^\phi)\prod_{m=1}^{l}(1-b_{k-l+m})P[Z_{k-l}=1]   \Bigg).
\end{align*}
We deduce that
\begin{align*}
\sup_{x,y}\E^{x,y}\Big[\big|\hat{f}(\tilde{\eta}_{-\infty}^n)-\hat{f}(\tilde{\omega}_{-\infty}^{n})  \big|\Big]&\leq
 \kappa \sum_{l=1}^k f_lP[Z_{k-l}=1]=\kappa P[Z_{k}=1]
\end{align*}
with
$$\kappa:=\textrm{var}_{n-M_{k}+1}(e^\phi)+\sup_{1\leq l\leq k}\frac{\textrm{var}_{n-M_{k-l}+1}(e^\phi)}{f_l}.
$$

Finally, since $\frac{\textrm{var}_{n-M_{k-l}+1}(e^\phi)}{b_{l-1}}\leq 2$ and $\prod_{j=0}^{\infty}(1-b_j)>0$ (using that $b_0 < 1$ and $\sum_{j=0}^\infty b_j<\infty$), we observe that
\begin{equation*}
\kappa\leq 1+\frac{\textrm{var}_{n-M_{k-l}+1}(e^\phi)}{b_{l-1}\prod_{j=0}^{\infty}(1-b_j)}\leq c_2
\end{equation*}
for some positive constant $c_2$ depending on $C$, $\delta$, and $\beta$.
\qed

\section{Applications}
\subsection{Functional central limit theorem (FCLT) for potentials with non-summable variations}
Let $\sigma > 0$. A function $h:A\to \R$ satisfies the FCLT, also called Weak Invariance Principle, when the process $\{\zeta_n(t), t \in [0,1], n\geq 1\}$ defined by
\begin{equation*}
\zeta_n(t) = \frac{1}{\sigma \sqrt{n}} \sum_{i = 0}^{\lfloor nt \rfloor} h(\eta_i),
\end{equation*}
converges weakly to a standard Brownian motion on $D[0,1]$. \cite[Section 4.3]{tyran2005invariance} showed that the FCLT holds when the potential $\phi$ has summable variations. A straightforward application of Theorem 1 in \cite{tyran2005invariance} and our Corollary~ \ref{Cor1} is the following FCLT for potentials with non-summable variations.
\begin{prop}
\label{fclt}
Assume that $\phi$ is regular and satisfies Assumption (A) with $\delta \in (1/2,1]$. Let $\mu$ be shift-invariant and compatible with $\phi$. Also, let $h:A \to \R$ be a function such that $\int h\circ \eta_0\; d\mu = 0$. If $\sigma^2 :=\int (h\circ \eta_0)^2\; d\mu > 0$, then $h$ satisfies the FCLT.
\end{prop}
\begin{proof}
Because $\phi$ is regular, the alphabet is finite. Hence $\sigma^2 \leq \|h\|^2_\infty < \infty$ as required by Theorem 1 in \cite{tyran2005invariance}. It remains to verify the condition on the mixing rate. Using the processes $\{\tilde{\eta}_i, i \in \Z\}$ and $\{\tilde{\omega}_i, i \in \Z\}$ introduced in the proof of Theorem \ref{Maintheo2}, it is sufficient to check that there exists $\gamma > 1/2$ such that
\begin{equation}\label{condFCLT}
\limsup_{n \rightarrow \infty} n^\gamma \sup_{x,y}\E^{x,y}[|h(\tilde{\eta}_n)-h(\tilde{\omega}_n)|] < \infty.
\end{equation}
We have from Corollary ~ \ref{Cor1} that, for all $x,y \in \X_-$,
\begin{align*}
\E^{x,y}[|h(\tilde{\eta}_n)-h(\tilde{\omega}_n)|] &\leq 2\|h\|_{\infty} \P^{x,y}[\tilde{\eta}_n\neq \tilde{\omega}_n]\\
&\leq \frac{c_1}{n^{\delta'}},
\end{align*}
where $\delta' < \delta$ and $c_1 $ is positive constant that depends on $C, \delta, \delta'$, and $h$. Taking $\delta'$ and $\gamma$ such that $1/2 < \gamma \leq \delta' <\delta$, we obtain \eqref{condFCLT}.
\end{proof}

\subsection{Chernoff-type inequality for potentials with non-summable variations}

Chernoff-type inequality gives finite sample bounds for deviations of additive functionals from their mean. When the variation rate of the potential is summable, we have an exponential inequality \cite{marton/1998, gallesco/gallo/takahashi/2014, chazottes2020optimal}. Nevertheless, the rate of concentration when the variation rate is not summable was an open question. Using our result, we can obtain the following stretched exponential inequality for sums of random variables.
\begin{prop}
\label{ConcIneq}
Assume that $\phi$ is regular and satisfies Assumption (A) with $\delta \in (1/2,1]$. Let $\mu$ be compatible with $\phi$. For all $\delta'<2\delta-1$, $n\geq 1$, $t\geq 0$ and all functions $h:A \to \R$ we have that
\begin{equation*}
\mu\Bigg[\Big|\frac{1}{n}\sum_{i=1}^n(h(\eta_i)-\E[h(\eta_i)])\Big|\geq t\Bigg]\leq 2\exp\Bigg\{-\frac{C_7n^{\delta'}t^2}{R(h)^2}\Bigg\}
\end{equation*}
where $R(h):=\max_{a\in A}h(a)-\min_{a\in A}h(a)$ and $C_7$  is a constant that depends on $c_2$, $\delta$ and~$\delta'$.
\end{prop}
\noindent
\begin{proof} This is a consequence of Theorem 1 of \cite{chazottes/collet/kulske/regig/2007} and Corollary~\ref{Cor1}. In order to apply Theorem 1 of \cite{chazottes/collet/kulske/regig/2007}, we need to estimate the terms $\|\overline{D}\|_{\ell^2(\N)}^2$ and $\|\delta f\|^2_{\ell^2(\N)}$ (for $f(x_1,\dots,x_n)=\frac{1}{n}\sum_{i=1}^n h(x_i)$) there.
We have that
\begin{equation*}
\|\overline{D}\|_{\ell^2(\N)}^2\leq \sup_{x,y\in \mathcal{X}_-} \Big(1+\sum_{i=1}^n\P^{x,y}[\hat{\eta}_i\neq \hat{\omega}_i]\Big)^2.
\end{equation*}
Now, for $\delta'<2\delta-1$, using Corollary~ \ref{Cor1}, we obtain that
\begin{equation}
\label{normD}
\|\overline{D}\|_{\ell^2(\N)}^2\leq c_1 n^{1-\delta'}
\end{equation}
for some positive constant $c_1$ depending on $C, \delta$ and $\delta'$.

For a given function $f: A^{n}\to \R$ we define the oscillation of $f$ at site $i\in \{1,\dots,n\}$, by
\[
\delta_if:= \sup_{x_j=x'_j, j \neq i}\big|f(x_1,\dots,x_n)-f(x'_1,\dots,x'_n)\big|.
\]
Now, taking $f(x_1,\dots, x_{n})= \frac{1}{n}\sum_{i=1}^{n}h(x_i)$ we have $\delta_i f=\frac{R(h)}{n}$ for $i\in \{1,\dots,n\}$. Thus, we obtain that
\begin{align}
\label{normf}
\|\delta f\|^2_{\ell^2(\N)}=\sum_{i=1}^{n}\left(\frac{R(h)}{n}\right)^2= \frac{R(h)^2}{n}.
\end{align}
Finally, using (\ref{normD}) and (\ref{normf}) in Theorem 1 of~\cite{chazottes/collet/kulske/regig/2007}, we obtain Proposition \ref{ConcIneq}.
\end{proof}

\subsection{Poisson autoregression model}
\label{Poisson}
As a second application of our results, we consider a model with countable infinite alphabet called Poisson autoregression, which is popular in applications \cite{kedem2005regression}. Only the Markovian case of these models were studied in the literature. We will show how we can choose the parameters of non-Markovian Possion autoregression models to satisfy the Assumption (A) and thus apply the results of Section~\ref{Results}.

Consider an absolutely converging sequence $(\beta_i)_{i\geq 1}$ and a sequence of non negative integers $(\gamma_i)_{i\geq 1}$ such that $S:=\sum_{i=1}^{\infty}|\beta_i|\gamma_i<\infty$.
Consider $A=\Z_+$ and the potential $\phi$ defined for all $x\in \mathcal{X}_-$ by
$$
\phi(x)=-\lambda(x_{-\infty}^{-1}) + x_0\log \lambda(x_{-\infty}^{-1}) - \sum_{k = 0}^{x_0}\log(k),
$$
where
$$
\lambda(x_{-\infty}^{-1})=\exp\Big\{\sum_{i=1}^{\infty}\beta_i (x_{-i}\wedge \gamma_i)\Big\}.
$$
For this model, we obtain that
\begin{equation}
\label{RTYUP}
\chi^2_k(\phi)=\sup_{a\in \mathcal{X}} \sup_{x,y\in  \mathcal{X}_-}\Bigg(e^{\lambda(a_{-k}^{-1}y)\Big(\frac{\lambda(a_{-k}^{-1}x)}{\lambda(a_{-k}^{-1}y)}-1\Big)^2}-1\Bigg).
\end{equation}
Now, since $e^{-S}  \leq \lambda(x_{-\infty}^{-1})\leq  e^S$ and  the exponential function is locally bilipschitz, using (\ref{RTYUP}), we have that
\begin{equation*}
c_1^{-1}  \Big(\sum_{i=k+1}^{\infty}|\beta_i|\gamma_i\Big)^2 \leq  \chi^2_k(\phi)\leq c_1 \Big(\sum_{i=k+1}^{\infty}|\beta_i|\gamma_i\Big)^2
\end{equation*}
where $c_1$ is a positive constant that depends only on $S$.

Finally, choosing the sequences $(\beta_i)_{i\geq 1}$ and $(\gamma_i)_{i\geq 1}$ such that
$$\frac{c_2^{-1}}{i^{\frac{3+\varepsilon}{2}}}\leq |\beta_i|\gamma_i \leq \frac{c_2}{i^{\frac{3+\varepsilon}{2}}}$$
 for some $\varepsilon>0$ and $c_2\geq 1$, we obtain that
\begin{equation*}
\frac{c_3^{-1}}{k^{1+\varepsilon}}\leq \chi^2_k(\phi)\leq \frac{c_3}{k^{1+\varepsilon}}.
\end{equation*}
where $c_3$ is a positive constant.

Finally, for this model, we mention that the existence of a shift-invariant probability measure compatible with $\phi$ is obtained by applying Theorem 5.1 of~\cite{johansson/oberg/pollicot/2007} with $K=e^{2\sinh S}$ and $\pi$ equal to the Poisson law with parameter $e^S$. Assumption (A) implies the square summability of the variation, which guarantees the uniqueness of the shift-invariant probability measure \cite[Corollary 4.2.]{johansson/oberg/pollicot/2007}.

\section*{Acknowledgements}
\noindent
C.~Gallesco and D.Y.~Takahashi would like to thank Sandro Gallo for several fruitful discussions that motivated this work. We also thank Leandro Cioletti for comments on a early version of this manuscript.
C.~Gallesco was partially supported by FAPESP (grant 2017/19876-4) and CNPq (grant 312181/2017-5). D.Y.~Takahashi thanks the support of FAPESP  Research, Innovation and Dissemination Center for Neuromathematics (grant 2013/ 07699-0).

\bibliographystyle{alpha}
\bibliography{MixingTime}

\end{document}


TITLE: Find all possible rational pairs of a parametric sextic and all cases where it is reducible for either parameter
QUESTION [0 upvotes]: Find all possible rational solutions pairs $(z,a)$ of the equation
$a^6 + a^4 (-18368 + 9184 z - 2912 z^2) + a^2 (61702144 - 61702144 z + 36814848 z^2 - 10694656 z^3 + 1748992 z^4) - z^2 (44281626624 - 44281626624 z + 13332971520 z^2 - 1131282432 z^3 + 28901376 z^4) = 0$.
I have the following solutions for $a = (0, \pm 24,\pm48,\pm72,\pm 96,\pm 144,\pm 192)$ and $z = (-2,0,1,4)$ which I think is the complete set.  A second question is to find all possible rational $a$ or $z$ such that this equation is reducible in either parameter.  I have the set $z =(-44,-2,0,1,44/23,67/19,4,134/29)$ which I also think is complete.  In this case for the $z$ values I get four pairs of cubics and one quartic all in parameter $a$.

REPLY [2 votes]: It seems hopeless to try to provably find all the rational points on the genus $7$ curve $C$ mentioned in the question using either Chabauty's method, or etale descent. I will simply explain how to find a map from the genus $7$ curve to the rank three elliptic curve mentioned in the comment.
First, I noticed that all of the exponents on $a$ are even, and so setting $b = a^{2}$, we get an equation $b^{3} + b^2 (-18368 + 9184z - 2912z^{2}) + \cdots$. This equation defined a curve of $X$ genus $4$. I computed the automorphism group of this curve modulo several primes, and noticed that the order was $2$ (for $p > 17$ or so). I then had Magma compute the automorphism group $G$ over $\mathbb{Q}$ and found it had order $2$. If we realize $X \subseteq \mathbb{P}^{3}$ in its canonical model, the automorphism can be given by a simple linear change of variables. This allows one to easily compute the curve quotient, which results in a plane cubic. Magma's built-in routines can then get one to the minimal Weierstrass model of $X/G$.
Searching for nice equations for the maps at each step results in the map from the original curve (in terms of $a$ and $z$), to the homogenous form of the elliptic curve ($y^{2} z = x^{3} - 159600xz^{2} + 31482500z^{3}$) given by the equations
$$
  x = a^{2} (470957175a^2 - 403726629600z^2 + 1240537082400z - 2481074164800)
$$
$$
  y = 6385126825a^4 - 6377874641600a^2z^2 + 33875545316800a^2z -
    67751090633600a^2 + 216129127244800z^4 - 8459911552153600z^3 +
    99706100436096000z^2 - 331145109327155200z + 331145109327155200
$$
$$
  z = a^{2} (1220184a^2 - 875686464z^2).
$$
If I hadn't tried several tricks to manually simplify the equations, they would have been several pages long.


TITLE: Second homotopy group of the wedge sum of $S^2$ with the presentation complex of a finitely generated group
QUESTION [0 upvotes]: I am reading a paper which makes the following claim:
let $G$ be a finitely presented group, and let $X$ be the presentation complex of $G$.
Let $X' = X \vee S^2$ be the wedge sum of $X$ with the sphere.
Then $X'$ induces a map $\alpha: S^2 \to X'$ which generates a free copy of $\mathbb{Z}[G] = \mathbb{Z}[\pi_1(X')]$ inside $\pi_2(X')$.
Here, $\mathbb{Z}[G]$ denotes the group ring of $G$ over the integers.
I don't follow with how $\alpha$ includes $\mathbb{Z}[G]$ in $\pi_2(X')$. What is the justification here?

REPLY [4 votes]: Let $\tilde Y$ denote the universal cover of $Y$. Then
$\tilde{X'}$ is formed from $\tilde X$ by gluing a copy of $S^2$ to each preimage of the basepoint of $X$. Then by the Hurewicz isomorphism, we have
$$ \pi_2(X') = \pi_2(\tilde{X'}) = H_2(\tilde X) \oplus \bigoplus_{g \in \pi_1(X)} \mathbb Z = \pi_2(X) \oplus \mathbb Z[\pi_1(X)],$$
since the classes of the glued spheres are independent and generate a summand of $H_2(\tilde X)$, and are freely permuted by $\pi_1(X)$.