Analysis of the Gibbs Sampler for Gaussian Hierarchical Models via multigrid decomposition
Giacomo Zanella\footnote{ Department of Decision Sciences
BIDSA and IGIER
Bocconi University
via Roentgen 1
20136 Milan
Italy. \href{mailto:giacomo.zanella@unibocconi.it
%\fotnote{Abstract to be re-written We study the convergence properties of the Gibbs Sampler in the context of posterior distributions arising from Bayesian analysis of Gaussian hierarchical models. We consider centred and non-centred parameterizations as well as their hybrids including the full family of {\em partially non-centred parameterizations. We develop a novel methodology based on multi-grid decompositions to derive analytic expressions for the convergence rates of the algorithm for an arbitrary number of layers in the hierarchy, while previous work was typically limited to the two-level case. Our work gives a complete understanding for the three-level symmetric case and this gives rise to approximations for the non-symmetric case. We also give analogous, if less explicit, results for models of arbitrary level. This theory gives rise to simple and easy-to-implement guidelines for the practical implementation of Gibbs samplers on conditionally Gaussian hierarchical models.
Gibbs_Sampler_for_Gaussian_hierarchal_models.tex
\documentclass[11pt,a4paper]{article} %{report} %

%\textwidth=18cm \hoffset=-1cm %
%\textheight=28cm \voffset=-4cm %
%\oddsidemargin 0cm
%\usepackage{harvard}
%usepackage[style=authoryear]{biblatex}
%\addbibresource{<your-bib-file.bib>} % note the .bib is required
%\bibliographystyle{imsart-nameyear}
\usepackage{natbib}
\bibliographystyle{plainnat}%unsrtnat
\usepackage{stmaryrd}
\usepackage{enumitem}
\usepackage{graphicx}
\usepackage{caption}
\usepackage{subcaption}
\usepackage{amssymb}
\usepackage{bbm}
\usepackage{verbatim}
%\usepackage{enumerate}
\usepackage[hyperindex]{hyperref}
\usepackage{amsmath, amsfonts, amsthm, amscd}       %formattazione matematica
\usepackage{fancybox}
\usepackage[usenames,dvipsnames,svgnames,table]{xcolor}
%\usepackage{bbold}
\usepackage[T1]{fontenc}
\setlength{\fboxsep}{1.5ex}
\newlength{\querylen}
\setlength{\querylen}{\textwidth}
\addtolength{\querylen}{-2\fboxsep}

\input{notation_Gibbs}

\usepackage{tikz}
\usepackage{tikz-qtree}
\usetikzlibrary{bayesnet}


\newcommand{\R}{\mathbb{R}}
\newcommand{\PR}{\textbf{P}}
\newcommand{\Z}{\mathbb{Z}}
\newcommand{\N}{\mathbb{N}}
\newcommand{\K}{\mathbb{K}}
\newcommand{\E}{\mathbb{E}}
\newcommand{\I}{\mathbb{I}}
\newcommand{\1}{\mathbbm{1}}
\newcommand{\sM}{\mathcal{M}}
\newcommand{\x}{\textbf{x}}
\newcommand{\X}{\textbf{X}}
\newcommand{\y}{\textbf{y}}
\newcommand{\e}{\textbf{e}}
\newcommand{\st}{\hspace{0.5mm}:\hspace{1.5mm}}
\newcommand{\query}[1]{\par\noindent%\hfill\smallskip\medskip
 \shadowbox{\begin{minipage}[t]{\querylen} {#1} \end{minipage}}\medskip}
\theoremstyle{plain}
\newtheorem{defi}{Definition}
\newtheorem*{defi*}{Definition}
\newtheorem{rmk}{Remark}
\newtheorem*{rmk*}{Remark}
\newtheorem{theorem}{Theorem}
\newtheorem{prop}{Proposition}
\newtheorem{condition}{Condition}
\newtheorem{ass}{Assumption}
\newtheorem{ex}{Example}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{coroll}[theorem]{Corollary}
\newtheorem{claim}{Claim}
\newtheorem{?}{Question}
\newtheorem{model}{Model}
\newtheorem{taggedmodelx}{Model}
\newenvironment{taggedmodel}[1]
 {\renewcommand\thetaggedmodelx{#1}\taggedmodelx}
 {\endtaggedmodelx}
\newtheorem{problem}{Problem}
\newtheorem{taggedalgorithmx}{Algorithm}
\newenvironment{taggedalgorithm}[1]
 {\renewcommand\thetaggedalgorithmx{#1}\taggedalgorithmx}
 {\endtaggedalgorithmx}
\newtheorem{taggedsamplerx}{Sampler}
\newenvironment{taggedsampler}[1]
 {\renewcommand\thetaggedsamplerx{#1}\taggedsamplerx}
 {\endtaggedsamplerx}

%\usepackage{titling}
%\newcommand{\subtitle}[1]{%
%  \posttitle{%
%    \par\end{center}
%    \begin{center}\large#1\end{center}
%    \vskip0.5em}%
%}
%%\titleformat{\section}[runin]
%%  {\normalfont\Large\bfseries}{\thesection}{1em}{}
%%\titleformat{\subsection}[runin]
%%  {\normalfont\large\bfseries}{\thesubsection}{1em}{}
%\usepackage{titlesec}
%\titleformat{\section}[runin]
%  {\normalfont\large\bfseries}{\thesection}{1em}{}
%\titleformat{\subsection}[runin]
%{\normalfont\itshape}{\thesubsection}{1em}{}
%\titleformat{\subsubsection}[runin]
%{\normalfont\itshape}{\thesubsubsection}{1em}{}

%\renewcommand\fotnote[2][]{\relax}
\begin{document}
\title{
Analysis of the Gibbs Sampler for Gaussian Hierarchical Models via multigrid decomposition
}
%\date{October 31, 2014}
\author{Giacomo Zanella\footnote{
Department of Decision Sciences, BIDSA and IGIER, Bocconi University, via Roentgen 1, 20136
Milan, Italy.
\href{mailto:giacomo.zanella@unibocconi.it}{giacomo.zanella@unibocconi.it}
}
\ and Gareth Roberts\footnote{
Department of Statistics, University of Warwick, Coventry, CV4 7AL, UK.
\href{mailto:gareth.o.roberts@warwick.ac.uk}{gareth.o.roberts@warwick.ac.uk}
}}
\maketitle
\begin{abstract}
%\fotnote{Abstract to be re-written}
We study the convergence properties of the Gibbs Sampler in the context of posterior distributions arising from Bayesian analysis of Gaussian hierarchical models. We consider
centred and non-centred parameterizations as well as their hybrids including the full family of
{\em partially non-centred} parameterizations.
We develop a novel methodology based on multi-grid decompositions to derive analytic expressions for the convergence rates of the algorithm for an arbitrary number of layers in the hierarchy, while previous work was typically limited to the two-level case.
Our work gives a complete understanding for the three-level symmetric case and this 
gives rise to approximations for the non-symmetric case. We also give analogous, if less explicit,
results for models of arbitrary level.
This theory gives rise to simple and easy-to-implement guidelines for the practical
implementation of Gibbs samplers on conditionally Gaussian hierarchical models.
\end{abstract}

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

MCMC is established as the computational workhorse of most Bayesian statistical analyses for complex models.  For conditionally conjugate
models, the Gibbs sampler remains the most natural algorithm of choice, particularly for hierarchical models which have proved
arguably the most popular general class of models for Bayesian analysis. However convergence can be a major issue and it has long been recognised that reparameterisation can drastically affect the convergence rate and hence efficiency of the sampler, see for example \cite{hills1992parameterization}.

%Hierarchical models are now widely adopted within Bayesian Statistical modelling for the componentwise construction of complex models. 
The popularity of hierarchical models
 stems partly from their flexibility, interpretability, and the fact that easily implementable MCMC algorithms such as the Gibbs sampler can be used to simulate their posterior distributions \citep{gelsmi90,smirob93}. Moreover, empirical experience with Gibbs samplers for hierarchical models
(see for example \cite{GelfandSahuCarlin1995}) suggests that the Gibbs sampler's convergence properties are aided substantially by  hierarchical structure, % of the model under investigation. 
suggesting a synbiosis between hierarchical models and the Gibbs sampler. While there is some  theoretical work to investigate this phenomenon (see for example \cite{RobertsSahu1997,mengvand,RobertsPapaspiliopoulosSkold2003,jones2004sufficient, PapaspiliopoulosRobertsSkold2007,mengyu,bass2016comparison})  comparatively little is known about the link between hierarchical models and MCMC convergence, and almost nothing
% is understood about this connection 
for models of hierarchical depth greater than two. The present paper offers the first steps towards such an understanding.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
We shall couch all our results in terms of $L^2$ rates of convergence. Specifically,
let $(\bbeta(s))_{s=1,2,\dots}$ be a Markov chain with stationary distribution $\pi$ and transition operator defined by $P^sf(\bbeta(0))=\E[f(\bbeta(s))|\bbeta(0)]$.
The \emph{rate of convergence} $\rho(\bbeta(s))$ associated to $(\bbeta(s))_{s=1,2,\dots}$ is defined as the smallest number $\rho$ such that for all $r>\rho$
\begin{equation}\label{eq:rates_of_conv}
\lim_{s\rightarrow\infty}
\frac{\|P^sf-\E_\pi[f]\|_{L^2(\pi)}}{r^s}
=0
\qquad
\forall 
f\in L^2(\pi)\,,
\end{equation}
where $L^2(\pi)$ denotes the space of square $\pi$-integrable functions, $\|\cdot\|_{L^2(\pi)}$ is its associated $L^2$-norm and $\E_\pi[f]=\int f\, d\pi$ is the expectation of $f$ with respect to $\pi$.
The rate of convergence $\rho(\bbeta(s))$ characterizes the speed at which $(\bbeta(s))_{s=1,2,\dots}$ converges to its stationary distribution $\pi$, with a simple argument giving that if
$$
T = \min \{s;\ 
\|P^sf-\E_\pi[f]\|_{L^2(\pi)} \le \epsilon \}
 $$ 
 then $T = O\left(\frac{1}{-\log(\rho)}\right)$.

%Although rates of convergence could in theory provide very informative and useful characterizations of convergence of MCMC algorithms, they can rarely be evaluated for practically used models. 
%The most notable exception is probably the case of Gibbs Sampler schemes with two components.
%In this work we derive simple, analytic expressions for the rates of convergence of the Gibbs Sampler for multilevel linear models with 3 or more components (such as model \ref{model:S3}), combining previous work \citep{.amit:1991:stochastic+convergence+distributions,.amit:1996:properties+convergence+perturbations,RobertsSahu1997} with a novel approach based on multi-grid analysis.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



The theoretical results given in this paper apply to the Gibbs sampler for Gaussian hierarchical models where we are able to extend and improve substantially on existing literature \citep{RobertsSahu1997,mengyu,bass2016comparison}
both in terms of generality of hierarchical structure and the availability of explicit rates. We shall study the connections between rates of convergence and parameterisation, for instance studying different options for centred and non-centred (also called {\em sufficient} and {\em ancillary} in \cite{mengyu}) parameterisations and their links to explicit rates of convergence of their corresponding Gibbs samplers. 

In general, the Gibbs sampler can be elegantly described in terms of orthogonal projections \citep{.amit:1991:stochastic+convergence+distributions,.amit:1996:properties+convergence+perturbations,diaconis2010},
although it is rarely possible to extract practical convergence information for
complex Gibbs samplers %emanating from complex Bayesian models 
from this theory. However one important family in which explicit and useful results concerning Gibbs sampler rates of converge are available is the family of multivariate Gaussian distributions. This paper will take advantage of this explicit
framework to study rates of convergence of the Gibbs sampler as they vary with
the hierarchical structure of the model. Our study goes substantially beyond existing studies of this problem \cite{RobertsSahu1997,mengyu,bass2016comparison,GaoOwen2016} and leads to simple and practical recommendations for parameterisation choice in Gibbs sampling for hierarchical models.

%The restriction of our work to Gaussian targets is clearly restrictive. However
Our results can be readily used for models which are conditionally Gaussian, for example for the very popular context (see for example \cite{gelman2013bayesian}) where there exist unknown variances at various levels of the hierarchy.
%\fotnote{Done Do we want to stress more the importance of this class of models (indeed one of the most widely used class of models by practitioners employing Bayesian techniques)?}
In that case our optimality results describe the optimal updating strategies for the hierarchical mean structure conditional on the variances.  For such Gaussian hierarchical models of sufficient complexity, block updating of the entire Gaussian component together is typically infeasible
($O(n^3)$ in the dimension ($n$) of the Gaussian to be updated).


Furthermore, the Gaussian linear models we study epitomise Bayesian multilevel modeling more generally, in the sense that many more elaborate models share the qualitative features of these more simple models.
Moreover, many non-Gaussian  hierarchical models can be well-approximated by Gaussian ones (for example for sufficiently large data sets), so that it is reasonable to conjecture that the qualitative conclusions (at least) of our study might remain valid when extrapolated to non-Gaussian models, rather like the analysis given in \cite{sahu1999convergence}. A detailed study of this question is left for future work.

Section \ref{sec:motiv} carefully introduces the 3-level hierarchical models we shall consider, and  provides motivating simulations. Then
in Section \ref{sec:multigrid} we shall give a  complete analysis for 3-level symmetric models (i.e.\ homogenous variances and symmetric data structure). At the heart of the analysis is a multigrid decomposition of the Gibbs sampler into completely independent Markov chains describing different levels of hierarchical granularity, Theorem \ref{thm:factorization_3}.
Although multigrid ideas have already been used in methodological contexts to design improved MCMC schemes \citep{GoodmanSokal1989,LiuSabatti2000}, to our knowledge they had never been used in theoretical contexts to study convergence rates.
%Indeed the methodology we present here is novel and provides new insight into the behavior of the Gibbs Sampler in hierarchical contexts.
We demonstrate that the slowest of these independent chains is that corresponding to the coarsest level and thus derive explicit expressions for the rates of convergence in symmetric contexts.
For 3-level non-symmetric models, we give bounds on rates based on comparisons with related symmetric models. Moereover we describe partial non-centering strategies and describe {\em bespoke} strategies
which optimise the degree of non-centering at each branch.

Section \ref{sec:k_levels} considers hierarchical models with arbitrary depth ($\ge 4$).
Using an appropriate auxiliary random walk, whose evolution through the hierarchical tree is governed by the parameters' squared partial correlations, we are able to extend the multigrid analysis to general tree structures and some non-symmetric cases.
We again demonstrate a fundamental multigrid decomposition in 
Theorem \ref{thm:factorization_general} where the coarsest level chain
converges the slowest, and we give explicit formulae for optimal partial non-centering strategies.
%\fotnote{Here we could say that the theory developed here provides an understanding of the convergence properties of Gibbs Sampler for linear hierarchical models that goes significantly beyond what currently available in terms of generality of hierarchical structure and availability of explicit rates. But maybe we already said this enough.}


%\tableofcontents
%\begin{itemize}[noitemsep,nolistsep]
%\cite{RobertsPapaspiliopoulosSkold2003,PapaspiliopoulosRobertsSkold2007} provided theoretical insight into the relationship between model parametrization and rate of convergence of the Gibbs Sampler for two level models.
%One of the contributions of this work is to provide a detailed analysis of the role of hierarchical parametrizations for models with an arbitrary number of levels, with particular attention to the three level case.
%\item The target distributions we will consider belong to the Gaussian family.
%Therefore one could use direct sampling methods rather than using Gibbs sampling strategies. 
%Nevertheless direct sampling requires $O(n^3)$ operations to invert the covariance matrix of a Gaussian distribution of dimension $n$. This makes it an impractical strategy in high-dimensional contexts typically arising in hierarchical models where the number of parameters is of the same order of the number of data-points.
%On the other hand a Gibbs sampler iteration requires $O(n)$ operations and thus is potentially more scalable.
%Indeed our results suggest that, when the optimal parametrization is chosen, the number of iterations to reach stationarity is $O(1)$ in the number of dimensions, making this kinds of strategies very appealing in big data contexts.
%\item Although this is not the main focus of this paper, the results we present can also be seen as a first step towards a rigorous characterization of the computational complexity of fitting multilevel hierarchical models (e.g. suggest linear complexity if variances are well identified).\fotnote{Maybe this bullet point is too vague and thus useless.}
%\item Overview of the organization of the paper and main results.
%\end{itemize}


\section{Three level hierarchical linear models}\label{sec:motiv}

We shall begin with a detailed study of the three-level case, giving a fairly complete analysis for the case of symmetric models,
which we shall follow up with some new results for the non-symmetric case. Therefore we shall begin by
considering the following three-level Gaussian hierarchical model.
\begin{taggedmodel}{S3}[Symmetric 3-levels hierarchical model]\label{model:S3}
Suppose
\begin{equation}\label{eq:NCP_3}
y_{ijk}=\mu+a_i+b_{ij}+\epsilon_{ijk},
\end{equation}
where $i$, $j$ and $k$ run from 1 to $I$, $J$ and $K$ respectively and $\epsilon_{ijk}$ are iid normal random variables with mean 0 and variance $\var{e}$.
We employ the standard bayesian model specification assuming 
$a_i\sim N(0,\var{a})$, $b_{ij}\sim N(0,\var{b})$ and a flat prior on $\mu$.
\end{taggedmodel}
Unless otherwise stated, we consider the variance terms $\var{a}$, $\var{b}$ and $\var{e}$ to be known.
Defining $\avec=(a_i)_{i}$, $\bvec=(b_{ij})_{i,j}$ and $\yvec=(y_{ijk})_{i,j,k}$, the Gibbs Sampler explores the posterior distribution $(\mu,\avec,\bvec)|\yvec$ by iteratively sampling from the full conditional distributions of $\mu$, $\avec$ and $\bvec$ as follows (see below for motivation of denoting such sampler as GS($1,1$)).
\begin{taggedsampler}{GS($1,1$)}\label{sampler:GSNN}
Initialize $\mu(0)$, $\avec(0)$ and $\bvec(0))$ and then iterate
\begin{enumerate}[noitemsep,nolistsep]
\item Sample $\mu(s+1)$ from $p(\mu|\avec(s),\bvec(s),\yvec)$;
\item Sample $a_i(s+1)$ from $p(a_i|\mu(s+1),\bvec(s),\yvec)$ for all $i$; 
\item Sample $b_{ij}(s+1)$ from $p(b_{ij}|\mu(s+1),\avec(s+1),\yvec)$ for all $i$ and $j$,
\end{enumerate}
where $p(\mu|\avec,\bvec,\yvec)$, $p(a_i|\mu,\bvec,\yvec)$ and $p(b_{ij}|\mu,\avec,\yvec)$ are the full conditionals of Model \ref{model:S3} (see Appendix \ref{appendix:full_conditionals} for explicit expressions).
\end{taggedsampler}

%%%%%% FIGURE: 4 parametrizations PLATE REPRESENTATION %%%%%%
\begin{figure}[h!]
\centering
\begin{subfigure}[b]{0.49\textwidth}
\centering
  \tikz{
    % nodes
     \node[obs] (y) {$y_{ijk}$};%
     \node[latent,above=of y,yshift=-0.5cm] (eta) {$\eta_{ij}$}; %
     \node[latent,above=of eta,yshift=-0.5cm] (gamma) {$\gamma_i$}; %
     \node[latent,above=of gamma,xshift=0cm,yshift=-0.5cm] (mu) {$\mu$}; %
    % plate
     \plate [inner sep=0.2cm,xshift=0cm,yshift=0.1cm] {plate1} {(y)} {$k=1,\dots,K$}; %
     \plate [inner sep=0.2cm,xshift=-0cm,yshift=0.1cm] {plate2} {(plate1) (eta)} {$j=1,\dots,J$};
     \plate [inner sep=0.2cm,xshift=-0cm,yshift=0.1cm] {plate3} {(plate2) (gamma)} {$i=1,\dots,I $}
    % edges
     \edge {eta} {y}
     \edge {gamma} {eta}
     \edge {mu} {gamma}
}   
        \caption{Fully centred parametrization}
        \label{fig:CC}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.49\textwidth}
\centering
  \tikz{
    % nodes
     \node[obs] (y) {$y_{ijk}$};%
     \node[latent,above=of y,xshift=0cm] (b) {$b_{ij}$}; %
     \node[latent,left=of b,xshift=-0.4cm] (a) {$a_i$}; %
     \node[latent,left=of a,xshift=0cm] (mu) {$\mu$}; %
    % plate
     \plate [inner sep=0.2cm,xshift=0.1cm,yshift=0.1cm] {plate1} {(y)} {$k=1,...,K$}; %
     \plate [inner sep=0.2cm,xshift=0cm,yshift=0.1cm] {plate2} {(plate1) (b)} {$j=1,...,J$};
     \plate [inner sep=0.2cm,xshift=0cm,yshift=0.1cm] {plate3} {(plate2) (a)} {$i=1,...,I $}
    % edges
     \edge {b} {y}
     \edge {a} {y}
     \edge {mu} {y}
}   
        \caption{Fully non-centred parametrization}
        \label{fig:NN}
\end{subfigure}
\begin{subfigure}[b]{0.49\textwidth}
\centering
  \tikz{
    % nodes
     \node[obs] (y) {$y_{ijk}$};%
     \node[latent,above=of y,xshift=0cm,yshift=-0.4cm] (eta) {$\eta_{ij}$}; %
     \node[latent,above=of eta,yshift=-0.4cm] (a) {$a_i$}; %
     \node[latent,left=of a,xshift=-0.7cm] (mu) {$\mu$}; %
    % plate
     \plate [inner sep=0.2cm,xshift=0cm,yshift=0.1cm] {plate1} {(y)} {$k=1,...,K$}; %
     \plate [inner sep=0.2cm,xshift=-0cm,yshift=0.1cm] {plate2} {(plate1) (eta)} {$j=1,...,J$};
     \plate [inner sep=0.2cm,xshift=-0cm,yshift=0.1cm] {plate3} {(plate2) (a)} {$i=1,...,I $}
    % edges
     \edge {eta} {y}
     \edge {a} {eta}
     \edge {mu} {eta}
}   
        \caption{Mixed parametrization: $(\mu,\avec,\etavec)$}
        \label{fig:NC}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.49\textwidth}
\centering
  \tikz{
    % nodes
     \node[obs] (y) {$y_{ijk}$};%
     \node[latent,above=of y,xshift=0cm,yshift=-0.4cm] (b) {$b_{ij}$}; %
     \node[latent,left=of b,xshift=-0.4cm] (gamma) {$\gamma_i$}; %
     \node[latent,above=of gamma,xshift=0cm,yshift=-0.2cm] (mu) {$\mu$}; %
    % plate
     \plate [inner sep=0.2cm,xshift=0.1cm,yshift=0.1cm] {plate1} {(y)} {$k=1,...,K$}; %
     \plate [inner sep=0.2cm,xshift=-0cm,yshift=0.1cm] {plate2} {(plate1) (eta)} {$j=1,...,J$};
     \plate [inner sep=0.2cm,xshift=-0cm,yshift=0.1cm] {plate3} {(plate2) (gamma)} {$i=1,...,I $}
    % edges
     \edge {b} {y}
     \edge {gamma} {y}
     \edge {mu} {gamma}
}   
        \caption{Mixed parametrization: $(\mu,\gavec,\bvec)$}
        \label{fig:CN}
\end{subfigure}
    \caption{Graphical representations of 3-levels hierarchical linear models under different parametrizations.}
    \label{fig:4parametrizations}
\end{figure}

The parametrization $(\mu,\avec,\bvec)$ induced by \eqref{eq:NCP_3} is often referred to as \emph{non-centred parametrization} (NCP) and it is contrasted with the \emph{centred para\-me\-tri\-za\-tion} (CP) obtained by replacing $a_i$ and $b_{ij}$ with $\gamma_i=\mu+a_i$ and $\eta_{ij}=\gamma_i+b_{ij}$ respectively.
Under the centred parametrization $(\mu,\gavec,\etavec)$ the model formulation becomes
\begin{equation}\label{eq:CP_3}
y_{ijk}\sim N(\eta_{ij},\var{e}),\qquad
\eta_{ij}\sim N(\gamma_i,\var{b}),\qquad
\gamma_i\sim N(\mu,\var{a}),\qquad
p(\mu)\propto1\,.
\end{equation}
Figures \ref{fig:NN} and \ref{fig:CC} provides a graphical representation of the two parametrizations.
In the $(\mu,\avec,\bvec)$ case $(1,1)$ refers to the fact that both levels 1 and 2 use a non-centred parametrization, while in the $(\mu,\gavec,\etavec)$ case $(0,0)$ indicates that both levels use a centred parametrization.
The resulting Gibbs sampler for the centred parametrization is as follows.
\begin{taggedsampler}{GS($0,0$)}\label{sampler:GSCC}
Initialize $\mu(0)$, $\gavec(0)$ and $\etavec(0))$ and then iterate
\begin{enumerate}[noitemsep,nolistsep]
\item Sample $\mu(s+1)$ from $p(\mu|\gavec(s),\etavec(s),\yvec)$;
\item Sample $\gamma_i(s+1)$ from $p(\gamma_i|\mu(s+1),\etavec(s),\yvec)$ for all $i$; 
\item Sample $\eta_{ij}(s+1)$ from $p(\eta_{ij}|\mu(s+1),\gavec(s+1),\yvec)$ for all $i$ and $j$,
\end{enumerate}
where $p(\mu|\gavec,\etavec,\yvec)$, $p(\gamma_i|\mu,\etavec,\yvec)$ and $p(\eta_{ij}|\mu,\gavec,\yvec)$ are the full conditionals induced by \eqref{eq:CP_3} (see Appendix \ref{appendix:full_conditionals} for explicit expressions).
\end{taggedsampler}
%In the $(\mu,\gavec,\etavec)$ case $(0,0)$ refers to the fact that both levels 1 and 2 use a centred parametrization, while in the $(\mu,\avec,\bvec)$ case $(1,1)$ indicates that both levels use a non-centred parametrization.
Together with the fully non-centred parametrization $(\mu,\avec,\bvec)$ and the fully centred parametrization $(\mu,\gavec,\etavec)$, one can also consider the mixed parametrizations given by $(\mu,\gavec,\bvec)$ and $(\mu,\avec,\etavec)$ and the corresponding Gibbs Sampler schemes $GS{(0,1)}$ and $GS(1,0)$.
See Figures \ref{fig:NC} and \ref{fig:CN} for graphical representations.

\subsection{Illustrative example}\label{sec:motivating_example}
As an illustrative example, we simulated data from Model \ref{model:S3} with $I=J=100$, $K=5$, $\mu=0$, $\sigma_{a}=\sigma_{e}=10$ and $\sigma_b=10^{-0.5}$.
This correspond to a scenario of high level of noise in the measurements.
We fit model \ref{model:S3} assuming the standard deviations $(\sigma_{a},\sigma_{b},\sigma_e)$ to be unknown and placing weakly informative priors, namely $\frac{1}{\var{a}}$, $\frac{1}{\var{b}}$ and $\frac{1}{\var{e}}$ a priori distributed according to an Inverse Gamma distribution with shape and rate parameters equal to $0.01$.
We compare the efficiency of the Gibbs sampling schemes corresponding to the four different parametrizations, denoting them by $GS(1,1)$, $GS(0,0)$, $GS(0,1)$ and $GS(1,0)$.
For this simple example we initialized the chains at true values of the parameters $(\mu,\avec,\bvec)$ and $(\sigma_a,\sigma_b,\sigma_e)$, which we know because we are in a simulated dataset example.
The more realistic case of starting the chains from randomly chosen states led to the same conclusions of this illustrative examples. %\fotnote{Do we need to add simulations on this in an appendix?}
Note that all the four schemes have the same starting states (modulo re-parametrization) to have a fair comparison.

Figure \ref{fig:level0_mixing} shows the mixing behaviour of the global mean $\mu$ and displays the potentially dramatic difference among mixing properties of the Gibbs Sampler under different parametrizations.
\begin{figure}[h!]
\centering
\includegraphics[width=\linewidth]{level0mixing.png}
\caption{
Mixing behaviour at level $0$ (i.e. $\mu$) under four different parametrizations.
The ranges of the y-axes are constant across parametrizations.
}
\label{fig:level0_mixing}
\end{figure}
Based on Figure \ref{fig:level0_mixing} one would certainly exclude using $GS(1,1)$
and
$GS(1,0)$ to fit the model under consideration and may be tempted to deduce that both
$GS(0,0)$ and 
$GS(0,1)$
lead to good mixing properties of the resulting chain.
However, as an additional check, a cautious practitioner may also explore the mixing of the parameters at the first level, namely $\avec$ and $\gavec$. 
Figure \ref{fig:level1_mixing} displays the behaviour of the global averages of such parameters, namely $\amean=\frac{\sum_{i} a_i}{I}$ and $\gamean=\frac{\sum_i \gamma_i}{I}$, in the first 3000 iterations.
\begin{figure}[h!]
\centering
\includegraphics[width=\linewidth]{level1mixing.png}
\caption{Mixing behaviour at level 1 (i.e.\ $\amean$ and $\gamean$) under four different parametrizations. 
The ranges of the y-axes are constant across parametrizations sharing the same parameters at level 1.}
\label{fig:level1_mixing}
\end{figure}
Again, we see a dramatic difference induced by different parametrizations and, somehow surprisingly, we see that, despite having good mixing behaviour at level 0 (i.e.\ $\mu$), $GS(0,0)$ displays very poor mixing behaviour at level 1 (i.e.\ $\gavec$).
It is then natural to explores also the mixing behaviour at level 2 and Figure \ref{fig:level2_mixing} does so again by plotting the global averages $\bmean=\frac{\sum_{ij}\beta_{ij}}{IJ}$ and $\etamean=\frac{\sum_{ij}\eta_{ij}}{IJ}$.
\begin{figure}[h!]
\centering
\includegraphics[width=\linewidth]{level2mixing.png}
\caption{
Mixing behaviour at level 2 (i.e.\ $\bmean$ and $\etamean$) under four different parametrizations. 
The ranges of the y-axes are constant across parametrizations sharing the same parameters at level 2.
}
\label{fig:level2_mixing}
\end{figure}
In this case $GS(1,1)$ and $GS(0,1)$ are the only one achieving good mixing.
Based on Figures \ref{fig:level0_mixing}, \ref{fig:level1_mixing} and \ref{fig:level2_mixing} it is natural to choose to fit the model using the sampler $GS(0,1)$ corresponding to the mixed parametrization $(\mu,\gavec,\bvec)$, as it is the only one providing a good mixing across all three levels.

This simple example shows many typical issues arising when fitting Bayesian multi-level models and raises many questions. 
For example, one would like to know what are good parameters to use to diagnose convergence, in order to avoid misleading conclusions like the one suggested by Figure \ref{fig:level0_mixing}.
In fact, while in two level model good mixing of the global hyperparameters such as $\mu$ typically indicates good global mixing, this is not true in other multi-level models. 
Indeed, it is legitimate to wonder whether diagnoses based only on the global means, like in Figures \ref{fig:level0_mixing}, \ref{fig:level1_mixing} and \ref{fig:level2_mixing}, are enough to deduce good mixing of the whole markov chain, which in our example has more than $10^4$ dimensions ($1+I+IJ$ mean components and $3$ precision components).
Below we will show that for Model \ref{model:S3}, mixing of the global means ensures mixing of the whole $(1+I+IJ)$-dimensional mean components of the chain given the variances (see e.g.\ Corollary \ref{thm:rate_eq_skeleton}).
Therefore it is enough to monitor the three global means and the three variances to ensure a reliable check of the chain mixing properties.

Even more crucially, it is desirable to have simple, robust and theoretically grounded guidance in choosing a computationally efficient parametrization, given the huge impact it can have on computational performances.
The centred parametrization is often a good choice but it can fail dramatically in case of highly noisy of weakly informative data and therefore it does not provide a "black-box" solution.
In particular, models with more than two levels clearly shows a more complex structure than the two level ones. 
For example the simple illustrative example of this section shows that in three-level models both the fully centred and fully non-centred parametrizations can fail to produce computationally efficient samplers.
The theoretical analysis developed in the next section will provide the required guidance.
For example Theorem \ref{thm:rates_3} below implies that, if 
$(I,J,K,\sigma_{a},\sigma_{b},\sigma_e)=(100,100,5,10,10^{-0.5},10)$, the $L^2$ rates of convergence (NB: numbers close to 1 mean slow convergence, see \eqref{eq:rates_of_conv} and discussion thereof) of the various Gibbs Samplers under consideration are 
$$
(\rho_{(1,1)},\rho_{(0,0)},\rho_{(0,1)},\rho_{(1,0)})
=
(0.995,0.998,0.005,0.999)\,.
$$
These numbers provide a quantitative and theoretically grounded description of the behaviour heuristically observed in this section and can be easily used to optimize performances (see e.g. Section \ref{sec:cond_optimal}).



\section{Multigrid decomposition}\label{sec:multigrid}
The basic ingredient of our analysis is the following multigrid decomposition.
Consider the four possible parametrization of Model \ref{model:S3}: 
$(\mu,\avec,\bvec)$,
$(\mu,\gavec,\etavec)$
and the mixed parametrizations $(\mu,\gavec,\bvec)$ and $(\mu,\avec,\etavec)$.
In order to provide a unified treatment, regardless of the chosen parametrization, we denote the parameters used by $(\bbeta^{(0)},\bbeta^{(1)},\bbeta^{(2)})$ and the resulting Gibbs Sampler by $GS(\bbeta)$.
For example, in the NCP case $\bbeta^{(0)}=\mu$, $\bbeta^{(1)}=\avec$, $\bbeta^{(2)}=\bvec$ and $GS(\bbeta)$ coincides with \ref{sampler:GSNN}.
First consider the map $\delta$ sending $\bbeta=(\bbeta^{(0)},\bbeta^{(1)},\bbeta^{(2)})$ to 
\begin{equation}\label{eq:multigrid_3}
%\bbeta=(\mu,\bbeta^{(1)},\bbeta^{(2)})
%\rightarrow
\delta(\bbeta)
=
\left( \begin{array}{c}
\delta^{(0)}\bbeta \\
\delta^{(1)}\bbeta \\
\delta^{(2)}\bbeta
\end{array}\right)
=
\left( \begin{array}{c}
\delta^{(0)}\bbeta^{(0)}\,,\,\delta^{(0)}\bbeta^{(1)}\,,\,\delta^{(0)}\bbeta^{(2)} \\
\delta^{(1)}\bbeta^{(1)}\,,\,\delta^{(1)}\bbeta^{(2)}\\
\delta^{(2)}\bbeta^{(2)}
\end{array}\right)\,,
\end{equation}
where, loosely speaking, $\delta^{(i)}\bbeta$ represent the increments of $\bbeta$ at the $i$-th coarseness level.
More precisely 
\begin{gather*}
\delta^{(0)}\beta^{(0)}
=
\beta^{(0)}\,,\qquad
\delta^{(0)}\bbeta^{(1)}
=
\beta^{(1)}_{\cdot}\,,\qquad
\delta^{(0)}\bbeta^{(2)}
=
\beta^{(2)}_{\cdot\cdot}\,,
%=\left(\mu,\frac{\sum_{i}\bbeta^{(1)}_i}{I},\frac{\sum_{i,j}\bbeta^{(2)}_{ij}}{IJ}\right)
\\
\delta^{(1)}\bbeta^{(1)}
=
\left(\beta^{(1)}_1-\beta^{(1)}_{\cdot},\dots,\beta^{(1)}_I-\beta^{(1)}_{\cdot}\right),\,
\delta^{(1)}\bbeta^{(2)}
=
\left(
\beta^{(2)}_{1\cdot}-\beta^{(2)}_{\cdot\cdot},\dots,
\beta^{(2)}_{I\cdot}-\beta^{(2)}_{\cdot\cdot}\right)\,,
\\
\delta^{(2)}\bbeta^{(2)}
=
\left(
\beta^{(2)}_{11}-\beta^{(2)}_{1\cdot},
\beta^{(2)}_{12}-\beta^{(2)}_{1\cdot},\dots,
\beta^{(2)}_{IJ}-\beta^{(2)}_{I\cdot},
\beta^{(2)}_{I(J-1)}-\beta^{(2)}_{I\cdot}\right)\,,
\end{gather*}
where
\begin{align*}
\beta^{(1)}_{\cdot}=\frac{\sum_{i}\beta^{(1)}_i}{I}
\,,\qquad
\beta^{(2)}_{\cdot\cdot}
=
\frac{\sum_{i,j}\beta^{(2)}_{ij}}{IJ}
\,,\qquad
\beta^{(2)}_{i\cdot}=\frac{\sum_{j}\beta^{(2)}_{ij}}{J}\,.
\end{align*}
It is easy to see that the map $\delta$ is a bijection between $\R^d$ and 
$\R^3\times (\R^{I})^*\times (\R^{I})^*\times_{i=1}^I(\R^{J})^*$, where $(\R^p)^*=\{(v_1,\dots,v_p)\in\R^p\,:\,\sum_{i=1}^pv_i=0\}$.
The dimensionality of $\delta\bbeta$ equals the one of $\bbeta$, which is $1+I+IJ$, because $\delta\bbeta$ has $3+2I+IJ$ parameters and $2+I$ constraints.
It is worth noting that the three subspaces spanned by $(\delta^{(0)}\bbeta^{(0)}\,,\,\delta^{(0)}\bbeta^{(1)}\,,\,\delta^{(0)}\bbeta^{(2)})$, 
$(\delta^{(1)}\bbeta^{(1)}\,,\,\delta^{(1)}\bbeta^{(2)})$ and 
$\delta^{(2)}\bbeta^{(2)}$
respectively do not depend on the choice of parametrization $\bbeta$.
Thus the multi-grid decomposition is intrinsic to the model, and not dependent on the particular parameterization being considered.
%\fotnote{Gareth: you said you wanted to add a sentence here on the fact that the multi-grid decomposition does not depend on the parametrization.}

%%%%%
\begin{theorem}\label{thm:factorization_3} \emph{($GS(\bbeta)$ factorizes under $\delta$)} Let $(\bbeta(s))_s$ be a Markov chain on $\R^d$ evolving according to $GS(\bbeta)$. 
Then $(\delta^{(0)}\bbeta(s))_s$, $(\delta^{(1)}\bbeta(s))_s$ and $(\delta^{(2)}\bbeta(s))_s$ are three independent Markov chains evolving according to blocked Gibbs Samplers with blocks given by
$(\delta^{(0)}\bbeta^{(0)}\,,\,\delta^{(0)}\bbeta^{(1)}\,,\,\delta^{(0)}\bbeta^{(2)})$, 
$(\delta^{(1)}\bbeta^{(1)}\,,\,\delta^{(1)}\bbeta^{(2)})$ and 
$\delta^{(2)}\bbeta^{(2)}$
respectively.
\end{theorem}

\tikzstyle{state}=[shape=ellipse,draw=blue!50,fill=blue!20]
\tikzstyle{observation}=[shape=rectangle,draw=orange!50,fill=orange!20]
\tikzstyle{lightedge}=[<-,dotted]
\tikzstyle{mainstate}=[state,thick]
\tikzstyle{mainedge}=[<-,solid]


%% GIBBS SAMPLER WITH FRACTIONAL TIMES
%\begin{figure}[htbp]
%\begin{center}
%\begin{tikzpicture}[]
%% 1st column
%\node               at (0,6) {$t$};
%\node[state] (s1_1) at (0,5) {$\bbeta^{(0)}$};
%\node[mainstate] (s2_1) at (0,4) {$\bbeta^{(1)}$};
%\node[state] (s3_1) at (0,3) {$\bbeta^{(2)}$};
%% 2nd column
%\node               at (2,6) {$t+\tfrac{1}{3}$};
%\node[mainstate] (s1_2) at (2,5) {$\bbeta^{(0)}$}
%    edge[mainedge] (s2_1)
%    edge[mainedge] (s3_1)
%;
%\node[state] (s2_2) at (2,4) {$\bbeta^{(1)}$}
%%    edge[lightedge] (s2_1)
%;
%\node[state] (s3_2) at (2,3) {$\bbeta^{(2)}$}
%%    edge[lightedge] (s3_1)
%;
%% 3rd column
%\node               at (4,6) {$t+\tfrac{2}{3}$};
%\node[mainstate] (s1_3) at (4,5) {$\bbeta^{(0)}$}
%%    edge[mainedge]  (s1_2)
%;
%\node[state] (s2_3) at (4,4) {$\bbeta^{(1)}$}
%    edge[mainedge] (s1_2)
%    edge[mainedge] (s3_2)
%;
%\node[state] (s3_3) at (4,3) {$\bbeta^{(2)}$}
%%    edge[lightedge] (s3_2)
%;
%% 4th column
%\node               at (6,6) {$t+1$};
%\node[mainstate] (s1_4) at (6,5) {$\bbeta^{(0)}$}
%%    edge[mainedge]  (s1_3)
%;
%\node[state] (s2_4) at (6,4) {$\bbeta^{(1)}$}
%%    edge[lightedge] (s2_3)
%;
%\node[state] (s3_4) at (6,3) {$\bbeta^{(2)}$}
%    edge[mainedge] (s1_3)
%    edge[mainedge] (s2_3)
%;
%\end{tikzpicture}
%\end{center}
%\caption{aa.}
%\end{figure}


%%%%PLOT EASIER 

%\begin{figure}[h!]
%\begin{center}
%\begin{tikzpicture}[]
%% 1st column
%\node               at (0,6) {$t$};
%\node[state] (s1_1) at (0,5) {$\bbeta^{(0)}$};
%\node[mainstate] (s2_1) at (0,3.5) {$\bbeta^{(1)}$};
%\node[state] (s3_1) at (0,2) {$\bbeta^{(2)}$};
%%2nd column
%\node               at (3,6) {$t+1$};
%\node[mainstate] (s1_2) at (3,5) {$\bbeta^{(0)}$}
%    edge[mainedge] (s2_1)
%    edge[mainedge] (s3_1)
%;
%\node[state] (s2_3) at (3,3.5) {$\bbeta^{(1)}$}
%    edge[mainedge] (s1_2)
%    edge[mainedge] (s3_1)
%;
%\node[state] (s3_4) at (3,2) {$\bbeta^{(2)}$}
%    edge[mainedge,bend right=45] (s1_2)
%    edge[mainedge] (s2_3)
%;
%% 3st column
%\node               at (7,6) {$t$};
%\node[state] (s1_3) at (7,5) {$\delta^{(0)}\bbeta$};
%\node[mainstate] (s2_3) at (7,3.5) {$\delta^{(1)}\bbeta$};
%\node[state] (s3_3) at (7,2) {$\delta^{(2)}\bbeta$};
%%4nd column
%\node               at (10,6) {$t+1$};
%\node[mainstate] (s1_4) at (10,5) {$\delta^{(0)}\bbeta$}
%    edge[mainedge] (s1_3)
%;
%\node[state] (s2_4) at (10,3.5) {$\delta^{(1)}\bbeta$}
%    edge[mainedge] (s2_3)
%;
%\node[state] (s3_4) at (10,2) {$\delta^{(2)}\bbeta$}
%    edge[mainedge] (s3_3)
%;
%\end{tikzpicture}
%\end{center}
%\caption{aa.}
%\end{figure}

%\begin{figure}[htbp]
%\begin{center}
%\begin{tikzpicture}[]
%% 1st column
%\node               at (0,6) {$t$};
%\node[state] (s1_1) at (0,5) {$\bbeta^{(0)}$};
%\node[mainstate] (s2_1) at (0,3.5) {$\bbeta^{(1)}$};
%\node[state] (s3_1) at (0,2) {$\bbeta^{(2)}$};
%%2nd column
%\node               at (3,6) {$t+1$};
%\node[mainstate] (s1_2) at (3,5) {$\bbeta^{(0)}$}
%    edge[mainedge] (s2_1)
%    edge[mainedge] (s3_1)
%;
%\node[state] (s2_3) at (3,3.5) {$\bbeta^{(1)}$}
%    edge[mainedge] (s1_2)
%    edge[mainedge] (s3_1)
%;
%\node[state] (s3_4) at (3,2) {$\bbeta^{(2)}$}
%    edge[mainedge,bend right=45] (s1_2)
%    edge[mainedge] (s2_3)
%;
%% 1st column
%\node               at (7,6+5) {$t$};
%\node[state] (s1_1) at (7,5+5) {$\delta^{(0)}\bbeta^{(0)}$};
%\node[mainstate] (s2_1) at (7,3.5+5) {$\delta^{(0)}\bbeta^{(1)}$};
%\node[state] (s3_1) at (7,2+5) {$\delta^{(0)}\bbeta^{(2)}$};
%%2nd column
%\node               at (10,6+5) {$t+1$};
%\node[mainstate] (s1_2) at (10,5+5) {$\delta^{(0)}\bbeta^{(0)}$}
%    edge[mainedge] (s2_1)
%    edge[mainedge] (s3_1)
%;
%\node[state] (s2_3) at (10,3.5+5) {$\delta^{(0)}\bbeta^{(1)}$}
%    edge[mainedge] (s1_2)
%    edge[mainedge] (s3_1)
%;
%\node[state] (s3_4) at (10,2+5) {$\delta^{(0)}\bbeta^{(2)}$}
%    edge[mainedge,bend right=45] (s1_2)
%    edge[mainedge] (s2_3)
%;
%% 1st column
%\node[state] (s1_1) at (7,5) {$\delta^{(1)}\bbeta^{(1)}$};
%\node[mainstate] (s2_1) at (7,3.5) {$\delta^{(1)}\bbeta^{(2)}$};
%%2nd column
%\node[mainstate] (s1_2) at (10,5) {$\delta^{(1)}\bbeta^{(1)}$}
%    edge[mainedge] (s2_1)
%;
%\node[state] (s2_3) at (10,3.5) {$\delta^{(1)}\bbeta^{(2)}$}
%    edge[mainedge] (s1_2)
%;
%% 3st column
%\node[state] (s1_3) at (7,2) {$\delta^{(2)}\bbeta^{(2)}$};
%%4nd column
%\node[mainstate] (s1_4) at (10,2) {$\delta^{(2)}\bbeta^{(2)}$}
%    edge[mainedge] (s1_3)
%;
%\end{tikzpicture}
%\end{center}
%\caption{aa.}
%\end{figure}


\newlength{\cola}
\setlength{\cola}{0cm}
\newlength{\colb}
\setlength{\colb}{1.5cm}
\newlength{\colc}
\setlength{\colc}{8.7cm}%8.2cm
\newlength{\cold}
\setlength{\cold}{11cm}%10.5cm
\newlength{\cole}
\setlength{\cole}{4cm}%4cm
\newlength{\colf}
\setlength{\colf}{6cm}%
\tikzstyle{state}=[shape=ellipse,draw=Apricot!50,fill=Apricot!20]
\tikzstyle{beta}=[state]
\tikzstyle{delta0}=[state,draw=-red!50!green!25,fill=-red!85!green!70]
\tikzstyle{delta1}=[state,draw=red!50,fill=red!25]
\tikzstyle{delta2}=[state,draw=green!40,fill=green!15]
\begin{figure}[h!]
\begin{center}
\begin{tikzpicture}[]
% 1 column
\node               at (\cola,6+2.5) {$s$};
\node[beta,inner sep=0.2pt] (s1_1) at (\cola,5+2.5) {$\bbeta^{(0)}$};
\node[beta,inner sep=0.2pt] (s2_1) at (\cola,3.5+2.5) {$\bbeta^{(1)}$};
\node[beta,inner sep=0.2pt] (s3_1) at (\cola,2+2.5) {$\bbeta^{(2)}$};
%2 column
\node               at (\colb,6+2.5) {$s+1$};
\node[beta,inner sep=0.2pt] (s1_2) at (\colb,5+2.5) {$\bbeta^{(0)}$}
    edge[mainedge] (s2_1)
    edge[mainedge] (s3_1)
;
\node[beta,inner sep=0.2pt] (s2_3) at (\colb,3.5+2.5) {$\bbeta^{(1)}$}
    edge[mainedge] (s1_2)
    edge[mainedge] (s3_1)
;
\node[beta,inner sep=0.2pt] (s3_4) at (\colb,2+2.5) {$\bbeta^{(2)}$}
    edge[mainedge,bend right=45] (s1_2)
    edge[mainedge] (s2_3)
;
% 5 column
\node               at (\cole,8.5) {$s$};
\node[delta0,inner sep=0.2pt] (s1_3) at (\cole,7.5) {$\delta^{(0)}\bbeta$};
\node[delta1,inner sep=0.2pt] (s2_3) at (\cole,6) {$\delta^{(1)}\bbeta$};
\node[delta2,inner sep=0.2pt] (s3_3) at (\cole,4.5) {$\delta^{(2)}\bbeta$};
%5 column
\node               at (\colf,8.5) {$s+1$};
\node[delta0,inner sep=0.2pt] (s1_4) at (\colf,7.5) {$\delta^{(0)}\bbeta$}
    edge[mainedge] (s1_3)
;
\node[delta1,inner sep=0.2pt] (s2_4) at (\colf,6) {$\delta^{(1)}\bbeta$}
    edge[mainedge] (s2_3)
;
\node[delta2,inner sep=0.2pt] (s3_4) at (\colf,4.5) {$\delta^{(2)}\bbeta$}
    edge[mainedge] (s3_3)
;
% 3.1 column
\node               at (\colc,5.3+4.5) {$s$};
\node[delta0,inner sep=0.2pt] (s1_1) at (\colc,4.5+4.5) {$\delta^{(0)}\bbeta^{(0)}$};
\node[delta0,inner sep=0.2pt] (s2_1) at (\colc,3.5+4.5) {$\delta^{(0)}\bbeta^{(1)}$};
\node[delta0,inner sep=0.2pt] (s3_1) at (\colc,2.5+4.5) {$\delta^{(0)}\bbeta^{(2)}$};
%4.1 column
\node               at (\cold,5.3+4.5) {$s+1$};
\node[delta0,inner sep=0.2pt] (s1_2) at (\cold,4.5+4.5) {$\delta^{(0)}\bbeta^{(0)}$}
    edge[mainedge] (s2_1)
    edge[mainedge] (s3_1)
;
\node[delta0,inner sep=0.2pt] (s2_3) at (\cold,3.5+4.5) {$\delta^{(0)}\bbeta^{(1)}$}
    edge[mainedge] (s1_2)
    edge[mainedge] (s3_1)
;
\node[delta0,inner sep=0.2pt] (s3_4) at (\cold,2.5+4.5) {$\delta^{(0)}\bbeta^{(2)}$}
    edge[mainedge,bend right=65] (s1_2)
    edge[mainedge] (s2_3)
;
% 3.2 column
\node[delta1,inner sep=0.2pt] (s1_1) at (\colc,6) {$\delta^{(1)}\bbeta^{(1)}$};
\node[delta1,inner sep=0.2pt] (s2_1) at (\colc,5) {$\delta^{(1)}\bbeta^{(2)}$};
%4.2 column
\node[delta1,inner sep=0.2pt] (s1_2) at (\cold,6) {$\delta^{(1)}\bbeta^{(1)}$}
    edge[mainedge] (s2_1)
;
\node[delta1,inner sep=0.2pt] (s2_3) at (\cold,5) {$\delta^{(1)}\bbeta^{(2)}$}
    edge[mainedge] (s1_2)
;
% 3.3 column
\node[delta2,inner sep=0.2pt] (s1_3) at (\colc,4) {$\delta^{(2)}\bbeta^{(2)}$};
%4.3 column
\node[delta2,inner sep=0.2pt] (s1_4) at (\cold,4) {$\delta^{(2)}\bbeta^{(2)}$}
    %edge[mainedge] (s1_3)
;
\end{tikzpicture}
\end{center}
\caption{
Illustration of Theorem \ref{thm:factorization_3}.
Left: transition from $\bbeta(s)$ to $\bbeta(s+1)$ in Sampler $GS(\bbeta)$. This is the common structure of a Gibbs Sampler with $3$ blocks.
Center: transition from $\delta\bbeta(s)$ to $\delta\bbeta(s+1)$ in Sampler $GS(\bbeta)$. 
There are three Markov chains, namely $(\delta^{(0)}\bbeta(s))_s$, $(\delta^{(1)}\bbeta(s))_s$ and $(\delta^{(2)}\bbeta(s))_s$, evolving independently. 
Right: more detailed structure of the evolution of $(\delta^{(0)}\bbeta(s))_s$, $(\delta^{(1)}\bbeta(s))_s$ and $(\delta^{(2)}\bbeta(s))_s$.
These are Gibbs Samplers with respectively $3$, $2$ and $1$ blocks.
%The first is a three stage Gibbs Sampler, the second a two-stage Gibbs Sampler and the last one produces independent samples of $\delta^{(2)}\bbeta^{(2)}$.
%Depending on the choice of parametrization some arrows may not be actually present (e.g. for the centred parametrization).
}
\end{figure}

\begin{rmk*}
The posterior independence of $\delta^{(0)}\bbeta$, $\delta^{(1)}\bbeta$ and $\delta^{(2)}\bbeta$ is a well-known fact following from properties of Gaussian distributions.
The remarkable fact following from Theorem \ref{thm:factorization_3} is that also the Markov chains induced by the Gibbs Sampler are independent (note that the independence of the random vector under the target measure is a necessary but not sufficient condition for the independence of a corresponding MCMC scheme).
\end{rmk*}

Despite its simplicity, Theorem \ref{thm:factorization_3} provides a powerful tool to analyze the Markov chain of interest $(\bbeta(s))_s$. 
In fact the factorization into independent Markov chains implies that the rate of convergence of $(\bbeta(s))_t$ is simply given by the worst rate of convergence among $\delta^{(0)}\bbeta(s)$, $\delta^{(1)}\bbeta(s)$ and $\delta^{(2)}\bbeta(s)$.
Interestingly, the slowest chain is always the chain at the highest level, namely $\delta^{(0)}\bbeta(s)$.
\begin{theorem}\label{thm:rates_ordering_3} (Hierarchical ordering of convergence rates) 
Let $\delta^{(0)}\bbeta(s)$, $\delta^{(1)}\bbeta(s)$ and $\delta^{(2)}\bbeta(s)$ be the Markov chains defined in Theorem \ref{thm:factorization_3}.
Then the associated convergence rates satisfy
$$
\rho(\delta^{(0)}\bbeta(s))\geq
\rho(\delta^{(1)}\bbeta(s))\geq
\rho(\delta^{(2)}\bbeta(s))=0\,.
$$
\end{theorem}

\subsection{Explicit rates of convergence}\label{sec:explicit_rates}
Theorems \ref{thm:factorization_3} and \ref{thm:rates_ordering_3} imply that the rate of convergence of the global chain $\bbeta(s)$ coincides with the one of the sub-chain $\delta^{(0)}\bbeta(s)$ sampling the global means $(\beta^{(0)},\beta^{(1)}_{\cdot},\beta^{(2)}_{\cdot\cdot})$. % is always the slowest subchain and thus it is the one deterining the speed of the global chain $\bbeta(s)$.
\begin{coroll}\label{thm:rate_eq_skeleton}(Rate of convergence of $GS(\bbeta) $)
Given the notation of Theorem \ref{thm:factorization_3},
$$
\rho(\bbeta(s))=
\rho(\delta^{(0)}\bbeta(s))\,.
$$
\end{coroll}
While the original Markov chain $(\bbeta(s))_{s\in\N}$ has $1+I+IJ$ components and thus is typically high-dimensional, the skeleton chain $(\delta^{(0)}\bbeta(s))_{s\in\N}$ is a low-dimensional chain (namely 3-dimensional) amenable to direct analysis, for example using the tools developed in \cite{RobertsSahu1997}.
Therefore using Corollary \ref{thm:rate_eq_skeleton} we can analytically evaluate the rates of convergence for the Gibbs Sampler under different parametrizations.
\begin{theorem}\label{thm:rates_3}
Given an instance of Model \ref{model:S3}, the rate of convergence of the four Gibbs Sampler schemes $GS(0,0)$, $GS(1,1)$, $GS(0,1)$ and $GS(1,0)$ are given by
\begin{align*}
\rho_{(0,0)}=&
1-\frac{\tvar{a}}{\tvar{a}+\tvar{b}}\,\frac{\tvar{b}}{\tvar{b}+\tvar{e}}\,,
\\
\rho_{(1,1)}=&
\max\left\{\frac{\tvar{a}}{\tvar{a}+\tvar{e}},\frac{\tvar{b}}{\tvar{b}+\tvar{e}}\right\}\,,
\\
\rho_{(0,1)}=&
1-\frac{\tvar{a}}{\tvar{a}+\tvar{e}}\,\frac{\tvar{e}}{\tvar{b}+\tvar{e}}\,,
\\
\rho_{(1,0)}=&
\max\left\{\frac{\tvar{a}}{\tvar{a}+\tvar{b}},\frac{\tvar{e}}{\tvar{b}+\tvar{e}}\right\}\,,
\end{align*}
where $\tvar{a}=\frac{\var{a}}{I}$, $\tvar{b}=\frac{\var{b}}{IJ}$ and $\tvar{e}=\frac{\var{e}}{IJK}$.
\end{theorem}
Theorem \ref{thm:rates_3} provides explicit and informative formulas regarding the interaction between choice of parametrization and resulting efficiency of the Gibbs Sampler for Model .
\begin{figure}[h!]
\centering
\includegraphics[width=\linewidth]{Rates_3_levels.png}
\caption{Plot of rates of convergence for three-levels Gaussian hierarchical models under different parametrizations.
 Color levels correspond to values of $\log(1-\rho)$, where $\rho$ is the rate of convergence, as a function of $\log(\tvar{a})$ and $\log(\tvar{b})$ for fixed $\log(\tvar{e})=0$.
}
\label{fig:rates_3}
\end{figure}
Figure \ref{fig:rates_3} summarizes graphically the dependence of the converge rates of different parametrizations from the values of the variances of various levels.

\subsection{Conditionally optimal parametrization}\label{sec:cond_optimal}
A natural and practically relevant question is what is the optimal parametrization (among the four possible choices $(\mu,\avec,\bvec)$, $(\mu,\gavec,\bvec)$,$(\mu,\avec,\etavec)$ and $(\mu,\gavec,\etavec)$) as a function of the normalized variance components $(\tvar{a},\tvar{b},\tvar{e})$.
Using the formulas of Theorem \ref{thm:rates_3} we can obtain a simple and easily interpretable answers.
\begin{coroll}[Optimal parametrization for Model \ref{model:S3}]\label{coroll:optimal_param_3levels}
The rate of convergence of the Gibbs Sampler targeting Model \ref{model:S3} is minimized by the following choice of parametrization:
\begin{itemize}
%\item use a centred parametrization $\etavec$ at the lowest level if and only if $\tvar{b}\geq\tvar{e}$,
\item if $\tvar{b}\geq\tvar{e}$ use a centred parametrization $\etavec$ at the lowest level, while if $\tvar{b}<\tvar{e}$ use a non-centred parametrization $\bvec$,
\item if $\tvar{a}\geq\tvar{b}+\tvar{e}$ use a centred parametrization $\gavec$ at the middle level, while if $\tvar{a}<\tvar{b}+\tvar{e}$ use a non-centred parametrization $\avec$.
\end{itemize}
\end{coroll}
The table in Figure \ref{table:optimal_param_3levels} provides a graphical representation of the decision process.
\begin{figure}[h!]
\begin{center}
\begin{tabular}{ l | c | c }
  & $\tvar{a}\geq \tvar{b}+\tvar{e}$ & $\tvar{a}< \tvar{b}+\tvar{e}$ \\
  \hline			
 $\tvar{b}\geq \tvar{e}$ & $(\mu,\gavec,\etavec)$ & $(\mu,\avec,\etavec)$ \\
  \hline			
 $\tvar{b}< \tvar{e}$  & $(\mu,\gavec,\bvec)$ & $(\mu,\avec,\bvec)$ \\
  \hline  
\end{tabular}\caption{Optimal parametrization for 3-levels hierarchical models as a function of the normalized variance components.}\label{table:optimal_param_3levels}
\end{center}
\end{figure}
This simple rule guarantees that the resulting Gibbs Sampler has a rate of converges smaller than $\frac{2}{3}$ (thus guaranteeing a high sampling efficiency for fixed variances).
The optimized convergence rate equals $\frac{2}{3}$ if and only if $\tvar{a}=\tvar{b}+\tvar{e}$ and $\tvar{b}=\tvar{e}$ (in which case all parametrizations are equivalent).

Figure \ref{table:optimal_param_3levels} implies that the choice of parametrization of a given level (i.e.\ whether it is computationally convenient to use a center or non-center para\-me\-tri\-zation) depends on the ratio between the normalized variance at the level under consideration and the sum of the normalized variances of the levels below.
This results extend previous intuition for the two-level case (e.g.\ \cite{RobertsPapaspiliopoulosSkold2003}) to deeper hierarchical levels (in this case three levels).
%\fotnote{Results needs to be discussed more carefully, for example relate to the heuristic in Gelfand et al}

Corollary \ref{coroll:optimal_param_3levels} allows for simple and effective strategies to ensure high sampling efficiency in practical implementations of Gibbs Sampling for Model \ref{model:S3} in the case of unknown variances.
Common implementations choose a parametrization of the gaussian component, say $\bbeta=(\bbeta^{(0)},\bbeta^{(1)},\bbeta^{(2)})$, and alternate sampling $\bbeta|(\sigma_a,\sigma_b,\sigma_e)$ with $GS(\bbeta^{(0)},\bbeta^{(1)},\bbeta^{(2)})$ and $(\sigma_a,\sigma_b,\sigma_e)|\bbeta$ (for which direct sampling is possible and computationally cheap). 
Indeed it is sufficient to include an ``if'' statement in the Gibbs Sampling implementation, after sampling $(\sigma_a,\sigma_b,\sigma_e)|\bbeta$, to choose the optimal parametrization $\bbeta$ given $(\sigma_a,\sigma_b,\sigma_e)$ according to Figure \ref{table:optimal_param_3levels}. 
This will ensure that the sampling step $\bbeta|(\sigma_a,\sigma_b,\sigma_e)$ will have a high efficiency.
This kind of strategies will be particularly robust and effective when there is little prior information on the values of $(\sigma_a,\sigma_b,\sigma_e)$, or during the transient phase of the MCMC when $(\sigma_a,\sigma_b,\sigma_e)$ may assume unexpected values, or when the variance components are not well identified in the posterior.
In all these situations it will be difficult to choose a priori a good parametrization and therefore optimizing it on-the-fly will be a much more robust and effective strategy.
Note that the process of choosing the optimal parametrization (Figure \ref{table:optimal_param_3levels}) has a computational cost which is negligible compared to a Gibbs Sampling sweep.


\subsection{Bounds for centred parametrization in the non-symmetric case}\label{sec:bounds}
The approach of Section \ref{sec:multigrid} provides an elegant decomposition into independent Markov chains that allows to evaluate analytically the convergence rates of the Gibbs Sampler for Model \ref{model:S3}.
However, Model \ref{model:S3} has two over-restrictive assumptions. First it assumes the variances terms $\var{b}$ and $\var{e}$ to be constant across the hierarchy (i.e. not to depend on $i$ and $j$) and secondly, more crucially, it assumes the study under consideration to be perfectly balanced (i.e. $J$ and $K$ not depending on $i$ and $j$).
In this section we consider the following more general case.
\begin{taggedmodel}{NS3}[Non-symmetric 3-levels hierarchical model]\label{model:3_nonsymmetric}
Consider a more general 3-levels model with centred parametrization
\begin{align*}
\qquad p(\mu)&\propto 1&\\
\qquad \gamma_i&\sim N(\mu,\var{a}) &&i=1,\dots,I,\\
\qquad \eta_{ij}&\sim N(\gamma_i,\var{b,i}) &&j=1,\dots,J_i,\\
\qquad y_{ijk}&\sim N(\eta_{ij},\var{e,ij}) &&k=1,\dots,K_{i,j},
\end{align*}
where variance components are assumed to be known.
\end{taggedmodel}
Although the multigrid factorization of Theorem \ref{thm:factorization_3} does not apply directly to Model \ref{model:3_nonsymmetric}, it can still be used to obtain upper and lower bounds on the rates of convergence.
%%%%%
%\begin{theorem}\label{thm:rates_3_uneven_CP}
%Given an instance of Model \ref{model:3_nonsymmetric} such that 
%\begin{equation}\label{eq:assumption_for_bound}
%\min_{i=1,\dots,I}
%\frac{J_i\precision{b,i}}{\precision{a}+J_i\precision{b,i}}
%\geq
%\max_{i=1,\dots,I}
%\frac{J_i\precision{b,i}}{\precision{a}+J_i\precision{b,i}}
%\frac{1}{J_i}\sum_{j=1}^{J_i}
%\frac{\precision{b,i}}{\precision{b,i}+K_{ij}\precision{e,ij}}
%\end{equation}
%the rate of convergence of the Gibbs Sampler with centred parametrization $(\mu,\gavec,\etavec)$ satisfies
%$$
%%1-
%%R_{a,b}R_{b,e}
%%\;\leq\;
%\rho_{(\mu,\gavec,\etavec)}
%\;\leq\;
%1-
%\frac{1}{I}
%\sum_{i=1}^{I}
%\frac{J_i\precision{b,i}}{\precision{a}+J_i\precision{b,i}}
%+
%\max_{i=1,\dots,I}
%\frac{J_i\precision{b,i}}{\precision{a}+J_i\precision{b,i}}
%\frac{1}{J_i}\sum_{j=1}^{J_i}
%\frac{\precision{b,i}}{\precision{b,i}+K_{ij}\precision{e,ij}}\,.
%$$
%\end{theorem}
\begin{theorem}\label{thm:rates_3_uneven_CP}
%\footnote{Giacomo: check whether somewhere in the paper we are jumping from variances to precisions without saying it!!}
Given an instance of Model \ref{model:3_nonsymmetric} we define
$$
\ratio{a}{b}^{(i)}
=
%\frac{J_i\precision{b,i}}{\precision{a}+J_i\precision{b,i}}
%\frac{\var{a}}{\var{a}+\frac{\var{b}}{J_i}},
\frac{\var{a}}{\var{a}+J_i^{-1}\var{b}},
\quad\hbox{and}\quad
\ratio{e}{b}^{(i)}=\frac{1}{J_i}\sum_{j=1}^{J_i}
%\frac{\precision{b,i}}{\precision{b,i}+K_{ij}\precision{e,ij}}\,.
\frac{K_{ij}^{-1}\var{e,ij}}{\var{b,i}+K_{ij}^{-1}\var{e,ij}}\,.
$$
Then if 
\begin{equation}\label{eq:assumption_for_bound}
\min_{i=1,\dots,I}
\ratio{a}{b}^{(i)}
\geq
\max_{i=1,\dots,I}
\ratio{a}{b}^{(i)}
\ratio{e}{b}^{(i)}
\end{equation}
the rate of convergence of the Gibbs Sampler with centred parametrization $(\mu,\gavec,\etavec)$ satisfies
$$
\rho%_{(\mu,\gavec,\etavec)}
\;\leq\;
1-
\frac{1}{I}
\sum_{i=1}^{I}
\ratio{a}{b}^{(i)}
+
\max_{i=1,\dots,I}
\ratio{a}{b}^{(i)}
\ratio{e}{b}^{(i)}
$$
\end{theorem}
%%%%
%\begin{theorem}\label{thm:rates_3_uneven_CP}
%Given an instance of Model \ref{model:3_nonsymmetric}, the rate of convergence of the Gibbs Sampler with centred parametrization $(\mu,\gavec,\etavec)$ satisfies
%$$
%1-
%R_{a,b}R_{b,e}
%\;\leq\;
%\rho_{(\mu,\gavec,\etavec)}
%\;\leq\;
%1-
%r_{a,b}r_{b,e},
%$$
%where
%$$
%r_{a,b}
%=
%\min_{i}\frac{J_{i}\precision{b,i}}{\precision{a}+J_{i}\precision{b,i}}
%\,,
%\quad
%R_{a,b}
%=
%\max_{i}\frac{J_{i}\precision{b,i}}{\precision{a}+J_{i}\precision{b,i}}
%$$
%and
%$$
%r_{b,e}
%=
%\min_{i}\frac{1}{J_i}\sum_{j=1}^{J_i}\frac{K_{ij}\precision{e,ij}}{\precision{b,i}+K_{ij}\precision{e,ij}}
%\,,
%\quad
%R_{b,e}
%=
%\max_{i}\frac{1}{J_i}\sum_{j=1}^{J_i}\frac{K_{ij}\precision{e,ij}}{\precision{b,i}+K_{ij}\precision{e,ij}}
%.
%$$
%\end{theorem}
%%%%
The results of Theorem \ref{thm:rates_3_uneven_CP} suggest that as the number of datapoints increase the efficiency of the Gibbs sampler with centred parametrization increases.
In fact, as $K_{ij}$ increases assumption \eqref{eq:assumption_for_bound} is eventually satisfied and the bound on the convergence rate goes to $0$ as $J_i$ and $K_{ij}$ increase.
Theorem \ref{eq:assumption_for_bound} provides rigorous theoretical support and characterization of the well known fact that the centred parametrization is to be preferred in contexts of large and informative datasets \citep{GelfandSahuCarlin1995,RobertsPapaspiliopoulosSkold2003}.
%%%%


\subsection{Handling heterogeneity with bespoke parametrizations}\label{sec:bespoke}
In non-symmetric cases such as Model \ref{model:3_nonsymmetric} all parametrizations previously considered, namely $(\mu,\avec,\bvec)$, $(\mu,\gavec,\etavec)$, $(\mu,\gavec,\bvec)$ and $(\mu,\avec,\etavec)$, may be computationally inefficient.
In this more general scenario, the computationally optimal strategy will involve centering some branches of the hierarchy and non-centering others. We will refer to these strategies as \emph{bespoke parametrizations}.
In order to design a bespoke parametrization, a question of particular interest is whether the optimal choice on each branch of the hierarchy can be taken independently of other branches or not. In the former case, the optimization process would have a local nature and the resulting algorithm could be implemented easily and efficiently, while in the latter case (optimal decision on each branch influenced by other branches) the global nature of the optimization problem would impose implementation challenges.

While the rate of convergence of the general three level case (Model \ref{model:3_nonsymmetric}) is intractable, the two-level case offers enough tractability and yet richness to provide insight in the choice of optimal heterogeneous parametrizations.
We therefore consider the following non-symmetric 2-levels model (which we describe in terms of precisions rather than variances for notational convenience).
\begin{taggedmodel}{NS2}[Non-symmetric 2-levels hierarchical model]\label{model:2_nonsymmetric}
Consider the following 2-levels model with centred parametrization
%\begin{align*}
%\qquad p(\mu)&\propto 1&\\
%\qquad \gamma_i&\sim N(\mu,\var{a}) &&i=1,\dots,I,\\
%\qquad y_{ij}&\sim N(\gamma_{i},\var{e,i}) &&j=1,\dots,J_{i},
%\end{align*}
\begin{align*}
\qquad p(\mu)&\propto 1&\\
\qquad \gamma_i&\sim N(\mu,1/\precision{a}) &&i=1,\dots,I,\\
\qquad y_{ij}&\sim N(\gamma_{i},1/\precision{e,i}) &&j=1,\dots,J_{i},
\end{align*}
where the precision components $(\precision{a},(\precision{e,i})_i)$ are assumed to be known.
\end{taggedmodel}
%\begin{model}\label{model:2_nonsymmetric}(Non-symmetric 2-levels hierarchical model) 
%Consider a more general 3-levels model with centred parametrization
%\begin{align*}
%\qquad p(\mu)&\propto 1&\\
%\qquad \gamma_i&\sim N(\mu,\var{a}) &&i=1,\dots,I,\\
%\qquad y_{ij}&\sim N(\gamma_{i},\var{e,i}) &&j=1,\dots,J_{i},
%\end{align*}
%where variance components are assumed to be known.
%\end{model}
%\begin{model}\label{model:S3}(Symmetric 3-levels hierarchical model) 
%Suppose
%\begin{equation}\label{eq:NCP_3}
%y_{ijk}=\mu+a_i+\epsilon_{ij},
%\end{equation}
%where $i$, $j$ and $k$ run from 1 to $I$, $J$ and $K$ respectively and $\epsilon_{ijk}$ are iid normal random variables with mean 0 and variance $\var{e}$.
%We employ the standard bayesian model specification assuming 
%$a_i\sim N(0,\var{a})$, $b_{ij}\sim N(0,\var{b})$ and a flat prior on $\mu$.
%\end{model}
\cite{RobertsPapaspiliopoulosSkold2003}
studied the symmetric version of Model \ref{model:2_nonsymmetric},
where $J_i=J$ and 
$\precision{e,i}=\precision{e}$ for all 
$i$ and some fixed $J$ and $\precision{e}$.
They showed that the rates of convergence induced by the centred and non-centred parametrizations are given respectively by 
\begin{equation}\label{eq:rate_2_symmetric}
\rho_{CP}=\frac{\precision{a}}{\precision{a}+\tprec{e}}
\quad\hbox{ and }\quad
\rho_{NCP}=\frac{\tprec{e}}{\precision{a}+\tprec{e}}
\,,
\end{equation}
where $\tprec{e}=J\precision{e}$.
%\begin{equation}\label{eq:rate_2_symmetric}
%\rho_{CP}=\frac{\tvar{e}}{\var{a}+\tvar{e}}
%\quad\hbox{ and }\quad
%\rho_{NCP}=\frac{\var{a}}{\var{a}+\tvar{a}}
%\,,
%\end{equation}
%where $\tvar{e}=\frac{\var{e}}{J}$.
The following Theorem provides an extension to the general non-symmetric case.
We consider Model \ref{model:2_nonsymmetric} with a bespoke parametrization 
$(\mu,\beta_1,\dots,\beta_I)$ defined by $I$ indicators $(\lambda_1,\dots,\lambda_I)\in\{0,1\}^I$ as $\beta_i=\gamma_i-\lambda_i\mu$, meaning that $\lambda_i$ equals 0 if component $i$ is centred and 1 if it is non-centred.
%%%%%
\begin{theorem}\label{thm:bespoke_parametrizations_2}
The rate of convergence of the Gibbs Sampler targeting Model \ref{model:2_nonsymmetric} with bespoke parametrization given by $(\lambda_1,\dots,\lambda_I)\in\{0,1\}^I$ is
\begin{equation}\label{eq:rate_2_heterogenous}
\rho_{\lambda_1\dots \lambda_I}
=
\frac{
\sum_{i\,:\,\lambda_i=1}\tprec{i}\frac{\tprec{i}}{\tprec{i}+\precision{a}}
+
\sum_{i\,:\,\lambda_i=0}\precision{a}\frac{\precision{a}}{\tprec{i}+\precision{a}}
}{
\sum_{i\,:\,\lambda_i=1}\tprec{i}+\sum_{i\,:\,\lambda_i=0}\precision{a}
}\,,
\end{equation}
where $\tprec{i}=J_i\precision{e,i}$.
\end{theorem}

Equation \eqref{eq:rate_2_heterogenous} shows that in the non-symmetric case, the GS rate of convergence is given by a weighted average of the precision ratios $\frac{\precision{a}}{\tprec{i}+\precision{a}}$ and $\frac{\tprec{i}}{\tprec{i}+\precision{a}}$ depending on whether each component is centred or not.
This has clear analogies with the symmetric case in \eqref{eq:rate_2_symmetric}.
The weights in the average of \eqref{eq:rate_2_heterogenous} are themselves a function of $(\lambda_1,\dots,\lambda_I)$, thus introducing dependence across components in terms of centering and the overall convergence rate.
However, starting from \eqref{eq:rate_2_heterogenous} we can see that for any $i\in\{1,\dots,I\}$
\begin{equation}\label{eq:centering_i}
\rho_{\lambda_1\dots \lambda_{i-1}0\lambda_{i+1}\dots \lambda_I}-\rho_{\lambda_1\dots \lambda_{i-1}1\lambda_{i+1}\dots \lambda_I}
=
\frac{\frac{\precision{a}\tprec{i}}{\precision{a}+\tprec{i}}(\precision{a}-\tprec{i})}{(\rho_{-i}+\precision{a})(\rho_{-i}+\tprec{i})}\,,
\end{equation}
where 
$
\rho_{-i}=\sum_{\ell\neq i \,:\,\lambda_\ell=0}\tprec{\ell}\frac{\tprec{\ell}}{\tprec{\ell}+\precision{a}}
+\sum_{\ell\neq i \,:\,\lambda_\ell=1}\tprec{a}\frac{\tprec{a}}{\tprec{\ell}+\precision{a}}\geq 0
$.
Equation \ref{eq:centering_i} implies that 
$\rho_{\lambda_1\dots \lambda_{i-1}0\lambda_{i+1}\dots \lambda_I}>\rho_{\lambda_1\dots \lambda_{i-1}1\lambda_{i+1}\dots \lambda_I}$ if and only if $\precision{a}>\tprec{i}$, which in turn implies the following corollary.
\begin{coroll}\label{coroll:optimal_param_2_levels}
Let $\bar{\lambda}_i=\1(\precision{a}>\tprec{i})$ for all $i$ from 1 to $I$. Then 
$$
\rho_{\bar{\lambda}_1\dots \bar{\lambda}_{I}}
\leq
\rho_{\lambda_1\dots \lambda_I}
\qquad\hbox{for any }(\lambda_1\dots \lambda_I)\in\{0,1\}^I\,.
$$
\end{coroll}
Corollary \ref{coroll:optimal_param_2_levels} implies that in Model \ref{model:2_nonsymmetric}, the choice of optimal parametrization in each branch of the three can be carried out independently following the simple rule:
for each $i$ in $\{1,\dots,I\}$ use centred parametrization $\gamma_i$ if and only if $\precision{a}\leq J_i\precision{e,i}$, otherwise use a non-centred parametrization $a_i=\gamma_i-\mu$.

As mentioned before, the non-symmetric three level case (Model \ref{model:3_nonsymmetric}) is not easily tractable and therefore we do not have neat results such as Corollary \ref{coroll:optimal_param_2_levels} regarding the design of bespoke parametrizations.
However, motivated by Corollary \ref{coroll:optimal_param_3levels}, Corollary \ref{coroll:optimal_param_2_levels} and Theorem \ref{eq:PNCP_3} below we can design an heuristically motivated criteria to design efficient bespoke parametrizations.
For each $i=1,\dots,I$, let $\lambda_{i}$ be the non-centering indicator of the $i$-th parameter (i.e. $\beta_i=\gamma_i-\lambda_i\mu$) and for each $j=1,\dots,J_i$ let $\lambda_{ij}$ be the non-centering indicator of the $ij$-th parameter (i.e.\ $\beta_{ij}=\eta_{ij}-\lambda_{ij}\gamma_i$).
The heuristic choice for $(\lambda_i)_{i}$ and $(\lambda_{ij})_{ij}$ is the following:
%\begin{equation}\label{eq:bespoke_3}
%\lambda_{ij}=\1(\precision{b,i}>K_{ij}\precision{e,ij})
%\quad\hbox{and}\quad
%\lambda_{i}=\1\left(\precision{a}>\sum_{j=1}^{J_i}\Big(\frac{1}{\precision{b,i}}+\frac{1}{K_{ij}\precision{e,ij}}\Big)^{-1}\right)\,.
%\end{equation}
\begin{equation}\label{eq:bespoke_3}
\lambda_{ij}=\1\left(\var{b,i}<\frac{\var{e,ij}}{K_{ij}}\right)%K_{ij}^{-1}\var{e,ij}
\quad\hbox{and}\quad
\lambda_{i}=\1\left(\frac{1}{\var{a}}>\sum_{j=1}^{J_i}
%\Big(\var{b,i}+K_{ij}^{-1}\var{e,ij}\Big)^{-1}
\frac{1}{\var{b,i}+K_{ij}^{-1}\var{e,ij}}
\right)\,.
\end{equation}
It is easy to see that in the symmetric case equation \eqref{eq:bespoke_3} coincides with the criteria derived in Section \ref{sec:cond_optimal} from the explicit rates of convergence analysis.
In the non-symmetric case, these heuristic criteria can be motivated by observing that the non-centering indicators $\lambda_i$ and $\lambda_{ij}$ defined by \eqref{eq:bespoke_3} equal 1 or 0 depending on whether the corresponding optimal partial non-centering coefficients $\lambda^{(1,0)}_i$ and $\lambda^{(2,1)}_{ij}$ (see Theorem \ref{thm:PNCP_3_levels} below) are above or below a half.
Note also that the definition of $\lambda_{ij}$ in \eqref{eq:bespoke_3} is consistent with the non-symmetric two level case of Corollary \ref{coroll:optimal_param_2_levels}.
%\fotnote{I'm not very satisfied with the motivation for these 3-level heuristic bespoke parametrizations.}




%%%%%
%%%%%
\subsection{Optimal partial non-centering}\label{sec:PNCP_3}
The centred and non-centred parametrizations can be seen as extrema of a continuum of parametrizations, namely partially non-centred parametrizations (PNCP) \cite{}.
PNCP can be obtained from the centred (or non-centred) parametrization by linear transformations which preserve the hierarchical structure.
For example, in the symmetric context of Model \ref{model:S3}, the family of PNCP can be obtained from the centred parametrization $(\mu,\gavec,\etavec)$ as
$\beta^{(0)}=\mu$,
$\beta^{(1)}_i=\gamma_i-\lambda^{(1,0)}\mu$,
 and 
$\beta^{(2)}_{ij}=\eta_{ij}-\lambda^{(2,1)}\gamma_i-\lambda^{(2,0)}\mu$,
for some real valued non-centering parameters $(\lambda^{(1,0)},\lambda^{(2,0)},\lambda^{(2,1)})$.
In the more general non-symmetric context of Model \ref{model:3_nonsymmetric} one may allow the non-centering parameters to depend on the location in the hierarchical tree, thus obtaining
\begin{align}\label{eq:PNCP_3}
\begin{split}
\beta^{(0)}&=\mu,\\
\beta^{(1)}_i&=\gamma_i-\lambda^{(1,0)}_i\mu,\\
\beta^{(2)}_{ij}&=\eta_{ij}-\lambda^{(2,1)}_{ij}\gamma_i-\lambda^{(2,0)}_{ij}\mu,
\end{split}
\end{align}
for some real valued $(\lambda^{(1,0)}_i)_i$ and $(\lambda^{(2,0)}_{ij},\lambda^{(2,1)}_{ij})_{ij}$.
\begin{theorem}\label{thm:PNCP_3_levels}
Let $(\mu,\gavec,\etavec)$ be given by Model \ref{model:3_nonsymmetric} conditioned on some value of $\textbf{y}$. If $(\lambda^{(1,0)}_i)_i$ and $(\lambda^{(2,0)}_{ij},\lambda^{(2,1)}_{ij})_{ij}$ are defined as 
\begin{gather*}
\lambda^{(1,0)}_i
=
\frac{\frac{1}{\var{a}}}{\frac{1}{\var{a}}+\sum_{j=1}^{J_i}\frac{1}{\var{b,i}+K_{ij}^{-1}\var{e,ij}}}
\\
\lambda^{(2,0)}_{ij}
=
0
\quad\hbox{and}\quad
\lambda^{(2,1)}_{ij}
=
\frac{K_{ij}^{-1}\var{e,ij}}{\var{b,i}+K_{ij}^{-1}\var{e,ij}}\,,
\end{gather*}
then the resulting $(\beta^{(0)},\beta^{(1)}_i,\beta^{(2)}_{ij})$ is a collection of independent random variables.
\end{theorem}
Theorem \ref{thm:PNCP_3_levels} implies that appropriate partial non-centering leads to i.i.d.\ samples even in 3 level models.
While any Gaussian random vector can be factorized into independent random variables via some linear transformation (e.g.\ a Cholesky decomposition), the interesting aspect of Theorem \ref{thm:PNCP_3_levels} is that for hierarchical models such factorization can be obtained within the class of hierarchical reparametrizations.
%\fotnote{Gareth: I added this sentence, check if that's fine. Feel free to remove it, or mention fast Cholesky decomposition for sparse matrices, or message passing, or whatever you feel appropriate...}

\section{Hierarchical linear models with arbitrary number of levels}\label{sec:k_levels}
In this section we consider Gaussian hierarchical models with an arbitrary number of levels, namely $k$ levels.
We refer to the highest level of the hierarchy (i.e. the one furthest away from the data) as level 0, down to level $k-1$ being the lowest level (i.e. closest to the data).
The 3 level model of Section \ref{sec:multigrid} is a special case of the theory developed here where $k=3$.
%Previous theorems (e.g.\ Theorem \ref{thm:factorization_3}) are special cases of the general results given in this section (Theorems \ref{thm:factorization_general}).

\subsection{Model formulation}
In order to allow for more generality and keep the notation concise, in this section we will use a graphical models notation.
In particular $T$ will denote a finite tree with $k$ levels and root $t_0\in T$. 
For each node $t\in T$ we will denote by $pa(t)$ the unique parent of $t$ and by $ch(t)$ the collection of children of $t$. 
Moreover we write $s\preceq t$ and $s\succeq t$ if $s$ is respectively an ancestor or a descendent of $t$ (with $s$ and $t$ possibly being equal) while $s\prec t$ and $s\succ t$ denote the same notions with the additional condition of $s\neq t$.
For each node $t\in T$ we denote by $\ell(t)$ the level of node $t$ (i.e. its distance from $t_0$).
For each $d\in\{0,\dots,k-1\}$ we denote by $T_d=\{t\in T\,:\,\ell(t)=d\}$ the collection of nodes at level $d$. For example we have $T_0=\{t_0\}$ and $T=\cup_{d=0}^{k-1} T_d$.
Noisy observations will occur only at leaf nodes.
The collection of leaf nodes is denoted as $T_L=\{t\in T\,:\,ch(t)=\emptyset\}$.
For simplicity we assume that all leaf nodes are at level $k-1$, i.e. $T_L=T_{k-1}$, although this assumption could be easily relaxed allowing some branches to be longer than others.

\begin{taggedmodel}{NSk}[$k$-levels hierarchical model]\label{model:NSk} 
Suppose that we have a hierarchy described by a tree $T$ with observations occuring at leaf nodes 
$T_L$.
We assume the following hierarchical model
\begin{align}%\label{eq:NCP_3}
y^{(i)}_{t}&\sim N(\gamma_t,
1/\precision{t}^{(e)}
%\var{\ell}
)
&t\in T_L\,,\\%,\, i\in\{1,\dots,n_t\}
\gamma_{t}&\sim N(\gamma_{pa(t)},
1/\precision{t}
)
&t\in T\backslash{t_0},
\end{align}
where $i\in\{1,\dots,n_t\}$ with $n_t$ being the number of observed data at leaf node $t$, $(\precision{t})_{t\in T\backslash t_0}$ and $(\precision{t}^{(e)})_{t\in T_L}$ are known precision components and all normal random variables are sampled independently.
Following the standard Bayesian model specification we assume a flat prior on $\gamma_{t_0}$.
%\footnote{Giacomo: make notation consistent with precisions being defined at the child level and not at the parent level.}
\end{taggedmodel}
We are interested in sampling from the posterior distribution of $\bgamma_T=(\gamma_t)_{t\in T}$ given some observations $\textbf{y}=(y_t)_{t\in T_L}$.
The centred parametrization $\bgamma_T$ of Model \ref{model:NSk} leads to the following Gibbs Sampler.
\begin{taggedsampler}{GS($\bgamma_T$)}\label{sampler:GSkCP}
Initialize $\bgamma_T(0)$ and then iterate the following kernel:\\
For $d=0,\dots,k-1$
\begin{itemize}[noitemsep,nolistsep]
\item Sample $\gamma_t(s+1)$ from 
$p(\gamma_t|\bgamma_{T_{d-1}}(s+1),\bgamma_{T_{d+1}}(s),\yvec)$
for all $t\in T_d$,
\end{itemize}
where $p(\gamma_t|\bgamma_{T_{d-1}},\bgamma_{T_{d+1}},\yvec)=p(\gamma_t|\bgamma_{T\backslash t},\yvec)$ is the full conditional distribution of $\gamma_t$ given by Model \ref{model:NSk}.
When $d$ equals 0 or $k-1$ the levels $\bgamma_{T_{d-1}}$ and $\bgamma_{T_{d+1}}$ have to be replaced by empty sets in the conditioning. 
\end{taggedsampler}

\subsection{Non centering and hierararchical reparametrizations}\label{sec:PNCP_k_levels}
Model \ref{model:NSk}  expresses Gaussian hierarchical models using a centred pa\-ra\-me\-tri\-za\-tion.
The corresponding non-centred version is given by the following example.
%%%%
\begin{ex}[Fully non-centred parametrization]\label{ex:NCPk}
Under the same setting as Model \ref{model:NSk}, define
\begin{align*}
y^{(i)}_{t}&\sim N\Big(\sum_{t\preceq \ell}\alpha_t,
1/\precision{t}^{(e)}
%\var{\ell}
\Big)
&t\in T_L,\\
\alpha_{t}&\sim N(0,
1/\precision{t}
%\var{pa(t)}
)
&t\in T\backslash{t_0},
\end{align*}
and assume a flat prior on $\alpha_{t_0}$.
\end{ex}
%%%%
The non-centred parametrization $\balpha_T$ can be obtained as a linear tranformation of the centred version $\bgamma_T$ of Model \ref{model:NSk}.
More generally, we will consider the class of parametrizations that can be obtained by reparametrizing  $\bgamma_T$ as follows.
%%%%
\begin{defi*}[Hierarchical reparametrizations]\label{defi:hierarchical_reparametrization}
Let $\bgamma_T=(\gamma_t)_{t\in T}$ be a random vector with elements indexed by a tree $T$.
We say that $\bbeta_T=(\beta_t)_{t\in T}$ is a hierarchical (linear) reparametrization of $\bgamma_T$ if
\begin{equation}\label{eq:hierarchical_reparametrization}
\beta_t=\sum_{r\preceq t}\lambda_{tr } \gamma_r\,\qquad t\in T,
\end{equation}
for some real-valued coefficients $\Lambda=(\lambda_{tr})_{r\preceq t,t\in T}$ satisfying $\lambda_{tt}\neq 0$ for all $t\in T$.
We denote \eqref{eq:hierarchical_reparametrization} by $\bbeta_T=\Lambda\bgamma_T$.
\end{defi*}
%%%%
Analogously to Section \ref{sec:PNCP_3}, we refer to the family of hierarchical reparametrizations of $\gavec_T=(\gamma_t)_{t\in T}$ as \emph{partially non-centred parametrizations} (PNCP) of Model \ref{model:NSk}.
Note that \eqref{eq:hierarchical_reparametrization} does not span the space of all linear transformations of $\bgamma_T$.
In fact $\Lambda=(\lambda_{tr})_{r\preceq t,t\in T}$ can be thought as a $|T|\times|T|$ matrix $\Lambda=(\lambda_{tr})_{r,t\in T}$ inducing a linear transformation of $\bgamma_T$ with the additional sparsity constraint of being zero on all elements $\lambda_{tr}$ such that $r\npreceq t$.
The following Lemma shows that the definition of the class of PNCP does not depend on the starting parametrization used to formulate Model \ref{model:NSk}.
For example, one could equivalently define the class of PNCP of Model \ref{model:NSk} as the collection of hierarchical reparametrization of the non-centred parametrization $\balpha_T$ of Example \ref{ex:NCPk}.
%%%%%
\begin{lemma}\label{lemma:hierarchical_reparametrization_invertible}
If $\bbeta_T$ is a hierarchical reparame\-tri\-za\-tion of $\bgamma_T$, then also $\bgamma_T$ is a hierarchical reparametrization of $\bbeta_T$.
\end{lemma}
%%%%%
%\begin{proof}
%See Appendix \ref{proof:hierarchical_reparametrization_invertible}.
%\end{proof}
%%%%%

As for the 3-levels case we are interested in assessing the computational efficiency of the different Gibbs Sampling schemes arising from different PNCP's.
For each PNCP $\bbeta_T$ the corresponding Gibbs Sampler scheme $GS(\bbeta_T)$ is defined analogously to $GS(\bgamma_T)$.
\begin{taggedsampler}{GS($\bbeta_T$)}\label{sampler:GSk}
Initialize $\bbeta_T(0)$ and then iterate the following kernel:\\
For $d=0,\dots,k-1$
\begin{itemize}[noitemsep,nolistsep]
\item Sample $\beta_t(s+1)$ from 
$
p(\beta_t|(\bbeta_{T_p}(s+1))_{0\leq p<d},(\bbeta_{T_p}(s))_{d<p\leq k-1},\yvec)
$
for all $t\in T_d$,
\end{itemize}
where $p(\beta_t|(\bbeta_{T_p})_{0\leq p<d},(\bbeta_{T_p})_{d<p\leq k-1},\yvec)=p(\beta_t|\bbeta_{T\backslash t},\yvec)$ is the full conditional distribution of $\beta_t$ given by Model \ref{model:NSk}.
\end{taggedsampler}

%%%%%%
%%%%%%
%\subsection{Conditional independence structure of hierarchical }
Sampler \ref{sampler:GSk} is easy to implement because the family of PNCP preserves the hierarchical structure of Model \ref{model:NSk}.
In fact, any PNCP of Model \ref{model:NSk} exhibits the following conditional independence structure:
\begin{equation}\tag{H}\label{eq:hierarchical_structure}
\beta_r\bot\beta_t|\bbeta_{T\backslash\{r,t\}}
\hbox{ unless }
r\preceq t
\hbox{ or }
t\preceq r\,.
\end{equation}
Note that this is a weaker condition than the one holding for the centred parametrization $\bgamma_T$. In the latter case, the conditional independence graph corresponds exactly to the tree $T$, meaning that if $r\neq t$
\begin{equation}\tag{T}\label{eq:tree_cond_indep}
\gamma_r\bot\gamma_t|\bgamma_{T\backslash\{r,t\}} 
\hbox{ unless }
r=pa(t)
\hbox{ or }
t=pa(r)\,.
\end{equation}
Despite being weaker than \eqref{eq:tree_cond_indep}, condition \eqref{eq:hierarchical_structure} still guarantees that all parameters at the same level are conditionally independent (thus simplifying the update of $\bbeta_{T_d}|\bbeta_{T\backslash T_d}$) and that the full conditional distribution of each $\beta_t$ depends only on the descendants or ancestors of $t$.
The following Lemma and Corollary provide a simple way to check that any PNCP of Model \ref{model:NSk} satisfies \eqref{eq:hierarchical_structure}.
%%%%%
\begin{lemma}\label{lemma:H_closed}
Property \eqref{eq:hierarchical_structure} is closed under hierarchical re-parametrizations, meaning that if $\bbeta_T$ satisfies \eqref{eq:hierarchical_structure} then any hierarchical re-parametrization of $\bbeta_T$ satisfies \eqref{eq:hierarchical_structure} too.
\end{lemma}
%%%%%
%\begin{proof}
%See Appendix \ref{proof:H_closed}
%\end{proof}
%%%%%


%%%%%
\begin{coroll}\label{coroll:PNCP}
Any PNCP $\bbeta_T$ of Model \ref{model:NSk} satisfies \eqref{eq:hierarchical_structure}.
\end{coroll}
%%%%%%
%\begin{proof}
%Follows from Lemma \ref{lemma:H_closed} and equation \eqref{eq:tree_cond_indep}.
%\end{proof}
%%%%%%

%%%%
%%%%
\subsection{Symmetry assumption}\label{sec:symmetry}
To provide a full analysis of the convergence properties of Sampler \ref{sampler:GSk} we need a symmetry assumption that we now define.
Let $\rho_{tr}$ denote the partial correlation $Corr\left(\beta_t,\beta_r\Big | \bbeta_{T\backslash\{t,r\}}\right)$, namely
%\footnote{Giacomo: check that the minus is consistently kept across proofs. MAybe we could avoid having the minus for $r\neq t$ and the plus for $r=t$.}
%\footnote{Giacomo: is this notation about correlations clashing with the notation about rates of convergence?}
\begin{align*}
\rho_{tr}&=-\frac{Q_{tr}}{\sqrt{Q_{tt}Q_{rr}}}&t\neq r\,,
\end{align*}
and $\rho_{tt}=1$ for all $t$. Here $Q$ is the precision matrix of $\bbeta_T$.
Let $\textbf{X}=(X_\ell)_{\ell=0}^{k-1}$ be a random walk going through $T$ from root to leaves as follows: $X_0=t_0$ almost surely and then, for $\ell\in\{0,\dots,k-2\}$
\begin{equation}\label{eq:non_symm_markov}
P(X_{\ell+1}=t\,|\,X_{\ell}=r)=\frac{\rho^2_{tr}}{\sum_{t'\in ch(r)}\rho^2_{t'r}}\1(t\in ch(r))\,.
\end{equation}
Equation \eqref{eq:non_symm_markov} implies that at each step the Markov chain $\textbf{X}$ jumps from the current state $r$ to one of its children $t\in ch(r)$ choosing $t$ proportionally to the squared partial correlation between $\beta_r$ and $\beta_t$.
Since $\ell(X_d)=d$ almost surely for all $d\in\{0,\dots,k-1\}$ we can use the following simplified notation: for any $t$ and $r$ in $T$ we use $P(t)$, $P(t|r)$ and $P(t\cap r)$ to denote respectively  $P(X_{\ell(t)}=t)$, $P(X_{\ell(t)}=t\,|\,X_{\ell(r)}=r)$ and $P(X_{\ell(t)}=t\,\cap\,X_{\ell(r)}=r)$.

Given the above definitions, we define the following symmetry condition: there exist a $k\times k$ symmetric matrix $C=(c_{dp})_{d,p=0}^{k-1}$ such that
\begin{align}\tag{S}\label{eq:k_symmetry}
\rho_{tr}
=
&
c_{\ell(r)\ell(t)}\sqrt{P(t|r)}
& r\preceq t\,,
\end{align}
and $\rho_{tr}=0$ if $r\npreceq t$ and $t\npreceq r$.
Note that $\rho_{tr}$ is invariant to coordinate-wise rescaling of $\bbeta_T$ and therefore both property \eqref{eq:k_symmetry} and the transition kernel of $\textbf{X}$ are invariant to rescalings.
Therefore we can consider, without loss of generality, the following rescaled version of $\bbeta_T$ defined by 
\begin{align}\label{eq:rescaling}
\tbeta_t&=\beta_t\sqrt{\frac{Q_{tt}}{P(t)}}&t\in T\,.
\end{align}
Given \eqref{eq:rescaling}, condition \eqref{eq:k_symmetry} can be written, in terms of the precision matrix of $\tbbeta_T=(\tbeta_t)_{t\in T}$ as
%\footnote{Giacomo: check that the this minus is consistently kept across paper!}
\begin{align}\tag{$\tilde{\hbox{S}}$}\label{eq:k_symmetry_rescaled}
\tilde{Q}_{tt}=P(t)
\quad\hbox{and}\quad
-\tilde{Q}_{tr}
=&
c_{\ell(t)\ell(r)}P(t \cap r)
\quad \hbox{for }t\neq r
\,.
\end{align}
The rescaled version $\tbbeta_T$ will be helpful later to design an appropriate multigrid decomposition of $\bbeta_T$.
Also, it can be seen that property \eqref{eq:k_symmetry_rescaled} is closed under symmetric hierarchical reparameterizations.
\begin{defi*}[Symmetric hierarchical reparametrizations]\label{defi:symmetric_reparametrization}
We say that $\bbeta_T=\Lambda\balpha_T$ is a symmetric hierarchical reparametrization of $\balpha_T$ if the coefficients of $\Lambda=(\lambda_{tr})_{r\preceq t,t\in T}$  depend only on the levels of $r$ and $t$ in the hierarchy $T$.
\end{defi*}
%%%%%
\begin{lemma}\label{lemma:S_closed} 
Property \eqref{eq:k_symmetry_rescaled} is closed under symmetric hierarchical repa\-ra\-me\-tri\-zations, meaning that if $\tbbeta_T$ satisfies \eqref{eq:k_symmetry_rescaled} then any symmetric hierarchical repara\-me\-trization of $\tbbeta_T$ satisfies \eqref{eq:k_symmetry_rescaled} too.
\end{lemma}
%%%%%
Various special cases of Model \ref{model:NSk} satisfy assumption \eqref{eq:k_symmetry}.
For example, we now consider three cases: a fully symmetric case (both the tree $T$ and the variances $(\tau_t)_{t\in T}$ are symmetric), a weakly symmetric case (non-symmetric tree and symmetric variances) and a non-symmetric case (both tree and variances non-symmetric).

%%%%%
\begin{taggedmodel}{Sk}[Symmetric $k$-levels hierarchical model]\label{model:Sk} 
Consider the $k$-level Gaussian Hierarchical model where the observed data are generated from
\begin{align*}
y_{i_1,\dots,i_{k-1},j}&\sim N(\gamma^{(k-1)}_{i_1,\dots,i_{k-1}},1/\tau_e)
&&(i_1,\dots,i_{k-1},j)\in[I_1]\times\dots\times[I_{k-1}]\times[J]\,,
\end{align*}
where $[N]=\{1,\dots,N\}$ for any positive integer $N$.
The parameters have the following hierarchical structure: for each level $d$ from $1$ to $k-1$
\begin{align*}
\gamma^{(d)}_{i_1,\dots,i_d}&\sim N(\gamma^{(d-1)}_{i_1,\dots,i_{d-1}},1/\tau_d) 
&&(i_1,\dots,i_d)\in[I_1]\times\dots\times[I_d]\,,
\end{align*}
Here $(\precision{1},\dots,\precision{k-1},\precision{e})$ are known precisions and the root parameter $\gamma^{(0)}$ is given a flat prior $p(\gamma^{(0)})\propto 1$.
For each $d\in\{1,\dots,k-1\}$ the positive integer $I_d$ represents the number of branches from level $d-1$ to level $d$.
\end{taggedmodel}
%[[\fotnote{The following two lines are probably useless}
%The vector of parameters $\bgamma$ induced by Model \ref{model:Sk} is structured as 
%%follows:
%%\begin{gather*}
%$\bgamma=\left(\bgamma^{(0)},\bgamma^{(1)},\dots,\bgamma^{(k)}\right)$
%%\in\R^N
%%\qquad\hbox{with }N=\sum_{d=0}^{k-1}I^{(d)}\,,
%%\end{gather*}
%where 
%%$I^{(d)}=\prod_{s=0}^dI_s$ and, for 
%each $\bgamma^{(d)}$ is the vector of parameters at level $d$
%$$
%\bgamma^{(d)}=(\gamma^{(d)}_{i_1,\dots,i_d})_{i_1=1,\dots,i_d=1}^{I_1,\dots,I_d}\in\R^{\prod_{s=0}^dI_s}\,.\quad]]
%$$
It is easy to see that the posterior distribution of Model \ref{model:Sk},
conditioned on the observed data $\yvec=(y_{i_1,\dots,i_{k-1},j})_{i_1,\dots,i_{k-1},j}$, 
satisfies \eqref{eq:k_symmetry}.
In this case the random walk $\textbf{X}$ defined by \eqref{eq:non_symm_markov} coincides with the natural random walk going through $T$.
%In this case the matrix $C=(c_{dp})_{d,p=0,\dots,k}$ is tridiagonal.


\begin{ex}[Weakly symmetric case]\label{ex:weakly_symm}
Another special case of Model \ref{model:NSk} satisfying \eqref{eq:k_symmetry} is given by the case of a general tree $T$ and precision terms defined as $\tau_t=\frac{\tau_{\ell(t)}}{\prod_{s\prec t}|ch(s)|}$ for all $t\in T$ and  $\tau^{(e)}_t=\frac{\tau_e}{n_t\prod_{s\prec t}|ch(s)|}$ , where $(\tau_1,\dots,\tau_{k},\tau_e)\in\R_+^{k+1}$ are level-specific precision terms.
This is an extension of Model \ref{model:Sk} where the lack of symmetry of $T$ is compensated by appropriate variance terms.
Condition \eqref{eq:k_symmetry} can be checked by evaluating the partial correlations $(\rho_{tr})_{t,r\in T}$ of the resulting vector $\bgamma_T$.
\end{ex}

\begin{ex}[Non-symmetric cases]\label{ex:non_symm}
In both cases previously considered (Model \ref{model:Sk} and Example \ref{ex:weakly_symm}) the auxiliary Markov chain $\textbf{X}$ defined in \eqref{eq:non_symm_markov} follows a natural random walk, in the sense that at each time the chain chooses the next state uniformly at random among children nodes.
However, condition \eqref{eq:k_symmetry} is also satisfied by non-symmetric cases where $\textbf{X}$ is not a natural random walk.
In particular any instance of Model \ref{model:NSk} such that 
\begin{align}\tag{S*}\label{eq:k_symmetry_appendix}
\sum_{r\in ch(t)}\rho_{tr}^2
\;=\;
c_{\ell(t)}
\qquad\hbox{for all }
t\in T\backslash T_L\,,
\end{align}
for some $(k-1)$-dimensional vector $(c_0,\dots,c_{k-2})$ induces a posterior distribution satisfying \eqref{eq:k_symmetry}.
In fact, in the context of Model \ref{model:NSk} conditions \eqref{eq:k_symmetry_appendix} and \eqref{eq:k_symmetry} are equivalent (this can be derived from \eqref{eq:tree_cond_indep} and \eqref{eq:non_symm_markov}).
\end{ex}

The cases previously considered are expressed in terms of centred pa\-ra\-me\-tri\-za\-tion, meaning that as all the  instances of Model \ref{model:NSk} they satisfy \eqref{eq:tree_cond_indep}.
Nevertheless Lemma \ref{lemma:S_closed} shows that any symmetric hierarchical reparametrization of a vector satisfying \eqref{eq:k_symmetry_rescaled} still satisfies \eqref{eq:k_symmetry_rescaled}.
This implies, for example, that the fully non-centred version of Model \ref{model:Sk} and any mixed strategy where some level is centred and some is not centred, still satisfies \eqref{eq:k_symmetry_rescaled} (after rescaling).
%\fotnote{Exposition of this last part should be improved}

Moreover, note that the exact analysis we will now provide for the Gibbs sampler on models satisfying \eqref{eq:k_symmetry_rescaled} can be used to provide bound on general cases that do not satisfy \eqref{eq:k_symmetry_rescaled} (see for example Theorem \ref{thm:rates_3_uneven_CP}).


%%%%%
%%%%%
\subsection{Multigrid decomposition}\label{sec:multigrid_k}
We now show how to use the multigrid decomposition to analyze the Gibbs Sampler for random vectors $\bbeta_T$ satisfying \eqref{eq:hierarchical_structure} and \eqref{eq:k_symmetry}.
Our aim is to provide a transformation of $\bbeta_T$
%\footnote{Giacomo: sometimes I use $\bbeta_T$, sometimes $\bbeta_T$...Should probably fix it.}
 that factorizes the Gibbs Sampler Markov Chain into independent and more tractable sub-chains.
Similarly to Section \ref{sec:multigrid} in the following we will often denote $\bbeta_{T_d}=(\beta_t)_{t\in T_d}$ by $\bbeta^{(d)}$. 
We proceed in two steps, first introducing the averaging operators $\phi^{(p)}$ and then the residual operators $\delta^{(p)}$. 
For any $p\leq d$ the averaging operator $\phi^{(p)}:\R^{T_d}\rightarrow\R^{T_p}$ is defined as
%\footnote{Giacomo: replace $\E[\beta_{X_d}|X_p=r]$ with $\E[\beta_{X_d}|\bbeta_T,X_p=r]$ across the paper}
\begin{align}\label{eq:phi_defi}
\phi^{(p)}_r\bbeta^{(d)}
=&
\E[\beta_{X_d}|\bbeta_T,X_p=r]&r\in T_p
\nonumber\\
=&
\sum_{t\in T_d}
\beta_{t}P(t|r)
&
\end{align}
where $\textbf{X}=(X_\ell)_{\ell=0}^{k-1}$
is the Markov chain defined by \eqref{eq:non_symm_markov}.
%is a Markov chain going through $T$ from root to leaves according to a natural random walk, meaning that $X_0=t_0$ almost surely and 
%$P(X_\ell,\cdot)=\hbox{Unif}(ch(X_\ell))$
%%, or equivalently $P(X_\ell,X_{\ell+1})=\frac{1}{n_{X_\ell}}\1(X_{\ell+1}\in ch(X_\ell))$
% for $\ell\in\{1,\dots,k-1\}$.
Loosely speaking $\phi^{(p)}\bbeta^{(d)}=\E[\beta_{X_d}|\bbeta_T,X_p]$
 can be interpreted as the averages of $\bbeta^{(d)}$ at the coarseness corresponding to the $p$-th level of the hierarchy.
In particular $\phi^{(d)}\bbeta^{(d)}=\bbeta^{(d)}$ and $\phi^{(0)}_{t_0}\bbeta^{(d)}=\E[\beta_{X_d}|\bbeta_T]$.
%, meaning that  ``$\phi^{(p)}\bbeta^{(d)}=\E[\bbeta^{(d)}|\mathcal{F}^{(p)}]$'' where $\mathcal{F}^{(p)}$ is the $\sigma$-algebra generated by the $p$-th level.
\begin{ex}[Averaging operators in the symmetric case]
Let $\bbeta_T=\bgamma_T$ be given by Model \ref{model:Sk}. Then 
\begin{align*}
\phi^{(p)}_r\bbeta^{(d)}
=&
\frac{1}{\prod_{\ell=p+1}^d I_\ell}
\left(\sum_{t\in T_d\,:\,t\succeq r}
\beta_{t}\right)
&r\in T_p\,.
\end{align*}
\end{ex}

Given $\phi$, we define the residual operators
$
\delta^{(p)}:\R^{T_d}
\rightarrow
%\times_{t\in T_{p-1}}(\R^{n_t})^*\subseteq
 \R^{T_p}
$
as
$\delta^{(p)}=(\delta^{(p)}_r)_{r\in T_p}$
with $\delta^{(p)}_r:\R^{T_d}
\rightarrow\R$ defined as
\begin{align}\label{eq:delta_defi}
\delta^{(p)}_r\bbeta^{(d)}=&
\phi^{(p)}_r\bbeta^{(d)}-\phi^{(p-1)}_{pa(r)}\bbeta^{(d)}
&r\in T_p
\end{align}
for $1\leq p\leq d\leq k-1$
and
$
\delta^{(0)}\bbeta^{(d)}=\phi^{(0)}\bbeta^{(d)}
$
for $0=p\leq d\leq k-1$.
Analogously to the 3-level case of Section \ref{sec:multigrid}, under suitable assumptions the residual operators $\delta^{(p)}$ decompose the Markov chain generated by the Gibbs Sampler into independent sub-chains.
To prove the result we first need the following lemma.
%%%%
\begin{lemma}[$p$-residuals interact only with $p$-residuals]\label{lemma:delta_expansion}
Let $\bbeta_T$ be a Gaussian random vector satisfying \eqref{eq:hierarchical_structure} and \eqref{eq:k_symmetry_rescaled}.
Then for any $p$ and $d$ with $0\leq p\leq d\leq k-1$, for all $t\in T_p$ we have
the identity
\begin{align*}
\E[\delta^{(p)}_r\bbeta^{(d)}|\bbeta\backslash\bbeta^{(d)}]
-\E[\delta^{(p)}_r\bbeta^{(d)}]
&=
\sum_{\ell\in\{p,\dots,k-1\}\backslash d}
c_{d\ell}%\frac{c_{d\ell}}{c_{dd}}
\left(\delta^{(p)}_r\bbeta^{(\ell)}-\E[\delta^{(p)}_r\bbeta^{(\ell)}]\right)
\,.
\end{align*}
\end{lemma}
%%%%
%\begin{proof}
%See Appendix \ref{proof:delta_expansion}
%\end{proof}
%%%%
Given Lemma \ref{lemma:delta_expansion} we can prove the following multigrid decomposition for hierarchical linear models.
%%%%
\begin{theorem}[Multigrid decomposition for $k$ levels]\label{thm:factorization_general}
Let $(\bbeta(s))_{s\in\N}$ be a Markov chain evolving according to \ref{sampler:GSk} with $\bbeta_T$ satisfying \eqref{eq:hierarchical_structure} and \eqref{eq:k_symmetry_rescaled}.
Then $(\delta^{(0)}\bbeta(s))_s$, \dots, $(\delta^{(k-1)}\bbeta(s))_s$ are $k$ independent Markov chains. 
Moreover, each chain $\delta^{(p)}\bbeta(s)=(\delta^{(p)}\bbeta^{(d)}(s))_{d=p}^{k-1}$ evolves according to the following blocked Gibbs sampler scheme with $(k-p)$ blocks:
for $d$ going from $p$ to $k-1$ sample
\begin{align}\label{eq:blocked_gibbs}
\delta^{(p)}\bbeta^{(d)}(s+1)
\;\sim\;
\mathcal{L}\left(\delta^{(p)}\bbeta^{(d)}|
(\delta^{(p)}\bbeta^{(\ell)}(s+1))_{p\leq \ell<d},(\delta^{(p)}\bbeta^{(\ell)}(s))_{d<\ell\leq k-1}
\right)\,,
\end{align}
where $\mathcal{L}(X|Y)$ denotes the conditional distribution of $X$ given $Y$.
\end{theorem}
%%%%
\newlength{\colA}
\setlength{\colA}{0cm}
\newlength{\colB}
\setlength{\colB}{1.55cm}
\newlength{\colE}
\setlength{\colE}{4.00cm}%4cm
\newlength{\colF}
\setlength{\colF}{6.1cm}%
\newlength{\colC}
\setlength{\colC}{8.55cm}%8.2cm
\newlength{\colD}
\setlength{\colD}{11cm}%10.5cm
\tikzstyle{state}=[shape=ellipse,draw=Apricot!50,fill=Apricot!20]
\tikzstyle{empty}=[]%shape=ellipse,draw=Apricot!50,fill=Apricot!20
\tikzstyle{beta}=[state]
%\tikzstyle{delta0}=[state,draw=blue!50,fill=blue!20]
%\tikzstyle{delta1}=[state,draw=red!50,fill=red!20]
%\tikzstyle{delta2}=[state,draw=green!50,fill=green!20]
\begin{figure}[h!]
\begin{center}
\begin{tikzpicture}[]
% 1 column
\node               at (\colA,6+2.5) {$s$};
\node[beta,inner sep=0.2pt] (s1_1) at (\colA,7.5) {$\bbeta^{(0)}$};
\node[empty] (s_1) at (\colA,6.75) {$\vdots$};
\node[beta,inner sep=0.2pt] (s2_1) at (\colA,5.85) {$\bbeta^{(1)}$};
\node[empty,inner sep=0.2pt] (s3_1) at (\colA,5.1) {$\vdots$};
\node[beta,inner sep=0.2pt] (s4_1) at (\colA,4.2) {$\bbeta^{(k-1)}$};
%2 column
\node               at (\colB,6+2.5) {$s+1$};
\node[beta,inner sep=0.2pt] (s1_2) at (\colB,7.5) {$\bbeta^{(0)}$}
    edge[mainedge] (s2_1)
    edge[mainedge] (s4_1)
;
\node[empty] (s_2) at (\colB,6.75) {$\vdots$}
;
\node[beta,inner sep=0.2pt] (s2_2) at (\colB,5.85) {$\bbeta^{(1)}$}
    edge[mainedge,bend right=45] (s1_2)
    edge[mainedge] (s4_1)
;
\node[empty,inner sep=0.2pt] (s3_2) at (\colB,5.1) {$\vdots$}
;
\node[beta,inner sep=0.2pt] (s4_2) at (\colB,4.2) {$\bbeta^{(k-1)}$}
    edge[mainedge,bend right=45] (s1_2)
    edge[mainedge,bend right=45] (s2_2)
%    edge[mainedge] (s3_2)
;
% 5 column
\node               at (\colE,8.5) {$s$};
\node[delta0,inner sep=0.2pt] (s1_3) at (\colE,7.5) {$\delta^{(0)}\bbeta$};
\node[empty] (s_3) at (\colE,6.75) {$\vdots$};
\node[delta1,inner sep=0.2pt] (s2_3) at (\colE,5.85) {$\delta^{(p)}\bbeta$};
\node[empty] (s3_3) at (\colE,5.1) {$\vdots$};
\node[delta2,inner sep=0.2pt] (s4_3) at (\colE,4.2) {$\delta^{(k-1)}\bbeta$};
%6 column
\node               at (\colF,8.5) {$s+1$};
\node[delta0,inner sep=0.2pt] (s1_4) at (\colF,7.5) {$\delta^{(0)}\bbeta$}
    edge[mainedge] (s1_3)
;
\node[empty] (s_4) at (\colF,6.75) {$\vdots$}
;
\node[delta1,inner sep=0.2pt] (s2_4) at (\colF,5.85) {$\delta^{(p)}\bbeta$}
    edge[mainedge] (s2_3)
;
\node[empty] (s3_4) at (\colF,5.1) {$\vdots$}
;
\node[delta2,inner sep=0.2pt] (s4_4) at (\colF,4.2) {$\delta^{(k-1)}\bbeta$}
    edge[mainedge] (s4_3)
;
% column C
\node               at (\colC,8.5) {$s$};
\node[delta1,inner sep=0.2pt] (s1_3) at (\colC,7.5) {$\delta^{(p)}\bbeta^{(p)}$};
\node[empty] (s_3) at (\colC,6.75) {$\vdots$};
\node[delta1,inner sep=0.2pt] (s2_3) at (\colC,5.85) {$\delta^{(p)}\bbeta^{(d)}$};
\node[empty,inner sep=0.2pt] (s3_3) at (\colC,5.1) {$\vdots$};
\node[delta1,inner sep=0.2pt] (s4_3) at (\colC,4.2) {$\delta^{(p)}\bbeta^{(k-1)}$};
% column D
\node               at (\colD,8.5) {$s+1$};
\node[delta1,inner sep=0.2pt] (s1_4) at (\colD,7.5) {$\delta^{(p)}\bbeta^{(p)}$}
    edge[mainedge] (s2_3)
    edge[mainedge,bend right=0] (s4_3)
;
\node[empty] (s_3) at (\colD,6.75) {$\vdots$}
;
\node[delta1,inner sep=0.2pt] (s2_4) at (\colD,5.85) {$\delta^{(p)}\bbeta^{(d)}$}
    edge[mainedge,bend right=55] (s1_4)
    edge[mainedge] (s4_3)
;
\node[empty,inner sep=0.2pt] (s3_4) at (\colD,5.1) {$\vdots$}
;
\node[delta1,inner sep=0.2pt] (s4_4) at (\colD,4.2) {$\delta^{(p)}\bbeta^{(k-1)}$}
    edge[mainedge,bend right=55] (s1_4)
    edge[mainedge,bend right=55] (s2_4)
;
\end{tikzpicture}
\end{center}
\caption{
Illustration of Theorem \ref{thm:factorization_general}.
Left: transition from $\bbeta(s)$ to $\bbeta(s+1)$ in Sampler \ref{sampler:GSk}. This is the common structure of a Gibbs Sampler with $k$ blocks.
Center: transition from $\delta\bbeta(s)$ to $\delta\bbeta(s+1)$ in Sampler \ref{sampler:GSk}. 
There are $k$ Markov chains, namely $(\delta^{(p)}\bbeta(s))_{s\in\N}$ for $p\in\{0,\dots,k-1\}$, evolving independently. 
Right: for each $p\in\{0,\dots,k-1\}$ the transition from $\delta^{(p)}\bbeta(s)=(\delta^{(p)}\bbeta^{(d)}(s))_{d=p}^{k-1}$ to $\delta^{(p)}\bbeta(s+1)=(\delta^{(p)}\bbeta^{(d)}(s+1))_{d=p}^{k-1}$ follows a Gibbs Sampler scheme with $(k-p)$ blocks.
}
\end{figure}
%\fotnote{Gareth: do we want to comment on the fact that single-site and blocked is equivalent for $\bbeta$ but not for $\delta\bbeta$?}
Theorem \ref{thm:factorization_general} implies that running a Gibbs sampler $(\bbeta(s))_s$ targeting distributions satisfying \eqref{eq:hierarchical_structure} is equivalent to running $k$ independent blocked Gibbs Samplers, one for each level of coarseness from $(\delta^{(0)}\bbeta(s))_s$ to $(\delta^{(k-1)}\bbeta(s))_s$.
%%%%%
\begin{coroll}\label{coroll:rate_k_levels}
Let $\bbeta_T$ satisfies \eqref{eq:hierarchical_structure} and \eqref{eq:k_symmetry_rescaled}.
Then the rate of convergence of \ref{sampler:GSk} is given by $\max\{\rho_0,\dots,\rho_{k-1}\}$ where for each $p\in\{0,\dots,{k-1}\}$, $\rho_p$ is the rate of convergence of $(\delta^{(p)}\bbeta(s))_s$.
\end{coroll}


%%%%%
%%%%%
\subsection{Hierarchical ordering of rates}\label{sec:ordering_k}
The convergence properties of these $k$ Markov chains are closely related, as the following theorem suggests.
\begin{theorem}\label{thm:kronecker}
The rate of convergence of $(\delta^{(p)}\bbeta(s))_s$ is given by the largest modulus eigenvalue of $(\I_{k-p}-L)^{-1}U$. Here $\I_{k-p}$ is the $(k-p)$ dimensional identity matrix, while $L$ and $U$ are, respectively, the strictly lower and strictly upper triangular part of $(c_{d\ell})_{d,\ell=p}^{k-1}$, with $c_{d\ell}$ given by \eqref{eq:k_symmetry_rescaled}.
%\fotnote{Giacomo: want to mention Kroenecker products here?}
\end{theorem}

Theorem \ref{thm:kronecker} implies that, from the rates of convergence point of view, the $k$ Markov chains updating $\delta^{(p)}\bbeta$ for $p=0,\dots,k-1$ behave as Gibbs samplers targeting a decreasing number of dimensions (from $k$ down to 1) of a single $k$-dimensional Gaussian distribution with precision matrix given by $-C$, where $C=(c_{d\ell})_{d,\ell=p}^{k-1}$ is given by \eqref{eq:k_symmetry_rescaled}.
This suggests that the convergence properties of the sub-chains will typically improve from that of $(\delta^{(0)}\bbeta(s))_s$ to that $(\delta^{(k-1)}\bbeta(s))_s$ and that the rate of convergence of $(\delta^{(0)}\bbeta(s))_s$ will typically determine the rate of the whole sampler \ref{sampler:GSk}.
In particular, in the centred parametrization case we can show that the rate of convergence is monotonically non-increasing from $(\delta^{(0)}\bbeta(s))_s$ to $(\delta^{(k-1)}\bbeta(s))_s$.
\begin{theorem}\label{thm:ordering_CP}\emph{(Ordering of rates for centred parametrization)}
Let $\bgamma$ be a Gaussian vector satisfying \eqref{eq:tree_cond_indep} and \eqref{eq:k_symmetry_rescaled} and 
let $(\bgamma(s))_{s\in\N}$ be the corresponding Markov chain evolving according to \ref{sampler:GSkCP}.
Then the convergence rates of the $k$ independent Markov chains $(\delta^{(0)}\bgamma(s))_s$, \dots, $(\delta^{(k-1)}\bgamma(s))_s$ satisfy
\begin{equation}\label{eq:ordering_k}
\rho(\delta^{(0)}\bgamma(s))\geq
\rho(\delta^{(1)}\bgamma(s))\geq
\cdots\geq
%\rho(\delta^{(k-2)}\bgamma(s))\geq
\rho(\delta^{(k-1)}\bgamma(s))=0\,.
\end{equation}
\end{theorem}
In Theorem \ref{thm:ordering_CP} we needed the additional assumption \eqref{eq:tree_cond_indep} to prove \eqref{eq:ordering_k}.
The reason is that, while in most cases the convergence rates of DSGS targeting a $n$-th dimensional Gaussian distribution improves if one of the coordinates is conditioned to a fixed value and the DSGS has to sample only from the remaining $(n-1)$ coordinates, this is not true in general.
Example 2 of \cite{RobertsSahu1997} provides a counter-example (see also \citealp[page 319]{Whittaker1990}).
In  \cite{RobertsSahu1997}, this example was used a counter-example regarding blocking strategies, it also works in the present context.
%\fotnote{Gareth: you mentioned you may try to make this sentence more clear. I CLARIFIED A BIT
%Some additional remarks (probably not worth being mentioned): 1- the result does not hold for DSGS but it holds for RSGS; 2- we didn't show that the Q-matrix of \cite[Example 2]{RobertsSahu1997} can be obtained from a hierarchical model and so strictly speaking the latter is not a counter-example for $\rho(\bbeta(s))=\rho(\delta^{(0)}\bbeta(s))$. I AGREE TO LEAVE OUT}
Theorem \ref{thm:ordering_CP} implies the following corollary.
\begin{coroll}\label{coroll:rate_k_levels_CP}
The rate of convergence of a DSGS targeting a Gaussian distribution satisfying \eqref{eq:tree_cond_indep} and \eqref{eq:k_symmetry_rescaled} is given by the largest squared eigenvalue of the $k$-dimensional matrix $C-\I_{k}$, where $C=(c_{d\ell})_{d,\ell=0}^{k-1}$ is defined by \eqref{eq:k_symmetry_rescaled} and $\I_{k}$ is the $k$-dimensional matrix.
\end{coroll}
In particular, considering the special case of Model \ref{model:Sk} it is easy to deduce the following Corollary.
\begin{coroll}\label{coroll:rate_k_levels_CP_symm}
The rate of convergence of $GS(\bgamma_T)$ targeting Model \ref{model:Sk} is given by the largest squared eigenvalue of the $k$-dimensional matrix
\begin{equation*}
\begin{pmatrix}
0 & r_{1}\\
(1-r_{2}) & 0 & r_{2}\\
 & \dots & \dots & \dots\\
 &  & (1-r_{k-2}) & 0 & r_{k-2}\\
 &  &  & (1-r_{k-1}) & 0
\end{pmatrix}
\end{equation*}
where 
$
r_\ell=
\frac{I_\ell\tau_{\ell}}{\tau_{\ell-1}+I_\ell\tau_{\ell}}
$
with $(\tau_1,\dots,\tau_{k-1})$ and $(I_1,\dots,I_{k-1})$ given by Model \ref{model:Sk}, $\tau_0=0$, $\tau_{k}=\tau_e$ and $I_{k}=J$.
%\begin{equation*}
%\begin{pmatrix}
%0 & c_1\\
%c_1 & 0 & c_2\\
% & \dots & \dots & \dots\\
% &  & c_{k-1} & 0 & c_{k}\\
% &  &  &c_{k} & 0
%\end{pmatrix}
%\end{equation*}
%where 
%$
%c_\ell=\sqrt{
%\frac{\tau_{\ell}}{\tau_{\ell}+I_{\ell+1}\tau_{\ell+1}}
%\frac{I_\ell\tau_{\ell}}{\tau_{\ell-1}+I_\ell\tau_{\ell}}
%}
%$
%with $(\tau_1,\dots,\tau_k)$ and $(I_1,\dots,I_k)$ given by Model \ref{model:Sk}, $\tau_0=0$ and $\tau_{k+1}=\tau_e$.
\end{coroll}

%%%%%
%%%%%
\subsection{Example: rates of convergence for 4-level models}\label{appendix:4_levels}
The results developed in Sections \ref{sec:multigrid_k} and \ref{sec:ordering_k} allow to analyze hierarchical models with an arbitrary number of levels.
For example we could consider 4-level extensions of Model \ref{model:S3}.
\begin{taggedmodel}{S4}\label{model:S4}(Symmetric 4-levels hierarchical model) 
Suppose
\begin{equation}\label{eq:NCP_4}
y_{ijk\ell}=\mu+a_i+b_{ij}+c_{ijk}+\epsilon_{ijk\ell},
\end{equation}
where $i$, $j$, $k$ and $\ell$ run from 1 to $I$, $J$, $K$ and $L$ respectively and $\epsilon_{ijk\ell}$ are iid normal random variables with mean 0 and variance $\var{e}$.
We employ a standard bayesian model specification assuming 
$a_i\sim N(0,\var{a})$, $b_{ij}\sim N(0,\var{b})$,  $c_{ijk}\sim N(0,\var{c})$ and a flat prior on $\mu$.
\end{taggedmodel}
In order to fit Model \ref{model:S4} with a Gibbs Sampler like \ref{sampler:GSk}, one could consider centering or non-centering each of the three levels $(a_i)_i$, $(b_{ij})_{ij}$ and $(c_{ijk})_{ijk}$.
Let $(\lambda_1,\lambda_2,\lambda_3)\in\{0,1\}^3$ be the non-centering indicators associated to the resulting in $8=2^3$ combinations. Here $\lambda_d=1$ indicates that the $d$-th level is non-centred while $\lambda_d=0$ indicates that it is centred.
The corresponding rates of convergence $\rho_{(\lambda_1,\lambda_2,\lambda_3)}$ can then be expressed in terms of the following normalized variance ratios
%Given two levels $i$ and $j$ we define the variance ratio $\ratio{i}{j}$ as
\begin{align*}
\ratio{i}{j}
&=
\frac{\tilde{\sigma}_i^2}{\tilde{\sigma}_i^2+\tilde{\sigma}_j^2}
&i,j\in\{1,2,3,4\}
\,,
\end{align*}
where $\tvar{1}=\frac{\var{a}}{I}$, $\tvar{2}=\frac{\var{b}}{IJ}$, $\tvar{3}=\frac{\var{c}}{IJK}$ and $\tvar{4}=\frac{\var{e}}{IJKL}$.
%where $\tvar{a}=\frac{\var{a}}{I}$, $\tvar{b}=\frac{\var{b}}{IJ}$, $\tvar{c}=\frac{\var{c}}{IJK}$ and $\tvar{e}=\frac{\var{e}}{IJKL}$.
If $\lambda_1=1$
%Using a non-centred parametrization at the highest level
 (i.e. using the non-centred parametrization $(a_i)_i$ at level $1$) the rates are 
\begin{align*}
%%%%%
\rho_{111}
&=
\max\{\ratio{1}{4},\ratio{2}{4},\ratio{3}{4}\}
\\
%%%%%
\rho_{110}
&=
\max\{\ratio{1}{3},\ratio{2}{3},\ratio{4}{3}\}
\\
%%%%%
\rho_{100}
&=
\max\{\ratio{1}{2},1-\ratio{2}{3}\ratio{3}{4}\}
\\
%%%%%
\rho_{101}
&=
\max\{\ratio{1}{2},1-\ratio{2}{4}\ratio{4}{3}\}
\end{align*}
If instead $\lambda_1=0$ (i.e. using $(\gamma_i)_i$ with $\gamma_i=\mu+a_i$ at level 1) we have four different parametrization with rates of convergence given by
\begin{align*}
%%%%%
\rho_{000}
&=
\frac{1}{2}\left(
1
+\ratio{2}{3}(\ratio{4}{3}-\ratio{1}{2})
+\sqrt{
(1+\ratio{2}{3}(\ratio{4}{3}-\ratio{1}{2}))^2
-4 \ratio{2}{1}\ratio{2}{3}\ratio{4}{3}}
\right)\,,
\\
%%%%%
\rho_{001}
&=
\frac{1}{2}\left(
1+\ratio{2}{4}(\ratio{3}{4}-\ratio{1}{2})
+\sqrt{
(1+\ratio{2}{4}(\ratio{3}{4}-\ratio{1}{2}))^2
-4 \ratio{2}{1}\ratio{2}{4}\ratio{3}{4}}
\right)\,,
\\
%%%%%
\rho_{011}
&=
\frac{1}{2}\left(
1-\ratio{1}{4}\ratio{4}{2}\ratio{4}{3}+
\ratio{2}{4}\ratio{3}{4}
+\sqrt{
(1-\ratio{1}{4}\ratio{4}{2}\ratio{4}{3}+
\ratio{2}{4}\ratio{3}{4})^2
-4 \ratio{2}{4}\ratio{3}{4}}
\right)\,,
\\
%%%%%
\rho_{010}
&=
\frac{1}{2}\left(
1-\ratio{1}{3}\ratio{3}{2}\ratio{3}{4}+
\ratio{2}{3}\ratio{4}{3}
+\sqrt{
(1-\ratio{1}{3}\ratio{3}{2}\ratio{3}{4}+
\ratio{2}{3}\ratio{4}{3})^2
-4 \ratio{2}{3}\ratio{4}{3}}
\right)\,.
%\\
%%%%%%
%2\rho_{010}
%&=
%1-\ratio{1}{3}\ratio{3}{2}\ratio{3}{4}+
%\ratio{2}{3}\ratio{4}{3}
%+\sqrt{
%(1-\ratio{1}{3}\ratio{3}{2}\ratio{3}{4}+
%\ratio{2}{3}\ratio{4}{3})^2
%-4 \ratio{2}{3}\ratio{4}{3}}
%\,.
\end{align*}
These rates can be easily derived from Corollary \ref{coroll:rate_k_levels} and Theorem \ref{thm:kronecker}.
Even if some of the rates obtained have more complicated expressions than in the two and three-levels case, especially when the highest level is centred, the theory developed in this section can be easily applied to models with an arbitrary number of levels.
It is worth noting that in this 4-level case the skeleton chain $\delta^{(0)}\bbeta$ is always the slowest chain for all mixed parametrizations (can be checked by inspection), even if we provided a general proof of this fact only for the centred parametrization.
The expressions given here can be easily used to derive conditionally optimal parametrizations for Model \ref{model:S4} given the rescaled variance components $(\tvar{i})_{i=1}^4$. For example, choosing whether to center or not each level by comparing the level-specific rescaled variances with the sum of the rescaled variances of the lower levels like in Section \ref{sec:cond_optimal} leads to rates of convergence upper bounded by $\frac{3}{4}$.
Note that, while this strategy is optimal for 3-level models (see Section \ref{sec:cond_optimal} and Corollary \ref{coroll:optimal_param_3levels}), the same strategy is not the optimal one for 4 levels. In fact the supremum over $(\tvar{i})_{i=1}^4$ of $\min_{(\lambda_1,\lambda_2,\lambda_3)\in\{0,1\}^3}\rho_{(\lambda_1,\lambda_2,\lambda_3)}$ is strictly smaller than $\frac{3}{4}$, but the improvement is small and does not seem sufficient to motivate a more complicated strategy.
%\fotnote{Gareth: I added a part to say that in 4 levels this strategy is not the optimal one. You may want to double-check it, iit's probably not great! OK!}

\subsection{Optimal partial non-centering}
In Section \ref{sec:PNCP_k_levels}, we defined the class of Partially Non-Centered Parametrizations (PNCP) $\bbeta_T$ of Model \ref{model:NSk}.
As for the 3-level case (see Section \ref{sec:PNCP_3}), even in this context one can find a optimal PNCP that factorizes the unknown parameters $\bbeta_T$ into independent random variables.
%\begin{theorem}\label{thm:PNCP_k_levels}
%Let $\bgamma_T$ be given by Model \ref{model:NSk} (conditioned on some value of $\textbf{y}$). Denote by $Q=(Q_{tr})_{t,r\in T}$ its precision matrix and define $\Lambda=(\Lambda_{tr})_{r\preceq t,t\in T}$ as 
%\begin{align*}
%\Lambda_{tt}&=1&t\in T,\\
%\Lambda_{t\,pa(t)}&=\frac{Q_{t\,pa(t)}}{D_{tt}}&t\in T\backslash \{t_0\},\\
%\Lambda_{tr}&=0&t\in t\in T,\;r\notin\{t,pa(t)\},
%\end{align*}
%where
%\begin{align*}
%D_{\ell\ell}&=1
%&\ell\in T_L
%\\
%D_{tt}&=Q_{tt}-\sum_{t'\in ch(t)}\Lambda_{t't}^2D_{t't'}
%&t\in T\backslash \{t_0\}\,.
%\end{align*}
%Then the resulting hierarchical reparametrization $\bbeta_T=\Lambda\bgamma_T$ is a collection of independent random variables.
%\end{theorem}
\begin{theorem}\label{thm:PNCP_k_levels}
Let $\bgamma_T$ be given by Model \ref{model:NSk} (conditioned on some value of $\textbf{y}$). Denote by $Q=(Q_{tr})_{t,r\in T}$ its precision matrix and define $\Lambda=(\lambda_{tr})_{r\preceq t,t\in T}$ as 
\begin{align*}
\lambda_{tt}&=1&t\in T,\\
\lambda_{t\,pa(t)}&=\frac{Q_{t\,pa(t)}}{D_{tt}}&t\in T\backslash \{t_0\},\\
\lambda_{tr}&=0&t\in T,\;r\notin\{t,pa(t)\},
\end{align*}
where
\begin{align}
D_{tt}&=Q_{tt}-\sum_{t'\in ch(t)}\lambda_{t't}^2D_{t't'}
&t\in T\backslash \{t_0\}\,.\label{eq:D_definition}
\end{align}
Then the resulting hierarchical reparametrization $\bbeta_T=\Lambda\bgamma_T$ is a collection of independent random variables.
\end{theorem}
Note that the re-parametrization coefficients $\Lambda=(\lambda_{tr})_{r\preceq t,t\in T}$ defined in can be easily calculated iteratively starting from the leaves going up to the root of $T$. 
In fact the optimal PNCP is related to the multiscale or change-of-resolution Kalman filter (see e.g.\ \citealp{Chou1994}) and more generally to belief propagation algorithms for trees. 

\section{Conclusions and future work}
In this work we studied the convergence properties of Gibbs Sampler algorithms in the context of Gaussian hierarchical models.
To do so we developed a novel analytic approach based on multigrid decompositions that allows to factorize the Markov chain of interest into independent and easier to analyze sub-chains.
This decomposition enables us to evaluate explicitly the $L^2$-rate of convergence in symmetric (e.g. Model \ref{model:Sk}) and weakly-symmetric cases (e.g. Examples \ref{ex:weakly_symm} and \ref{ex:non_symm}). These results can then be used to provide upper bounds in the more general non-symmetric case (e.g. Theorem \ref{thm:rates_3_uneven_CP}).
%The general results are proven for models with an arbitrary number of levels and the 3-level case is considered in more depth to provide practically useful indications (e.g. Sections \ref{sec:explicit_rates}, \ref{sec:cond_optimal}, \ref{sec:bounds} and \ref{sec:bespoke}).


We believe that the results presented in this work can provide a novel and significant contribution to understanding the behavior of MCMC algorithms (and in particular the Gibbs Sampler) in the extremely popular context of hierarchical linear models.
%From a more theoretical point of view, we believe that the results presented in this work provide a significant and novel contribution to understanding the behavior of MCMC algorithms (and in particular the Gibbs Sampler) for hierarchical linear models.
%e believe that the use of a multigrid decomposition to analyze the Gibbs Sampler on hierarchical models offers an original and valuable viewpoint on 
%At the same time, our resu
In fact, together with explicit formulas for $L^2$-rates of convergences, the multigrid decomposition we presented here provides a simple and intuitive theoretical characterizations of practical behaviors commonly observed in practice when fitting hierarchical models with MCMC, such as slower mixing for hyper-parameters at higher levels (see Theorems \ref{thm:rates_ordering_3} and \ref{thm:ordering_CP}), algorithmic scalability with width of the hierarchy but not with height (e.g. Theorem \ref{thm:rates_3} and Corollary \ref{coroll:rate_k_levels_CP_symm}) and good performances of centred parametrization in data-rich contexts (Theorem \ref{thm:rates_3_uneven_CP}).
Also, it is worth noting that these results are not limited to hierarchical models with two levels but rather can be applied to models with an arbitrary number of layers.

The results presented here have also various practically useful implications in terms of algorithmic implementations of MCMC for hierarchical linear models. 
For example, our analysis provides easy and theoretically grounded indications to choose the computationally optimal parametrization on-the-fly (see Sections \ref{sec:explicit_rates}, \ref{sec:cond_optimal}), how to handle heterogeneity (see Section \ref{sec:bespoke}) and indicates which parameters to monitor in the convergence diagnostic process (see Theorem \ref{thm:rates_ordering_3} and discussion at the end of Section \ref{sec:motivating_example}).
Also, the bounds in Theorem \ref{thm:rates_3_uneven_CP} can be used to have a theoretical guarantee of good convergence properties with large datasets.

The present work could be extended in many directions. 
First one could extend the results for non-symmetric cases by improving and generalizing the bounds of Theorem \ref{thm:rates_3_uneven_CP}, by providing upper bounds on the rates of convergence for the optimal bespoke parametrization (Section \ref{sec:bespoke}) and ideally by weakening the symmetry assumption in \eqref{eq:k_symmetry}.
Also it would be worth exploring more deeply the use of partially-non centred parametrization and their computational feasibility.
Finally an ambitious aim, involving non-trivial generalizations of the approach of \cite{.amit:1996:properties+convergence+perturbations}, would be to extend to other tractable distributions within the exponential family beyond the Gaussian.
%\fotnote{Gareth: added a short sentence n non-Gaussian case, maybe it's too vague, what do you think? HAVE EXTENDED SLIGHTLY}

In terms of classes of models considered, a natural and important extension would be to consider the multivariate case (where each parameter $\gamma_t$ is a multivariate random vector) and the regression case. 
This would make the results presented here applicable to a large class of widely used models.
We expect many results developed in this work to extend to the multivariate and regression case, even if in that context the role played by non-symmetric cases will be more crucial.
%\fotnote{Mention that regression cases are somehow less symmetric than the univariate case considered here?}
Another important class of models that would be worth considering and approaching with methodologies analogous to the ones developed here are models based on Gaussian processes commonly used, for example, in spatial statistics (see e.g.\ \citealp{BassSahu2016}).
Also, it would be interesting to compare nested structures like the one considered here to crossed-effect structures (see e.g. \citealp{GaoOwen2016}).


In terms of algorithms considered it would be relatively easy to extend the machinery developed in this work to other Gibbs Sampling schemes, such as the random scan Gibbs Sampling case, where the component to be updated is chosen uniformly at random at each step.
In this case the resulting sub-chains would not evolve completely independently because they would be coupled by the random choice of which component to update. 
Nevertheless each sub-chain would still be a Markov chain marginally and the overall rate of the vector of coupled Markov chains would still be amenable to direct analysis.
Also, in the random scan case the reversibility of the induced Markov chains would allow to prove orderings results under weaker assumptions than Theorem \ref{thm:ordering_CP}. 

While this work is focused on $L^2$-rates of convergence, the same approach could be used to derive bounds on the distance (e.g.\ total variation or Wasserstein) between the distribution of the Markov chain at a given iteration and the target distribution.
%Finally, it would be interesting to use the results of this work to derive bounds on the distance between the distribution of the Markov chain at a given iteration and the target distribution.
%Although this is not the point of view we took in this work, it is often the case that statements on $L^2$-rates of convergence can be turned into similar statements on the total variation or the Wasserstein distance from the target (see for example \cite{RajaratnamSparks2015}).
Such a formulation would be interesting to extend the recent growth in literature on providing rigorous characterizations of the computational complexity of Bayesian hierarchical linear models, see for example  \cite{RajaratnamSparks2015,roberts2016complexity,johndrow2015approximations}.
%\fotnote{Gareth: do you think it's worth citing the paper I sent you by \cite{RajaratnamSparks2015}? Could also cite the paper by you and jeff on complexity bounds via diffusion limits, although the context is quite different. I HAVE ADDED THIS AND A COUPLE MORE ALSO}.
In order to provide full characterizations, however, the case of unknown variances should be considered (see e.g.\ \cite{jones2004sufficient} for the two level case).
%\fotnote{Gareth: I wrote this final paragraph on distributional distance and computational complexity. What do you think? Feel free to change it or remove it.}

\subsection*{Acknowledgments}
The authors are grateful for stimulating discussions with Omiros Papaspiliopoulos
and Art Owen.
GZ supported in part% by EPSRC Research Grant EP/D002060, 
by an EPSRC Doctoral Prize fellowship and by the European Research Council (ERC) through StG ``N-BNP'' 306406. GOR acknowledges support from
EPSRC through grants EP/K014463/1 (i-Like) and EP/D002060/1 (CRiSM).

\bibliography{GibbsSampler_bibliography}

%\newpage
\appendix




%%%%%
%%%%%
\section{Full conditional distributions of \ref{sampler:GSNN} and \ref{sampler:GSCC}}\label{appendix:full_conditionals}
\begin{taggedsampler}{GS($1,1$)}
Initialize $\mu(0)$, $\avec(0)$ and $\bvec(0)$ and then iterate
\begin{align*}
\mu(s+1)
&\sim N\left(\ymean-\amean(s)-\bmean(s),\frac{\var{e}}{n}\right)\\
a_i(s+1)
&\sim N\left(\frac{\var{a}(\yimean-\mu(s+1)-\bimean(s))}{\var{a}+\frac{\var{e}}{JK}},\frac{\var{a}\frac{\var{e}}{JK}}{\var{a}+\frac{\var{e}}{JK}}\right)
&\forall\, i
\\
b_{ij}(s+1)
&\sim N\left(\frac{\var{b}(\yijmean-\mu(s+1)-a_i(s+1))}{\var{b}+\frac{\var{e}}{K}},\frac{\var{b}\frac{\var{e}}{K}}{\var{b}+ \frac{\var{e}}{K} }\right)
&\forall\, i,j
\end{align*}
where we use the dot subscript to indicate averaging over one dimension, meaning that
$\amean=\sum_{i}\frac{a_i}{I}$,
$\bmean=\sum_{i,j}\frac{b_{ij}}{IJ}$,
$\ymean=\sum_{i,j,k}\frac{y_{ijk}}{IJK}$,
$\bimean=\sum_{j}\frac{b_{ij}}{J}$,
$\yimean=\sum_{j,k}\frac{y_{ijk}}{JK}$
and
$\yijmean=\sum_{k}\frac{y_{ijk}}{K}$.
\end{taggedsampler}
\begin{taggedsampler}{GS($0,0$)}
Initialize $\mu(0)$, $\gavec(0)$ and $\etavec(0))$ and then iterate
\begin{align*}
\mu(s+1)
&\sim N\left(\gamean(s),\frac{\var{a}}{I}\right)
&\\
\gamma_i(s+1)
&\sim N\left(\frac{\frac{\var{b}}{J}\mu(s+1)+\var{a}\etai(s)}{\frac{\var{b}}{J}+\var{a}},\frac{\var{a}\frac{\var{b}}{J}}{\var{a}+\frac{\var{b}}{J}}\right)
&\forall\, i\\
\eta_{ij}(s+1)
&\sim N\left(\frac{\frac{\var{e}}{K}\gamma_i(s+1)+\var{b}\yijmean}{\var{b}+\frac{\var{e}}{K}},\frac{\var{b}\frac{\var{e}}{K}}{\var{b}+\frac{\var{e}}{K}}\right)
&\forall\, i,j
\end{align*}
where as before the dot subscript indicates averaging over indices.
\end{taggedsampler}

%%%%%
%%%%%
\section{Theorem proofs}\label{appendix:proofs}
In this section we list proofs following the order of appearance in the paper. 
This is different from the mathematical chronology. 
To follow the latter, the reader should start from the results (lemmas, theorems and corollaries) for $k$-level models, namely the results from Lemma \ref{lemma:hierarchical_reparametrization_invertible} to Theorem \ref{thm:PNCP_k_levels} and then move to the ones for $3$-level models, namely the results from Theorem \ref{thm:factorization_3} to Theorem \ref{thm:PNCP_3_levels}.

%%%%%
%%%%%
\subsection*{Proof of Theorems \ref{thm:factorization_3}, \ref{thm:rates_ordering_3} and \ref{thm:rates_3}}\label{proof:factorization_3}
%%%%%
Theorems \ref{thm:factorization_3}, \ref{thm:rates_ordering_3} and \ref{thm:rates_3} are substantially special cases of the Theorems of Section \ref{sec:k_levels}.
In particular Theorem \ref{thm:factorization_3} is a special case of Theorem \ref{thm:factorization_general} for $k=2$ and Theorems \ref{thm:rates_ordering_3} and \ref{thm:rates_3} can be directly verified as follows.
Using Corollary \ref{coroll:rate_k_levels} we can evaluate the rates of convergence of the three subchains
$
\rho(\delta^{(0)}\bbeta(s))
$, 
$\rho(\delta^{(1)}\bbeta(s))$
and
$\rho(\delta^{(2)}\bbeta(s))$
for the four parametrzations under consideration in Section \ref{sec:multigrid} and check by inspection that 
$$
\rho(\delta^{(0)}\bbeta(s))
\rho(\delta^{(1)}\bbeta(s))\geq
\rho(\delta^{(2)}\bbeta(s))=0
$$
and that the rates of convergence $\rho(\delta^{(0)}\bbeta(s))$ are the ones given by Theorem \ref{thm:rates_3}.
In particular the rates of convergence of $\delta^{(0)}\bbeta(s))$, 
$\delta^{(1)}\bbeta(s)$
and
$\delta^{(2)}\bbeta(s)$
under $GS(\bbeta)$ are given by the following Table.
%\begin{figure}[h!]
%%\begin{center}
%\begin{tabular}{ l |c| c | c | c}
%  & $(\mu,\avec,\bvec)$ & $(\mu,\gavec,\etavec)$ & $(\mu,\gavec,\bvec)$ & $(\mu,\avec,\etavec)$ \\
%  \hline			
% $\rho(\delta^{(0)}\bbeta(s))$ & 
%$\frac{\tvar{a}}{\tvar{a}+\tvar{e}}\vee\frac{\tvar{b}}{\tvar{b}+\tvar{e}}$ & 
%$1-\frac{\tvar{a}}{\tvar{a}+\tvar{b}}\frac{\tvar{b}}{\tvar{b}+\tvar{e}}$ &
%$1-\frac{\tvar{a}}{\tvar{a}+\tvar{e}}\frac{\tvar{e}}{\tvar{b}+\tvar{e}}$ &
%$\frac{\tvar{a}}{\tvar{a}+\tvar{b}}\vee\frac{\tvar{e}}{\tvar{b}+\tvar{e}}$
%\\
%  \hline			
% $\rho(\delta^{(1)}\bbeta(s))$ & 
%$\frac{\tvar{a}}{\tvar{a}+\tvar{e}}\frac{\tvar{b}}{\tvar{b}+\tvar{e}}$ & 
%$\frac{\tvar{a}}{\tvar{a}+\tvar{b}}\left(1-\frac{\tvar{b}}{\tvar{b}+\tvar{e}}\right)$ &
%$\frac{\tvar{a}}{\tvar{a}+\tvar{e}}\left(1-\frac{\tvar{e}}{\tvar{b}+\tvar{e}}\right)$ &
%$\frac{\tvar{a}}{\tvar{a}+\tvar{b}}\frac{\tvar{e}}{\tvar{b}+\tvar{e}}$
%\\
%  \hline  
% $\rho(\delta^{(2)}\bbeta(s))$  & 0 & 0 & 0 & 0 \\
%  \hline  
%\end{tabular}\caption{Rates of convergence of $\delta^{(0)}\bbeta(s))$, 
%$\delta^{(1)}\bbeta(s)$
%and
%$\delta^{(2)}\bbeta(s)$ for $GS(\bbeta)$ under various parametrizations.}\label{table:rates_3levels}
%%\end{center}
%\end{figure}
\begin{figure}[h!]
%\begin{center}
\centering
\begin{tabular}{ l |c| c | c }
& $\rho(\delta^{(0)}\bbeta(s))$ &  $\rho(\delta^{(1)}\bbeta(s))$ &  $\rho(\delta^{(2)}\bbeta(s))$ 
\\
\hline			
$(\mu,\avec,\bvec)$ & 
$\frac{\tvar{a}}{\tvar{a}+\tvar{e}}\vee\frac{\tvar{b}}{\tvar{b}+\tvar{e}}$ & 
$\frac{\tvar{a}}{\tvar{a}+\tvar{e}}\frac{\tvar{b}}{\tvar{b}+\tvar{e}}$ &
0
\\
\hline			
 $(\mu,\gavec,\etavec)$ &
$1-\frac{\tvar{a}}{\tvar{a}+\tvar{b}}\frac{\tvar{b}}{\tvar{b}+\tvar{e}}$ &
$\frac{\tvar{a}}{\tvar{a}+\tvar{b}}\left(1-\frac{\tvar{b}}{\tvar{b}+\tvar{e}}\right)$ &
0
\\
\hline
$(\mu,\gavec,\bvec)$ & 
$1-\frac{\tvar{a}}{\tvar{a}+\tvar{e}}\frac{\tvar{e}}{\tvar{b}+\tvar{e}}$ &
$\frac{\tvar{a}}{\tvar{a}+\tvar{e}}\left(1-\frac{\tvar{e}}{\tvar{b}+\tvar{e}}\right)$ &
0
\\
\hline
$(\mu,\avec,\etavec)$ &
$\frac{\tvar{a}}{\tvar{a}+\tvar{b}}\vee\frac{\tvar{e}}{\tvar{b}+\tvar{e}}$ &
$\frac{\tvar{a}}{\tvar{a}+\tvar{b}}\frac{\tvar{e}}{\tvar{b}+\tvar{e}}$ &
0
\\
\hline
\end{tabular}\caption{Rates of convergence of $\delta^{(0)}\bbeta(s))$, 
$\delta^{(1)}\bbeta(s)$
and
$\delta^{(2)}\bbeta(s)$ for $GS(\bbeta)$ under various parametrizations.}\label{table:rates_3levels}
%\end{center}
\end{figure}
%%%%%%
%\begin{proof}[Proof or Theorem \ref{thm:rates_ordering_3}]
%[[At the moment we simply wrote down the rates of convergence for the four parametrizations and checked]]
%From Corollary \ref{coroll:rate_k_levels} it follows that
%\end{proof}
%
%%%%%%
%\begin{proof}[Proof or Theorem \ref{thm:rates_3}]
%[[
%To be written. Simple calculations, now written in\\ GS$\_$for$\_$three$\_$levels$\_$Gaussian$\_$hierarch$\_$model.pdf
%]]
%\end{proof}
%%%%%
\subsection*{Proof of Corollary \ref{thm:rate_eq_skeleton}}
By Theorem \ref{thm:factorization_3} the rate of convergence of the whole chain $(\bbeta(s))_{s\in \N}$ coincides with the maximum of the rates of the subchains, meaning that 
$\rho(\bbeta(s))=\max\{
\rho(\delta^{(0)}\bbeta(s)),
\rho(\delta^{(1)}\bbeta(s)),
\rho(\delta^{(2)}\bbeta(s))
\}$.
By Theorem \ref{thm:rates_ordering_3} the latter maximum equals $\rho(\delta^{(0)}\bbeta(s))$.

%%%%%
\subsection*{Proof of Corollary \ref{coroll:optimal_param_3levels}}
Follows from Theorem \ref{thm:rates_3} by checking that for both 
$\balpha=\gavec$ and $\balpha=\avec$, the inequality
$\rho_{(\mu,\balpha,\etavec)}\leq \rho_{(\mu,\balpha,\bvec)}$
holds if and only if 
$\tvar{b}\geq\tvar{e}$;
and for both $\balpha=\gavec$ and $\balpha=\avec$ the inequality
$\rho_{(\mu,\gavec,\balpha)}\leq \rho_{(\mu,\avec,\balpha)}$
holds if and only if 
$\tvar{a}\geq\tvar{b}+\tvar{e}$.

%%%%%
\subsection*{Proof of Theorem \ref{thm:rates_3_uneven_CP}}\label{appendix:proof_bound}
Given an instance of Model \ref{model:3_nonsymmetric} with variance terms $(\var{a},(\var{b,i})_i,(\var{e,ij})_{ij})$ satisfying \eqref{eq:assumption_for_bound}, the proof will proceed by comparing the original Gibbs Sampler with an auxiliary Gibbs Sampler targeting a different instance of Model \ref{model:3_nonsymmetric} with variance terms $(\var{a},(\var{b,i})_i,(\bar{\sigma}_{e,ij}^2)_{ij})$ satisfying \eqref{eq:k_symmetry_appendix} and thus allowing direct analysis using Corollary \ref{coroll:rate_k_levels_CP}.
%We will provide an upper bound on the convergence properties of the Gibbs Sampler targeting such model by comparison with a Gibbs Sampler targeting a different instance of Model \ref{model:3_nonsymmetric} which satisfies \eqref{eq:k_symmetry_appendix} and thus is amenable to direct analysis using Corollary \ref{coroll:rate_k_levels_CP}.
In the context of Model \ref{model:3_nonsymmetric}, \eqref{eq:k_symmetry_appendix} reduces to requiring $\sum_{j=1}^{J_i}\rho_{\gamma_i \eta_{ij}}^2$ to be constant over $i$, where $\rho_{\gamma_i \eta_{ij}}$ is the partial correlation $Corr(\gamma_i,\eta_{ij}|\mu,(\gamma_\ell)_{\ell\neq i},(\eta_{\ell s})_{(\ell s)\neq (ij)})$ as in Section \ref{sec:symmetry}.
By computing the partial correlations of Model \ref{model:3_nonsymmetric} it can be checked that $\sum_{j=1}^{J_i}\rho_{\gamma_i \eta_{ij}}^2=r^{(i)}_{a,b}r^{(i)}_{e,b}$, where $r^{(i)}_{a,b}$ and $r^{(i)}_{e,b}$ are defined in Theorem \ref{thm:rates_3_uneven_CP}.
For each $i=1,\dots,I$ we define auxiliary variance terms $(\bar{\sigma}_{e,ij}^2)_{j=1}^{J_i}$ such that $\bar{\sigma}_{e,ij}^2\geq\var{e,ij}$ for all $j=1,\dots,J_{i}$ and 
\begin{equation}\label{eq:proof_bounds}
r^{(i)}_{a,b}\,%\frac{J_i\precision{b,i}}{\precision{a}+J_i\precision{b,i}}
\frac{1}{J_i}\sum_{j=1}^{J_i}
\frac{K_{ij}^{-1}\bar{\sigma}_{e,ij}^2}{\var{b,i}+K_{ij}^{-1}\bar{\sigma}_{e,ij}^2}
=
\max_{\ell=1,\dots,I}
r^{(\ell)}_{a,b}r^{(\ell)}_{e,b}\,.%\frac{J_i\precision{b,i}}{\precision{a}+J_i\precision{b,i}}
%\frac{\precision{b,\ell}}{\precision{b,\ell}+K_{\ell j}\precision{e,\ell j}}\,.
\end{equation}
Such $(\bar{\sigma}_{e,ij}^2)_{j=1}^{J_i}$ exist because $r^{(i)}_{a,b}\geq \max_{\ell=1,\dots,I}r^{(\ell)}_{a,b}r^{(\ell)}_{e,b}$ by \eqref{eq:assumption_for_bound} and the left hand side of \eqref{eq:proof_bounds} can take any value in $(0,r^{(i)}_{a,b}]$ for $(\bar{\sigma}_{e,ij}^2)_{j=1}^{J_i}$ belonging to $[0,\infty)$.
\eqref{eq:proof_bounds} implies that the instance of Model \ref{model:3_nonsymmetric} with variance terms $(\var{a},(\var{b,i})_i,(\bar{\sigma}_{e,ij}^2)_{ij})$ satisfies \eqref{eq:k_symmetry_appendix} 
with $c_0=\sum_{i=1}^I\rho_{\mu\gamma_i}^2=1-\frac{1}{I}\sum_{i=1}^Ir^{(i)}_{a,b}$
and
$c_1=\sum_{j=1}^{J_i}\rho_{\gamma_i\eta_{ij}}^2=\max_{\ell=1,\dots,I}
r^{(\ell)}_{a,b}r^{(\ell)}_{e,b}$.
As discussed in Example \ref{ex:non_symm}, for models with centred parametrization like Model \ref{model:3_nonsymmetric}, \eqref{eq:k_symmetry_appendix} implies \eqref{eq:k_symmetry} and, after rescaling, \eqref{eq:k_symmetry_rescaled}.
In this case the matrix $C=(c_{dp})_{d,p=0}^2$ is given by
\begin{equation*}
C=
\begin{pmatrix}
1 & \sqrt{c_0} & 0\\
\sqrt{c_0} & 1 & \sqrt{c_1}\\
0 & \sqrt{c_1} & 1
\end{pmatrix}\,.
\end{equation*}
Therefore, by Corollary \ref{coroll:rate_k_levels_CP}, the rate of convergence of the Gibbs Sampler targeting the posterior distribution of Model \ref{model:3_nonsymmetric} with variance terms $(\var{a},(\var{b,i})_i,(\bar{\sigma}_{e,ij}^2)_{ij})$ is given by $c_0+c_1$, which is the largest squared eigenvalue of $C-\I_3$, where $\I_3$ is the 3-dimensional identity matrix.

Finally we show that the Gibbs Sampler rate of convergence induced by the auxiliary variance terms $(\var{a},\var{b,i},\bar{\sigma}_{e,ij}^2)$ is greater or equal than the original one given by $(\var{a},\var{b,i},\var{e,ij})$.
Denote by $Q$ the precision matrix of the original posterior distribution and by $\bar{Q}$ the auxiliary one.
By deriving $Q$ and $\bar{Q}$ from the definition of Model \ref{model:3_nonsymmetric}, it is easy to see that the only terms of $Q$ and $\bar{Q}$ affected by replacing $\var{e,ij}$ with $\bar{\sigma}_{e,ij}^2$ are $(Q_{\eta_{ij}\eta_{ij}})_{ij}$ and  $(\bar{Q}_{\eta_{ij}\eta_{ij}})_{ij}$.
Moreover $\bar{\sigma}_{e,ij}^2\geq\var{e,ij}$ implies $Q_{\eta_{ij}\eta_{ij}}=\frac{1}{\var{b,i}}+\frac{K_{ij}}{\var{e,ij}}\geq \frac{1}{\var{b,i}}+\frac{K_{ij}}{\bar{\sigma}_{e,ij}^2}=\bar{Q}_{\eta_{ij}\eta_{ij}}$.
The result then follows from Theorem 7 of \cite{RobertsSahu1997}, which states that the convergence rate of a deterministic scan Gibbs Sampler with single-site update is a non-increasing function of the diagonal elements of the target precision matrix (when the off-diagonal terms are kept constant).
%
%two matrices are equal differ only in the terms $(Q_{\eta_{ij}\eta_{ij}})_{ij}$ and  $(\tilde{Q}_{\eta_{ij}\eta_{ij}})_{ij}$
%
%This follows from \cite[Thm.7]{RobertsSahu1997}, which states that the converges rate of the Gibbs Sampler is 
%
%Finally note that the only elements of the precision matrix under consideration affected by  replacing $\var{e,ij}$ with $\bar{\sigma}_{e,ij}^2$ are 
%
%Finally it is easy to check that the $A$-matrix $(A_{tr})_{rt}=(\frac{Q_{tr}}{Q_tt})$ given by $(\precision{a},\precision{b,i},\tprec{e,ij})$ dominates from above the one given by $(\precision{a},\precision{b,i},\precision{e,ij})$ because the $A$ matrix is monotone with respect to the terms $\precision{e,ij}$.
%Therefore the corresponding rates of convergence have the same domination and we have the thesis of Theorem \ref{thm:rates_3_uneven_CP}.
%\begin{equation*}
%\begin{pmatrix}
%1 & \sqrt{c_0} & 0\\
%\sqrt{c_0} & 1 & \sqrt{c_1}\\
%0 & \sqrt{c_1} & 1
%\end{pmatrix}
%\end{equation*}
%where in this case the vector $(c_0,c_1)$, defined in \eqref{eq:k_symmetry_appendix},
%is given by
%$c_0=\frac{1}{I}\sum_{i=1}^I\frac{J_i\precision{b,i}}{\precision{a}}$
%and
%$c_1=\max_{i=1,\dots,I}\frac{1}{J_i}\sum_{j=1}^{J_i}
%\frac{J_i\precision{b,i}}{\precision{a}+J_i\precision{b,i}}
%\frac{\precision{b,i}}{\precision{b,i}+K_{ij}\precision{e,ij}}$.
%It follows that the rate of the model with precisions $(\precision{a},\precision{b,i},\tprec{e,ij})$ is given by $c_0+c_1$.
%Finally it is easy to check that the $A$-matrix given by $(\precision{a},\precision{b,i},\tprec{e,ij})$ dominates from above the one given by $(\precision{a},\precision{b,i},\precision{e,ij})$ because the $A$ matrix is monotone with respect to the terms $\precision{e,ij}$.
%Therefore the corresponding rates of convergence have the same domination and we have the thesis of Theorem \ref{thm:rates_3_uneven_CP}.

%%%%%
\subsection*{Proof of Theorem \ref{thm:bespoke_parametrizations_2}}\label{proof:bespoke_parametrizations_2}
The Markov chain under consideration is a Gibbs Sampler sweeping through $(\mu,\beta_1,\dots,\beta_I)|\textbf{y}$ for some observed data $\textbf{y}=((y_{ij})_{j=1}^{J_i})_{i=1}^I$ assuming Model \ref{model:2_nonsymmetric} and $\beta_i=\gamma_i-\lambda_i\mu$ with $\lambda_i\in\{0,1\}$.
%We can assume, without loss of generality, that $y_{ij}=0$ for all $i$ and $j$ by replacing $\mu$ and $\beta_i$ with $\mu-\E[\mu|\textbf{y}]$ and $\beta_i-\E[\beta_i|\textbf{y}]$ respectively.
To compute the Gibbs Sampler rate of convergence we first need to compute the $(I+1)\times (I+1)$ matrix $A$ indexed by $(\alpha_1,\alpha_2)\in\{\mu,\beta_1,\dots,\beta_I\}\times \{\mu,\beta_1,\dots,\beta_I\}$ and defined as $A_{\alpha_1\alpha_2}=-\frac{Q_{\alpha_1\alpha_2}}{Q_{\alpha_1\alpha_1}}$
for $\alpha_1\neq \alpha_2$ and $0$ for $\alpha_1=\alpha_2$, where $Q$ is the precision matrix of $(\mu,\beta_1,\dots,\beta_I)|\textbf{y}$.
See \cite[Sec.2.2]{RobertsSahu1997} for more details on the derivation of such $A$-matrices.
By computing the precision matrix of $(\mu,\beta_1,\dots,\beta_I)|\textbf{y}$ it is easy to see that $A$ is given by 
\begin{equation*}
\begin{pmatrix}
0 & A_{\mu\beta_1} & \cdots & A_{\mu\beta_I}\\
A_{\beta_1\mu} & 0 & \cdots &0\\
\vdots & \vdots & \ddots & \vdots\\
A_{\beta_I\mu} & 0 & \cdots & 0 
\end{pmatrix}
\end{equation*}
where for all $i$
$$
A_{\mu\beta_i}=\frac{-\lambda_i\tprec{i}+(1-\lambda_i)\precision{a}}{\sum_{\ell=1}^I\lambda_\ell\tprec{\ell}+(1-\lambda_\ell)\precision{a}}
\,,\quad
A_{\beta_i\mu}=\frac{-\lambda_i\tprec{i}+(1-\lambda_i)\precision{a}}{\tprec{i}+\precision{a}}\,.
$$
From \citet[Thm.1]{RobertsSahu1997} the rate of convergence of the Gibbs sampler of interest equals the largest modulus eigenvalue of $B=(\I_{I+1}-L)^{-1}U$, where $\I_{I+1}$ is the $(I+1)$-dimensional identity matrix, $L$ is the lower triangular part of $A$ and $U=A-L$.
Simple calculations show that in this case
\begin{equation*}
B=
\begin{pmatrix}
0 &A_{\mu\beta_1}&\cdots&A_{\mu\beta_I}
\end{pmatrix}\varotimes
\begin{pmatrix}
1 \\A_{\beta_1\mu}\\\vdots\\A_{\beta_I\mu}
\end{pmatrix}
=
\begin{pmatrix}
0 & A_{\mu\beta_1} & \cdots & A_{\mu\beta_I}\\
0 &  A_{\mu\beta_1}A_{\beta_1\mu} & 0 & A_{\mu\beta_I}A_{\beta_1\mu}\\
\vdots &  \vdots & \ddots & \vdots\\
0 &  A_{\mu\beta_1}A_{\beta_I\mu}0 & \cdots & A_{\mu\beta_I}A_{\beta_I\mu}\,.
\end{pmatrix}
\end{equation*}
Finally note that $B$ has $I$ eigenvalues equal to 0 and one equal to 
$$
\sum_{i=1}^I A_{\mu\beta_i}A_{\beta_i\mu}
=
\frac{
\sum_{i\,:\,\lambda_i=1}\tprec{i}\frac{\tprec{i}}{\tprec{i}+\precision{a}}
+
\sum_{i\,:\,\lambda_i=0}\precision{a}\frac{\precision{a}}{\tprec{i}+\precision{a}}
}{
\sum_{i\,:\,\lambda_i=1}\tprec{i}+\sum_{i\,:\,\lambda_i=0}\precision{a}
}\,.
$$

%%%%%
\subsection*{Proofs of Lemmas \ref{lemma:hierarchical_reparametrization_invertible} and  \ref{lemma:H_closed}}
In order to prove Lemmas \ref{lemma:hierarchical_reparametrization_invertible} and  \ref{lemma:H_closed} we need some preliminary results on matrices $M=(M_{tr})_{t,r\in T}$ indexed by elements of a tree $T$.
%For any matrix $M=(M_{tr})_{t,r\in T}$ we define the following conditions
%\begin{align}\tag{H$_m$}\label{eq:H}
%M_{tr}=&0
%\hbox{ unless }
%t\preceq r
%\hbox{ or }
%r\preceq t\,,
%\\
%\tag{L$_m$}\label{eq:L}
%M_{tr}=&0
%\hbox{ unless }
%t\succeq r
%\hbox{ and }
%M_{tt}\neq 0 \quad\forall t\in T
%\,,
%\\
%\tag{L$^*_m$}\label{eq:LL}
%M_{tr}=&0
%\hbox{ unless }
%t\succ r
%\,,
%\\
%\tag{U$_m$}\label{eq:U}
%M_{tr}=&0
%\hbox{ unless }
%t\preceq r
%\hbox{ and }
%M_{tt}\neq 0 \quad\forall t\in T
%\,,
%\\
%\tag{U$^*_m$}\label{eq:UU}
%M_{tr}=&0
%\hbox{ unless }
%t\preceq r
%\,.
%\end{align}
%Here \eqref{eq:H} refers to the matrix analogous of \eqref{eq:hierarchical_structure}, 
%while 
%\eqref{eq:L},
%\eqref{eq:LL},
%\eqref{eq:U}
%and
%\eqref{eq:UU}
%refers to matrices that are, respectively, lower-triangular, strictly lower-triangular, upper-triangular and strictly upper-triangular with respect to the tree $T$, with additional conditions on the diagonals for  \eqref{eq:L}
%and
%\eqref{eq:U}.
%First we show that if a matrix $M$ satisfies respectively \eqref{eq:L} or \eqref{eq:U}, then it is invertible and its inverse still satisfies respectively \eqref{eq:L} or \eqref{eq:U}.
\begin{lemma}[Triangular matrices on trees]\label{lemma:triangular}
Suppose that a matrix $L=(L_{tr})_{t,r\in T}$ satisfies the following lower-triangularity condition
\begin{align}
\tag{L$_m$}\label{eq:L}
L_{tr}=&0
\hbox{ unless }
t\succeq r
\hbox{ and }
L_{tt}\neq 0 \quad\forall t\in T
\,.
\end{align}
Then $L$ is invertible and its inverse satisfies \eqref{eq:L}.
Similarly, if $U=(U_{tr})_{t,r\in T}$ satisfies the following upper-triangularity condition
\begin{align}
\tag{U$_m$}\label{eq:U}
U_{tr}=&0
\hbox{ unless }
t\preceq r
\hbox{ and }
U_{tt}\neq 0 \quad\forall t\in T
\,,
\end{align}
then $U$ is invertible and its inverse still satisfies \eqref{eq:U}.
\end{lemma}
\begin{proof}
Suppose that $L$ satisfies \eqref{eq:L}.
Also without loss of generality suppose $L_{tt}=1$ for all $t\in T$ by rescaling. 
Then write $L=(\mathbb{I}_T+N)$ where $\mathbb{I}_T$ is the $|T|\times |T|$ identity matrix and $N$ satisfies the following strict lower-triangularity condition
\begin{align}
\tag{L$^*_m$}\label{eq:LL}
N_{tr}=&0
\hbox{ unless }
t\succ r
\,.
\end{align}
Consider 
$N^2_{tr}=\sum_{s\in T}N_{ts}N_{sr}=\sum_{s\in T\,:\,r\prec s\prec t}N_{ts}N_{sr}$.
From the last expression it follows that $N^2_{tr}\neq 0$ implies
% there exist an $s\in T$ such that $r\prec s\prec t$. Therefore to have $(N)^2_{tr}\neq 0$ it must hold 
$r\prec t$ and $|\ell(r)-\ell(t)|\geq 2$, where $\ell(t)$ denote the level of $t$ in the tree $T$ .
Iterating the same argument we have that $N^p_{tr}\neq 0$ implies $r\prec t$ and $|\ell(r)-\ell(t)|\geq p$.
It follows that $N^p$ satisfies \eqref{eq:LL} for all $p\geq1$ and that $N^p=\boldsymbol{0}_T$ for all  $p\geq k$ where $\boldsymbol{0}_T$ is the $|T|\times |T|$ zero matrix.
Here $k$ indicates the number of levels of $T$, as in Section \ref{sec:k_levels}.
From $N^p=\boldsymbol{0}_T$ for $p\geq k$ it follows that 
$$
(\mathbb{I}_T+N)(\mathbb{I}_T+\sum_{p=1}^{k-1} (-1)^pN^p)
%=
%\mathbb{I}_T+\sum_{p=1}^{k-1} (-1)^pN^p+\sum_{p=0}^{k-1} (-1)^{p}N^{p+1}
=
\mathbb{I}_T+(-1)^{k-1}N^{k}
=
\mathbb{I}_T
$$ 
and therefore $L^{-1}=(\mathbb{I}_T+N)^{-1}=\mathbb{I}_T+\sum_{p=1}^{k-1} (-1)^pN^p$.
Since $N^p$ satisfies \eqref{eq:LL} for all $p\geq1$ it follows that 
$L^{-1}=\mathbb{I}_T+\sum_{p=1}^{k-1} (-1)^pN^p$ satisfies \eqref{eq:L}.

The analogous statement for \eqref{eq:U} can be deduced by observing that $U$ satisfies 
\eqref{eq:U} if and only if its transpose satisfies \eqref{eq:L}.
%, then $L^{-1}$ satisfies \eqref{eq:U} too.
%
%An analogous argument can be used to show that if $L$ satisfies \eqref{eq:U}, then $L^{-1}$ satisfies \eqref{eq:U} too.
\end{proof}

Lemma \ref{lemma:triangular} can be used to deduce Lemma \ref{lemma:hierarchical_reparametrization_invertible} in a straightforward way.
\begin{proof}[Proof of Lemma \ref{lemma:hierarchical_reparametrization_invertible}]
Suppose that $\bbeta_T=\Lambda\bgamma_T$ is a hierarchical reparametrization of $\bgamma_T$.
This is equivalent to say that $\Lambda=(\Lambda_{tr})_{t,r,\in T}$ is a matrix satisfying \eqref{eq:L} and that for all $t\in T$, $\bbeta_t=\sum_{s\in T}\Lambda_{ts}\gamma_s$. 
Lemma \ref{lemma:triangular} implies that $\Lambda$ is invertible and that $Z=\Lambda^{-1}$ satisfies \eqref{eq:L}.
Therefore $\bgamma_T=Z\bbeta_T$ is a hierarchical reparametrization of $\bbeta_T$.
\end{proof}

To prove Lemma \ref{lemma:H_closed} we need an additional preliminary result.

%%%%
\begin{lemma}[Closure of \ref{eq:H}]\label{lemma:closure_Hm}
Suppose that $M=(M_{tr})_{t,r\in T}$ satisfies the following condition
\begin{align}\tag{H$_m$}\label{eq:H}
M_{tr}=&0
\hbox{ unless }
t\preceq r
\hbox{ or }
r\preceq t\,.
\end{align}
Then if we multiply $M$ from the right with a martix $L$ satisfying \eqref{eq:L}, the product $ML$ still satisfy \eqref{eq:H}.
Similarly if we multiply $M$ from the left  with a martix $U$ satisfying \eqref{eq:U}, the product $UM$ still satisfy \eqref{eq:H}.
\end{lemma}
\begin{proof}
Consider
$(UM)_{tr}=\sum_{s\in T}U_{ts}M_{sr}$. From \eqref{eq:H} and \eqref{eq:U}, the elements $U_{ts}M_{sr}$ in the latter sum are non-zero only when $s$ belongs to the intersection of $\{s\in T\,:\,s\succeq t\}$ and $\{s\in T\,:\,s\succeq r\hbox{ or }s\preceq r\}$.
It is easy to see that the intersection of these two sets is non-empty only if $t\succeq r$ or $t\preceq r$.
Therefore $UM$ satisfies \eqref{eq:H}.
The argument to show that $ML$ satisfies \eqref{eq:H} is analogous.
\end{proof}
%%%%
We now combine Lemmas \ref{lemma:triangular} and \ref{lemma:closure_Hm} to prove Lemma \ref{lemma:H_closed}.
%%%%
\begin{proof}[Proof of Lemma \ref{lemma:H_closed}]
Suppose that $\bbeta_T$ satisfies \eqref{eq:hierarchical_structure}.
This is equivalent to saying that its precision matrix $Q^{(\bbeta)}=(Q^{(\bbeta)}_{tr})_{t,r\in T}$ satisfies \eqref{eq:H}.
Consider a hierarchical reparametrization of $\bbeta_T$ denoted by $\Lambda\bbeta_T$.
Then its precision matrix is given by 
$Q^{(\Lambda\bbeta)}=(\Lambda^T)^{-1}Q^{(\bbeta)}\Lambda^{-1}$ where $\Lambda^T$ denote the transpose of  $\Lambda$.
By definition of hierarchical reparametrizations, $\Lambda$ satisfies \eqref{eq:L}.
Therefore Lemma \ref{lemma:triangular} implies that $\Lambda^{-1}$ satisfies \eqref{eq:L} and, consequently, Lemma \ref{lemma:closure_Hm} implies that $Q^{(\bbeta)}\Lambda^{-1}$ satisfies \eqref{eq:H}.
Since $\Lambda$ satisfies \eqref{eq:L}, then $\Lambda^T$ satisfies \eqref{eq:U} and thus Lemma \ref{lemma:triangular} implies that $(\Lambda^T)^{-1}$ satisfies \eqref{eq:U}.
We can then apply Lemma \ref{lemma:closure_Hm} to $(\Lambda^T)^{-1}$ and $Q^{(\bbeta)}\Lambda^{-1}$ to deduce that $(\Lambda^T)^{-1}Q^{(\bbeta)}\Lambda^{-1}$ satisfies \eqref{eq:H} and thus $\Lambda\bbeta_T$ satisfies \eqref{eq:hierarchical_structure}.
\end{proof}
Corollary \ref{coroll:PNCP} follows easily from equation \eqref{eq:tree_cond_indep} and Lemma \ref{lemma:H_closed}.

%%%%%%
%\subsection*{Proof of Lemma \ref{lemma:hierarchical_reparametrization_invertible}}\label{proof:hierarchical_reparametrization_invertible}
%Fix a branch $B$ of the hierarchy $T$ and order the elements of $\bbeta_B=(\bbeta_t)_{t\in B\subseteq T}$ as $\bbeta_B=(\beta_0,\dots,\beta_k)$ going from root to leaf.
%Let $\Lambda_B=(\lambda_{d\ell})_{d,\ell=0}^{k-1}$ be the corresponding $k\times k$ matrix with $\lambda_{\beta_d\beta_\ell}=\lambda_{d\ell}$ given by \eqref{eq:hierarchical_reparametrization} for $\ell\leq d$ and $\lambda_{\ell d}=0$ for $\ell> d$.
%Since $\Lambda_B$ is lower-triangular and $\lambda_{d d}\neq 0$ for all $d=0,\dots,k$, then $\Lambda^{-1}=Z=(\zeta_{d\ell})_{\ell,d=1}^{k-1}$ exist and is lower-triangular. 
%From \eqref{eq:hierarchical_reparametrization} and $\Lambda^{-1}=Z$ it follows
%\begin{equation}\label{eq:reparametrization_inverse}
%\alpha_t=
%\sum_{r\preceq t}\zeta_{d\ell(r)}\, \beta_r\,,
%\end{equation}
%meaning that hierarchical reparametrizations are invertible and that their inverse is still a hierarchical reparametrizations.
%
%%%%%%
%\subsection*{Proof of Lemma \ref{lemma:H_closed}}\label{proof:H_closed}
%Let $\bbeta$ be a hierarchical reparametrization of $\balpha$.
%Suppose $p\leq d$, $t\in T_d$ and $r\in T_p$. Consider
%\begin{multline*}
%Cov\left(\beta_t,\beta_r
%\big|\bbeta_{T\backslash\{t,r\}}\right)
%=\\
%Cov\left(
%\lambda_{d d} \alpha_t+\sum_{s\prec t}\lambda_{d\ell(s)} \alpha_s,
%\lambda_{p p}\alpha_r+\sum_{s\prec r}\lambda_{p\ell(s)} \alpha_s
%\big|\bbeta_{T\backslash\{t,r\}}\right)\,.
%\end{multline*}
%For any $s\prec t$, equation \eqref{eq:reparametrization_inverse} implies that $\alpha_s$ is a deterministic function of $\bbeta_{T\backslash\{t,r\}}$. 
%Similarly for $\alpha_s$ and $s\prec r$.
%Therefore
%\begin{align*}
%Cov\left(\beta_t,\beta_r
%\big|\bbeta_{T\backslash\{t,r\}}\right)
%&=
%Cov\left(
%\lambda_{d d} \alpha_t,\lambda_{p p}\alpha_r
%\big|\bbeta_{T\backslash\{t,r\}}\right)\\
%&=
%\lambda_{d d}\lambda_{p p}Cov\left(
% \alpha_t,\alpha_r
%\big|\bbeta_{T\backslash\{t,r\}}\right)\,.
%\end{align*}
%which is 0 because $\boldsymbol{\alpha}$ satisfies \eqref{eq:hierarchical_structure}.
%\fotnote{MISSING: the last statement is not justified because conditioning in on $\bbeta$ and not $\balpha$.}


%%%%%
%%%%%
\subsection*{Proof of Lemma \ref{lemma:S_closed}}
The strategy to prove Lemma \ref{lemma:S_closed} is similar to the one used to prove Lemmas \ref{lemma:hierarchical_reparametrization_invertible} and  \ref{lemma:H_closed} above, with the difference that we have to check that the symmetry condition is preserved under the operations considered. To do so we first prove two auxiliary lemmas.

\begin{lemma}[Symmetric triangular matrices on trees]\label{lemma:symmetric_triangular}
Suppose that a matrix $L=(L_{tr})_{t,r\in T}$ satisfies the following symmetric lower-triangularity condition
\begin{align}
\tag{SL$_m$}\label{eq:SL}
L_{tr}=&l_{\ell(t)\ell(r)}\1(t\succeq r)
\hbox{ and }
l_{dd}\neq 0
\hbox{ for all }
 d\in\{0,\dots,k-1\}
\,,
\end{align}
where $(l_{pd})_{p,d=0}^{k-1}$ is a $k\times k$ real valued matrix.
Then $L$ is invertible and its inverse satisfies \eqref{eq:SL}.
Similarly, if $U=(U_{tr})_{t,r\in T}$ satisfies the following symmetric upper-triangularity condition
\begin{align}
\tag{SU$_m$}\label{eq:SU}
U_{tr}=&u_{\ell(t)\ell(r)}\1(t\preceq r)
\hbox{ and }
u_{dd}\neq 0
\hbox{ for all }
 d\in\{0,\dots,k-1\}
\,,
\end{align}
with $(u_{pd})_{p,d=0}^{k-1}$ being a $k\times k$ real valued matrix,
then $U$ is invertible and its inverse still satisfies \eqref{eq:SU}.
\end{lemma}
\begin{proof}
Suppose that $L$ satisfies \eqref{eq:SL}.
Also without loss of generality suppose $L_{tt}=1$ for all $t\in T$ by rescaling.
Since $L$ also satisfies \eqref{eq:L}, arguing as in the proof of Lemma \ref{lemma:triangular} we can write 
$L^{-1}=\mathbb{I}_T+\sum_{p=1}^{k-1} (-1)^pN^p$ 
where
$N=L-\mathbb{I}_T$ 
and $N$ satisfies the following symmetric strict lower-triangularity condition
\begin{align}
\tag{SL$^*_m$}\label{eq:SSL}
N_{tr}=&l_{\ell(t)\ell(r)}\1(t\succ r)
\hbox{ and }
l_{dd}\neq 0
\hbox{ for all }
 d\in\{0,\dots,k-1\}
\,.
\end{align}
%where $(l_{pd})_{p,d=0}^{k-1}$ is a $k\times k$ real valued matrix.
From 
$N^2_{tr}=
\sum_{s\in T}N_{ts}N_{sr}=
\sum_{s\in T\,:\,r\prec s\prec t}N_{ts}N_{sr}=
\sum_{\ell'\,:\,\ell(r)\prec \ell' \prec \ell(t)}l_{\ell(t)\ell'}l_{\ell'\ell(r)}$
we deduce that $N^2$ still satisfies \eqref{eq:SSL} for some different $(l_{pd})_{p,d=0}^{k-1}$.
Iterating the same argument we have that $N^p$ satisfies \eqref{eq:SSL} for all $p\geq 1$.
Since $L^{-1}=\mathbb{I}_T+\sum_{p=1}^{k-1} (-1)^pN^p$ it follows that $L^{-1}$ satisfies \eqref{eq:SL}.

%An analogous argument can be used to show that if $U$ satisfies \eqref{eq:SU}, then $U^{-1}$ satisfies \eqref{eq:SU} too.
The analogous statement for \eqref{eq:SU} can be deduced by observing that $U$ satisfies \eqref{eq:SU} if and only if its transpose satisfies \eqref{eq:SL}.
\end{proof}
%o prove Lemma \ref{lemma:S_closed} we need an additional preliminary result.
%%%%
%\begin{lemma}[Closure of \ref{eq:S}]\label{lemma:closure_Sm}
%Suppose that $M=(M_{tr})_{t,r\in T}$ satisfies the following condition
%\fotnote{This condition is identical to S tilde}
%\begin{align}\tag{S$_m$}\label{eq:S}
%M_{tr}=
%&
%m_{\ell(t)\ell(r)}
%P(X_{\ell(r)}=r\cap X_{\ell(t)}=t)
%\,,
%\end{align}
%where $(X_\ell)_{\ell=0}^{k-1}$ is a Markov chain going through $T$ from root to leaves\footnote{Define better here.}
%Then if we multiply $M$ from the right with a martix $L$ satisfying \eqref{eq:SL}, the product $ML$ still satisfy \eqref{eq:S}.
%Similarly if we multiply $M$ from the left  with a martix $U$ satisfying \eqref{eq:SU}, the product $UM$ still satisfy \eqref{eq:S}.
%\end{lemma}
\begin{lemma}\label{lemma:closure_Sm}%[Closure of \eqref{eq:k_symmetry_rescaled}]
Let $M=(M_{tr})_{t,r\in T}$ and $L=(L_{tr})_{t,r\in T}$ satisfy respectively \eqref{eq:k_symmetry_rescaled} and \eqref{eq:SL}.
Then the product $ML$ satisfies \eqref{eq:k_symmetry_rescaled} for some different matrix $(c_{pd})_{p,d=1}^{k-1}$.
Similarly, if $U=(U_{tr})_{t,r\in T}$ satisfies \eqref{eq:SU} then the product $UM$ satisfy \eqref{eq:k_symmetry_rescaled}.
\end{lemma}
\begin{proof}
First consider the product $ML$.
Lemma \ref{lemma:closure_Hm} implies that $ML$ satisfies \eqref{eq:H} and therefore $(ML)_{tr}=0$ unless $t\preceq r$ or $t\succeq r$.
Given $t\preceq r$ or $t\succeq r$ and using \eqref{eq:k_symmetry_rescaled} and \eqref{eq:SL} we have
\begin{align*}
&(ML)_{tr}
=
\sum_{s\in T}M_{ts}L_{sr}
=
\sum_{s\in T\,:\,s\succeq r}M_{ts}L_{sr}
\\&=
\sum_{s\in T\,:\,s\succeq (t\vee r) }M_{ts}L_{sr}
+
\sum_{s\in T\,:\,r\preceq s\prec t}M_{ts}L_{sr}
\\&=
\sum_{s\in T\,:\,s\succeq (t\vee r)}
m_{\ell(t)\ell(s)}l_{\ell(s)\ell(r)}P(s)
+
\sum_{s\in T\,:\,r\preceq s\prec t}
m_{\ell(t)\ell(s)}l_{\ell(s)\ell(r)}P(t)
\\&=
P(t\vee r)
\sum_{\ell'=\ell(t\vee r)}^{k-1}m_{\ell(t)\ell'}l_{\ell'\ell(r)}
+
P(t)
\sum_{\ell'=\ell(r)}^{\ell(t)-1}m_{\ell(t)\ell'}l_{\ell'\ell(r)}
\\&=
P(t\vee r)
\sum_{\ell'=\ell(r)}^{k-1}m_{\ell(t)\ell'}l_{\ell'\ell(r)}\,,
\end{align*}
where the sum from $\ell(r)$ to $\ell(t)$ equals 0 if $\ell(r)\geq \ell(t)$.
The latter equation implies that $ML$ satisfies \eqref{eq:k_symmetry_rescaled}.

To prove that the product $UM$ satisfies \eqref{eq:k_symmetry_rescaled} note that $M$ and $U$ satisfying \eqref{eq:k_symmetry_rescaled} and \eqref{eq:SU} respectively is equivalent to $M^T$ and $U^T$ satisfying \eqref{eq:k_symmetry_rescaled} and \eqref{eq:SL} respectively.
Therefore, by the first part of this Lemma, $(UM)^T=M^TU^T$ satisfies \eqref{eq:SL} and thus $((UM)^T)^T=UM$ satisfies \eqref{eq:k_symmetry_rescaled}.
%\fotnote{Can't we just say that lower triangular is the transpose of upper triangular?}
%Similarly, given $t\preceq r$ or $t\succeq r$ and using \eqref{eq:S} and \eqref{eq:SU} we have
%\begin{align*}
%&(UM)_{tr}
%=
%\sum_{s\in T}U_{ts}M_{sr}
%\\&=
%\sum_{s\in T\,:\,s\succeq (t\vee r) }U_{ts}M_{sr}
%+
%\sum_{s\in T\,:\,t\preceq s\prec r}U_{ts}M_{sr}
%\\&=
%\sum_{s\in T\,:\,s\succeq (t\vee r)}
%u_{\ell(t)\ell(s)}m_{\ell(s)\ell(r)}P(s)
%+
%\sum_{s\in T\,:\,t\preceq s\prec r}
%u_{\ell(t)\ell(s)}m_{\ell(s)\ell(r)}P(r)
%\\&=
%P(t\vee r)
%\sum_{\ell'=\ell(t\vee r)}^{k-1}u_{\ell(t)\ell'}m_{\ell'\ell(r)}
%+
%\frac{1}{\prod_{v\in T\,:\,v\prec t}n_v}
%\sum_{\ell'=\ell(r)}^{\ell(t)-1}u_{\ell(t)\ell'}m_{\ell'\ell(r)}
%\\&=
%\frac{1}{\prod_{v\in T\,:\,v\prec (t\vee r)}n_v}
%\sum_{\ell'=\ell(r)}^{k}u_{\ell(t)\ell'}m_{\ell'\ell(r)}\,,
%\end{align*}
%which implies that $UM$ satisfies \eqref{eq:S}.
\end{proof}
%%%%
We now combine Lemmas \ref{lemma:symmetric_triangular} and \ref{lemma:closure_Sm} to prove Lemma \ref{lemma:S_closed}.
%%%%
\begin{proof}[Proof of Lemma \ref{lemma:S_closed}]
Suppose that the precision matrix $Q^{(\tbbeta)}=(Q^{(\tbbeta)}_{tr})_{t,r\in T}$ of $\tbbeta_T$ satisfies \eqref{eq:k_symmetry_rescaled}.
Consider a symmetric hierarchical reparametrization of $\tbbeta_T$ denoted by $\Lambda\tbbeta_T$.
Then its precision matrix is given by 
$Q^{(\Lambda\tbbeta)}=(\Lambda^T)^{-1}Q^{(\tbbeta)}\Lambda^{-1}$ where $\Lambda^T$ denote the transpose of  $\Lambda$.
By definition of symmetric hierarchical reparametrizations, $\Lambda$ satisfies \eqref{eq:SL}.
Therefore Lemma \ref{lemma:symmetric_triangular} implies that $\Lambda^{-1}$ satisfies \eqref{eq:SL} and, consequently, Lemma \ref{lemma:closure_Sm} implies that $Q^{(\tbbeta)}\Lambda^{-1}$ satisfies \eqref{eq:k_symmetry_rescaled}.
Since $\Lambda$ satisfies \eqref{eq:SL}, then $\Lambda^T$ satisfies \eqref{eq:SU} and thus Lemma \ref{lemma:triangular} implies that $(\Lambda^T)^{-1}$ satisfies \eqref{eq:SU}.
We can then apply Lemma \ref{lemma:closure_Sm} to $(\Lambda^T)^{-1}$ and $Q^{(\tbbeta)}\Lambda^{-1}$ to deduce that $Q^{(\Lambda\tbbeta)}=(\Lambda^T)^{-1}Q^{(\tbbeta)}\Lambda^{-1}$ satisfies \eqref{eq:k_symmetry_rescaled}.
\end{proof}




%%%%%
\subsection*{Proof of Lemma \ref{lemma:delta_expansion}}\label{proof:delta_expansion}
Suppose $\bbeta_T$ has zero mean (otherwise replace $\bbeta_T$ by $\bbeta_T-\mathbb{E}[\bbeta]$).
As in Section \ref{sec:multigrid_k}, given any $r$ and $t$ in $T$ we denote $P(X_{\ell(t)}=t|X_{\ell(r)}=r)$ by $P(t|r)$.
Using \eqref{eq:hierarchical_structure}, for any $d\in\{0,\dots,k-1\}$ and $t\in T_d$, we can write the full conditional expectation of $\beta_t$ as
\begin{align*}
\E[\beta_t|\bbeta_{T\backslash t}]
%&=
%\sum_{r\in T\backslash t}
%\frac{
%Cov\left(\beta_t,\beta_r|\bbeta\backslash\{\beta_t,\beta_r\}\right)
%}{
%Var\left(\beta_t|\bbeta_{T\backslash t}\right)
%}
%\beta_r
%\\
&=
\sum_{r\prec t}
A_{tr}
\beta_r
+
\sum_{r\succ t}
A_{tr}
\beta_r\,,
\end{align*}
where $A_{tr}=-\frac{Q_{tr}}{Q_{tt}}$ for any $r\neq t$.
Note that \eqref{eq:k_symmetry_rescaled} implies
\begin{equation}\label{eq:symm_A}
A_{tr}=
\frac{c_{\ell(t)\ell(r)}P(r\cap t)}{P(t)}
=
c_{\ell(t)\ell(r)}P(r|t)\,.
\end{equation}
It follows 
\begin{align*}
\E[\beta_t|\bbeta_{T\backslash t}]
&=
\sum_{r\prec t}
c_{d\ell(r)}P(r|t)
\beta_r
+
\sum_{r\succ t}
c_{d\ell(r)}P(r|t)
\beta_r
\\
%&=
%\sum_{r\prec t}
%c_{d\ell(r)}
%\beta_r
%+
%\sum_{\ell=d+1}^{k}
%c_{d\ell}
%\phi_t^{(d)}\bbeta^{(\ell)}\\
&=
\sum_{\ell\neq d}
c_{d\ell}
\E[\beta_{X_\ell}|\bbeta_T,X_d=t]\,.
\end{align*}
Since the last equation does not depend on $\bbeta^{(d)}$ we have
$
\E[\beta_t|\bbeta_{T\backslash t}]=
\E[\beta_t|\bbeta\backslash\bbeta^{(d)}].$
For any $r\in T_p$, by definition of $\phi^{(p)}_r\bbeta^{(d)}$, we have
\begin{align*}%\label{eq:phi_expansion}
\E[\phi^{(p)}_r\bbeta^{(d)}|\bbeta\backslash\bbeta^{(d)}]
&=
\sum_{t\in T_d}
\Pr(t|r)
\E[\beta_{t}|\bbeta\backslash\bbeta^{(d)}]
\\
&=
\sum_{t\in T_d}
\Pr(t|r)
\sum_{\ell\neq d}
c_{d\ell}
\E[\beta_{X_\ell}|\bbeta_T,X_d=t]
\\
&=
\sum_{\ell\neq d}
c_{d\ell}
\sum_{t\in T_d}
\Pr(t|r)
\E[\beta_{X_\ell}|\bbeta_T,X_d=t]
\\
&=
\sum_{\ell\neq d}
c_{d\ell}
\E[\beta_{X_\ell}|\bbeta_T,X_p=r]
\,.
\end{align*}
From the latter equation and the definition of $\delta^{(p)}_r\bbeta^{(d)}$ it follows
\begin{align*}
\E[\delta^{(p)}_r\bbeta^{(d)}|\bbeta\backslash\bbeta^{(d)}]
&=
\E[\phi^{(p)}_r\bbeta^{(d)}|\bbeta\backslash\bbeta^{(d)}]-\E[\phi^{(p-1)}_{pa(r)}\bbeta^{(d)}|\bbeta\backslash\bbeta^{(d)}]
\\
&=
\sum_{\ell\neq d}
c_{d\ell}
\left(
\E[\beta_{X_\ell}|\bbeta_T,X_p=r]
-
\E[\beta_{X_\ell}|\bbeta_T,X_{p-1}=pa(r)]
\right)
\\
&=
\sum_{\ell\in\{p,\dots,k-1\}\backslash d}
c_{d\ell}
\delta^{(p)}_r\bbeta^{(\ell)}
\,.
\end{align*}

%%%%%
\subsection*{Proof or Theorem \ref{thm:factorization_general}}\label{proof:factorization_general}
To prove Theorem \ref{thm:factorization_general} we first need the following lemma.
\begin{lemma}\label{lemma:residual_independence}
Given $d\in\{0,\dots,k-1\}$ and $p,p'\in\{0,\dots,d\}$ with $p\neq p'$, %itholds the following conditional independence
\begin{align*}
\delta^{(p)}\bbeta^{(d)}
\bot
\delta^{(p')}\bbeta^{(d)}
\,|\,
\bbeta\backslash\bbeta^{(d)}
\,.
\end{align*}
\end{lemma}
\begin{proof}
Let $d\in\{0,\dots,k-1\}$, $p,p'\in\{0,\dots,d\}$ with $p< p'$ and $\bbeta\backslash\bbeta^{(d)}$ be fixed.
To make the notation more compact we denote $\E[\cdot|\bbeta\backslash\bbeta^{(d)}]$ by $\tE[\cdot]$, $\phi^{(p)}_{r}\bbeta^{(d)}$ by $\tilde{\phi}^{(p)}_{r}$ and $P(X_{p'}=r'|X_p=r)$ by $P(r'|r)$ for all $r\in T_p$ and $r'\in T_{p'}$.
By replacing $\beta_t$ with $\beta_t-\tE[\beta_t]$ we can suppose without loss of generality that $\tE[\beta_t]=0$ for all $t\in T_d$ and therefore 
$\tE[\tilde{\phi}^{(p)}_{r}]=\tE[\tilde{\phi}^{(p')}_{r'}]=0$ and 
$\tE[\delta^{(p)}_r\bbeta^{(d)}]=\tE[\delta^{(p')}_{r'}\bbeta^{(d)}]=0$
for all $r\in T_{p}$ and $r'\in T_{p'}$.
By definition $\tilde{\phi}^{(p)}_r=\sum_{s\in T_{p'}}
P(s|r)
\tilde{\phi}^{(p')}_s$ and therefore
\begin{equation}
\tE[\tilde{\phi}^{(p)}_r\tilde{\phi}^{(p')}_{r'}]
=
\sum_{s\in T_{p'}}
P(s|r)
\tE[\tilde{\phi}^{(p')}_s\tilde{\phi}^{(p')}_{r'}]
=
P(r'|r)
\tE[(\tilde{\phi}^{(p')}_{r'})^2]\,,\label{eq:phi_dot_product}
\end{equation}
where we used $\tE[\tilde{\phi}^{(p')}_s\tilde{\phi}^{(p')}_{r'}]=0$  for $r'\neq s$ and $r',s\in T_{p'}$, which follows from the conditional independence of $(\beta_t)_{t\in T_d\,,\,t\succeq r'}$ and $(\beta_t)_{t\in T_d\,,\,t\succeq s}$ given $\bbeta\backslash\bbeta^{(d)}$.
Note that $P(r'|r)$ in \eqref{eq:phi_dot_product} could be 0.
From $\tE[\beta_t]=0$ for any $t\in T_d$ and \eqref{eq:k_symmetry_rescaled} we have $\tE[\beta_t^2]=P(t)^{-1}$ and therefore
\begin{align}
\tE[(\tilde{\phi}^{(p')}_{r'})^2]
=&
\tE[(\sum_{t\in T_d}P(t|r')\beta_t)^2]
=
\sum_{t\in T_d}P(t|r')^2\tE[\beta_t^2]
\nonumber\\
=&
\sum_{t\in T_d}\frac{P(t|r')^2}{P(t)}
=
\frac{1}{P(r')}\sum_{t\in T_d}P(t|r')
=
\frac{1}{P(r')}\,.
\label{eq:phi_squared}
\end{align}
Combining \eqref{eq:phi_dot_product} and \eqref{eq:phi_squared} we have 
$\tE[\tilde{\phi}^{(p)}_r\tilde{\phi}^{(p')}_{r'}]=0$ if $r'\nsucc r$ and
\begin{align}\label{eq:zero_simplification}
\tE[\tilde{\phi}^{(p)}_r\tilde{\phi}^{(p')}_{r'}]
=&
\frac{P(r'|r)}{P(r')}
=
\frac{1}{P(r)}
&
\hbox{if }r'\succ r\,.
\end{align}
From the last equality and the definition of $\delta^{(p)}_r\bbeta^{(d)}$ in \eqref{eq:delta_defi} we have
\begin{align*}
\tE&[\delta^{(p)}_r\bbeta^{(d)}\delta^{(p')}_{r'}\bbeta^{(d)}]
=\\&
\tE[(\tilde{\phi}^{(p)}_r-\tilde{\phi}^{(p-1)}_{pa(r)})(\tilde{\phi}^{(p')}_{r'}-\tilde{\phi}^{(p'-1)}_{pa(r')})]
=\\&
\tE[
\tilde{\phi}^{(p)}_r\tilde{\phi}^{(p')}_{r'}
-\tilde{\phi}^{(p-1)}_{pa(r)}\tilde{\phi}^{(p')}_{r'}
-\tilde{\phi}^{(p)}_r\tilde{\phi}^{(p'-1)}_{pa(r')}
+\tilde{\phi}^{(p-1)}_{pa(r)}\tilde{\phi}^{(p'-1)}_{pa(r')}
]
=\\&
\frac{P(r'|r)}{P(r')}
-
\frac{P(r'|pa(r))}{P(r')}
-
\frac{P(pa(r')|r)}{P(pa(r'))}
+
\frac{P(pa(r')|pa(r))}{P(pa(r'))}
=0,
\end{align*}
where the last equality is trivial if $pa(r)\nprec r'$ and can be deduced from \eqref{eq:zero_simplification} otherwise.
The desired conditional independence follows from
$\tE[\delta^{(p)}_r\bbeta^{(d)}\delta^{(p')}_{r'}\bbeta^{(d)}]=0=\tE[\delta^{(p)}_r\bbeta^{(d)}]\tE[\delta^{(p')}_{r'}\bbeta^{(d)}]$
for all $r\in T_p$ and $r'\in T_{p'}$.
\end{proof}

\begin{proof}[Proof or Theorem \ref{thm:factorization_general}]
Theorem \ref{thm:factorization_general} follows easily from Lemmas \ref{lemma:delta_expansion} and \ref{lemma:residual_independence} as follows.
For each $d\in\{0,\dots,k-1\}$
the sampling step 
\begin{align}\label{eq:sample_original}
\bbeta^{(d)}(s+1)
\;\sim\;
\mathcal{L}\left(\bbeta^{(d)}|
\big(\bbeta^{(\ell)}(s+1)\big)_{0\leq \ell<d}\,,
\big(\bbeta^{(\ell)}(s)\big)_{d<\ell\leq k-1}
\right)
\end{align}
in Sampler \ref{sampler:GSk} is equal in distribution to sampling jointly the (d+1) residuals 
\begin{align*}
(\delta^{(p)}\bbeta^{(d)}(s+1))_{p=0}^d
\;\sim\;
\mathcal{L}\left((\delta^{(p)}\bbeta^{(d)})_{p=0}^d|
\big(\bbeta^{(\ell)}(s+1)\big)_{0\leq \ell<d}\,,
\big(\bbeta^{(\ell)}(s)\big)_{d<\ell\leq k-1}
\right)\,.
\end{align*}
From the conditional independence statement in Lemma  \ref{lemma:residual_independence} the latter is equivalent to sampling independently each residual $\delta^{(p)}\bbeta^{(d)}(s+1)$ from
$$
\mathcal{L}\left(\delta^{(p)}\bbeta^{(d)}|
\big(\bbeta^{(\ell)}(s+1)\big)_{0\leq \ell<d}\,,
\big(\bbeta^{(\ell)}(s)\big)_{d<\ell\leq k-1}
\right)\,.
$$
Moreover, from Lemma \ref{lemma:delta_expansion}
\begin{align*}
\mathcal{L}\left(\delta^{(p)}\bbeta^{(d)}|
\bbeta\backslash\bbeta^{(d)}
\right)
=
\mathcal{L}\left(\delta^{(p)}\bbeta^{(d)}|
(\delta^{(p)}\bbeta^{(\ell)})_{\ell\in\{p,\dots,k-1\}\backslash d}
\right)\,.
\end{align*}
Therefore the original sampling step in \eqref{eq:sample_original} is equivalent to sampling independently
\begin{align}\label{eq:sample_reparametrized}
\delta^{(p)}\bbeta^{(d)}(s+1)
\;\sim\;
\mathcal{L}\left(\delta^{(p)}\bbeta^{(d)}|
\big(\delta^{(p)}\bbeta^{(\ell)}(s+1)\big)_{p\leq \ell<d}\,,
\big(\delta^{(p)}\bbeta^{(\ell)}(s)\big)_{d<\ell\leq k-1}
\right)\,,
\end{align}
for $p=0,\dots,d$.
The thesis follows from the equivalence between \eqref{eq:sample_original} and \eqref{eq:sample_reparametrized}.
\end{proof}

\subsection*{Proof of Corollary \ref{coroll:rate_k_levels}}
\begin{proof}[Proof of Corollary \ref{coroll:rate_k_levels}]
The map from $\bbeta_T$ to $(\delta^{(0)}\bbeta,\dots,\delta^{(k-1)}\bbeta)$ is an injective linear transformation. 
The injectivity holds because for any $d\in\{0,\dots,k-1\}$ and $t\in T_d$ we can reconstruct $\beta_t$ from $(\delta^{(0)}\bbeta^{(d)},\dots,\delta^{(d)}\bbeta^{(d)})$
$$
\sum_{r\preceq t}\delta_r^{(\ell(r))}\bbeta^{(d)}
=
\phi_{t_0}^{(0)}\bbeta^{(d)}+
\sum_{t_0\prec r\preceq t}\left(\phi_r^{(\ell(r))}\bbeta^{(d)}-\phi_{pa(r)}^{(\ell(r)-1)}\bbeta^{(d)}\right)
=
\phi_{t}^{(d)}\bbeta^{(d)}
=
\beta_t\,.
$$
It follows that 
$(\delta\bbeta(s))_{s\in \N}=(\delta^{(0)}\bbeta(s),\dots,\delta^{(k-1)}\bbeta(s))_{s\in \N}$ is a Markov chain with the same rate of convergence of the original chain $(\bbeta(s))_{s\in \N}$.
%Moreover, from Theorem \ref{thm:factorization_general} it follows that $(\delta^{(0)}\bbeta(s))_s$, \dots, $(\delta^{(k)}\bbeta(s))_s$ are $k$ independent Markov chains.
Then the thesis follows from Theorem \ref{thm:factorization_general} and the fact that the rate of convergence of a collection of independent Markov chains equals the supremum of the rates of convergence of the single chains.
\end{proof}

%%%%
\subsection*{Proof of Theorem \ref{thm:kronecker}}
\begin{proof}[Proof of Theorem \ref{thm:kronecker}]
We are interested in the rate of convergence of 
the blocked sampler targeting 
$\delta^{(p)}\bbeta=(\delta^{(p)}\bbeta^{(p)},\dots,\delta^{(p)}\bbeta^{(k-1)})$ 
and evolving according to
\eqref{eq:blocked_gibbs}.
Consider first the case $p\in\{1,\dots,k-1\}$.
Note that $\delta^{(p)}\bbeta$ has a singular variance-covariance matrix because for each $t\in T_{p-1}$ and $d\in\{p,\dots,k-1\}$ it follows from \eqref{eq:phi_defi} and \eqref{eq:delta_defi} that
\begin{align}\label{eq:delta_sum_zero}
\sum_{r\in ch(t)}P(r|t)\delta^{(p)}_r\bbeta^{(d)}=0
\end{align}
and therefore some elements of $(\delta^{(p)}\bbeta)$ are linear combinations of the others.
In order to use standard tools it is more convenient to work with non-singular Gaussian random vectors.
To do so it is sufficient to consider a sub-vector of $\delta^{(p)}\bbeta$ obtained by removing from $T_{p}$ one children node for each parent node in $T_{p-1}$.
More formally, let $f$ be an arbitrary map from $T_{p-1}$ to $T_{p}$ such that $f(t)\in ch(t)$ for all $t\in T_{p-1}$ and then define the subset $T'_p\subseteq T_p$ as $T'_p= T_p\backslash f(T_{p-1})$.
It is then easy to see that the resulting sub-vector
$\delta_{T'_p}^{(p)}\bbeta
%=(\delta_r^{(p)}\bbeta)_{r\in T'_p}
=(\delta_{T'_p}^{(p)}\bbeta^{(p)},\dots,\delta_{T'_p}^{(p)}\bbeta^{(k-1)})
$
with $\delta_{T'_p}^{(p)}\bbeta^{(d)}=(\delta_r^{(p)}\bbeta^{(d)})_{r\in T'_p}$ for all $d\in\{p,\dots,k-1\}$ has an invertible variance-covariance matrix.
Moreover, since each $\delta^{(p)}\bbeta^{(d)}$ is a function of the corresponding $\delta_{T'_p}^{(p)}\bbeta^{(d)}$ via \eqref{eq:delta_sum_zero} it follows that the blocked sampler targeting 
$\delta^{(p)}\bbeta$ 
and evolving according to
\eqref{eq:blocked_gibbs} is equivalent in distribution to a blocked Gibbs Sampler targeting $\delta_{T'_p}^{(p)}\bbeta$ and evolving according to
\begin{align}\label{eq:blocked_gibbs_reduced}
\delta^{(p)}_{T'_p}\bbeta^{(d)}(s+1)
\;\sim\;
\mathcal{L}\left(\delta^{(p)}_{T'_p}\bbeta^{(d)}|
(\delta^{(p)}_{T'_p}\bbeta^{(\ell)}(s+1))_{p\leq \ell<d},(\delta^{(p)}_{T'_p}\bbeta^{(\ell)}(s))_{d<\ell\leq k-1}
\right)\,,
\end{align}
for $d\in\{p,\dots,k-1\}$.
Let 
%$Q_{\delta^{(p)}_{T'_p}\bbeta}$%=(Q_{\delta_r^{(p)}\bbeta^{(d)}\,\delta_{r'}^{(p)}\bbeta^{(d')}})_{r,r'\in T'_p\,,\, d,d'\in\{p,\dots,k-1\}}
% be the precision matrix of $\delta^{(p)}_{T'_p}\bbeta$ and
 $A_{\delta^{(p)}_{T'_p}\bbeta}=(A_{\delta_r^{(p)}\bbeta^{(d)}\,\delta_{r'}^{(p)}\bbeta^{(d')}})_{r,r'\in T'_p\,,\, d,d'\in\{p,\dots,k-1\}}$ be the $A$-matrix associated to the Gibbs sampler in \eqref{eq:blocked_gibbs_reduced}, defined by
\begin{multline}\label{eq:A_mat_defi}
\E[\delta^{(p)}_r\bbeta^{(d)}|\delta^{(p)}_{T'_p}\bbeta\backslash \delta^{(p)}_r\bbeta^{(d)}]
-\E[\delta^{(p)}_r\bbeta^{(d)}]
=\\
\sum_{\delta^{(p)}_{r'}\bbeta^{(d')}\in \delta^{(p)}_{T'_p}\bbeta\backslash \delta^{(p)}_r\bbeta^{(d)}}
%_{d'\in\{p,\dots,k-1\}\backslash d}
A_{\delta_r^{(p)}\bbeta^{(d)}\,\delta_{r'}^{(p)}\bbeta^{(d')}}
\left(\delta^{(p)}_{r'}\bbeta^{(d')}-\E[\delta^{(p')}_{r'}\bbeta^{(d')}]\right)
\,.
\end{multline}
See \cite[Sec.2.2]{RobertsSahu1997} for more details on the properties of such $A$-matrices.
%\fotnote{Gareth: check this proof, in particular this definition of $A$-matrix.}
%The latter matrix is defined as zero on the diagonal and $A_{ij}=-\frac{Q_{ij}}{Q_{ii}}$ on the off-diagonal entries $(i,j)$.
%See \cite[Sec.2.2]{RobertsSahu1997} for more details on the properties of such $A$-matrices.
% and satisfies 
%$
%\E[\delta^{(p)}_{r}\bbeta^{(d)}|\delta^{(p)}_{T'_p}\bbeta\backslash \delta^{(p)}_{r}\bbeta^{(d)}]-\E[\delta^{(p)}_{r}\bbeta^{(d)}]
%=
%\sum_{r'\in T'_p\,,\, d'\in\{p,\dots,k-1\}}\left(\E[\delta^{(p)}_{r}\bbeta^{(d)}|\delta^{(p)}_{T'_p}\bbeta\backslash \delta^{(p)}_{r}\bbeta^{(d)}]-\E[\delta^{(p)}_{r}\bbeta^{(d)}|\delta^{(p)}_{T'_p}\bbeta\backslash \delta^{(p)}_{r}\bbeta^{(d)}]\right)\,.
%$
Then Lemma \ref{lemma:delta_expansion} implies that 
$A_{\delta_r^{(p)}\bbeta^{(d)}\,\delta_{r'}^{(p)}\bbeta^{(d')}}=c_{dd'}$
if $r=r'$ and $d'\in\{p,\dots,k-1\}\backslash d$ and 0 otherwise.
The latter is equivalent to the equation
\begin{equation}\label{eq:A_matrix_delta}
A_{\delta^{(p)}_{T'_p}\bbeta}=\left(C^{(p)}-\I_{k-p}\right)\varotimes \I_{|T'_p|}
\end{equation}
where $C^{(p)}$ is the $(k-p)\times (k-p)$ square matrix $C^{(p)}=(c_{dd'})_{d,d'=p}^{k-1}$, $\I_n$ denotes the $n$ dimensional identity matrix, $\varotimes$ denotes the Kronecker product of matrices and $|T'_p|=|T_p|-|T_{p-1}|$ is the cardinality of $T'_p$.
Then note that from \citet[Thm.1]{RobertsSahu1997} the rate of convergence of the Gibbs sampler in \eqref{eq:blocked_gibbs_reduced} equals the largest modulus eigenvalue of $B=(\I_{(k-p)|T'_p|}-L)^{-1}U$, where $L$ is the lower triangular part of $A_{\delta^{(p)}_{T'_p}\bbeta}$ and $U=A_{\delta^{(p)}_{T'_p}\bbeta}-L$.
Using basic properties of the Kroenecker product we can see that 
\begin{align*}
B=&
(\I_{(k-p)}\varotimes\I_{|T'_p|}-\tilde{L}\varotimes\I_{|T'_p|})^{-1}\left(\tilde{U}\varotimes \I_{|T'_p|}\right)
\\=&
\left((\I_{(k-p)}-\tilde{L})\varotimes \I_{|T'_p|}\right)^{-1}\left(\tilde{U}\varotimes \I_{|T'_p|}\right)
\\=&
\left((\I_{(k-p)}-\tilde{L})^{-1}\varotimes \I_{|T'_p|}\right)\left(\tilde{U}\varotimes \I_{|T'_p|}\right)
\\=&
\left((\I_{(k-p)}-\tilde{L})^{-1} \tilde{U}\right)\varotimes \I_{|T'_p|}
\end{align*}
where $\tilde{L}$ is the lower triangular part of $(C^{(p)}-\I_{k-p})$ and $\tilde{U}=(C^{(p)}-\I_{k-p})-\tilde{L}$.
From
$B=\tilde{B}\varotimes \I_{|T'_p|}$, where $\tilde{B}=(\I_{(k-p)}-\tilde{L})^{-1} \tilde{U}$, it follows that the unique eigenvalues of $B$ are the same as the unique eigenvalues of $\tilde{B}$ and thus the largest modulus eigenvalue of $B$ equals the one of $\tilde{B}$.
The case $p=0$ is analogous (with no need to consider a sub-vector of $\delta^{(0)}\bbeta$ and $T'_0$ being equal to $T_0$ itself) and trivial to check.
\end{proof}

%%%%%
\subsection*{Proof of Theorem \ref{thm:ordering_CP}}
Theorem \ref{thm:factorization_general} shows that $(\delta^{(p)}\bgamma(s))_s$ for $p\in\{0,\dots,k-1\}$ are $k$ independent Markov chains.
Arguing as in the proof of Theorem \ref{thm:kronecker} above for each $p$ we consider $(\delta_{T'_p}^{(p)}\bgamma(s))_s$ rather than $(\delta^{(p)}\bgamma(s))_s$ to avoid working with singular Gaussian random vectors.
For any $p$ from 0 to $k-1$, \eqref{eq:tree_cond_indep} implies that the $Q$-matrix of $\delta_{T'_p}^{(p)}\bgamma$ is tridiagonal
(with $k-p$ blocks corresponding to $\delta_{T'_p}^{(p)}\bgamma^{(p)}$ up to $\delta_{T'_p}^{(p)}\bgamma^{(k-1)}$).
Therefore Theorem 5 of \citep{RobertsSahu1997} implies that the rate of convergence of the Gibbs Sampler in \eqref{eq:blocked_gibbs_reduced} is given by $\rho(\delta_{T'_p}^{(p)}\bgamma(s))=\lambda(A_{\delta_{T'_p}^{(p)}\bgamma(s)})^2$ where $A_{\delta_{T'_p}^{(p)}\bgamma(s)}$ is the $A$-matrix defined as in \eqref{eq:A_mat_defi} and $\lambda(A_{\delta_{T'_p}^{(p)}\bgamma(s)})$ is its largest eigenvalue.
Then, using \eqref{eq:A_matrix_delta} from the proof of Theorem \ref{thm:kronecker}, it follows $\lambda(A_{\delta_{T'_p}^{(p)}\bgamma(s)})^2=\lambda(C^{(p)}-\I_{k-p})^2$ and thus the rate of convergence of $(\delta^{(p)}\bgamma(s))_s$ is given by $\rho(\delta^{(p)}\bgamma(s))=\lambda(C^{(p)}-\I_{k-p})^2$. 
Noting that $C^{(p+1)}-\I_{k-(p+1)}$ is obtained from $C^{(p)}-\I_{k-p}$ by removing the first row and column, the desired inequality $\rho(\delta^{(p)}\bgamma(s))=\lambda(C^{(p)}-\I_{k-p})^2 \geq \lambda(C^{(p+1)}-\I_{k-(p+1)})^2=\rho(\delta^{(p+1)}\bgamma(s))$ follows by applying the Cauchy interlacing theorem (see e.g. \citealp{bhatia2013}), which states that the eigenvalues of a principal submatrix of a symmetric matrix interlace the original eigenvalues.

%%%%%
\subsection*{Proof of Theorems \ref{thm:PNCP_3_levels} and \ref{thm:PNCP_k_levels}}\label{proof:PNCP_k_levels}
Since Theorem \ref{thm:PNCP_3_levels} is a special case of Theorem \ref{thm:PNCP_k_levels} we only need to prove the latter.
Let $\Lambda$ be defined as in the statement of Theorem \ref{thm:PNCP_k_levels}.
The precision matrix of the hierarchical reparametrization $\Lambda\bbeta_T$ is given by $D=(\Lambda^T)^{-1}Q\Lambda^{-1}$ where $Q$ is the precision matrix of $\bbeta_T$.
To show that $\Lambda\bbeta_T$ is made of independent random variables it suffices to show that $D$ is diagonal. In particular we show that  $\Lambda^T D \Lambda=Q$ with $D$ being a diagonal matrix with elements defined by \eqref{eq:D_definition}.
Since $\Lambda^T D\Lambda$ is symmetric, we can consider $(\Lambda^T D\Lambda)_{rt}$ assuming $t\succeq r$ without loss of generality. Since $D$ is diagonal and $\Lambda=(\lambda_{ss'})_{s,s'\in T}$ is non-zero only for $s'\in\{s,pa(s)\}$ we have
\begin{align*}
(\Lambda^T D\Lambda)_{rt}
&=\sum_{r'\in T}\sum_{s\in T}\lambda_{sr}D_{sr'}\lambda_{r't} \\
&=\sum_{r'\in T}\lambda_{r'r} D_{r'r'}\lambda_{r't}
=D_{rr}\lambda_{rt}+\sum_{r'\in ch(r)}\lambda_{r'r} D_{r'r'}\lambda_{r't}\,.
\end{align*}
Therefore $(\Lambda^T D\Lambda)_{rt}=0$ unless $r=t$ or $t=pa(r)$.
In the latter case 
$(\Lambda^T D\Lambda)_{r\,pa(r)}=
D_{rr}\lambda_{rpa(r)}=
D_{rr}\frac{Q_{r\,pa(r)}}{D_{rr}}=Q_{r\,pa(r)}$.
In the former case (i.e.\ $r=t$), 
$
(\Lambda^T D\Lambda)_{rr}=
D_{rr}+\sum_{r'\in ch(r)}\lambda_{r'r}^2 D_{r'r'}=
Q_{rr},
$
where we used \eqref{eq:D_definition} for $D_{rr}$.
It follows $\Lambda^T D \Lambda=Q$ as desired.
%
%
%%%%%%
%%%%%%
%\section{Multigrid decomposition for $k$ levels under weaker symmetry assumptions and proof of Theorem \ref{thm:rates_3_uneven_CP}}\label{appendix:weaker_conditions}
%We now show how to develop a multigrid decomposition under weaker symmetry conditions than \eqref{thm:rates_3_uneven_CP} and use the results to prove the bound in Theorem \ref{thm:rates_3_uneven_CP}.
%Since it is enought to prove Theorem \ref{thm:rates_3_uneven_CP}, in this section we consider only the centred parametrization case (i.e.\ vectors $\bgamma_T$ satisfying \eqref{eq:tree_cond_indep}).
%%\fotnote{easy to extend to all parametrizations (not only centred one). However only centred is then used to provide upper bounds on convergence rates. Worth to write it in more general version?}
%The weaker symmetry requirement replacing \eqref{eq:k_symmetry} is the following: there exist a $k$-dimensional vector $(c_0,\dots,c_{k-1})$ such that the $A$-matrix $(A_{tr})_{t,r\in T}$ of the random vector $\bgamma_T$ under consideration satisfies
%\begin{align}\tag{S*}\label{eq:k_symmetry_appendix}
%\sum_{r\in ch(t)}A_{rt}A_{tr}
%\;=\;
%c_{\ell}
%\qquad\hbox{for all }
%t\in T_\ell\,.
%\end{align}
%Recalling that $A_{tr}=-\frac{Q_{tr}}{Q_{tt}}$ by definition it is easy to check that if a vector $\bgamma_T$ satisfies \eqref{eq:tree_cond_indep} and \eqref{eq:k_symmetry} it automatically satisfies \eqref{eq:k_symmetry_appendix}.
%the converse, however is not true.
%In this section we will prove the following theorem.
%\begin{theorem}\label{coroll:rate_k_levels_CP_appendix}
%The rate of convergence of a DSGS targeting a Gaussian distribution satisfying \eqref{eq:tree_cond_indep} and \eqref{eq:k_symmetry_appendix} is given by the largest squared eigenvalue of the $k$-dimensional matrix
%\begin{equation*}
%\begin{pmatrix}0 & 1\\
%c_0 & 0 & 1\\
% & \dots & \dots & \dots\\
% &  & c_{k-2} & 0 & 1\\
% &  &  & c_{k-1} & 0
%\end{pmatrix}
%\end{equation*}
%where the vector $(c_0,\dots,c_{k-1})$ is given by \eqref{eq:k_symmetry_appendix}.
%\end{theorem}
%To prove Theorem \ref{coroll:rate_k_levels_CP_appendix} we will employ a multigrid decomposition approach analogous to Section \ref{sec:multigrid_k}.
%Unlike Section \ref{sec:multigrid_k}, in this case the auxiliary Markov chain $\textbf{X}=(X_\ell)_{\ell=0}^{k-1}$ going through the tree $T$ evolves in a non-symmetric manner.
%In particular $X_0=t_0$ almost surely and for $\ell\in\{1,\dots,k-1\}$
%\begin{equation}\label{eq:non_symm_markov}
%P(X_{\ell+1}=t\,|\,X_{\ell}=r)=\frac{A_{tr}A_{rt}}{\sum_{t'\in ch(r)}A_{t'r}A_{rt'}}\1(t\in ch(r))\,.
%\end{equation}
%Equation \eqref{eq:non_symm_markov} implies that at each step the Markov chain jumps from the current state $r$ to one of its children $t\in ch(r)$ choosing $t$ proportionally to the squared correlation between $r$ and $t$, namely $A_{tr}A_{rt}=\frac{Q_{tr}^2}{Q_{tt}Q_{rr}}$.
%\fotnote{Actually this is the inverse correlation, isn't it?}
%To obtain the correct multigrid decomposition of $\bgamma$ we first need to define a zero-mean and rescaled version, defining
%\begin{equation}\label{eq:rescale}
%\tgamma_t=\sqrt{\frac{Q_{tt}}{P(X_{\ell(t)}=t)}}\left(\gamma_t-\E[\gamma_t]\right)\,.
%\end{equation}
%\fotnote{We could probably avoid the centering}
%Note that the Gibbs Sampler is invariant to translation and rescaling, so the Gibbs Sampler with parametrization $\btgamma=(\tgamma_t)_{t\in T}$ is equivalent to the one with parametrization  $\bgamma=(\gamma_t)_{t\in T}$.
%Given \eqref{eq:non_symm_markov} and \eqref{eq:rescale} we define the operators $\phi$ and $\delta$ as in Section \ref{sec:multigrid_k}.
%In particular for any $p\leq d$ the averaging operator $\phi^{(p)}:\R^{T_d}\rightarrow\R^{T_p}$ is defined as
%\begin{align}\label{eq:phi_defi_tgamma}
%\phi^{(p)}_r\btgamma^{(d)}
%=&
%\E[\tgamma_{X_d}|X_p=r]
%%=&
%%\sum_{t\in T_d}
%%\tgamma_{t}P(X_d=t|X_p=r)
%&r\in T_p
%\end{align}
%where $\textbf{X}=(X_\ell)_{\ell=0}^{k-1}$ is defined by \eqref{eq:non_symm_markov}.
%Given $\phi$, we define the residual operators
%$
%\delta^{(p)}:\R^{T_d}
%\rightarrow
%\times_{t\in T_{p-1}}(\R^{n_t})^*\subseteq \R^{T_p}
%$
%as
%$\delta^{(p)}=(\delta^{(p)}_r)_{r\in T_p}$
%with $\delta^{(p)}_r:\R^{T_d}
%\rightarrow\R$ defined as
%\begin{align}\label{eq:delta_defi}
%\delta^{(p)}_r\btgamma^{(d)}=&
%\phi^{(p)}_r\btgamma^{(d)}-\phi^{(p-1)}_{pa(r)}\btgamma^{(d)}
%&r\in T_p
%\end{align}
%for $1\leq p\leq d\leq k-1$
%and
%$
%\delta^{(0)}\btgamma^{(d)}=\phi^{(0)}\btgamma^{(d)}
%$
%for $0=p\leq d\leq k-1$.
%
%\fotnote{ATTENTION: All the proof following are exactly identical to the other versions. The only thing needed is $\tilde{Q}_{tr}=c_{\ell(t)\ell(r)}P(t\vee r)$. Should remove the double version!}
%
%
%\vspace{10mm}
%
%
%%%%%
%\begin{theorem}[Multigrid decomposition under weaker symmetry]\label{thm:factorization_general_appendix}
%Let $(\bgamma(s))_{s\in\N}$ be a Markov chain evolving according to \ref{sampler:GSkCP} with $\bgamma_T$ satisfying \eqref{eq:tree_cond_indep} and \eqref{eq:k_symmetry_appendix}.
%Then $(\delta^{(0)}\btgamma(s))_s$, \dots, $(\delta^{(k)}\btgamma(s))_s$ are $k$ independent Markov chains. 
%Moreover, each chain $\delta^{(p)}\btgamma(s)=(\delta^{(p)}\btgamma^{(d)}(s))_{d=p}^{k-1}$ evolves according to the following blocked Gibbs sampler scheme with $(k-p+1)$ blocks:
%for $d$ going from $p$ to $k$ sample
%\begin{align*}
%\delta^{(p)}\btgamma^{(d)}(s+1)
%\;\sim\;
%\mathcal{L}\left(\delta^{(p)}\bgamma^{(d)}|
%(\delta^{(p)}\btgamma^{(\ell)}(s+1))_{p\leq \ell<d},(\delta^{(p)}\btgamma^{(\ell)}(s))_{d<\ell\leq k-1}
%\right)\,.
%\end{align*}
%\end{theorem}
%%%%%
%Similarly to Section \ref{sec:multigrid_k}, to prove Theorem \ref{thm:factorization_general_appendix} we need the following lemma.
%%%%%
%\begin{lemma}[$p$-residuals interact only with $p$-residuals]\label{lemma:delta_expansion_appendix}
%Let $\bgamma$ be a Gaussian random vector satisfying \eqref{eq:tree_cond_indep} and \eqref{eq:k_symmetry_appendix}.
%Then for any $p$ and $d$ with $1\leq p\leq d\leq k-1$ itholds
%\begin{align*}
%\E[\delta^{(p)}_t\bgamma^{(d)}|\bgamma\backslash\bgamma^{(d)}]
%&=
%c_{d-1}
%%\frac{Z_{\ell-1}}{Z_\ell}
%\,\delta^{(p)}_t\bgamma^{(d-1)}
%+
%%\frac{Z_{\ell+1}}{Z_\ell}\,
%\delta^{(p)}_t\bgamma^{(d+1)}
%&
%t\in T_p
%\,,
%\end{align*}
%with the $(d-1)$th and $(d+1)$th term to be removed in case $d$ equals 0 or $k$ respectively. 
%\end{lemma}
%%%%%
%\begin{proof}
%Suppose $\bgamma$ has zero mean (otherwise replace $\bgamma$ by $\bgamma-\mathbb{E}[\bgamma]$).
%Using \eqref{eq:tree_cond_indep}, for any $d\in\{1,\dots,k-1\}$ and $t\in T_d$, we can write the full conditional expectation of $\gamma_t$ as
%\begin{align*}
%\E[\gamma_t|\bgamma_{T\backslash t}]
%&=
%A_{tpa(t)}
%\gamma_{pa(t)}
%+
%\sum_{t'\in ch(t)}
%A_{tt'}
%\gamma_r\,.
%\end{align*}
%Using the definition of $\phi^{(d)} $ in \eqref{eq:phi_defi_appendix} and the fact that the last equation does not depend on $\bgamma^{(d)}$ we have
%\begin{align*}
%\E[w^{(d)}_t\gamma_t|\bgamma\backslash\bgamma^{(d)}]
%&
%=
%A_{tpa(t)}
%A_{pa(t)t}
%w^{(d-1)}_{pa(t)}\gamma_{pa(t)}
%+
%\sum_{t'\in ch(t)}
%w^{(d+1)}_{t'}
%\gamma_{t'}\,.
%\end{align*}
%Again using \eqref{eq:phi_defi_appendix} and \eqref{eq:k_symmetry_appendix} we have, for any $r\in T_p$
%\begin{align*}%\label{eq:phi_expansion}
%&\E[\phi^{(p)}_r\bgamma^{(d)}|\bgamma\backslash\bgamma^{(d)}]
%\\
%&=
%\sum_{t\in T_d\,:\,t\succeq r}
%\E[w^{(d)}_t\gamma_{t}|\bgamma\backslash\bgamma^{(d)}]
%\\
%&=
%\sum_{t\in T_d\,:\,t\succeq r}
%\left(
%A_{tpa(t)}
%A_{pa(t)t}
%w^{(d-1)}_{pa(t)}\gamma_{pa(t)}
%+
%\sum_{t'\in ch(t)}
%w^{(d+1)}_{t'}
%\gamma_{t'}
%\right)
%\\
%&=
%\sum_{t'\in T_{d-1}\,:\,t'\succeq r}
%w^{(d-1)}_{t'}\gamma_{t'}
%\left(
%\sum_{t\in ch(t')}
%A_{t\,t'}
%A_{t't}
%\right) 
%+
%\sum_{t'\in T_{d+1}\,:\,t'\succeq r}
%w^{(d+1)}_{t'}
%\gamma_{t'}
%%\\
%%&=
%%\sum_{t\in T_{d-1}\,:\,t\succeq r}
%%w^{(d-1)}_{t}\gamma_{t}
%%\left(
%%\sum_{r\in ch(t)}
%%A_{r\,pa(r)}
%%A_{pa(r)r}
%%\right) 
%%+
%%\sum_{t\in T_{d+1}\,:\,t\succeq r}
%%w^{(d+1)}_{t}
%%\gamma_t
%\\
%&=
%c_{d-1}\phi^{(p)}_r\bgamma^{(d-1)}
%+
%\phi^{(p)}_r\bgamma^{(d+1)}
%\,.
%\end{align*}
%From the latter equation and \eqref{eq:delta_defi} it follows
%\begin{align*}
%\E[\delta_r^{(p)}\bgamma^{(d)}|\bgamma\backslash\bgamma^{(d)}]
%=&
%\E[\phi^{(p)}_r\bgamma^{(d)}|\bgamma\backslash\bgamma^{(d)}]-\E[\phi^{(p-1)}_{pa(r)}\bgamma^{(d)}|\bgamma\backslash\bgamma^{(d)}]
%\\
%=&
%c_{d-1}\phi^{(p)}_r\bgamma^{(d-1)}
%+
%\phi^{(p)}_r\bgamma^{(d+1)}-
%\\
%&
%c_{d-1}\phi^{(p-1)}_{pa(r)}\bgamma^{(d-1)}
%-
%\phi^{(p-1)}_{pa(r)}\bgamma^{(d+1)}
%\\
%=&
%c_{d-1}\delta^{(p)}_r\bgamma^{(d-1)}
%+
%\delta^{(p)}_{r}\bgamma^{(d+1)}
%\,.
%\end{align*}
%\end{proof}
%%%%%
%\begin{lemma}\label{lemma:residual_independence_gamma}
%For any $d\in\{0,\dots,k-1\}$ and $p,p'\in\{0,\dots,d\}$ with $p\neq p'$ itholds the following conditional independence
%\begin{align*}
%\delta^{(p)}\bbeta^{(d)}
%\bot
%\delta^{(p')}\bbeta^{(d)}
%\,|\,
%\bbeta\backslash\bbeta^{(d)}
%\,.
%\end{align*}
%\end{lemma}
%\begin{proof}
%\fotnote{THE PROOF NEEDS TO BE UPDATED!!}
%Let $d\in\{0,\dots,k-1\}$ and $p,p'\in\{0,\dots,d\}$ with $p\neq p'$, say $p< p'$.
%To make the notation more compact we will denote $\E[\cdot|\bbeta\backslash\bbeta^{(d)}]$ by $\tE[\cdot]$, $\phi^{(p)}_{r}\bbeta^{(d)}$ by $\tilde{\phi}^{(p)}_{r}$ and $P(X_{p'}=r'|X_p=r)$ by $P(r'|r)$ for all $r\in T_p$ and $r'\in T_{p'}$.
%We can suppose without loss of generality that $\tE[\beta_t]=0$ for all $t\in T_d$ by replacing $\beta_t$ with $\beta_t-\tE[\beta_t]$.
%Using the definition of $\phi^{(p)}_r\bbeta^{(d)}$ in \eqref{eq:phi_defi} we have
%\begin{align}
%\tE[\tilde{\phi}^{(p)}_r\tilde{\phi}^{(p')}_{r'}]
%=
%\sum_{s\in T_{p'}}
%P(s|r)
%\tE[\tilde{\phi}^{(p')}_s\tilde{\phi}^{(p')}_{r'}]
%=
%P(r'|r)
%\tE[(\tilde{\phi}^{(p')}_{r'})^2]\,,\label{eq:phi_dot_product_gamma}
%\end{align}
%where we used $\tE[\tilde{\phi}^{(p')}_s\tilde{\phi}^{(p')}_{r'}]=0$  for $r'\neq s$ and $r',s\in T_{p'}$, which follows from the conditional independence of $(\beta_t)_{t\in T_d\,,\,t\succeq r'}$ and $(\beta_t)_{t\in T_d\,,\,t\succeq s}$ given $\bbeta\backslash\bbeta^{(d)}$.
%Note that $P(r'|r)$ in \eqref{eq:phi_dot_product} could be 0.
%Then, using \eqref{eq:k_symmetry} we have
%\begin{align}
%\tE[(\tilde{\phi}^{(p')}_{r'})^2]
%=&
%\tE[(\sum_{t\in T_d\,:\,t\succeq r'}P(t|r')\beta_t)^2]
%=
%\sum_{t\in T_d\,:\,t\succeq r'}P(t|r')^2\tE[\beta_t^2]
%\nonumber\\
%=&
%\sum_{t\in T_d\,:\,t\succeq r'}\frac{P(t|r')^2}{c_{dd}P(t)}
%=
%\frac{1}{c_{dd}P(r')}\sum_{t\in T_d\,:\,t\succeq r'}P(t|r')
%=
%\frac{1}{c_{dd}P(r')}\,.
%\label{eq:phi_squared_gamma}
%\end{align}
%Combining \eqref{eq:phi_dot_product} and \eqref{eq:phi_squared} we have 
%$\tE[\tilde{\phi}^{(p)}_r\tilde{\phi}^{(p')}_{r'}]=0$ if $r'\nsucc r$ and
%\begin{align}\label{eq:zero_simplification_gamma}
%\tE[\tilde{\phi}^{(p)}_r\tilde{\phi}^{(p')}_{r'}]
%=&
%\frac{P(r'|r)}{c_{dd}P(r')}
%=
%\frac{1}{c_{dd}P(r)}
%&
%\hbox{if }r'\succ r\,.
%\end{align}
%From the last equality and the definition of $\delta^{(p)}_r\bbeta^{(d)}$ in \eqref{eq:delta_defi} we have
%\begin{align*}
%\tE&[\delta^{(p)}_r\bbeta^{(d)}\delta^{(p')}_{r'}\bbeta^{(d)}]
%=\\&
%\tE[(\tilde{\phi}^{(p)}_r-\tilde{\phi}^{(p-1)}_{pa(r)})(\tilde{\phi}^{(p')}_{r'}-\tilde{\phi}^{(p'-1)}_{pa(r')})]
%=\\&
%\tE[
%\tilde{\phi}^{(p)}_r\tilde{\phi}^{(p')}_{r'}
%-\tilde{\phi}^{(p-1)}_{pa(r)}\tilde{\phi}^{(p')}_{r'}
%-\tilde{\phi}^{(p)}_r\tilde{\phi}^{(p'-1)}_{pa(r')}
%+\tilde{\phi}^{(p-1)}_{pa(r)}\tilde{\phi}^{(p'-1)}_{pa(r')}
%]
%=\\&
%\frac{1}{c_{dd}}\left(
%\frac{P(r'|r)}{P(r')}
%-
%\frac{P(r'|pa(r))}{P(r')}
%-
%\frac{P(pa(r')|r)}{P(pa(r'))}
%+
%\frac{P(pa(r')|pa(r))}{P(pa(r'))}
%\right)
%=0,
%\end{align*}
%where the last equality is trivial if $pa(r)\nprec r'$ and can be deduced from \eqref{eq:zero_simplification} otherwise.
%The desired conditional independence follows from 
%$\E[\delta^{(p)}_r\bbeta^{(d)}\delta^{(p')}_{r'}\bbeta^{(d)}|\bbeta\backslash\bbeta^{(d)}]=0$
%for all $r\in T_p$ and $r'\in T_{p'}$.
%\end{proof}
%
%%%%%
%\begin{proof}[Proof or Theorem \ref{thm:factorization_general_appendix}]\ 
%Follows from Lemma \ref{lemma:delta_expansion_appendix} and Lemma \ref{lemma:residual_independence_gamma} in an analogous way of the proof of Theorem \ref{thm:factorization_general}.
%\end{proof}
%%%%%
%
%We can deduce a statement analogous to Theorem \ref{thm:kronecker}.
%%%%%%
%\begin{lemma}\label{lemma:kronecker_appendix}
%From Lemma \ref{lemma:delta_expansion_appendix} it follows that blocked Gibbs Sampler updating $\delta^{(p)}\bbeta=(\delta^{(p)}\bbeta^{(p)},\dots,\delta^{(p)}\bbeta^{(k)})\in \R^{(k-p+1)I^{(p)}}$ has $A$-matrix given by
%$$
%A_{\delta^{(p)}\bbeta}=\tilde{C}_{(p:k),(p:k)}\varotimes \I_{I^{(p)}}\,,
%$$
%where $\varotimes$ denotes the Kronecker product between matrices and $\tilde{C}=(\tilde{c}_{d\ell})_{d,\ell=0}^{k-1}$ is a $k\times k$ square matrix defined by 
%\dots%$\tilde{c}_{d\ell}=\frac{c_{d\ell}}{c_{dd}}$ for $d\neq\ell$ and $\tilde{c}_{dd}=0$.
%It follows that
%$\lambda(A_{\delta^{(p)}\bbeta})=\lambda(\tilde{C}_{(p:k),(p:k)})$.
%\end{lemma}
%%%%%%
%\begin{proof}\ 
%[[Follows from Lemma \ref{lemma:delta_expansion_appendix} and maybe some reference on Kronecker product]]
%\end{proof}
%
%\begin{theorem}\label{thm:ordering_CP_appendix}\emph{(Ordering of rates for centred parametrization)}
%Let $\bbeta$ be a Gaussian random vector satisfying \eqref{eq:tree_cond_indep} and \eqref{eq:k_symmetry_appendix}. 
%Let $(\bbeta(s))_{t\in\N}$ be the corresponding Markov chain evolving according to the Deterministic Scan Gibbs sampler (DSGS) scheme (e.g.,\ see Theorem \ref{thm:factorization_general}).
%Then the convergence rates of the $k+1$ independent Markov chains $(\delta^{(0)}\bbeta(s))_t$, \dots, $(\delta^{(k)}\bbeta(s))_t$ satisfy
%\begin{equation}
%\rho(\delta^{(0)}\bbeta(s))\geq
%\rho(\delta^{(1)}\bbeta(s))\geq
%\cdots\geq
%\rho(\delta^{(k-1)}\bbeta(s))\geq
%\rho(\delta^{(k)}\bbeta(s))=0\,.
%\end{equation}
%\end{theorem}
%\begin{proof}
%Follows from Theorem \ref{thm:factorization_general_appendix} and Lemma \ref{lemma:kronecker_appendix} in an analogous way to the proof of Theorem \ref{thm:ordering_CP}.
%\end{proof}
%
%Theorem \ref{thm:ordering_CP_appendix} implies the following corollary.
%\begin{coroll}\label{coroll:rate_k_levels_CP_appendix}
%The rate of convergence of a DSGS targeting a Gaussian distribution satisfying \eqref{eq:tree_cond_indep} and \eqref{eq:k_symmetry_appendix} is given by the largest squared eigenvalue of the $k$-dimensional matrix
%\begin{equation*}
%\begin{pmatrix}0 & 1\\
%c_0 & 0 & 1\\
% & \dots & \dots & \dots\\
% &  & c_{k-2} & 0 & 1\\
% &  &  & c_{k-1} & 0
%\end{pmatrix}
%\end{equation*}
%where the vector $(c_0,\dots,c_{k-1})$ is given by \eqref{eq:k_symmetry_appendix}.
%\end{coroll}


\end{document}


notation_Gibbs.tex
% DEFINING NOTATION AND SYMBOL FOR THE 
% ``Analysis of Gibbs Sampler for Gaussian hierarchical models''
%  .tex file

\newcommand{\sderr}{\sigma_e}
\newcommand{\varerr}{\sigma_e^2}
\newcommand{\sda}{\sigma_a}
\newcommand{\vara}{\sigma_a^2}
\newcommand{\sdb}{\sigma_b}
\newcommand{\varb}{\sigma_b^2}
\newcommand{\tvara}{\tilde{\sigma}_a^2}
\newcommand{\tvarb}{\tilde{\sigma}_b^2}
\newcommand{\tvarerr}{\tilde{\sigma}_e^2}
\newcommand{\avec}{\textbf{a}}
\newcommand{\bvec}{\textbf{b}}
\newcommand{\yvec}{\textbf{y}}
\newcommand{\amean}{a_{\cdot}}
\newcommand{\bmean}{b_{\cdot\cdot}}
\newcommand{\bimean}{{b}_{i\cdot}}
\newcommand{\ymean}{y_{\cdot\cdot\cdot}}
\newcommand{\yimean}{y_{i\cdot\cdot}}
\newcommand{\yijmean}{y_{ij\cdot}}
\newcommand{\bivec}{\textbf{b}_{1:I,\cdot}}
\renewcommand{\H}{\mathbb{H}}
\newcommand{\trho}{\tilde{\rho}}
\newcommand{\tc}{\tilde{c}}

\newcommand{\PNCP}{P_{NCP}}
\newcommand{\PCP}{P_{CP}}

\newcommand{\preca}{\tau_a}
\newcommand{\precb}{\tau_b}
\newcommand{\precerr}{\tau_e}
\newcommand{\tpreca}{\tilde{\tau}_a}
\newcommand{\tprecb}{\tilde{\tau}_b}
\newcommand{\tprecerr}{\tilde{\tau}_e}


\newcommand{\gavec}{\boldsymbol{\gamma}}
\newcommand{\etavec}{\boldsymbol{\eta}}
\newcommand{\gamean}{\gamma_{\cdot}}
\newcommand{\etamean}{\eta_{\cdot\cdot}}
\newcommand{\etai}{\eta_{i\cdot}}
\newcommand{\etaivec}{\underline{\bar{\eta}}_{i\cdot}}

\newcommand{\ratio}[2]{r_{#1,#2 }}
\newcommand{\var}[1]{\sigma_{#1}^2}
\newcommand{\tvar}[1]{\tilde{\sigma}_{#1}^2}
\newcommand{\precision}[1]{\tau_{#1}}
\newcommand{\tprec}[1]{\tilde{\tau}_{#1}}

\newcommand{\ttau}{\tilde{\tau}}
\newcommand{\tE}{\tilde{\E}}

\newcommand{\balpha}{{\boldsymbol{\alpha}}}
\newcommand{\btheta}{{\boldsymbol{\theta}}}
\newcommand{\bbeta}{{\boldsymbol{\beta}}}
\newcommand{\bgamma}{{\boldsymbol{\gamma}}}
\newcommand{\btgamma}{{\boldsymbol{\tilde{\gamma}}}}
\newcommand{\tgamma}{{\tilde{\gamma}}}

\newcommand{\mbeta}{{\bar{\beta}}}
\newcommand{\mbbeta}{{\bar{\bbeta}}}
\newcommand{\tbeta}{\tilde{\beta}}
\newcommand{\tbbeta}{\tilde{\bbeta}}
Gibbs_Sampler_for_Gaussian_hierarchal_models.bbl
\begin{thebibliography}{26}
\providecommand{\natexlab}[1]{#1}
\providecommand{\url}[1]{\texttt{#1}}
\expandafter\ifx\csname urlstyle\endcsname\relax
  \providecommand{\doi}[1]{doi: #1}\else
  \providecommand{\doi}{doi: \begingroup \urlstyle{rm}\Url}\fi

\bibitem[Amit(1991)]{.amit:1991:stochastic+convergence+distributions}
Yali Amit.
\newblock On rates of convergence of stochastic relaxation for {G}aussian and
  non-{G}aussian distributions.
\newblock \emph{J. Multivariate Anal.}, 38\penalty0 (1):\penalty0 82--99, 1991.

\bibitem[Amit(1996)]{.amit:1996:properties+convergence+perturbations}
Yali Amit.
\newblock Convergence properties of the {G}ibbs sampler for perturbations of
  {G}aussians.
\newblock \emph{Ann. Statist.}, 24\penalty0 (1):\penalty0 122--140, 1996.

\bibitem[Bass and Sahu(2016{\natexlab{a}})]{BassSahu2016}
Mark~R Bass and Sujit~K Sahu.
\newblock A comparison of centring parameterisations of gaussian process-based
  models for bayesian computation using mcmc.
\newblock \emph{Statistics and Computing}, pages 1--22, 2016{\natexlab{a}}.

\bibitem[Bass and Sahu(2016{\natexlab{b}})]{bass2016comparison}
Mark~R Bass and Sujit~K Sahu.
\newblock A comparison of centring parameterisations of gaussian process-based
  models for bayesian computation using mcmc.
\newblock \emph{Statistics and Computing}, pages 1--22, 2016{\natexlab{b}}.

\bibitem[Bhatia(2013)]{bhatia2013}
Rajendra Bhatia.
\newblock \emph{Matrix analysis}, volume 169.
\newblock Springer Science \& Business Media, 2013.

\bibitem[Chou et~al.(1994)Chou, Willsky, and Benveniste]{Chou1994}
Kenneth Chien-ko Chou, Alan~S Willsky, and Albert Benveniste.
\newblock Multiscale recursive estimation, data fusion, and regularization.
\newblock \emph{IEEE Transactions on Automatic Control}, 39\penalty0
  (3):\penalty0 464--478, 1994.

\bibitem[Diaconis et~al.(2010)Diaconis, Khare, and Saloff-Coste]{diaconis2010}
Persi Diaconis, Kshitij Khare, and Laurent Saloff-Coste.
\newblock Stochastic alternating projections.
\newblock \emph{Illinois J. Math.}, 54\penalty0 (3):\penalty0 963--979, 2010.

\bibitem[Gao and Owen(2016)]{GaoOwen2016}
Katelyn Gao and Art~B Owen.
\newblock Efficient moment calculations for variance components in large
  unbalanced crossed random effects models.
\newblock \emph{arXiv preprint arXiv:1602.00346}, 2016.

\bibitem[Gelfand and Smith(1990)]{gelsmi90}
Alan~E Gelfand and Adrian~FM Smith.
\newblock Sampling-based approaches to calculating marginal densities.
\newblock \emph{Journal of the American statistical association}, 85\penalty0
  (410):\penalty0 398--409, 1990.

\bibitem[Gelfand et~al.(1995)Gelfand, Sahu, and Carlin]{GelfandSahuCarlin1995}
Alan~E Gelfand, Sujit~K Sahu, and Bradley~P Carlin.
\newblock Efficient parametrisations for normal linear mixed models.
\newblock \emph{Biometrika}, 82\penalty0 (3):\penalty0 479--488, 1995.

\bibitem[Gelman et~al.(2013)Gelman, Carlin, Stern, Dunson, Vehtari, and
  Rubin]{gelman2013bayesian}
Andrew Gelman, John~B Carlin, Hal~S Stern, David~B Dunson, Aki Vehtari, and
  Donald~B Rubin.
\newblock Bayesian data analysis, 3rd edition, 2013.

\bibitem[Goodman and Sokal(1989)]{GoodmanSokal1989}
Jonathan Goodman and Alan~D Sokal.
\newblock Multigrid monte carlo method. conceptual foundations.
\newblock \emph{Physical Review D}, 40\penalty0 (6):\penalty0 2035, 1989.

\bibitem[Hills and Smith(1992)]{hills1992parameterization}
Susan~E Hills and Adrian~FM Smith.
\newblock Parameterization issues in bayesian inference.
\newblock \emph{Bayesian statistics}, 4:\penalty0 227--246, 1992.

\bibitem[Johndrow et~al.(2015)Johndrow, Mattingly, Mukherjee, and
  Dunson]{johndrow2015approximations}
James~E Johndrow, Jonathan~C Mattingly, Sayan Mukherjee, and David Dunson.
\newblock {Approximations of Markov Chains and Bayesian Inference}.
\newblock \emph{arXiv preprint arXiv:1508.03387}, 2015.

\bibitem[Jones and Hobert(2004)]{jones2004sufficient}
Galin~L Jones and James~P Hobert.
\newblock Sufficient burn-in for gibbs samplers for a hierarchical random
  effects model.
\newblock \emph{The Annals of Statistics}, 32\penalty0 (2):\penalty0 784--817,
  2004.

\bibitem[Liu and Sabatti(2000)]{LiuSabatti2000}
Jun~S Liu and Chiara Sabatti.
\newblock {Generalised Gibbs sampler and multigrid Monte Carlo for Bayesian
  computation}.
\newblock \emph{Biometrika}, pages 353--369, 2000.

\bibitem[Meng and Van~Dyk(1997)]{mengvand}
Xiao-Li Meng and David Van~Dyk.
\newblock The em algorithm -- an old folk-song sung to a fast new tune.
\newblock \emph{Journal of the Royal Statistical Society: Series B (Statistical
  Methodology)}, 59\penalty0 (3):\penalty0 511--567, 1997.

\bibitem[Papaspiliopoulos et~al.(2003)Papaspiliopoulos, Roberts, and
  Skold]{RobertsPapaspiliopoulosSkold2003}
O~Papaspiliopoulos, Gareth~O Roberts, and M~Skold.
\newblock Non-centered parameterizations for hierarchical models and data
  augmentation (with discussion).
\newblock In JM~Bernardo, MJ~Bayarri, JO~Berger, AP~Dawid, D~Heckerman, AFM
  Smith, and M~West, editors, \emph{Bayesian Statistics 7}, pages 307--326.
  Oxford University Press, New York, 2003.

\bibitem[Papaspiliopoulos et~al.(2007)Papaspiliopoulos, Roberts, and
  Sk{\"o}ld]{PapaspiliopoulosRobertsSkold2007}
Omiros Papaspiliopoulos, Gareth~O Roberts, and Martin Sk{\"o}ld.
\newblock A general framework for the parametrization of hierarchical models.
\newblock \emph{Statistical Science}, pages 59--73, 2007.

\bibitem[Rajaratnam and Sparks(2015)]{RajaratnamSparks2015}
Bala Rajaratnam and Doug Sparks.
\newblock {MCMC-based inference in the era of big data: A fundamental analysis
  of the convergence complexity of high-dimensional chains}.
\newblock \emph{arXiv preprint arXiv:1508.00947}, 2015.

\bibitem[Roberts and Rosenthal(2016)]{roberts2016complexity}
Gareth~O Roberts and Jeffrey~S Rosenthal.
\newblock Complexity bounds for markov chain monte carlo algorithms via
  diffusion limits.
\newblock \emph{Journal of Applied Probability}, 53\penalty0 (02):\penalty0
  410--420, 2016.

\bibitem[Roberts and Sahu(1997)]{RobertsSahu1997}
Gareth~O Roberts and Sujit~K Sahu.
\newblock Updating schemes, correlation structure, blocking and
  parameterization for the gibbs sampler.
\newblock \emph{Journal of the Royal Statistical Society: Series B (Statistical
  Methodology)}, 59\penalty0 (2):\penalty0 291--317, 1997.

\bibitem[Sahu and Roberts(1999)]{sahu1999convergence}
Sujit~K Sahu and Gareth~O Roberts.
\newblock On convergence of the em algorithmand the gibbs sampler.
\newblock \emph{Statistics and Computing}, 9\penalty0 (1):\penalty0 55--64,
  1999.

\bibitem[Smith and Roberts(1993)]{smirob93}
Adrian~FM Smith and Gareth~O Roberts.
\newblock Bayesian computation via the gibbs sampler and related markov chain
  monte carlo methods.
\newblock \emph{Journal of the Royal Statistical Society. Series B
  (Methodological)}, pages 3--23, 1993.

\bibitem[Whittaker(1990)]{Whittaker1990}
Joe Whittaker.
\newblock Graphical models in applied multivariate analysis, 1990.

\bibitem[Yu and Meng(2011)]{mengyu}
Yaming Yu and Xiao-Li Meng.
\newblock {To center or not to center: That is not the question -- An
  Ancillarity-Sufficiency Interweaving Strategy (ASIS) for boosting MCMC
  efficiency}.
\newblock \emph{Journal of Computational and Graphical Statistics}, 20\penalty0
  (3):\penalty0 531--570, 2011.

\end{thebibliography}
