Multi-Level-REML-Avignon-Arxiv.tex
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Latex file for:
%
% ``Multi Level Restricted Maximum Likelihood Covariance 
%   Estimation and Kriging for Large Spatial Datasets''
% 
% 
% % 
%   by Julio Castrillon, Marc G. Genton, Rio Yokota
%   
%   
%   Last changed: 13 November 2013
% 
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
\documentclass[12pt]{article}
\usepackage{epsfig, epsf, graphicx, float, color}
\usepackage{pstricks, psfrag}
\usepackage{amssymb,amsmath,mathrsfs}
\usepackage{verbatim,enumerate}
\usepackage{setspace}
\usepackage{natbib}

\usepackage{tikz-qtree,tikz-qtree-compat}

\usepackage{tikz}
\usepackage{mathpazo}
\usetikzlibrary{positioning}
% \usetikzlibrary{calc,fadings,d

% \usepackage[figlist,nofighead]{endfloat}

% \usepackage[nomarkers,tablesonly]{endfloat}

\setlength{\oddsidemargin}{-0.125in} \setlength{\topmargin}{-0.5in}
\setlength{\textwidth}{6.5in} \setlength{\textheight}{9in}

\setlength{\textheight}{9in} \setlength{\textwidth}{6.5in}
\setlength{\topmargin}{-40pt} \setlength{\oddsidemargin}{0pt}
\setlength{\evensidemargin}{0pt}

\setlength{\textheight}{9.4in} \setlength{\textwidth}{6.8in}
\setlength{\topmargin}{-71pt} \setlength{\oddsidemargin}{0pt}
\setlength{\evensidemargin}{-6pt} \tolerance=500
%\input psfig.tex
\setlength{\topmargin}{-56pt} \setlength{\oddsidemargin}{-6pt}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\wt{\widetilde}
\def\diag{\hbox{diag}}
\def\wh{\widehat}
\def\AIC{\hbox{AIC}}
\def\BIC{\hbox{BIC}}
%- Makes the section title start with Appendix in the appendix environment
\newcommand{\Appendix}
{%\appendix
\def\thesection{Appendix~\Alph{section}}
%\def\thesubsection{\Alph{section}.\arabic{subsection}}
\def\thesubsection{A.\arabic{subsection}}
}
\def\diag{\hbox{diag}}
%\def\log{\hbox{log}}
\def\bias{\hbox{bias}}
\def\Siuu{\boldSigma_{i,uu}}
\def\ANNALS{{\it Annals of Statistics}}
\def\BIOK{{\it Biometrika}}
\def\whT{\widehat{\Theta}}
\def\STATMED{{\it Statistics in Medicine}}
\def\STATSCI{{\it Statistical Science}}
\def\JSPI{{\it Journal of Statistical Planning \&amp; Inference}}
\def\JRSSB{{\it Journal of the Royal Statistical Society, Series B}}
\def\BMCS{{\it Biometrics}}
\def\COMMS{{\it Communications in Statistics, Theory \& Methods}}
\def\JQT{{\it Journal of Quality Technology}}
\def\STIM{{\it Statistics in Medicine}}
\def\TECH{{\it Technometrics}}
\def\AJE{{\it American Journal of Epidemiology}}
\def\JASA{{\it Journal of the American Statistical Association}}
\def\CDA{{\it Computational Statistics \& Data Analysis}}
\def\JCGS{{\it Journal of Computational and Graphical Statistics}}
\def\JCB{{\it Journal of Computational Biology}}
\def\EXT{{\it Extremes}}
\def\BIOINF{{\it Bioinformatics}}
\def\JAMA{{\it Journal of the American Medical Association}}
\def\JNUTR{{\it Journal of Nutrition}}
\def\JCGS{{\it Journal of Computational and Graphical Statistics}}
\def\LETTERS{{\it Letters in Probability and Statistics}}
\def\JABES{{\it Journal of Agricultural and
                      Environmental Statistics}}
\def\JASA{{\it Journal of the American Statistical Association}}
\def\ANNALS{{\it Annals of Statistics}}
\def\JSPI{{\it Journal of Statistical Planning \& Inference}}
\def\TECH{{\it Technometrics}}
\def\BIOK{{\it Bio\-me\-tri\-ka}}
\def\JRSSB{{\it Journal of the Royal Statistical Society, Series B}}
\def\BMCS{{\it Biometrics}}
\def\COMMS{{\it Communications in Statistics, Series A}}
\def\JQT{{\it Journal of Quality Technology}}
\def\SCAN{{\it Scandinavian Journal of Statistics}}
\def\AJE{{\it American Journal of Epidemiology}}
\def\STIM{{\it Statistics in Medicine}}
\def\ANNALS{{\it Annals of Statistics}}
\def\whT{\widehat{\Theta}}
\def\STATMED{{\it Statistics in Medicine}}
\def\STATSCI{{\it Statistical Science}}
\def\JSPI{{\it Journal of Statistical Planning \& Inference}}
\def\JRSSB{{\it Journal of the Royal Statistical Society, Series B}}
\def\BMCS{{\it Biometrics}}
\def\COMMS{{\it Communications in Statistics, Theory \& Methods}}
\def\JQT{{\it Journal of Quality Technology}}
\def\STIM{{\it Statistics in Medicine}}
\def\TECH{{\it Technometrics}}
\def\AJE{{\it American Journal of Epidemiology}}
\def\JASA{{\it Journal of the American Statistical Association}}
\def\CDA{{\it Computational Statistics \& Data Analysis}}
\def\dfrac#1#2{{\displaystyle{#1\over#2}}}
\def\VS{{\vskip 3mm\noindent}}
\def\boxit#1{\vbox{\hrule\hbox{\vrule\kern6pt
          \vbox{\kern6pt#1\kern6pt}\kern6pt\vrule}\hrule}}
\def\refhg{\hangindent=20pt\hangafter=1}
\def\refmark{\par\vskip 2mm\noindent\refhg}
\def\naive{\hbox{naive}}
\def\itemitem{\par\indent \hangindent2\pahttprindent \textindent}
\def\var{\hbox{var}}
\def\cov{\hbox{cov}}
\def\corr{\hbox{corr}}
\def\trace{\hbox{trace}}
\def\refhg{\hangindent=20pt\hangafter=1}
\def\refmark{\par\vskip 2mm\noindent\refhg}
\def\Normal{\hbox{Normal}}
\def\povr{\buildrel p\over\longrightarrow}
\def\ccdot{{\bullet}}
\def\bse{\begin{eqnarray*}}
\def\ese{\end{eqnarray*}}
\def\be{\begin{eqnarray}}
\def\ee{\end{eqnarray}}
\def\bq{\begin{equation}}
\def\eq{\end{equation}}
\def\bse{\begin{eqnarray*}}
\def\ese{\end{eqnarray*}}
\def\pr{\hbox{pr}}
\def\CV{\hbox{CV}}
\def\wh{\widehat}
\def\T{^{\rm T}}
\def\myalpha{{\cal A}}
\def\th{^{th}}

% Color corrections in text
\newcommand{\corb}[1]{\textcolor{black}{#1}}
\newcommand{\corred}[1]{\textcolor{black}{#1}}
\newcommand{\corblue}[1]{\textcolor{black}{#1}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Marc Definitions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\renewcommand{\baselinestretch}{1.2} % Change this 1.5 or whatever
\newcommand{\qed}{\hfill\hfill\vbox{\hrule\hbox{\vrule\squarebox
   {.667em}\vrule}\hrule}\smallskip}
\newcommand{\bbR}{\mathbb{R}}
\newcommand{\bbX}{\mathbb{X}}
\newcommand{\bbN}{\mathbb{N}}
\newcommand{\bbE}{\mathbb{E}}
\newcommand{\bbS}{\mathbb{S}}
\newcommand{\bD}{\mathbf{D}}
\newcommand{\bI}{\mathbf{I}}
\newcommand{\bL}{\mathbf{L}}
\newcommand{\bG}{\mathbf{G}}
\newcommand{\bW}{\mathbf{W}}
\newcommand{\bP}{\mathbf{P}}
\newcommand{\bU}{\mathbf{U}}
\newcommand{\bC}{\mathbf{C}}
\newcommand{\bK}{\mathbf{K}}
\newcommand{\bM}{\mathbf{M}}
\newcommand{\bJ}{\mathbf{J}}
\newcommand{\bH}{\mathbf{H}}
\newcommand{\bS}{\mathbf{S}}
\newcommand{\bV}{\mathbf{V}}
\newcommand{\bX}{\mathbf{X}}
\newcommand{\bY}{\mathbf{Y}}
\newcommand{\bZ}{\mathbf{Z}}
\newcommand{\bh}{\mathbf{h}}
\newcommand{\bx}{\mathbf{x}}
\newcommand{\by}{\mathbf{y}}
\newcommand{\bv}{\mathbf{v}}
\newcommand{\bz}{\mathbf{z}}
\newcommand{\bs}{\mathbf{s}}
\newcommand{\ba}{\mathbf{a}}
\newcommand{\bb}{\mathbf{b}}
\newcommand{\bo}{\mathbf{o}}
\newcommand{\bc}{\mathbf{c}}
\newcommand{\bve}{\mathbf{e}}
\newcommand{\bu}{\mathbf{u}}
\newcommand{\bw}{\mathbf{w}}
\newcommand{\bm}{\mathbf{m}}
\newcommand{\bg}{\mathbf{g}}
\newcommand{\cS}{\cal{S}}
\newcommand{\bmu}{\boldsymbol{\mu}}
\newcommand{\bxi}{\boldsymbol{\xi}}
\newcommand{\bsigma}{\boldsymbol{\sigma}}
\newcommand{\bgamma}{\boldsymbol{\gamma}}
\newcommand{\btau}{\boldsymbol{\tau}}
\newcommand{\brho}{\boldsymbol{\rho}}
\newcommand{\blambda}{\boldsymbol{\lambda}}
\newcommand{\bdelta}{\mathbf{\delta}}
\newcommand{\btheta}{\boldsymbol{\theta}}
\newcommand{\bvartheta}{\boldsymbol{\vartheta}}
\newcommand{\bpsi}{\boldsymbol{\psi}}
\newcommand{\bphi}{\boldsymbol{\phi}}
\newcommand{\bepsilon}{\boldsymbol{\epsilon}}
\newcommand{\balpha}{\boldsymbol{\alpha}}
\newcommand{\bbeta}{\boldsymbol{\beta}}
\newcommand{\bSigma}{\boldsymbol{\Sigma}}
\newcommand{\bLambda}{\boldsymbol{\Lambda}}
\newcommand{\bOmega}{\boldsymbol{\Omega}}
\newcommand{\br}{\boldsymbol{r}}
\newcommand{\0}{\mathbf{0}}
\newcommand{\1}{\mathbf{1}}
\newcommand{\binfty}{\boldsymbol{\infty}}
\newcommand{\E}{\mbox{E}}
\def\sD{\mathcal{D}}
\def\sN{\mathcal{N}}
\def\R{\Bbb{R}}
\newcommand{\verteq}[0]{\begin{turn}{90} $=$\end{turn}}
%\newcommand{\Pr}{\mbox{Pr}}
% \renewcommand{\baselinestretch}{1.25}
\pagenumbering{arabic}

\newcommand*{\ackname}{Acknowledgements}
% Definitions, Lemmas, Propositions and Theorems:
%\theoremstyle{definition}
\newtheorem{defi}{Definition}
%\theoremstyle{plain}
\newtheorem{theo}{Theorem}
\newtheorem{prop}{Proposition}
\newtheorem{lem}{Lemma}
%\theoremstyle{remark}
\newtheorem{rem}{Remark}
\newtheorem{asum}{Assumption}
\newtheorem{exam}{Example}

% \makeatletter
% \let\orig@makecaption\@makecaption
% \AtBeginTables{\renewcommand{\@makecaption}[2]{#1}}
% \AtBeginFigures{\renewcommand{\@makecaption}[2]{#1}}
% \makeatother

% \underset{x}{\operatorname{argmax}} 

% \DeclareMathOperator*{\argmax}{\arg\!\max} 

\newcommand{\argmax}{\operatornamewithlimits{argmax}}
\newcommand{\argmin}{\operatornamewithlimits{argmin}}

\begin{document}
\thispagestyle{empty} \baselineskip=28pt \vskip 5mm
\begin{center} {\Huge{\bf Multi-Level Restricted Maximum Likelihood
      Covariance Estimation and Kriging for Large Non-Gridded
      Spatial Datasets}}
\end{center}

\baselineskip=12pt \vskip 10mm

\begin{center}\large
Julio E. Castrill\'on-Cand\'as$^{1}$, Marc G.~Genton$^{2}$, and Rio
Yokota$^{3}$ \makeatletter{ \renewcommand*{\@makefnmark}{}
  \footnotetext{SRI Center for Uncertainty Quantification in
    Computational Science and Engineering$^{1}$; Computer, Electrical
    and Mathematical Sciences and Engineering$^{2}$, King Abdullah
    University of Science and Technology, Thuwal 23955-6900, Saudi
    Arabia; Tokyo Institute of Technology Global Scientific and
    Computing Center$^{3}$, 2-12-1 i7-2 O-okayama Meguro-ku, 152-8550,
    Tokyo, Japan.  E-mails: uvel@alum.mit.edu,
    marc.genton@kaust.edu.sa, rioyokota@gsic.titech.ac.jp}
  \makeatother}

\bigskip

{\it To appear in Spatial Statistics, (CC BY-NC-ND),
  doi:10.1016/j.spasta.2015.10.006.




}
\end{center}

\baselineskip=17pt \vskip 10mm \centerline{\today} \vskip 10mm

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{center}
{\large{\bf Abstract}}
\end{center}

\baselineskip=12pt

\begin{doublespace}
  We develop a multi-level restricted Gaussian maximum likelihood
  method for estimating the covariance function parameters and
  computing the best unbiased predictor. Our approach produces a new
  set of multi-level contrasts where the deterministic parameters of
  the model are filtered out thus enabling the estimation of the
  covariance parameters to be decoupled from the deterministic
  component. Moreover, the multi-level covariance matrix of the
  contrasts exhibit fast decay that is dependent on the smoothness of
  the covariance function.  Due to the fast decay of the multi-level
  covariance matrix coefficients only a small set is computed with a
  level dependent criterion.  We demonstrate our approach on problems
  of up to 512,000 observations with a Mat\'{e}rn covariance function
  and highly irregular placements of the observations. In addition,
  these problems are numerically unstable and hard to solve with
  traditional methods.

% Our results show that we are able to solve the parameter estimation
% problem for 100,000 observation locations and the Kriging estimate for
% 512,000 observations with good accuracy ($10^{-3}$ relative error) for
% the Mat\'{e}rn covariance function.
\end{doublespace}

\par\vfill\noindent {\bf KEY WORDS:} Fast Multipole Method; Hierarchical Basis; High
Performance Computing; Sparsification of Covariance Matrices

\par\medskip\noindent {\bf Short title}:
Multi-Level Restricted Maximum Likelihood and Kriging


\clearpage\pagebreak\newpage \pagenumbering{arabic}
\baselineskip=25.9pt


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
\label{Introduction}

Consider the following model for a Gaussian spatial random field $Z$:
\begin{equation}
Z(\bs) = \bm(\bs)\T \bbeta+\epsilon(\bs), \qquad \bs \in \R^d,
\label{Introduction:noisemodel}
\end{equation}
where $\bm\in\R^p$ is a known function of the spatial location $\bs$,
$\bbeta\in\R^p$ is an unknown vector of coefficients, and
$\epsilon$ is a stationary mean zero Gaussian random field with
parametric covariance function
$C(\bs,\bs';\btheta)=\cov\{\epsilon(\bs),\epsilon(\bs')\}$
having an unknown vector $\btheta\in\R^w$ of parameters. We observe
the data vector $\bZ=(Z(\bs_1),\ldots,Z(\bs_n))\T$ at locations
$\mathbb{S}:=\{ \bs_{1},\dots,\bs_{n}\}$, where $\bs_{1} \neq
\bs_{2} \neq \bs_{3} \neq$ $\dots \neq \bs_{n-1} \neq \bs_{n}$, and
wish to: 1) estimate the unknown vectors $\bbeta$ and $\btheta$; and
2) predict $Z(\bs_0)$, where $\bs_0$ is a new spatial location. These
two tasks are particularly challenging when the sample size $n$ is
large.

To address the estimation part, let $\bC(\btheta)=\cov(\bZ,\bZ\T)\in
\R^{n \times n}$ be the covariance matrix of $\bZ$ and assume it is
nonsingular for all $\btheta\in\R^w$.  Define $\bM=\big( \bm(\bs_1)
\ldots \bm(\bs_n)\big)\T\in \R^{n\times p}$ and assume it is of full
rank, $p$. The model \eqref{Introduction:noisemodel} leads to the
vectorial formulation
\begin{equation}
{\bf Z} = {\bf M \bbeta}+{\boldsymbol \epsilon},
\label{Introduction:vectormodel}
\end{equation}
where $\boldsymbol \epsilon$ is a Gaussian random vector, ${\boldsymbol \epsilon}
\sim {\cal N}_{n}(\0,\bC(\btheta))$. Then the log-likelihood function
is
\begin{equation}
\ell(\bbeta,\btheta)=-\frac{n}{2}\log(2\pi)-\frac{1}{2}\log
\det\{\bC(\btheta)\}-\frac{1}{2}(\bZ-\bM\bbeta)\T\bC(\btheta)^{-1}(\bZ-\bM\bbeta),
\label{Introduction:loglikelihood}
\end{equation}
which can be profiled by generalized least squares with
\begin{equation}
\hat \bbeta(\btheta)=\{\bM\T \bC(\btheta)^{-1} \bM\}^{-1}\bM\T \bC(\btheta)^{-1}\bZ. \label{GLSbeta}
\end{equation}
A consequence of profiling is that the maximum likelihood estimator
(MLE) of $\btheta$ then tends to be biased.  A solution to this
problem is to use restricted maximum likelihood (REML) estimation
which consists in calculating the log-likelihood of $n-p$ linearly
independent contrasts, that is, linear combinations of observations
whose joint distribution does not depend on $\bbeta$, from the set
$\bY=\{\bI_n-\bM(\bM\T\bM)^{-1}\bM\T\}\bZ$.  In this paper, we propose
% a judicious choice of the functions in $\bm$ defining the matrix
% $\bM$, and
a new set of contrasts that lead to significant
computational benefits (with good accuracy) when computing the REML
estimator of $\btheta$ for large sample size~$n$.

%We propose a method that
%  produces a new set of contrasts with fast decay of the covariance
%  matrix regardless of the decay rate of the covariance
%  function (but subject to the smoothness).

To address the prediction part, consider the best unbiased
predictor $\hat Z(\bs_0)=\lambda_0+\blambda\T\bZ$ where
$\blambda=(\lambda_1,\ldots,\lambda_n)\T$. The unbiasedness constraint
implies $\lambda_0=0$ and $\bM\T\blambda=\bm(\bs_0)$.  The
minimization of the mean squared prediction error
E$[\{Z(\bs_0)-\blambda\T\bZ\}^2]$ under the constraint
$\bM\T\blambda=\bm(\bs_0)$ yields
\begin{equation}
\hat Z(\bs_0)=\bm(\bs_0)\T\hat \bbeta+\bc(\btheta)\T
\bC(\btheta)^{-1}(\bZ-\bM\hat \bbeta), \label{KrigBLUP}
\end{equation}
where $\bc(\btheta)=\cov\{\bZ,Z(\bs_0)\}\in \R^{n}$ and $\hat \bbeta$ is
defined in (\ref{GLSbeta}).  In this paper, we propose a new
transformation of the data vector $\bZ$ leading to a decoupled
multi-level description of the model \eqref{Introduction:noisemodel}
without any loss of structure. This multi-level representation leads
to significant computational benefits when computing the kriging
predictor $\hat Z(\bs_0)$ in (\ref{KrigBLUP}) for large sample size
$n$.

Previous work has been performed to maximize
\eqref{Introduction:loglikelihood}.  The classical technique is to
compute a Cholesky factorization of $\bC$. However, this requires
${\cal O}(n^2)$ memory and ${\cal O}(n^3)$ computational steps, thus
impractical for large scale problems.

Under special structures of the covariance matrix, i.e., fast decay of
the covariance function, a tapering technique can be used to sparsify
the covariance matrix and thus increase memory and computational
efficiency (\cite{Furrer2006,Kaufman2008}). These techniques
are good when applicable but tend to be restrictive. For a review of various
approaches to spatial statistics for large datasets, see \cite{Sun2012}.

Recently we have seen the advent of solving the optimization problem
\eqref{Introduction:loglikelihood} from a computational numerical
perspective.  \cite{Anitescu2012} developed a matrix-free approach for
computing the maximum of the log-likelihood
\eqref{Introduction:loglikelihood} based on a stochastic programming
reformulation.  This method relies on Monte Carlo approximation of the
derivative of the score function with respect to the covariance
parameters $\btheta$ to compute the maximization
\eqref{Introduction:loglikelihood}. The authors show promising results
for a grid geometry of the placement of the observations.  However for
a non-grid geometry the cost of computing the preconditioner becomes
${\cal O}(n^2)$ and it is not clear how many iterations for
convergence are needed as the geometry deviates from a
grid. \corb{Moreover, due to the slow convergence rate of the Monte
  Carlo method ($\eta^{-1/2}$ convergence rate where $\eta$ is the
  number of realizations) many samples might be required before a
  suitable estimate is obtained.} The previous work was extended in
\cite{Stein2013}.  Although the results are impressive (1,000,000 +
size problems), the approach is restricted to regular grid geometries
with partially occluded areas.

\cite{Stein2012} presented a difference filter preconditioning for
large covariance matrices not unlike our multi-level method. By
constructing a preconditioner based on the difference filter the
number of iterations of a Preconditioned Conjugate Gradient (PCG)
drops significantly. However, the authors can only construct a
preconditioner for irregularly placed observations in 1D and for a
regular grid in higher dimension. Moreover, the authors point out that
the restrictions on the spectral density of the random field $Z$
are strong.


In \cite{Stein2004} the authors proposed a REML method in combination
with an approximation of the likelihood. This approach uses a
truncation method to compute an approximation of the likelihood
function. It appears to be effective if the truncated terms have small
correlations. \corb{However, if the covariance function has a slow
  decay then we expect that this approximation will not be accurate
  unless a large neighborhood is incorporated.} Moreover, this paper
does not include an analysis of the error with respect to the
truncation.

In \cite{Sun2015} the authors proposed new unbiased estimating
equations based on score equation approximations. The inverse
covariance matrix is approximated with a sparse inverse Cholesky
decomposition.  \corb{As in \cite{Stein2004} the approximation is
  expected to be fast and accurate for locally correlated observations
  but will suffer from slow decay of the covariance function.}
Moreover, the results are limited to grid-like geometries.

In the next section we present the basic ideas behind our approach.
In Section \ref{MultiLevelREML} we show the construction of a
multi-level basis from the observations points. In Section
\ref{MultiLevelCovarianceMatrix} we describe how to efficiently construct
a multi-level covariance matrix that arises from the new basis. In
Section \ref{multilevelestimator} a multi-level estimator is proposed. In
Section \ref{multilevelkriging} the multi-level kriging approach is
described. In Section \ref{numericalstudy} hard to solve numerical
examples are provided and compared with traditional methods. In Section
\ref{multilevelconclusions} we give concluding remarks.  Proofs are
relegated to the Appendix A and a notation summary can be found in
Appendix B.  We also include computational and mathematical details in
the remarks. However, these may be skipped on a first reading except
for the more mathematically oriented reader.

\section{Multi-Level REML and Kriging Basic Approach}

We now present the main ideas of our proposal. Denote by ${\cal
  P}^{p}(\mathbb{S})$ the span of the columns of the design matrix
${\bf M}$. Let $\bL \in \R^{ p \times n} $ be an orthogonal projection
from $\R^n$ to ${\cal P}^{p}(\mathbb{S})$ and $\bW \in \R^{(n-p)
  \times n}$ be an orthogonal projection from $\R^n$ to ${\cal
  P}^{p}(\mathbb{S})^{\perp}$, the orthogonal complement of ${\cal
  P}^{p}(\mathbb{S})$.  Moreover we assume that the operator $\left[
\begin{array}{c}
\bW \\
\bL
\end{array}
\right ]$ is orthonormal.

By applying the operator ${\bf W}$ to
\eqref{Introduction:vectormodel} we obtain $\bZ_{W}={\bf W Z} = {\bf
  W} ({\bf M \bbeta}+ {\boldsymbol \epsilon}) = {\bf W{\boldsymbol
    \epsilon}}$. Our first observation is that the trend
contribution ${\bf M} \bbeta$ is filtered out from the data ${\bf
  Z}$. We can now formulate the estimation of the covariance
parameters $\btheta$ without the trend. The new log-likelihood
function becomes
\begin{equation}
\ell_{W}(\btheta)
=\corb{-\frac{n-p}{2}\log(2\pi)-\frac{1}{2}\log
\det\{\bC_{W}(\btheta)\}
-\frac{1}{2}\bZ_{W}\T\bC_{W}(\btheta)^{-1}\bZ_{W},}
\label{Introduction:multilevelloglikelihood}
\end{equation}
where $\bC_{W}(\btheta) = \bW \bC(\btheta) \bW \T$ and $\bZ_{W}\sim
{\cal N}_{n-p}(\0, \bW \bC(\btheta) \bW\T)$. As shown in Section
\ref{multilevelestimator} , to estimate the coefficients $\btheta$ it is
not necessary to compute $\ell_W(\btheta)$ but a multi-resolution
version.

A consequence of the filtering is that we obtain an unbiased
estimator.  Moreover, a further consequence is \corred{that if $\bv
  \neq 0$ then}
\begin{equation}
0 < \min_{\bv \in \mathbb{R}^{n}} \frac{\bv\T\bC(\btheta)\bv}{ \|\bv\|^2 } 
\leq 
\min_{\bv \in \mathbb{R}^{n} \backslash {\cal P}^{p}(\bbS)}
 \frac{\bv\T\bC(\btheta)\bv}{ \|\bv\|^2 } 
\leq \max_{\bv \in \mathbb{R}^{n} \backslash {\cal P}^{p}(\bbS)} \frac{\bv\T\bC(\btheta)\bv}
{ \|\bv\|^2 }  
\leq \max_{\bv \in \mathbb{R}^{n}} \frac{\bv\T\bC(\btheta)\bv}{ \|\bv\|^2 }.
\label{Introduction:eqn1}
\end{equation}
This implies that the condition number of $\bC_{W}(\btheta)$ is less
than or equal to the condition number of $\bC(\btheta)$. Thus
computing the inverse of $\bC_{W}(\btheta)$ will be in general more
stable than for $\bC(\btheta)$. In practice, computing the inverse of
$\bC_{W}(\btheta)$ will be much more stable than $\bC(\btheta)$ (See
the results in Tables \ref{SRBF:table1} and \ref{SRBF:table2}).
\corb{High condition number are very bad for numerical methods in
  general. In general any numerical method will suffer if the
  condition number is high.}  Finally, the uncertainties in the
parameter estimates obtained from
(\ref{Introduction:multilevelloglikelihood}) can be quantified using
the Godambe information matrix as described in Sect. 2 and Appendix B
of \cite{Stein2004}.



% \begin{equation}
% \hat{\ell}_{W}(\btheta)
% =-\frac{\hat{n}}{2}\log(2\pi)-\frac{1}{2}\log
% \det\{\hat{\bC}_{W}(\btheta)\}
% -\frac{1}{2}\hat{\bZ}_{W}\T \bC_{W}(\btheta)^{-1}\hat{\bZ}_{W},
% \label{Introduction:multilevelloglikelihoodreduced}
% \end{equation}
% where $\hat{n} \leq n$, $\hat{\bZ}_W := \hat{\bZ}_W(1 : \hat{n})$ and
% $\bC_{W}(\btheta):=\mathbb{E}[\hat{\bZ}_W \hat{\bZ}_W^T]$. In
% section \ref{multilevelestimator} we indicate how to choose $\hat{n}$,
% which will corresponds to the different levels of resolution.

% \begin{rem}
As shown in Section \ref{MultiLevelCovarianceMatrix}, for covariance
functions that are differentiable up to a degree $\tilde f + 1$
(except at the origin), such as the Mat\'{e}rn, our approach leads
to covariance matrices $\bC_{W}$ where most of the coefficients are
small and thus can be safely eliminated.  We construct a level
dependent criterion approach to determine which entries are computed
and the rest are set to zero. With this approach we can now construct
a sparse covariance matrix $\tilde{\bC}_{W}$ that is close to
$\bC_{W}$ in a matrix norm sense even if the observations are highly
correlated with distance. 

The sparsity of $\tilde{\bC}_{W}$ will depend on the following: i) a
positive integer $\tau$, which is a multi-level distance criterion;
ii) a positive integer $\tilde f$, which is the degree of the
multi-level basis and associated accuracy parameters $\tilde p$; and
iii) the smoothness of the covariance function. The accuracy of
$\tilde{\bC}_{W}$ will depend monotonically on these parameters, i.e.,
as we increase $\tau$ and $\tilde f$ (and respectively $\tilde p$) the
matrix $\tilde{\bC}_{W}$ will be closer to ${\bC}_{W}$ in a norm
sense. This is explained in detail in Section
\ref{MultiLevelCovarianceMatrix}.

The choice of the projectors $\bL$ and $\bW$ will determine how
efficiently each likelihood function
\eqref{Introduction:multilevelloglikelihood} evaluation is
solved. {\it Indeed, we desire the transformation to have the following
properties:}
i) {\bf Stability:} The matrices $\bL$ and $\bW$ have orthogonal
  rows and the stacked matrix $[\bL; \bW]$ is orthonormal;
ii) {\bf Fast computation:} The computational cost of applying the
  matrix $[\bL; \bW]$ to a vector is ${\cal
    O}(n(\log{ n})^{\xi})$ for some small integer $\xi$;
iii) {\bf Fast log determinant computation:} The computational cost
  of computing $\log \det\{\tilde{\bC}_{W}(\btheta)\}$ to be bounded 
   by ${\cal
    O}(n^{3/2})$ in 2D and ${\cal
    O}(n^{2})$ in 3D. We also want to restrict the memory storage to ${\cal
    O}(n (\log{n})^{\xi})$;
iv) {\bf Fast inversion:} The computational cost of computing
  $\bC_{W}(\btheta)^{-1} \bZ_{W}$ to a desired accuracy $\varepsilon$ 
is better than ${\cal
  O}(n^{2})$. Memory storage is also desirable to be restricted to
  ${\cal O}(n (\log{n})^{\xi})$;
v) {\bf Accuracy:} Determinant computation and inversion are
also required to be accurate. \corred{We achieve 
the properties i) - v) in this paper.}

In Section \ref{MultiLevelREML} we describe how to construct
multi-level matrices $\bL$ and $\bW$ that satisfy properties i) and
ii) for most practical observation location placements (random for
example). \corred{We apply $\bL$ and $\bW$ to construct the sparse
  multi-level covariance matrix $\tilde \bC_W(\btheta)$.} The
determinant of the multi-level sparse covariance matrix
\corred{$\tilde \bC_{W}(\btheta)$} and the term \corred{$\tilde
  \bC_{W}(\btheta)^{-1}\bZ_{W}$} are computed by exploiting an
accurate sparse \corred{Cholesky} representation \corred{of} $\bC_{W}$
(properties iii) and v) ). The term $\bC_{W}(\btheta)^{-1} \bZ_{W}$
can also be computed by applying a Preconditioned Conjugate Gradient
(PCG) to a desired accuracy (properties iv) and v) ).  In Section
\ref{multilevelkriging} the multi-level kriging method is described.
In Section \ref{numericalstudy} we demonstrate the efficiency of our
method for numerous covariances and irregularly placed observations.
We are able to solve the fast inversion for up to 512,000 observations
to a relative accuracy of $10^{-3}$ with respect to the {\it
  unpreconditioned} system. It is important to note that the achieved
accuracy of preconditioned system will not necessarily imply
accuracy of unpreconditioned system if the condition number of the
preconditioner is high.  Furthermore, we test our approach to estimate
the covariance parameters of problems of up to 128,000
observations. In addition, the accuracy of the kriging estimates are
tabulated for different size problems. 



% {\red {\bf Notes to guide research (Not to be included in current
% form in the paper):}

% \begin{itemize}
% \item {\bf Hard problem:} With the advent of large scale sensor
%   networks current methods to perform covariance estimation and
%   Kriging employ either direct and/or dimensionality reduction or make
%   very restrictive assumptions on the placement of the nodes.  Direct
%   methods are accurate but slow, ${\cal O}(N^3)$, making it unfeasible
%   to perform estimation for large scale problem. Dimensionality
%   reduction method are either restrictive or ad hoc.

% \item {\bf Current Technical Idea:} In \cite{Castrillon2013} we
%   developed a Multi-Level method to transform the Kriging problem into
%   a Multi-Level representation that is much easier to solve. This
%   method is very effective even for highly irregular placement of the
%   nodes.  We propose to use the Multi-Level representation to estimate
%   directly the covariance matrix.

% \item {\bf Impact:} If successful, large scale problem (500,000+
%   source points) will be possible to be computed on a single
%   workstation instead of needing a supercomputer. 

% \item {\bf Objective:} To develop a fast and reliable restricted
%   maximum likelihood Estimator based on a Multi-Level representation.
% \end{itemize}
% }
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Multi-Level Basis}
\label{MultiLevelREML}


In this section we establish the general structure of the Multi-Level
Basis (MB) that is used solve the estimation and prediction problem.
We refer the reader to \cite{Castrillon2013} for a detailed
description.  The MB can then be used to: (i) form the multi-level
REML function \eqref{Introduction:multilevelloglikelihood}; (ii)
sparsify the covariance matrix $\bC_{W}(\btheta)$; and (iii) improve
the conditioning over the covariance matrix $\bC(\btheta)$. But first,
we establish some notation and definitions:
\begin{itemize}


\item Let $\alpha:= (\alpha_{1},\dots,\alpha_{d}) \in \mathbb{Z}{^d}$,
  $|\alpha| := \alpha_{1}+\dots+\alpha_{d}$, $\bx : = [x_1,\dots,x_d]\T$
  and $D^\alpha_{\bx} := \frac{\partial^{\alpha_{1} + \dots +
      \alpha_{d}}} {\partial x_{1}^{\alpha_{1}} \dots \partial
    x_{d}^{\alpha_{d}}}$. For any $h \in \bbN_0$ (where $\mathbb{N}_0
  := \mathbb{N} \cup \{0\}$) let ${\cal Q}^d_h$ be the set of
  monomials $\{x_1^{\alpha_1} \dots x_d^{\alpha_d}\,\,\,|\,\,\,
  |\alpha| \leq h\}$. Furthermore, let $\bM_h$ be the design matrix
  with respect to all the monomials in ${\cal Q}^d_h$. The
  number of monomials of degree $h$ with dimension $d$ is
  $\begin{pmatrix} d + h \\ h \end{pmatrix}$. 


\item We shall restrict the Gaussian spatial random field
  \eqref{Introduction:noisemodel} design matrix $\bM$ to $\bM_f$,
  where $f$ is the degree of the model. Thus $p$ will be equal to the
  number of monomials in ${\cal Q}^d_f$, which is $p
  := \begin{pmatrix} d + f \\ f \end{pmatrix}$.


\item Let $\tilde f \geq f$ be the degree of the multi-level basis and
  $\bM_{\tilde f}$ the associated design matrix.  The number of
  monomials in ${\cal Q}^d_{\tilde f}$ shall be referred as the
  accuracy parameter $\tilde p :=\begin{pmatrix} d + \tilde f
  \\ \tilde f \end{pmatrix}$.  These parameters are chosen by the
  user and are used to construct the multi-level basis.

\item Let $\bC(\btheta) : = \{ \phi(r_{i,j}; \btheta) \}$ where $\phi$
  is the covariance function, $r_{i,j} := \|\bs_i - \bs_j\|_2$ and
  $\bs_i,\bs_j \in \R^{d}$ for $i,j = 1, \dots, n$. Alternatively we
  refer to $\phi(\bx,\by; \btheta)$ as $\phi(r; \btheta)$, where
  $r:=\|\bx - \by\|_2$ and $\bx,\by \in \R^d$. Suppose ${\cal
    P}^{p}(\mathbb{S})$ is the span of the columns of the design
  matrix $\bM_f$.  We now assume that $\phi(r; \btheta)$ is a positive
  definite function and $C^{\tilde{f}+1}(\R)$ for all $r \in \R$
  except at the origin.

\item For any index $i,j \in \mathbb{N}_{0}$, $1 \leq i \leq n$, $1 \leq j
\leq n$, let $\bve_{i}[j] = \delta[i-j]$, where
$\delta[\cdot]$ is the discrete Kronecker delta function.

\end{itemize}

\begin{rem}
In practice instead of using the set of monomials ${\cal Q}^d_f$ we
use the set of Chebyshev polynomials of the first kind as these lead
to a more stable numerical scheme. However, for simplicity of the
presentation we keep it to monomials.
\end{rem}

% We now show how to construct a multilevel basis for $\R^n$.

The first step is to decompose the locations into a
series of multi-level cubes of dimension $d$. Without loss of
generality we assume that all the locations are contained in a unit
cube $B^{0}_{0}$ at level 0 and index 0. If the number of locations
inside $B^{0}_{0}$ is more than $p$ then equally subdivide $B^{0}_{0}$
into $2^{d}$ cubes ($B^{1}_{0}$, \dots, $B^{1}_{2^d - 1}$), where $d$ is the
number of dimensions. If the number of locations is $p$ or less then
stop, associate every location with the cube $B^{0}_{0}$ and denote
this as a leaf cube.  Otherwise, for each non-empty cube $B^{q}_{k}$ at
level $q = 1$ and index $k$ if the number of locations is more than
$p$ then subdivide, otherwise associate all the locations to
$B^{q}_{k}$ and denote this as a leaf cube.  This process is repeated for
all the subdivided cubes at levels $q = 2, \dots,$ until no
subdivisions are possible. The result is a tree structure with $0,
\ldots, t$ levels (See Algorithm~1 in \cite{Castrillon2013} for more
details).  We denote the leaf cubes as all the non-empty cubes that
contain at most $p$ locations, i.e., they will correspond to the leafs
of the tree structure.

\begin{rem}
For practical cases, $t$ increases proportionally to
$\log{n}$. If the inter location spacing collapses as $n^{-q}$ , where
$q$ is independent of $n$ , then $q\log{n}$ levels are needed, see
Section 4 in \cite{Beatson1997} for details.
\end{rem}

Suppose that there is a one-to-one mapping between the set of unit
vectors ${\cal E}:=\{\bve_{1},\dots,\bve_{n}\}$, which we denote as
leaf unit vectors, and the set of locations $\{
\bs_{1},\dots,\bs_{n}\}$, i.e. $\bs_{i} \longleftrightarrow \bve_{i}$
for all $i = 1, \dots, n$. It is clear that the space of
$\{\bve_{1},\dots,\bve_{n}\}$ is $\bbR^{n}$. The next step is to replace
${\cal E}$ with a new basis of $\R^{n}$ that is multi-level,
orthonormal and gives us the desired properties i) and ii) from Section
\ref{Introduction}. In \cite{Castrillon2013} the reader can find full
details on such a construction. However, for the sake of clarity for
the rest of the paper we associate the multi-level domain decomposition
to the multi-level basis:

\begin{figure}

\vspace{-5cm}

\begin{tabular}{c c}

\begin{tikzpicture}[scale=.75] 
  \begin{scope} 
    [place/.style={circle,draw=blue!50,fill=blue!20,thick,
      inner sep=0pt,minimum size=1.5mm}]
    \draw[step=8,gray,very thin] (0, 0) grid (8, 8);

    % Anchors
    \node[anchor=center] at (0, -4.5) {$$};
    \node[anchor=center] at (0,   10) {$$};
    % \node[anchor=center] at (16, -4.5) {$$};
    % \node[anchor=center] at (16,   10) {$$};

    \node at (0.5,7.5) [place] {};
    \node at (0,8.0) [] {($\bs_1$,$\bve_1)$};

    \node at (0.75,7.2) [place] {};
    \node at (0.3,6.3) [place] {};
    \node at (0.6,6.7) [place] {};
    \node at (1.3,6.3) [place] {};
    \node at (1.5,6.8) [place] {};

    \node at (2.5,4.5) [place] {};
    % \node at (3.5,7.5) [place] {};
    \node at (0.3,4.5) [place] {};
    \node at (1.6,5.7) [place] {};
    \node at (2.6,6.3) [place] {};
    \node at (3.2,5.2) [place] {};

    \node at (0.5,3.5) [place] {};
    \node at (1.5,2.5) [place] {};
    \node at (2.3,2.2) [place] {};
    \node at (3.3,2.5) [place] {};


    \node at (2.3,0.3) [place] {};
    \node at (2.7,0.5) [place] {};
    \node at (2.2,1.4) [place] {};
    \node at (2.6,1.4) [place] {};
    \node at (3.3,0.5) [place] {};
    \node at (3.7,0.7) [place] {};
    \node at (3.2,1.3) [place] {};
    \node at (3.6,1.5) [place] {};


    \node at (4.3,2.3) [place] {}; \node at (5.7,3.5) [place] {};
    \node at (6.2,3.4) [place] {}; \node at (7.3,2.4) [place] {};
    \node at (7.5,0.5) [place] {}; \node at (7.75,1) [] {($\bs_n$,$\bve_n)$};
    \node at (8.5,8.5) [] {$B^{0}_0$};
  \end{scope}

\end{tikzpicture}
\vspace{-8cm}
& 

\begin{tikzpicture}[scale=1]

    \node[anchor=center] at (0, -4.5) {$$};
    \node[anchor=center] at (0,   10) {$$};

\begin{scope}[xshift=5cm, yshift=4cm,
place/.style={circle,draw=blue!50,fill=blue!20,thick,
      inner sep=0pt,minimum size=1.5mm},
placer/.style={circle,draw=red!50,fill=red!20,thick,
      inner sep=0pt,minimum size=1.5mm},
]
\Tree [.\node[placer]{$B^{0}_{0}$}; 
             [.\node[place]{$B^{1}_{0}$};
                    [.\node[place]{$B^{2}_{0}$}; 
                            \node[place]{$B^{3}_{0}$}; 
                            \node[place]{$B^{3}_{1}$}; 
                            \node[place]{$B^{3}_{2}$}; 
                    ] 
                    [.\node[place]{$B^{2}_{1}$}; ] 
                    [.\node[place]{$B^{2}_{2}$}; ] 
                    [.\node[place]{$B^{2}_{3}$}; ] 
             ]                                         
             [.\node[place]{$B^{1}_{1}$}; 
                    [      [. $\vdots$ ] ]
             ]   
             [.\node[place]{$B^{1}_{2}$};
                    [      [.\node[place]{$B^{2}_{8}$}; ]
                           [.\node[place]{$B^{2}_{10}$}; ] 
                           [.\node[place]{$B^{2}_{11}$}; ]
                                          ]]] 
\end{scope}
\end{tikzpicture} 

%\Tree [.$B^{0}_{0}$ {$B^{1}_{0}$} {$B^{2}_{0}$} {$B^{3}_{0}$} ]
\\

\begin{tikzpicture}[scale=.75] 
  \begin{scope} 
    [place/.style={circle,draw=blue!50,fill=blue!20,thick,
      inner sep=0pt,minimum size=1.5mm}]
    \draw[step=4,gray,very thin] (0, 0) grid (8, 8);

    % Anchors
    \node[anchor=center] at (0, -4.5) {$$};
    \node[anchor=center] at (0,   10) {$$};

    \node at (0.5,7.5) [place] {};
    \node at (0,8.0) [] {($\bs_1$,$\bve_1)$};

    \node at (0.75,7.2) [place] {};
    \node at (0.3,6.3) [place] {};
    \node at (0.6,6.7) [place] {};
    \node at (1.3,6.3) [place] {};
    \node at (1.5,6.8) [place] {};

    \node at (2.5,4.5) [place] {};
    % \node at (3.5,7.5) [place] {};
    \node at (0.3,4.5) [place] {};
    \node at (1.6,5.7) [place] {};
    \node at (2.6,6.3) [place] {};
    \node at (3.2,5.2) [place] {};

    \node at (0.5,3.5) [place] {};
    \node at (1.5,2.5) [place] {};
    \node at (2.3,2.2) [place] {};
    \node at (3.3,2.5) [place] {};


    \node at (2.3,0.3) [place] {};
    \node at (2.7,0.5) [place] {};
    \node at (2.2,1.4) [place] {};
    \node at (2.6,1.4) [place] {};
    \node at (3.3,0.5) [place] {};
    \node at (3.7,0.7) [place] {};
    \node at (3.2,1.3) [place] {};
    \node at (3.6,1.5) [place] {};


    \node at (4.3,2.3) [place] {}; \node at (5.7,3.5) [place] {};
    \node at (6.2,3.4) [place] {}; \node at (7.3,2.4) [place] {};
    \node at (7.5,0.5) [place] {}; \node at (7.75,1) [] {($\bs_n$,$\bve_n)$};
    \node at (8.5,8.5) [] {$B^{0}_0$};
    \node at (3.5,3.5) [scale=0.75] {$B^1_1$};
    \node at (7.5,3.5) [scale=0.75] {$B^1_2$};
    \node at (3.5,7.5) [scale=0.75] {$B^1_0$};

 \end{scope}
\end{tikzpicture} 
\vspace{-8cm}

& 

\begin{tikzpicture}[scale=1]

    \node[anchor=center] at (0, -4.5) {$$};
    \node[anchor=center] at (0,   10) {$$};

\begin{scope}[xshift=5cm, yshift=4cm,
place/.style={circle,draw=blue!50,fill=blue!20,thick,
      inner sep=0pt,minimum size=1.5mm},
placer/.style={circle,draw=red!50,fill=red!20,thick,
      inner sep=0pt,minimum size=1.5mm}
]
\Tree [.\node[placer]{$B^{0}_{0}$}; 
             [.\node[placer]{$B^{1}_{0}$};
                    [.\node[place]{$B^{2}_{0}$}; 
                            \node[place]{$B^{3}_{0}$}; 
                            \node[place]{$B^{3}_{1}$}; 
                            \node[place]{$B^{3}_{2}$}; 
                    ] 
                    [.\node[place]{$B^{2}_{1}$}; ] 
                    [.\node[place]{$B^{2}_{2}$}; ] 
                    [.\node[place]{$B^{2}_{3}$}; ] 
             ]                                         
             [.\node[placer]{$B^{1}_{1}$}; 
                    [      [. $\vdots$ ] ]
             ]   
             [.\node[placer]{$B^{1}_{2}$};
                    [      [.\node[place]{$B^{2}_{8}$}; ]
                           [.\node[place]{$B^{2}_{10}$}; ] 
                           [.\node[place]{$B^{2}_{11}$}; ]
                                          ]]] 
\end{scope}
\end{tikzpicture} 

\\


\begin{tikzpicture}[scale=.75] 
  \begin{scope} 
    [place/.style={circle,draw=blue!50,fill=blue!20,thick,
      inner sep=0pt,minimum size=1.5mm}]

    \draw[step=4,gray,very thin] (0, 0) grid (8, 8);
    \draw[step=2,gray,very thin] (0, 0) grid (4, 4);
    \draw[step=2,gray,very thin] (0, 4) grid (4, 8);
    \draw[step=2,gray,very thin] (4, 0) grid (8, 4);

    \draw[step=1,gray,very thin] (0, 6) grid (2, 8);
    \draw[step=1,gray,very thin] (2, 0) grid (4, 2);
    

    % Anchors
    \node[anchor=center] at (0, -4.5) {$$};
    \node[anchor=center] at (0,   10) {$$};

    \node at (0.5,7.5) [place] {};
    \node at (0,8.0) [] {($\bs_1$,$\bve_1)$};

    \node at (0.75,7.2) [place] {};
    \node at (0.3,6.3) [place] {};
    \node at (0.6,6.7) [place] {};
    \node at (1.3,6.3) [place] {};
    \node at (1.5,6.8) [place] {};

    \node at (2.5,4.5) [place] {};
    % \node at (3.5,7.5) [place] {};
    \node at (0.3,4.5) [place] {};
    \node at (1.6,5.7) [place] {};
    \node at (2.6,6.3) [place] {};
    \node at (3.2,5.2) [place] {};

    \node at (0.5,3.5) [place] {};
    \node at (1.5,2.5) [place] {};
    \node at (2.3,2.2) [place] {};
    \node at (3.3,2.5) [place] {};


    \node at (2.3,0.3) [place] {};
    \node at (2.7,0.5) [place] {};
    \node at (2.2,1.4) [place] {};
    \node at (2.6,1.4) [place] {};
    \node at (3.3,0.5) [place] {};
    \node at (3.7,0.7) [place] {};
    \node at (3.2,1.3) [place] {};
    \node at (3.6,1.5) [place] {};


    \node at (4.3,2.3) [place] {}; \node at (5.7,3.5) [place] {};
    \node at (6.2,3.4) [place] {}; \node at (7.3,2.4) [place] {};
    \node at (7.5,0.5) [place] {}; \node at (7.75,1) [] {($\bs_n$,$\bve_n)$};
    \node at (8.5,8.5) [] {$B^{0}_0$};
    \node at (3.5,7.5) [scale=0.75] {$B^2_3$};
    \node at (7.5,3.5) [scale=0.75] {$B^2_{11}$};
    \node at (1.5,3.5) [scale=0.75] {$B^2_4$};

 \end{scope}
\end{tikzpicture} 
&

\begin{tikzpicture}[scale=1]

    \node[anchor=center] at (0, -4.5) {$$};
    \node[anchor=center] at (0,   10) {$$};

\begin{scope}[xshift=5cm, yshift=4cm,
placer/.style={circle,draw=red!50,fill=red!20,thick,
      inner sep=0pt,minimum size=1.5mm}]
\Tree [.\node[placer]{$B^{0}_{0}$}; 
             [.\node[placer]{$B^{1}_{0}$};
                    [.\node[placer]{$B^{2}_{0}$}; 
                            \node[placer]{$B^{3}_{0}$}; 
                            \node[placer]{$B^{3}_{1}$}; 
                            \node[placer]{$B^{3}_{2}$}; 
                    ] 
                    [.\node[placer]{$B^{2}_{1}$}; ] 
                    [.\node[placer]{$B^{2}_{2}$}; ] 
                    [.\node[placer]{$B^{2}_{3}$}; ] 
             ]                                         
             [.\node[placer]{$B^{1}_{1}$}; 
                    [      [. $\vdots$ ] ]
             ]   
             [.\node[placer]{$B^{1}_{2}$};
                    [      [.\node[placer]{$B^{2}_{8}$}; ]
                           [.\node[placer]{$B^{2}_{10}$}; ] 
                           [.\node[placer]{$B^{2}_{11}$}; ]
                                          ]]] 
\end{scope}
\end{tikzpicture} 

\\

\end{tabular}

\vspace{-3cm}

\caption{Multi-level domain decomposition of the location of
  observations for $d = p =2$. All the observation locations $\bs_i$
  are assumed to be contained in the unit cube $B^{0}_{0}$ (colored red
  node in tree). If the number of observations are greater than $p =
  2$ we subdivide into equal boxes. At the end we obtain a multi-level
  decomposition of the observation points.}
\label{MLRLE:fig1}
\end{figure}


\begin{itemize}
\item 
For each non empty cube $B^i_k$, for $i = 0, \dots, t$, associate a
series of multi-level basis vectors $\{ \bpsi^{i,k}_{\tilde k_1},
\bpsi^{i,k}_{\tilde k_2}, \dots \}$ that have the following property:
\begin{equation}
\bg^T  \bpsi^{i,k}_{\tilde k_j} = \sum_{a=1}^{n} \bg[a] \bpsi^{i,k}_{\tilde k_j}[a] = 0,
\label{hbconstruction:eqn2}
\end{equation}
for $j = 0,\dots$ and for all the vectors $\bg$ that are columns of
$\bM_{\tilde f}$. Furthermore, let $\bW^{i,k}$ be a matrix
$[ \bpsi^{i,k}_{\tilde k_1}, \bpsi^{i,k}_{\tilde k_2}, \dots ]$.

\item For $i = 0, \dots, t$ let $\bW_{i} := [\bW^{i,0}, \bW^{i,1},
  \dots]$

\item If $\tilde p > p$ we will have an extra $\tilde p - p$
  vectors corresponding to a level $-1$ for the initial cube
  $B^{0}_{0}$. Now, associate $\tilde p - p$ multi-level vectors $\{
  \bpsi^{-1,0}_{0},\bpsi^{-1,0}_{1}, \dots, \bpsi^{-1,0}_{\tilde p - p}
  \}$ that have the following property:
\begin{equation}
\bg\T \bpsi^{-1,0}_{j} =
\sum_{a=1}^{n} \bg[a] \bpsi^{-1,0}_{j}[a] = 0,
\label{hbconstruction:eqn3}
\end{equation}
for $j = 1,\dots, \tilde p - p$ and for all the vectors $\bg$ that are
columns of $\bM_{\tilde f}$. Similarly as above, let $\bW_{-1}:= [
  \bpsi^{-1,0}_{0},\bpsi^{-1,0}_{1}, \dots, \bpsi^{-1,0}_{\tilde p - p}
]$.

\item In total we will have $n - p$ multi-level vectors
and the transform matrix $\bW \in \R^{(p - n) \times n}$ is built as
$\bW : = [\bW_t, \dots, \bW_0, \bW_{-1}]^T$.

\item Now, it is clear that $\bW \bg = {\bf 0}$ for any $\bg \in
  \bM_{f}$. To complete the basis to span $\R^{n}$ we need $p$ more
  orthonormal vectors. In \cite{Castrillon2013} it is shown how to
  compute such a basis and stack the vectors as rows in the matrix
  $\bL \in \R^{p \times n}$.

\end{itemize}

With the construction of $\bW$ and $\bL$ we will have the following
properties: a) the matrix $ \bP := \left[
\begin{array}{c}
\bW \\
\bL
\end{array}
\right
]$ 
is orthonormal, i.e., $\bP\bP\T= \bI_n$;
b) any vector $\bv \in \R^{n}$ can be written as $\bv = \bL\T\bv_{L} +
  \bW\T\bv_{W}$ where $\bv_{L} \in \R^{p} $ and $\bv_{W} \in \R^{n-p}$ are
  unique; c) the matrix $\bW$ contains at most ${\cal O}(nt)$ non-zero
  entries and $\bL$ contains at most ${\cal O}(np)$ non-zero
  entries. This implies that for any vector $\bv \in \R^{n}$ the
  computational cost of applying $\bW \bv$ is at most ${\cal
    O}(nt)$ and $\bL \bv$ is at most ${\cal O}(np)$.
 
\section{Multi-Level Covariance Matrix}
\label{MultiLevelCovarianceMatrix}

In this section we show how we can use the matrix $\bW$ to produce a
highly sparse representation of $\bC_W(\btheta)$ with a
level-dependent tapering technique. 
% We refer the reader to Appendix A
% for a detailed description of the multi-level basis.

With the MB we can transform the observation data vector $\bZ$ by
applying the matrix $\bW$. This leads to the multi-level
log-likelihood function
\eqref{Introduction:multilevelloglikelihood}. The covariance matrix
$\bC(\btheta)$ is now transformed into $\bC_{W}(\btheta)$ with the
structure shown in Figure \ref{multilevelcov:fig1} where each of the
blocks $\bC^{i,j}_W = \bW_i \bC(\btheta) \bW_j \T$ for all $i,j = 0,
\dots, t$. This implies that %the $l^{th}$ row and $k^{th}$ column of
the entries of the matrix $\bC^{i,j}_W$ are formed from all the
interactions of the MB vectors between level $i$ and $j$. Thus for any
$\bpsi^{i,k}_{\tilde{l}}$ and $\bpsi^{j,l}_{\tilde{k}}$ vectors there
is a unique entry of $\bC^{i,j}_W$ of the form $
(\bpsi^{i,k}_{\tilde{k}})\T \bC(\btheta) \bpsi^{j,l}_{\tilde{l}}$.
The blocks $\bC_W^{i,j}$, where $i=-1$ or $j=-1$, correspond to the
case where the accuracy term $\tilde{p} > p$.
\begin{figure}[h!] 
\psfrag{a}[l]{$\bC^{t,t}_{W}$} 
\psfrag{b}[c]{$\dots$}
  \psfrag{c}[c]{\hspace{1mm} {$\bC^{t,i}_{W}$} } 
  \psfrag{d}[c]{ \hspace{5mm}  $\vdots$}
  \psfrag{e}[l]{$\ddots$} 
  \psfrag{f}[c]{$\vdots$}
  \psfrag{g}[c]{\hspace{5mm} $\bC^{-1,t}_{W}$} 
  \psfrag{h}[c]{$\dots$}
  \psfrag{i}[c]{\small $\bC^{i,i}_{W}$}
\begin{center}
\includegraphics[width=2.7in,height=2.10in,trim=0 0 110 0,clip]{./NewKmatrix}
\end{center}
\caption{Organization of the multi-level covariance matrix
  $\bC_W(\btheta)$. For this figure $i = -1$.}
\label{multilevelcov:fig1}
\end{figure}


The following lemma relates the covariance function $\phi$, the degree
$\tilde{f}$ (corresponding to the accuracy parameter $\tilde{p}$) of
the design matrix $\bM_{\tilde{f}}$ to the decay of the entries of the
matrix $\bC_W(\btheta)$.
\begin{lem} Let $B_{\ba}$ be the smallest ball in ${\mathbb R}^{d}$
with radii $r_{\ba}$ centered around the midpoint $\ba \in {\mathbb
  R}^{d}$ of the cube $B^{i}_{l}$ such that $B^{i}_{l} \subset
B_{\ba}$. Similarly, let $B_{\bb}$ be the smallest ball in ${\mathbb
  R}^{d}$ with radii $r_{\bb} \in {\mathbb R}^{d}$ centered around the
midpoint $\bb$ of the cube $B^{j}_{k}$ such that $B^{j}_{k} \subset
B_{\bb}$.  Now, since $\bpsi^{i,l}_{\tilde{l}}$ and
$\bpsi^{j,k}_{\tilde{k}}$ satisfy the moment
orthogonality condition from equations \eqref{hbconstruction:eqn2} 
and \eqref{hbconstruction:eqn3} for
all $\bg \in {\cal P}^{\tilde{p}}(\bbS)$ then the following bound holds:
\begin{equation}
|(\bpsi^{i,k}_{\tilde{k}})\T \bC(\btheta) \bpsi^{j,l}_{\tilde{l}} |
\leq
\sum_{|\alpha|=\tilde{f}+1} \sum_{|\beta|=\tilde{f}+1}
\frac{r_{\ba}^{\alpha}}{\alpha!} \frac{r_{\bb}^{\beta}}{\beta!}
\sup_{\bx \in B_{\ba}, \by \in B_{\bb}}
 |D^{\alpha}_{\bx} D^{\beta}_{\by} \phi(\bx,\by;\btheta)|,
\label{multilevelcov:eqn1}
\end{equation}
\label{multilevelcov:lemma1}
for $i,j = 0,\dots,t$.
\end{lem}
From Lemma \ref{multilevelcov:lemma1} we observe that the decay of the
entries of $\bC_W(\btheta)$ is dependent on the magnitude of the
derivatives of the covariance function $\phi(r;\btheta)$, the size of
$B_{\ba}$ and $B_{\bb}$ and the degree of ${\cal Q}^d_{\tilde{f}}$. Thus if
$\phi(r;\btheta)$ is smooth on $B_{\ba}$ and $B_{\bb}$ the entries of
$\bC_W(\btheta)$ will be small.



\begin{exam}
In Figure \ref{multilevelcov:fig2} we show a comparison between (a)
the covariance matrix $\bC(\btheta)$ and (b) the multi-level
covariance matrix $\bC_W(\btheta)$ for the following example: 1)
$\phi(r;\btheta) : = \exp(-r)$ and $d = 3$.  2) The observation
locations ($n = 8000$) are sampled from a uniform distribution on the
unit cube $[0,1]^3$. The actual values of the observations are not
necessary for this example.  3) $f = 3$ (leading to $p = 20$
monomials).  4) We sort the $x_1$ direction location from 0 to 1. This
is done for visualization reasons so that we may observe the decay
in the matrix $\bC(\btheta)$.
\begin{figure}[h!]
\begin{center}
\begin{tabular}{cc}
   \includegraphics[ width=3.0in, height=3in
   ]{./CovarianceExp1} & \includegraphics[ width=3.0in,
     height=3in ]{./HBCovarianceExp1} \\ (a)
   $\log_{10}{\mbox{abs}(\bC(\btheta))}$ & (b) $\log_{10}{\mbox{abs}(\bC_W(\btheta))}$ \\
\end{tabular} 
\end{center}
  \caption{Covariance matrix comparison between covariance matrices
    (a) $\log_{10}{\mbox{abs}(\bC(\btheta))}$ and (b)
    $\log_{10}{\mbox{abs}(\bC_W(\btheta))}$ for the exponential covariance
    function $\phi(r) = \exp(-r)$ with $n = 8000$, $p = 20$ and $d =
    3$.}
\label{multilevelcov:fig2}
\end{figure}
Notice that the decay of $\bC(\btheta)$ is dependent on the covariance
function $\phi(r;\btheta)$. It is clear that for this case a tapering
technique would not be very effective as most of the entries are
comparable in magnitude. In contrast a few of the entries of
$\bC_W(\btheta)$ with high magnitude are concentrated around
particular regions while most of the entries have very small
magnitudes making a hierarchical tapering technique to sparsify the
matrix a viable option.
\label{multilevelcov:example1}
\end{exam}

To produce a sparse matrix from $\bC_W(\btheta)$ we execute the
following multi-level tapering technique: 

\begin{itemize}
\item For all cubes $B^i_k$ at level $i$ let $L^{i,0}_k := B^i_k$ and
  $L^{i,j}_k := L^{i,j-1}_k \cup \{$union of all cubes at level $i$
  that share a face or corner with $L^{i,j-1}_k \}$ for $j =
  0,1,\dots$. A construction example is shown in Figure
  \ref{multilevelcov:fig3} for level $i$. Now, perform this
  construction for $i = 0,\dots,t$.

\begin{figure}[!h]
\begin{center}
\begin{tikzpicture}[scale=0.7] 
 \begin{scope} 
   [place/.style={circle,draw=blue!50,fill=blue!20,thick,
       inner sep=0pt,minimum size=0.5mm}]

    \draw[step=1, very thick] (1, 0) grid (2, 1);
    \node at (1.5,0.5) [] {$B^i_k$};
    \node at (1.5,-1) [] {$L^{i,0}_k$};

    \draw[step=1, dashed] (4, 0) grid (7, 3);
    \draw[step=1,blue!50, very thick ] (5, 1) grid (6, 2);

    \draw[very thick] (4, 0) -- (7, 0);
    \draw[very thick] (7, 0) -- (7, 3);
    \draw[very thick] (7, 3) -- (4, 3);
    \draw[very thick] (4, 3) -- (4, 0);

    \node at (5.5,1.5) [] {$L^{i,0}_k$};
    \node at (5.5,-1) [] {$L^{i,1}_k$};
    
    \draw[step=3,gray, dashed] (9, 0) grid (18, 9);
    \draw[step=3,blue!50, very thick] (12, 3) grid (15, 6);
    
    \draw[very thick] (9,  0) -- (18, 0);
    \draw[very thick] (18, 0) -- (18, 9);
    \draw[very thick] (18, 9) -- (9, 9);
    \draw[very thick] (9, 9) -- (9, 0);


    \node at (13.5,4.5) [] {$L^{i,1}_k$};
    \node at (13.5,-1) [] {$L^{i,2}_k$};
   
\end{scope}
\end{tikzpicture}
\end{center}

\caption{Construction of expanded cubes $L^{i,\tau}_{k}$, $\tau = 0,
  1, \dots$ from initial cube $B^{i}_{k}$.}
\label{multilevelcov:fig3}
\end{figure}

\item Set a user given constant $\tau \in \mathbb{N}_0$

\item The entry of $\bC_W(\btheta)$ corresponding to
  $(\bpsi^{i,k}_{\tilde k})\T \bC(\btheta) \bpsi^{j,l}_{\tilde l}$,
  for $i,j = 0,\dots,\tau$ is computed if the following level
  dependent criterion is true: If $\mbox{($j \geq i$ and $B^{j}_{l}
    \subset L^{i,\tau}_{k}$) or ($j < i$ and $B^{i}_{k} \subset
    L^{j,\tau}_{l}$)}$ is true for the given $\tau \in \mathbb{N}_0$
  then compute the entry $(\bpsi^{i,k}_{\tilde k})\T \bC(\btheta)
  \bpsi^{j,l}_{\tilde l}$.


\item For the case that $i = -1$ or $j = -1$ the entry
  corresponding to $(\bpsi^{i,k}_{\tilde k})\T \bC(\btheta)
  \bpsi^{j,l}_{\tilde l}$ is always computed.

%% \item For the case that $i = -1$ and $j = 0, \dots, t$ the entry
%%   corresponding to $(\bpsi^{-1,k}_{k})\T \bC(\btheta)
%%   \bpsi^{j,l}_{l}$ is computed if $\mbox{( $B^{j}_{l} \subset
%%     L^{0,\tau}_{k}$)}$ is true (i.e. always).

%% \item For the case that $i = 0,\dots,t$ and $j = -1$ the entry
%%   corresponding to $(\bpsi^{i,k}_{k})\T \bC(\btheta)
%%   \bpsi^{-1,l}_{l}$ is computed if ($B^{i}_{k} \subset
%%   L^{0,\tau}_{l}$) is true (i.e. always).

%% \item For the case that $i = -1$ and $j = -1$ the entry
%%   corresponding to $(\bpsi^{-1,k}_{k})\T \bC(\btheta)
%%   \bpsi^{-1,l}_{l}$ is always computed.

\end{itemize}

From this distance criterion we can apriori determine which entries of
the sparse matrix $\tilde \bC_W(\btheta)$ are to be computed. For any
given row of the matrix $\tilde \bC_W(\btheta)$ corresponding to level
$i$ and index $k$ construct the expanded cube $L^{i,\tau}_{k}$. Now,
for $j = i, \dots, t$ find all the cubes $B^j_l$ with the
corresponding index $l$ that are contained in $L^{i,\tau}_{k}$ (See
Figure \ref{multilevelcov:fig5}). For $j = 0, \dots, i - 1$ find all
the extended cubes $L^{j,\tau}_{l}$ such that $B^i_k \subset
L^{j,\tau}_{l}$. For $i,j = 0,\dots,t$ this action can be performed
efficiently by using the tree shown in Figure \ref{MLRLE:fig1}. 

\begin{figure}[!h]
\begin{center}
\begin{tikzpicture}[scale=0.7] 
 \begin{scope} 
   [place/.style={circle,draw=blue!50,fill=blue!20,thick,
       inner sep=0pt,minimum size=0.5mm}]

    %% \draw[step=1, very thick] (1, 0) grid (2, 1);
    %% \node at (1.5,0.5) [] {$B^i_k$};
    %% \node at (1.5,-1) [] {$L^{i,0}_k$};

    %% \draw[step=1, dashed] (4, 0) grid (7, 3);
    %% \draw[step=1,blue!50, very thick ] (5, 1) grid (6, 2);

    %% \draw[very thick] (4, 0) -- (7, 0);
    %% \draw[very thick] (7, 0) -- (7, 3);
    %% \draw[very thick] (7, 3) -- (4, 3);
    %% \draw[very thick] (4, 3) -- (4, 0);

    %% \node at (5.5,1.5) [] {$L^{i,0}_k$};
    %% \node at (5.5,-1) [] {$L^{i,1}_k$};
    
    \draw[step=1, gray,  dashed] (9, 0) grid (18, 9);

    \draw[step=1, blue!50, very thick] (16, 6) grid (17, 7);
    \node at (16.5,6.5) [ ] {$B^{j}_{l_1}$};

    \draw[step=1, blue!50, very thick] (11, 7) grid (12, 8);
    \node at (11.5,7.5) [ ] {$B^{j}_{l_2}$};

    \draw[step=1, blue!50, very thick] (10, 1) grid (11, 2);
    \node at (10.5,1.5) [ ] {$B^{j}_{l_3}$};

    \draw[step=1, blue!50, very thick] (14, 2) grid (15, 3);
    \node at (14.5,2.5) [ ] {$B^{j}_{l_4}$};

    \draw[step=1, blue!50, very thick] (16, 4) grid (17, 5);
    \node at (16.5,4.5) [ ] {$B^{j}_{l_5}$};

    \draw[step=1, very thick] (13, 4) grid (14, 5);
    \node at (13.5,4.5) [ ] {$B^{i}_k$};

    \draw[very thick] (9,  0) -- (18, 0);
    \draw[very thick] (18, 0) -- (18, 9);
    \draw[very thick] (18, 9) -- (9, 9);
    \draw[very thick] (9, 9) -- (9, 0);

    \node at (13.5,-1) [] {$B^j_{l_1}, \dots, B^j_{l_5} \subset
      L^{i,\tau}_k$; $\tau = 2$.};

    \node at (9.5,8.5) [] {$L^{i,\tau}_k$};
   
\end{scope}
\end{tikzpicture}
\end{center}

\caption{Example of finding all the boxes $B^j_{l}$ that are
contained in $L^{i,\tau}_k$ for $\tau = 2$.}
\label{multilevelcov:fig5}
\end{figure}

With this criterion we can produce a highly sparse matrix
$\tilde{\bC}_{W}(\btheta)$ that is close to $\bC_{W}(\btheta)$ in the
matrix 2-norm sense.  

\begin{rem}
The error $\| \tilde{\bC}_{W}(\btheta) - \bC_{W}(\btheta)\|_{2}$ will
be monotonically decreasing with respect to the smoothness of the
covariance function, the size of the degree of the multi-level basis
$\tilde f \geq f$ (accuracy parameter $\tilde{p} \geq p$) and the size
of $\tau$.  For a sufficiently large $\tau$ and $\tilde f$ the error
$\| \tilde{\bC}_{W}(\btheta) - \bC_{W}(\btheta)\|_{2}$ will be small
and the matrix $\tilde{\bC}_{W}(\btheta)$ becomes positive definite.

Now, the number of nonzeros of $\tilde{\bC}_{W}(\btheta)$ will
increase as we increase $\tau$ and $\tilde p$. To be able to
determine the size for $\tau$ and $\tilde p \geq p$ it is
helpful to derive an expression for the error $\|
\tilde{\bC}_{W}(\btheta) - \bC_{W}(\btheta) \|_{2}$ vs the number of
non zeros of $\tilde{\bC}_{W}(\btheta)$.

Error estimates can be derived for $\| \tilde{\bC}_{W}(\btheta) -
\bC_{W}(\btheta) \|_{2}$ with respect to the smoothness of the
covariance function, $\tilde{p}$ and $\tau$, but this is beyond the
scope of the present paper. In practice for the polynomial based model
${\cal Q}^d_{\tilde{f}}$ we set the level dependent criterion
parameter $\tau := 1$ and increase $\tilde{f}$ (and $\tilde{p}$) until
at least $\tilde{\bC}_{W}(\btheta)$ is positive definite. Moreover,
the sparse Cholesky factorization code in the Suite Sparse package
(\cite{Chen2008,Davis2009,Davis2005,Davis2001,Davis1999}) that is used
in this paper informs the user if the matrix is not positive definite.
\end{rem}

In \cite{Castrillon2013} the authors described how to apply a Kernel
Independent Fast Multipole Method (KIFMM) by \cite{ying2004} to compute
all the diagonal blocks $\tilde{\bC}^{i,i}_W(\btheta)$ for $i =
0,\dots,t$ in ${\cal O}(nt)$ computational steps to a fixed accuracy
$\varepsilon_{FM} > 0$. This approach can be easily extended to
compute all the blocks $\tilde{\bC}^{i,j}_W(\btheta)$ for $i,j =
-1,\dots,t$ in ${\cal O}(n(t + 1)^2)$.

\begin{rem}
The KIFMM by \cite{ying2004} is very flexible as it allows a large
class of covariance functions to be used including the exponential,
Gaussian and Mat\'{e}rn.  However, the computational efficiency is
mostly dependent on the implementation of the covariance function
(since the KIFMM computational cost is ${\cal O}(n)$) and the accuracy
parameter of the solver. For all the numerical experiments in this
paper the accuracy parameter is set to medium ($10^{-6}$ to $10^{-8}$)
or high ($10^{-8}$ or higher).

Due to the lack of a fast math C++ library for the Mat\'{e}rn
covariance function, we create a Hermite cubic spline interpolant of
the covariance function with the multithreaded Intel Math Kernel
Library (MKL) data fitting package. To generate a compact
representation of the interpolant we implement an $h$-adaptive mesh
generator in 1D such that the absolute error over the range $(0, 2.5]$
  is less than TOL. From \cite{Elden2004} given that the covariance
  function $\phi(r;\btheta) \in C^{4}(\R)$, $r \in \R$, on each mesh
  element (starting at $x_0 \in \R$) with length $h$ we can guarantee
  that the absolute error for the cubic Hermite interpolant is less
  than TOL if $\frac{h^4}{384} \max_{x \in
    [x_{0},x_{0}+h]}{\phi^{(4)}(x;\btheta)} < TOL$, where $\phi^{(4)}$
  refers to the fourth derivative with respect to $x$.  In this work
  we set $TOL = 5 \times 10^{-9}$.  Numerical test confirmed $TOL$
  accuracy for the Mat\'{e}rn covariance function with less than 200
  adaptive mesh nodes. This is sufficient for the numerical examples
  in this paper.
\label{multilevelcov:rem1}
\end{rem}


In Figure \ref{multilevelcov:fig4} we show an example of a sparse
matrix produced for $\tau=1$ for $n = 8,000$ observation
locations sampled from a uniform distribution on the unit cube.
Notice that the entries of the matrix that are not covered by the
sparsity pattern are around $10^{-7}$ times smaller, implying the
hierarchical sparsity technique will lead to good accuracy.

The total sparsity for this example is 46\% (23\% since the matrix
$\tilde{\bC}_W(\btheta)$ is symmetric), however, the sparsity density
improves significantly as $n$ increases as we expect the number of
non-zero entries of $\tilde{\bC}_W(\btheta)$ to increase at most as
${\cal O}((t+2)^2n)$ with the number of observations $n$ (See 
\cite{Castrillon2013}). 

In Figure \ref{multilevelcov:fig6} the sparsity pattern of the matrix
$\tilde{\bC}_W(\btheta)$ is shown for $n = 64,000$ observation
locations sampled from a uniform distribution on the unit cube.  For
this case the design matrix $\bM_f$ is constructed from $p = 20$
monomials (i.e. up to cubic polynomials) and $\tau = 1$.  The sparsity
of this example is $8.2\%$ (4.1 \% since the matrix
$\tilde{\bC}_W(\btheta)$ is symmetric).


\begin{figure}[b!]
\begin{center}
\begin{tabular}{cc}
   \includegraphics[ width=3.0in,
     height=3in ]{./HBCovarianceExp1-1} &
   \includegraphics[ width=3.0in, height=3in
   ]{./CoverKw}\\ (a) $\log_{10}{ abs( \bC_W(\btheta) )}$ & (b)
   Overlayed sparsity pattern in blue  \\
\end{tabular} 
\end{center}
  \caption{Sparsity pattern overlayed on $\bC_W(\btheta)$. Notice that
    most of the entries that are not covered by the blue boxes are around
    $10^{-7}$ times smaller in magnitude than the covered
    entries.}
\label{multilevelcov:fig4}
\end{figure}


\begin{figure}[b!]
\begin{center}
   \includegraphics[ width=2.7in,
     height=3.3in ]{./Sparse64p3}
\end{center}
  \caption{Sparsity pattern (8.2 \% non zeros) for $\tilde
    \bC_W$ with $\tau = 1$ and $n = 64,000$.}
\label{multilevelcov:fig6}
\end{figure}

%% \begin{rem}  Suppose that we are interested in obtaining a higher
%% degree decay $\tilde{f} > f$, for example, when $f$ is small. This
%% is obtained by constructing the design matrix $\bM_{\tilde{f}}$ and
%% assigning ${\cal P}^{\tilde{p}}(\mathbb{S})$ (where $\tilde{p} > p$)
%% as the span of its columns. From Remark
%% \ref{MultiLevel:Construction:remark} we are able to construct a
%% multi-level basis such that \eqref{multilevelcov:eqn1} is satisfied
%% with respect to the degree $\tilde{f}$.
%% \end{rem}

\section{Multi-Level Estimator}
\label{multilevelestimator}

As the result section shows it is not necessary to compute the entire
sparse matrix $\bC_W(\btheta)$ to obtain a good estimate of the
covariance parameter $\btheta$. Due to the multi-resolution properties
of the MB we can construct a partial multi-resolution likelihood
function that is effective.
%% By constructing the projector operators $\tilde{\bL}_i:\R^n
%% \rightarrow {\cal P}^{p}(\mathbb{S}) \oplus W_{-1}(\mathbb{S}) \oplus
%% W_{0}(\mathbb{S}) \oplus \dots W_{i-1}(\mathbb{S})$ and
%% $\tilde{\bW}_{i }: \R^n \rightarrow W_{i }(\mathbb{S}) \oplus\dots
%% \oplus W_{t}(\mathbb{S})$, $i = -1, \dots, t$, by using the procedure
%% in Section \ref{MultiLevel:Construction}

We can produce a series of multi-resolution likelihood functions
$\tilde{\ell}^i_{W}(\btheta)$, $i = -1,\dots t$ by applying the
partial transform $[\bW_t \T, \dots, \bW_i \T ]$ to the data $\bZ$,
thus
\begin{equation}
\tilde{\ell}^i_{W}(\btheta)
=-\frac{\tilde{n}}{2}\log(2\pi)-\frac{1}{2}\log
\det\{\tilde{\bC}_{W}^i(\btheta)\}
-\frac{1}{2}(\tilde{\bZ}^i_{W})\T\tilde{\bC}^i_{W}(\btheta)^{-1}\tilde{\bZ}^i_{W},
\label{Introduction:multilevelloglikelihoodreduced2}
\end{equation}
where $\tilde{\bZ}^i_W :=[\bW_t \T, \dots, \bW_i \T ] \T \bZ$,
$\tilde{n}$ is the length of $\tilde{\bZ}^i_W$ and
$\tilde{\bC}^i_{W}(\btheta)$ is the $\tilde{n} \times \tilde{n}$
upper-left submatrix of $\tilde{\bC}_{W}(\btheta)$.

% $\tilde{\bC}^i_{W}(\btheta):=\mathbb{E}[\tilde{\bZ}^i_W (\tilde{\bZ}_W^i)
%   \T]$.

\subsection{Computation of $\log{\det\{ \tilde{\bC}^i_{W}(\btheta)\}}$}

Since $\bC(\btheta)$ is symmetric positive definite from
\eqref{Introduction:eqn1} it can be shown that $\bC_{W}(\btheta)$ is
also symmetric positive definite.  It can also be shown that for a
sufficiently large $\tau$ and/or $\tilde{p}$ the matrix
$\tilde{\bC}^i_{W}(\btheta)$ will also be symmetric positive definite.
An approach to computing the determinant of
$\tilde{\bC}^i_{W}(\btheta)$ is to apply a sparse Cholesky
factorization technique such that $\bG\bG\T =
\tilde{\bC}^i_{W}(\btheta)$ where $\bG$ is a lower triangular
matrix. Since the eigenvalues of $\bG$ are located on the diagonal we
have that $\log \det \{\tilde{\bC}^i_{W}(\btheta)\} = 2 \sum_{i =
  1}^{\tilde{n}} \log{\bG_{ii}}$.

To reduce the fill-in of the factorization matrix $\bG$ we
apply the matrix reordering technique in Suite Sparse 4.2.1 package
(\cite{Chen2008,Davis2009,Davis2005,Davis2001,Davis1999}) 
% (\cite{Davis2009}, \cite{Chen2008}
% \cite{Davis2005},\cite{Davis2001},\cite{Davis1999})
with the Nested Dissection (NESDIS) function
% and/or the METIS 4.0.1 (\cite{Karpis1998})
package. The sparse Cholesky factorization is performed with the {\it lchol}
command from Suite Sparse 4.2.1 package.

Although in practice the combined NESDIS and sparse Cholesky
factorization is highly efficient, as shown by our numerical results,
a worse case complexity bound can be obtained.  For example, it can be
shown that by using the planar graph separation theorem (see
\cite{George1973}, \cite{Gilbert1987}) a worse case complexity of ${\cal
  O}(n^{3/2})$ and ${\cal O}(n (\log{n})^2 )$ storage is achieved in
2D. Similarly, the worse case complexity in 3D is ${\cal O}(n^{2})$.


\begin{exam}
Continuing Example \ref{multilevelcov:example1} we compute $\log{\det
  \{ \bC_{W}(\btheta)\}}$ and the approximation \linebreak
$\log{\det\{{\tilde{\bC}}_{W}(\btheta)\}}$ for $\tau = 0,1,2,\infty$
by applying the sparse Cholesky factorization. In
Table~\ref{multilevelcov:table1} we tabulated the absolute
$\varepsilon_{abs}:=|\log{\det\{{\tilde{\bC}}_{W}(\btheta)\}} -
\log{\det\{{\bC}_{W}(\btheta)\}}| $ and relative $\varepsilon_{rel}
:=\frac{\varepsilon_{abs}} { |\log \det \{\bC_{W}(\btheta)\}|} $
errors.  For $\tau = 0$ we obtain a very sparse matrix (4\% density),
but leads to a non-positive definite matrix, which is not valid for
the computation of the determinant. For $\tau = 1$ the matrix
$\tilde{\bC}_W(\btheta)$ becomes positive definite. As we increase
$\tau$ the approximation $\log{\det\{{\tilde{\bC}}_{W}(\btheta)\}}$
becomes more accurate. However, the density of
$\tilde{\bC}_{W}(\btheta)$ also increases.

\begin{table}[htp]
\caption{Log determinant errors comparisons between
   $\tilde{\bC}_{W}(\btheta)$ and ${\bC}_{W}(\btheta)$.}
\begin{center}
\begin{tabular}{c c c c c}
$\tau$ & density (\%) & $\log{\det\{{\tilde{\bC}}_{W}(\btheta)\}}$ &
$\varepsilon_{abs}$ 
& 
$\varepsilon_{rel}
$
\\ 
\hline
0 & $4$   & not positive definite & -- & -- \\
1 & $23$ & $-5.184354 \times 10^3$ & $2.23 \times 10^{-2}$ & $4.31
\times 10^{-6}$ \\
2 & $38$ & $-5.184332 \times 10^3$ & $3.78 \times 10^{-4}$ & 
$7.29 \times 10^{-8}$\\
$\infty$ &  $50$  & $-5.184331 \times 10^3$ & 0 & 0
\end{tabular}
\end{center}
\label{multilevelcov:table1}
\end{table}
\end{exam}

\subsection{Computation of $(\tilde{\bZ}^i_W)\T \tilde{\bC}^i_{W}(\btheta)^{-1}\tilde{\bZ}^i_{W}$}
\label{mutlilevelcomputation}

We have two choices for the computation of $(\tilde{\bZ}^i_W) \T \tilde
\bC^i_{W}(\btheta)^{-1} \tilde \bZ^i_{W}$. We can either use a Cholesky
factorization of $\tilde{\bC}^i_{W}(\btheta)$ or Preconditioned
Conjugate Gradient (PCG) coupled with a KIFMM.  The PCG choice
requires significantly less memory and allows more control of the
error. \corb{However, we already computed the sparse Cholesky factorization
of $\tilde{\bC}^i_{W}(\btheta)$ for the computation of the
determinant}. Thus we can use the same factors to compute $
(\tilde \bZ_{W}^i) \T\tilde{\bC}^i_{W}(\btheta)^{-1} \tilde \bZ^i_{W}$.

%% The alternative PCG choice requires significantly less memory and
%% allows more control of the error. The numerical computation of
%% $\bZ\T_{W}\tilde{\bC}^i_{W}(\btheta)^{-1}\bZ_{W}$ can be performed by
%% solving the following system of equations
%% \begin{equation}
%% \bar{\bC}_{W}(\btheta) \bar{\balpha}_W (\btheta) = \bar{\bZ}_{W},
%% \label{multilevelcov:system}
%% \end{equation}
%% where $\bar{\bC}_{W}(\btheta) := \bD_W^{-1} (\btheta) \tilde{\bC}^i_{W}(\btheta)
%% \bD_W^{-1}(\btheta)$, $\bar{\balpha}_W (\btheta):= \bD_W(\btheta)
%% \balpha_W (\btheta)$, $\bar{\bZ}_{W} := \bD_W^{-1} (\btheta) \bZ_{W}$,
%% and $\bD_W(\btheta) := \diag\{\tilde{\bC}^i_{W}(\btheta)\}$ is a diagonal
%% preconditioner. The final computation is then simply
%% $\bZ\T_{W}\tilde{\bC}^i_{W}(\btheta)^{-1}\bZ_{W} = \bZ\T_W \balpha_W$.  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Multi-Level Kriging}
\label{multilevelkriging}

An alternative formulation for obtaining the estimate $\hat \bZ(\bs_0)$ is
by solving the following problem
\begin{eqnarray}
\left( {{\begin{array}{*{20}c}
 \bC(\btheta) \hfill & \bM_f \hfill \\
 \bM_f\T \hfill & \0 \hfill \\
\end{array} }} \right)\left( {{\begin{array}{*{20}c}
 \hat \bgamma \hfill \\
 \hat \bbeta \hfill \\
\end{array} }} \right)=\left( {{\begin{array}{*{20}c}
 \bZ \hfill \\
 \0 \hfill \\
\end{array} }} \right).
\label{kriging:problem}
\end{eqnarray}
\corb{It is not hard to show that the solution of this problem
leads to equation \eqref{GLSbeta} and $\hat \bgamma(\btheta) =
\bC^{-1}(\btheta)\{\bZ - \bM_f \hat \hat \bbeta(\btheta)\}$ or
alternatively $\hat{\bbeta} = (\bM\T_f\bM_f)^{-1}\bM_f(\bZ -
\bC(\btheta)\hat{\bgamma})$.} The best unbiased predictor is
evaluated as
\begin{equation}
  \hat Z(\bs_0)
  =\bm(\bs_0)\T\hat \bbeta(\btheta)+\bc(\btheta)\T
  \hat \bgamma(\btheta)
\label{Kriging}
\end{equation}
and the Mean Squared Error (MSE) at the target point $\bs_0$ is given by
\begin{equation}
1 + 
\tilde{\bu}
\T
(\bM_f \T 
\bC(\btheta)^{-1} \bM_f )^{-1}
\tilde{\bu}
-\bc(\btheta)\T\bC^{-1}(\btheta)\bc(\btheta)
\label{Kriging:MSE}
\end{equation}
where $\tilde{\bu}\T := (\bM_f \bC^{-1}(\btheta)
\bc(\btheta) - \bm(\bs_0))$.

The computational cost for computing $\hat \bbeta(\btheta)$, $\hat
\bgamma(\btheta)$ and the MSE accurately using a direct method is
${\cal O}(n^3)$, which is unfeasible for large size problems. We
propose a much faster approach.

From \eqref{kriging:problem} we observe that $\bM_f\T \hat
\bgamma(\btheta) = \0$. This implies that $\hat{\bgamma} \in \R^{n}
\backslash {\cal P}^{p}(\mathbb{S})$ and can be uniquely rewritten as
$\hat{\bgamma} = \bW\T \bgamma_W$ for some $\bgamma_W \in
\R^{n-p}$. We can rewrite $\bC(\btheta) \hat \bgamma + \bM_f \hat
\bbeta = \bZ$ as  
\begin{equation}
\bC(\btheta) \bW\T \hat \bgamma_{W} + \bM_f \hat \bbeta =
       \bZ.
\label{Kriging:eqn1}
\end{equation}
Now apply the matrix $\bW$ to equation \eqref{Kriging:eqn1} and we
obtain $ \bW \{\bC(\btheta) \bW\T \hat \bgamma_{W} + \bM_f \hat
\bbeta\} = \bW \bZ.$ Since $\bW \bM_f = \0$ then $\bC_{W}(\btheta)
\hat \bgamma_{W} = \bZ_{W}$. Applying the preconditioner $\bD^{-1}_W
(\btheta)$, where $\bD_W (\btheta) := \diag(\bC_W (\btheta))$,
we have the system of equations $\bar{\bC}_{W}(\btheta)
\bar{\bgamma}_W (\btheta) = \bar{\bZ}_{W}$ where $\bar{\bgamma}_W
(\btheta) := \bD_W(\btheta) \bgamma_W (\btheta)$,
$\bar{\bC}_{W}(\btheta) := \bD^{-1}_W(\btheta) \bC_{W}(\btheta)
\bD^{-1}_W(\btheta)$ and $\bar{\bZ}_{W} := \bD^{-1}_W(\btheta)
\bZ_{W}$.


This system of equations is solved by a combination of a KIFMM and
PCG.  If $\bC_W(\btheta)$ and $\bD_W(\btheta)$ are symmetric positive
definite then an effective method to solve $\bar{\bC}_{W}(\btheta)
\bar{\bgamma}_W (\btheta) = \bar{\bZ}_{W}$ is the PCG method
implemented in PETSc by 
\cite{petsc-web-page,petsc-user-ref,petsc-efficient}. 

\begin{lem}
If the covariance function $\phi$ is positive definite, then the
matrix $\bD_W(\btheta)$ is always symmetric positive definite.
\label{multilevelkriging:lemma2}
\end{lem}

The matrix-vector products $\bC_W(\btheta) \bv$, where $\bv \in
\R^{n-p}$, are computed in ${\cal O}(n)$ computational steps to a
fixed accuracy $\varepsilon_{FM} > 0$.  The total computational cost
is ${\cal O}(kn(t+2))$, where $k$ is the number of iterations needed
to solve $\bar{\bC}_{W}(\btheta) \bar{\bgamma}_W (\btheta) =
\bar{\bZ}_{W}$ to a predetermined accuracy $\varepsilon_{PCG} > 0$.

It is important to point out that the introduction of a
preconditioner can degrade the performance of the PCG, in particular,
if the preconditioner is ill-conditioned. The accuracy of the PCG
method $\varepsilon_{PCG}$ has to be set such that the accuracy of the
{\it unpreconditioned} system $\bC_{W}(\btheta) \balpha_W (\btheta) =
\bZ_{W}$ is below a user given tolerance $\varepsilon > 0$.

\corb{We compute $\hat \bgamma = \bW\T \bgamma_W$ and $\hat{\bbeta} =
  (\bM\T_f\bM_f)^{-1}\bM_f \T (\bZ - \bC(\btheta)\hat{\bgamma})$ in at
  most ${\cal O}(np^2 + p^{3})$ computational steps. The matrix vector product
  $\bc(\btheta)\T \hat{\bgamma}(\btheta)$ is computed in ${\cal O}(n)$
  steps. Finally, the total cost for computing the estimate
  $\hat{\bZ}(\bs_0)$ from \eqref{Kriging} is ${\cal O}(np^2 + p^{3} +  kn(t+2))$}.

In Appendix C we show a procedure to compute the MSE fast.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Numerical Study and Statistical Examples}
\label{numericalstudy}

In this section we test the numerical efficiency and accuracy of our
solver for computing the terms $\log{\det
  \{\tilde{\bC}^{i}_{W}(\btheta)\}}$ and $\bar{\bgamma}_W
(\hat{\btheta}) = \bar{\bC}^{-1}_{W}(\hat{\btheta}) \bar{\bZ}_{W}$ for
Mat\'ern covariances. Our results show that we are able to solve
problems of up to 128,000 observations and kriging up to 512,000 size
problems with good accuracy. Our approach is not limited to 128,000
for parameter estimation. This was the maximum we could test due to
the memory limitation on our workstation in creating observations
larger than 128,000. We now describe the data sets.

{\bf Data set \#1 and \#2:} The sets of observation locations
$\bS^d_{1}, \dots, \bS^d_{10}$ vary from 1,000 to 512,000 and we
assume that $\bS^d_{l} \subset \bS^d_{l+1}$ for $l = 1, \dots, 9$ for
$d = 2$ and $d = 3$.  The observations locations are sampled from a
uniform distribution over the unit square $[0,1]^2$ for $d = 2$ (data
set \#1) and for $[0,1]^3$ for $d = 3$ (data set \#2), as shown in
Figure \ref{SRBF}. The target points $\bs_0$ are set to 1000 random
points across the domain $[0,1]^2$ (data set \#1) and $[0,1]^3$(data
set \#2). We shall refer to $\bZ^d_1, \dots \bZ^d_{10}$ as the observation
values associated with $\bS^d_1, \dots \bS^d_{10}$.

{\bf Data set \#3:} We take the data set generated by $\bS^d_{9}$ for
$d = 2$ (256,000 observation points) and carve out two disks located
at (1/4,1/4) and (3/4,3/4) with radii 1/4. This generates 100,637
observation points; see Figure \ref{SRBF}(c) for an example with 1,562
observation points randomly extracted from the data set.

We now test our approach on the Mat\'{e}rn covariance function
$\phi(r;\btheta):=\frac{1}{\Gamma(\nu)2^{\nu-1}} \left(
\sqrt{2\nu}\frac{r}{\rho} \right)^{\nu} $ $K_{\nu} \left(
\sqrt{2\nu}\frac{r}{\rho} \right)$, where $\Gamma$ is the gamma
function and $K_{\nu}$ is the modified Bessel function of the second
kind. All results are executed on a single CPU (4 core Intel i7-3770
CPU @ 3.40GHz.) with Linux Ubuntu 13.04.

\begin{figure}[htpb]
\begin{center}
\begin{tabular}{c c c}
\includegraphics[width=2in, height=2in]{./DataSet1.eps} &
\includegraphics[width=2in, height=2in]{./DataSet2.eps} &
\includegraphics[width=2in, height=2in]{./DataSet3.eps} 
\\
(a) & (b) & (c)
\end{tabular}
\end{center}
\caption{Data set examples: (a) Data set \#1: One thousand observation
  points randomly generated from a uniform distribution on $[0,1]^2$.
  (b) Data set \#2: One thousand observation points randomly generated
  from a uniform distribution on $[0,1]^3$.  The color of the
  observation locations represent the distance along the right axis
  coordinate. (c) Data set \#3: Two disks of 1562 randomly generated
  observation locations. The disks are contained in a square
  but are represented in a multi-level representation.}
\label{SRBF}
\end{figure}

\vspace{-.4cm}

\subsection{Parameter Estimation}

In this section we present the results for data set \#1 and \#3
for the Mat\'{e}rn parameter estimation.

Suppose we have two realizations of the Gaussian spatial random field
$Z(\bs)$ with the Mat\'{e}rn covariance for data set \#1 (2D).  We set
$f = 3$ ($p = 10$) and $\tilde{f} = 4,5,6$ (corresponding to
$\tilde{p} = 15,21,28$) and fix the covariance parameters to $\btheta=(\nu,\rho)
= (3/4,1/6)$.  Two realizations $\bZ^2_6$ ($n = 64,000$) and $\bZ^2_7$
($n = 128,000$) are generated from these parameters. For each
observation values $\bZ^2_{6}$ and $\bZ^2_{7}$ (and locations) apply
the transformation $\bW$ to compute $\bZ^2_{W,6}$ and $\bZ^2_{W,7}$
and solve the optimization problems $\hat{\btheta}_{i} : =
  \argmax_{\btheta \in (0,\infty) \times (0,\infty)}
  \hat{\ell}_{W}(\bZ^2_{W,i},\btheta)$ for $i = 6$ and $7$.

\begin{table}[b!]
\caption{Estimation results for data sets \#1 and \#3.  The
  observation data is generated with covariance parameters $\nu=0.75$
  and $\rho=1/6$.  The degree of the model is $f = 3$ (cubic), which
  gives $p = 10$ monomials and we set $\tau :=1$.  \corb{ For all the
    experimental runs the number of MB levels is $t = 6$.}  The first
  column is the size of the problem. Columns 2 and 3 are the
  parameters that are chosen to build the MB and the multi-level
  estimator.  \corb{Column 4 is the corresponding $\tilde p$ for
    $\tilde f$ and the maximum level.  Columns 5 and 6 are the errors
    of $\hat \nu - \nu$ and $\hat \rho - \rho$.  Column 7 is the
    percentage of non-zeros of the Cholesky factors. Column 8 is
    self-explanatory.  Column 9 and 10 are the approximate wall clock
    computational time (in seconds) needed to compute
    $\tilde{\bC}^{i}_W(\btheta)$ and to perform Cholesky factorization
    respectively.} The total time for each Newton iteration is about
  $t_{cons} (s) + t_{chol} (s)$.  For each problem it takes about 50
  Newton iterations to converge.  \corb{Estimation results for data
    set \#3 with observation data generated with covariance parameters
    $\nu=1.25$ and $\rho=1/6$.}  }
\begin{center}
Estimation results for data set \#1 (2D). 
\\
\corb{
\begin{tabular}{ r r r | r || r r r r r r r}
\multicolumn{1}{c}{$n$} &
\multicolumn{1}{c}{$\tilde f$} &
\multicolumn{1}{c}{$i$} &
\multicolumn{1}{c}{$\tilde p$} & 
% \multicolumn{1}{c}{$t$} &
\multicolumn{1}{c}{$\hat{\nu} - \nu$} &
\multicolumn{1}{c}{$\hat{\rho} - \rho$} &
\multicolumn{1}{c}{$nz(\bG) (\%)$} &
\multicolumn{1}{c}{size($\tilde{\bC}^i_W$)} &
\multicolumn{1}{c}{$t_{cons} (s)$} &
\multicolumn{1}{c}{$t_{chol} (s)$} &
 \\ 
%32,000 &  5 &  5    & 15 &  0.7354 & 0.1809  &  1.6 & 16640 &    8 &   0 \\
%32,000 &  5 &  4    & 15 &  0.7490 & 0.1689  &  2.4 & 28160 &   45 & 5 \\
%32,000 &  5 &  3    & 15 &  0.7514 & 0.1672  &  6.4 & 31040 &  155 & 5 \\
%32,000 &  - &  2    & 21 &  0.7501 & 0.1682  & 18.1 & 31664 &  660 & 24  \\
\hline \hline
64,000 &  6 &  6 & 28 &  -0.0759 &  0.0333 & 8.9 &   23  &  14 & 0  \\
64,000 &  6 &  5 & 28 &   0.0182 & -0.0132 & 1.7 & 35328 &  40 & 1  \\
64,000 &  6 &  4 & 28 &  -0.0043 &  0.0046 & 4.5 & 56832 & 230 &11  \\
64,000 &  6 &  3 & 28 &  -0.0049 &  0.0048 & 10.7 & 62208 & 961 &65  \\
\hline
64,000 &  5 &  6 & 21 &  0.0071 &  -0.0146 & 0.4 &   810 &  13 & 0  \\
64,000 &  5 &  5 & 21 &  0.0037 &  -0.0027 & 1.7 & 42496 &  43 & 2  \\
64,000 &  5 &  4 & 21 & -0.0030 &   0.0048 & 3.7 & 58624 & 220 &12  \\
64,000 &  5 &  3 & 21 & -0.0048 &   0.0046 & 7.0 & 62656 & 750 &32  \\
\hline
64,000 &  4 &  6 & 15 & -0.0080 & 0.0098 & 0.2 &  7749 &  13 & 0  \\
64,000 &  4 &  5 & 15 & -0.0047 & 0.0043 & 2.0 & 48640 &  53 & 3  \\
64,000 &  4 &  4 & 15 & -0.0068 & 0.0062 & 3.4 & 60160 & 161 & 9  \\
64,000 &  4 &  3 & 15 & -0.0051 & 0.0048 & 4.4 & 63040 & 550 &15  \\
%\hline
\hline \hline 
128,000 &  6 &  6 & 28 &  0.0010 & -0.0011 & 0.3 & 17179 & 75 & 0 \\
%Total  8600
128,000 &  6 &  5 & 28 &  0.0025 & -0.0020 & 2.1 & 99328 & 350 &13 \\
%Total 30000
128,000 &  6 &  4 & 28 & -0.0002 &  0.0005 & 4.0 &  120832 & 1200 & 70 \\
%Total 122000
\hline
128,000 &  5 &  6 & 21 & -0.0010 &  0.0015 & 0.5 &   42154 &   80 &  0 \\
%Total 9595
128,000 &  5 &  5 & 21 &  0.0004 & -0.0002 & 1.9 &  106496 &  300 & 14 \\
128,000 &  5 &  4 & 21 & -0.0016 &  0.0017 & 3.3 &  122624 & 1000 & 50 \\
%Total 100281
% \hline
\end{tabular} 
}
\\
\vspace{5mm}
Estimation results for data set \#3 (2D). \\
\corb{
\begin{tabular}{ r r r | r || r r r r r r r}
\multicolumn{1}{c}{$n$} & 
\multicolumn{1}{c}{$\tilde f$} &
\multicolumn{1}{c}{$i$} &
\multicolumn{1}{c}{$\tilde{p}$} &
% \multicolumn{1}{c}{$t$} &
\multicolumn{1}{c}{$\hat{\nu} - \nu$} &
\multicolumn{1}{c}{$\hat{\rho} - \rho $} &
\multicolumn{1}{c}{$nz(\bG) (\%)$} &
\multicolumn{1}{c}{size($\tilde{\bC}^i_W$)} &
\multicolumn{1}{c}{$t_{cons} (s)$} &
\multicolumn{1}{c}{$t_{chol} (s)$} &
 \\ 
\hline
100,637 &  6 &  6  & 66 &  0.0548 & -0.0237 & 0.5 &  2613 &  60 &  0 \\
%Total time = 4337.433602 
100,637 &  6 &  5  & 66 & -0.0031 & 0.0020 & 3.1 & 72231 & 600 & 12 \\
%Total  83452
\end{tabular}
}
\vspace{5mm}
\\


\end{center}
\label{SRBF:table6}
\end{table}
The optimization problem from the log-likelihood function
\eqref{Introduction:multilevelloglikelihood} is solved using {\it
  fminsearch} from the optimization toolbox in MATLAB with the local
minimizer search for $\nu$ in the interval $[1/2,5/4]$ and $\rho$ in
the interval $[1/7,1/5]$.  We set the parameter criterion to $\tau :=
1$, and the {\it fminsearch} tolerance is set to $10^{-3}$.  In Table
\ref{SRBF:table6} we tabulate the results for the parameter estimates
$\hat \nu$ and $\hat \rho$ for different problem sizes of data sets
\#1 and \#3 for the user defined parameters $\tilde f$ for the
construction of the MB and $i$ for the construction of the matrix
$\tilde \bC^{i}_W$. We notice that the estimates of $(\nu,\rho)$ seem to
approach the actual values as we increase the number of
observations. In particular, for $n = 128,000$ the estimate
$(\hat{\nu},\hat{\rho}) := (0.7498,0.1672)$ is very close to the
actual noise model parameters $(\nu, \rho) = (0.75,1/6)$ of the
covariance function. We observe that as we increase the number of levels
(i.e. decrease $i$) in the covariance matrix the absolute error decays
until it stagnates, usually by the time that the covariance matrix is
for two levels.  We also report the wall clock times (i.e. actual time
it took the executable to run, not to be confused with CPU time that
is unreliable as a measure) for computing each Newton iteration.  The
total number of Newton iterations is approximately~50.

In Table \ref{SRBF:table6} the results for parameter estimation
with data set \#3 are tabulated.  The realization is obtained from the
Gaussian random field $Z(\bs)$ with $n = 100,637$, $\nu = 1.25$ and
$\rho = 1/6$. For this case the absolute error is 0.25\% for the
estimate $\hat{\nu}$ and 0.04\% for $\hat{\rho}$.

\begin{table}[b!]
\caption{ (a) Statistical results for data set \#1 with observation
  data generated with covariance parameters $\nu=0.75$ and $\rho=1/6$.
  \corb{The parameter $M = 100$ is the number of realizations of the
    stochastic model, $\bbE_M [\hat{\nu} - \nu]$ is the bias of $M$
    estimates of $\nu$ and $std_M [\hat{\nu}]$ is the standard
    deviation of $M$ estimates of $\nu$.}  As $i$ is reduced from $t$
  to $t-1$ there is a significant drop in $std_M
  [\hat{\nu}]$. However, for $i < t - 1$ the standard deviation of the
  estimates does not improve significantly. This indicates that a good
  estimate can be obtained for $i = t-1$ and there is not much gain in
  reducing $i$, which increases the computational cost in computing
  $\tilde{\bC}^{i}_W$. }

\begin{center}
Statistical results for data set \#1 with multiple realizations.
\\
\corb{
\begin{tabular}{ r r r | r r || r r r r r r}
\multicolumn{1}{c}{$n$} & 
\multicolumn{1}{c}{$\tilde f$} &
\multicolumn{1}{c}{$i$} &
\multicolumn{1}{c}{$\tilde{p}$} &
\multicolumn{1}{c}{$t$} &
% \multicolumn{1}{c}{$M$} &
\multicolumn{1}{c}{$\bbE_M [\hat{\nu} - \nu]$} &
\multicolumn{1}{c}{$\bbE_M [\hat{\rho} - \rho]$} &
\multicolumn{1}{c}{$std_M [\hat{\nu}]$} &
\multicolumn{1}{c}{$std_M [\hat{\rho}]$} &
 \\ 
\hline
\hline
32,000 &  4 &  5  & 15 &  5 & $-6.0 \times 10^{-4}$  
& 
$  1.0 \times 10^{-3}$  
& $1.3 \times 10^{-2}$ 
& $1.0 \times 10^{-2}$ \\
% 0.749396757600804   0.167667787161979
% 0.012810986650296   0.010443579627680
32,000 &  4 &  4  & 15 & 5 & $-7.2 \times 10^{-4}$  
& $7.0 \times 10^{-4}$ 
& $5.9  \times 10^{-3}$ & 
$4.5 \times 10^{-3}$
      \\
%   0.749282066112105   0.167362223586580
%   0.005949633428934   0.004481192755942
32,000 &  4 &  3  & 15 & 5 & $-7.0 \times 10^{-4}$  
& $6.0 \times 10^{-4}$ 
&
$5.6  \times 10^{-3}$ &
$4.0   \times 10^{-3}$ \\
% 0.749304921099859   0.167269740204625
% 0.005616655940988   0.004026141142878
\hline
64,000 &  4 &  6  & 15 & 6 & $6.8 \times 10^{-4}$ 
& $1.1 \times 10^{-3}$ 
& $2.0  \times 10^{-2}$ &
$1.9   \times 10^{-2}$ \\
% 0.750677755620135   0.167804754810837
% 0.019937400717640   0.018966124613043
64,000 &  4 &  5  & 15 & 6 & 
$7.4 \times 10^{-4}$ 
& $-5.6 \times 10^{-4}$ 
& $5.4  \times 10^{-3}$ &
 $4.6   \times 10^{-3}$ \\
% 0.750736889798586   0.166108160045356
% 0.005444589164041   0.004580786955817
64,000 &  4 &  4  & 15 & 6 & 
$2.5 \times 10^{-4}$ 
& $-1.7 \times 10^{-4}$ 
& $3.9  \times 10^{-3}$ &
 $3.3   \times 10^{-3}$ \\
% 0.750246306660986   0.166491721001726
% 0.003929097723537   0.003252657421859
\hline
128,000 &  6 &  6  & 28 & 6 &
$-1.3 \times 10^{-3}$ 
& $1.5 \times 10^{-3}$ 
& $8.3  \times 10^{-3}$ &
 $7.7   \times 10^{-3}$ \\
% 0.748733890355902   0.168171138625029
% 0.008331258397086   0.007695950309419
128,000 &  6 &  5  & 28 & 6 & 
$-6.2 \times 10^{-4}$ 
& $6.5 \times 10^{-4}$ 
& $3.7  \times 10^{-3}$ &
 $3.3   \times 10^{-3}$ \\
% 0.749379081773050   0.167320684155065
% 0.003667656022770   0.003306102738729
\end{tabular}
}
\end{center}
\label{SRBF:table7}
\end{table}

In Table \ref{SRBF:table7} we generate $M=100$ realizations of the
stochastic model for data set \#1, to analyze the mean and standard
deviation of the Mat\'{e}rn covariance parameter estimates. The mean
estimate $\bbE_M [\hat{\nu}]$ refers to the mean of $M$ estimates
$\hat{\nu}$ for the the $M$ realizations of the stochastic
model. Similarly, $std_M [\hat{\nu}]$ refers to the the standard
deviation of the $M$ realizations. We observe that the mean appears to
approach the covariance parameters as we decrease $i$. As $i$ is
reduced from $t$ to $t-1$ there is a significant drop in the term
$std_M [\hat{\nu}]$. However, for $i < t - 1$ the standard deviation
$std_M [\hat{\nu}]$ does not improve significantly. Therefore, there
is not much gain in improving the estimate by decreasing $i$, which
increases the computational cost in computing $\tilde{\bC}^{i}_W$.

\subsection{Numerical examples for computing $\bC_{W}(\btheta)^{-1}\bZ_{W}$
and Kriging}

We test our approach for solving the system of equations
$\bar{\bC}_{W}(\btheta)\bar{\bgamma}_{W} = \bar \bZ_{W}$ (that we have
to solve to obtain the kriging predictor) on the data sets \#1 and
\#2. We also include results showing the kriging prediction error
between the multi-level and direct methods.

We first test the PCG method with data set \#2 (3D) on three test
cases: (a) $\btheta_a = (\nu,\rho)= (3/4,1/6)$, (b) $\btheta_b =(1,1/6)$ and
(c) $\btheta_c = (5/4,1/6)$. The value $\rho = 1/6$ gives us an
approximate decay of $5 \%$ (which is reasonable in practice) from the
center of the cube along each dimensional axis.  The PCG relative
error tolerance $epsilon_{PCG} > 0$ is set to a value that leads
to a relative error $\epsilon = 10^{-3}$ of the {\it unpreconditioned}
system $\bC_{W}(\btheta) \hat \bgamma_{W} = \bZ_{W}$.

In Table \ref{SRBF:table1} we report the total wall clock times and
iterations for computing $\hat{\bbeta}$, $\hat{\bgamma}$ and the
target $\hat \bZ(\bs_0)$ for data set \#2 (3D) with the Mat\'{e}rn
covariance function . The polynomial accuracy of the model is set to
cubic ($f = 3$, $p = 20$) and the accuracy parameter $\tilde{p}$ is
set to 20 (which corresponds to $\tilde{f} = 3$).  We look at three
cases: For (a) ($\btheta_a = (3/4,1/6)$) we set the KIFMM accuracy
to medium and the number of iterations increase as ${\cal
  O}(n^{0.58})$. For (b) ($\btheta_b = (1,1/6)$) we set the KIFMM
accuracy to medium and the number of iterations increases as ${\cal
  O}(n^{0.58})$. For (c) ($\btheta_c = (5/4,1/6)$) we set the KIFMM
accuracy to high and the number of iterations increases as ${\cal
  O}(n^{0.77})$.

\corb{In Table \ref{SRBF:table1} we also report the number of
  iterations needed for solving $\bC^{-1}(\btheta)\bZ$ with $10^{-3}$
  accuracy with a CG method.} In this case the number of iterations is
about 10 times larger than the multi-level version.  Moreover, for
solving the kriging problem (e.g. equation \eqref{GLSbeta}), $p$ such
problems have to be solved. Thus, it is at least about 200 times
faster since $p = 20$ for this case. An alternative is to solve
\eqref{kriging:problem}. However, in general it is not positive
definite. The matrix is highly ill-conditioned also making it
difficult to solve with an iterative solver such as generalized
minimal residual method (see \cite{Castrillon2013}).

\begin{table}[t!]
\caption{Diagonal pre-conditioned results for computing
  $\bar{\bC}_{W}(\btheta) \bar{\bgamma}_{W} = \bar{\bZ}_{W}$ for data
  set \#2 (3D) with the Mat\'{e}rn covariance $\btheta =
  (\nu,\rho)$. We look at three cases: (a) $\btheta_a = (3/4,1/6)$ (b)
  $\btheta_b = (1,1/6)$ and (c) $\btheta_c = (5/4,1/6)$.  The relative
  error of the residual of the unpreconditioned system is set to
  $\varepsilon = 10^{-3}$. The KIFMM is set to medium accuracy for (a)
  and (b), and set to high accuracy for (c). The second column is the
  number of iterations needed to obtain $10^{-3}$ relative error of
  the unpreconditioned system with $\bC_W(\btheta)$. We denote as
  itr($\bC_W$) as the number of CG iterations needed for convergence
  until the desired accuracy is achieved.  The third column is the
  number of iterations for solving $\bC^{-1}(\btheta)\bZ$ with $10^{-3}$
  accuracy.  The fourth column is the residual tolerance needed for
  convergence of $10^{-3}$ relative error for the unpreconditioned
  system. The fifth column presents the wall clock times for
  initialization (basis construction and preconditioner
  computation). The PCG iteration wall clock times for $\bC_W$ are
  given in the sixth column. The last column presents the total wall
  clock time to compute $\bar{\bgamma}_{W} =
  \bar{\bC}_{W}(\btheta)^{-1} \bar \bZ_{W}$.}
  % The second column is the number of PCG iterations. The
  % third column is the residual tolerance needed for convergence of the
  % system $\bar{\bC}_{W}(\btheta)\bar{\bgamma}_{W} = \bar{\bZ}_{W}$.
  % The fourth column is the wall clock time for initialization and
  % computation of the diagonal pre-conditioner. The fifth column is the
  % wall clock time of executing the iteration. The last column is the
  % total wall clock time for solving $\hat{\bgamma}$, $\hat{\bbeta}$
  % and $\hat{\bZ}(\bs_0)$.}
\begin{center}
\corb{(a) $\btheta_a = (3/4,1/6)$, $d = 3$, $f = 3$ ($p = 20$), $\tilde{f} = 3$ ($\tilde{p} = 20$)} \\
\begin{tabular} { r r r r r r r }
  \multicolumn{1}{c}{$n$} & itr($\bC_W$) & itr($\bC$) & \multicolumn{1}{c}{$\varepsilon_{PCG}$}  & Diag. (s) & Itr (s) & Total (s) \\
  \hline
  16,000  &    166 & 1296 & $1.02 \times 10^{-4}$   &      80 &     113 &     193 \\
  32,000  &    247 & 3065 & $9.88 \times 10^{-5}$   &     215 &     321 &     536 \\
  64,000  &    372 & 5517 & $1.00 \times 10^{-4}$   &     665 &    1226 &    1891 \\
  128,000 &    547 &    - & $4.84 \times 10^{-5}$   &    2060 &    3237 &    5397 \\
  256,000 &    847 &    - & $5.00 \times 10^{-5}$   &    5775 &    9885 &   15660 \\
  512,000 &   1129 &    - & $3.74 \times 10^{-5}$   &   17896 &   33116 &   51012 \\
\end{tabular}\\
\corb{(b) $\btheta_b = (1,1/6)$, $d = 3$, $f = 3$ ($p = 20$), $\tilde{f} = 3$ ($\tilde{p} = 20$)} \\
\begin{tabular} { r r r r r r r }
  \multicolumn{1}{c}{$n$} & itr($\bC_W$) & itr($\bC$) & \multicolumn{1}{c}{$\varepsilon_{PCG}$}  & Diag. (s) & Itr (s) & Total (s) \\
  \hline
  16,000  &    293 &  2970 &$8.23 \times 10^{-5}$   &     79  &     198 &     277 \\
  32,000  &    470 &  7786 &$8.18 \times 10^{-5}$   &    213  &     607 &     820 \\
  64,000  &    760 & 15808 &$7.09 \times 10^{-5}$   &    662  &    2495 &    3157 \\
  128,000 &   1167 &     - &$3.00 \times 10^{-5}$   &   2050  &    7109 &    9159 \\
  256,000 &   1961 &     - &$3.27 \times 10^{-5}$   &   5789  &   22878 &   28667 \\  
% 512,000 &        &     - &$     \times 10^{-5}$   &         &         &         \\
\end{tabular}\\
% \begin {tabular} { r r r r r r }
%   \multicolumn{1}{c}{$n$} & itr & \multicolumn{1}{c}{$\varepsilon_{PCG}$}  & Diag. (s) & Itr (s) & Total (s) \\
%   \hline
%   16,000  &    500 &  $6.28 \times 10^{-5}$   &         &         &         \\
%   32,000  &    861 &  $7.52 \times 10^{-5}$   &         &         &         \\
%   64,000  &   1556 &  $4.65 \times 10^{-5}$   &         &         &         \\
%   128,000 &   2928 &  $2.11 \times 10^{-5}$   &         &         &         \\
%   256,000 &        &  $     \times 10^{-5}$   &         &         &         \\
% \end {tabular}\\
% (c) $\btheta = (5/4,1/6)$\\
\corb{(c) $\btheta_c = (5/4,1/6)$, $d = 3$,
$f = 3$ ($p = 20$), $\tilde{f} = 3$ ($\tilde{p} = 20$)} \\
\begin{tabular} { r r r r r r r }
  \multicolumn{1}{c}{$n$} & itr($\bC_W$) & itr($\bC$) & \multicolumn{1}{c}{$\varepsilon_{PCG}$}  & Diag. (s) & Itr (s) & Total (s) \\
  \hline
  16,000  &    500 &  5953 & $6.27 \times 10^{-5}$   &     138 &     580 &    718   \\
  32,000  &    827 & 17029 & $7.29 \times 10^{-5}$   &     346 &    1574 &   1920  \\
  64,000  &   1567 & 37018 & $4.45 \times 10^{-5}$   &     910 &    6474 &   7384  \\
  128,000 &   2381 &     - & $2.23 \times 10^{-5}$   &    3974 &   25052 &  29026  \\
  256,000 &   4299 &     - & $2.61 \times 10^{-5}$   &   10322 &   72374 &  82696   \\
% 512,000 &        &  $     \times 10^{-5}$   &         &         &         \\
\end{tabular}\\
\end{center}
\label{SRBF:table1}
\end{table}

In Table \ref{SRBF:table2} the results for computing $\hat{\bbeta}$,
$\hat{\bgamma}$ and $\hat Z(\bs_0)$ for 1000 target points $\bs_0$ for
data set \#1 (2D) with the Mat\'{e}rn covariance function are
tabulated. We have three test cases: (a) $\btheta_a =(\nu,\rho)=
(1/2,1/6)$ (note for this case we obtain an exponential covariance
function), (b) $\btheta_b = (3/4,1/6)$, and (c) $\btheta_c = (1,1/6)$.
For (a) and (b) the KIFMM accuracy is set to medium. For (c) the KIFMM
accuracy is set to high. For this case the relative residual accuracy
for the unpreconditioned system is fixed at $10^{-2}$.

\corb{We note that for this case the results are even more impressive
  than for the 3D case}. For Table \ref{SRBF:table2} (c) the CG solver
stagnated and we terminated the iteration after 100,000. At this point
the matrix $\bC(\btheta)$ is highly ill-conditioned. In contrast, with
$\bC_W(\btheta)$ we are still able to solve the problem even for
128,000 size problem.


\begin{table}[t!]
\caption{Diagonal pre-conditioned results for computing
  $\bar{\bC}_{W}(\btheta)\bar{\bgamma}_{W} = \bar{\bZ}_{W}$ for the
  Mat\'{e}rn covariance with $\btheta = (\nu,\rho)$ for data set \#1 (2D). We look at three
  cases: (a) $\btheta_a = (3/4,1/6)$, (b) $\btheta_b = (1,1/6)$ and
  (c) $\btheta_c = (5/4,1/6)$.  The relative error of the residual of
  the unpreconditioned system is set to $\varepsilon = 10^{-2}$. The
  KIFMM is set to medium accuracy for (a), (b) and set to high
  accuracy for (c).}
\begin{center}
%% \begin{tabular} { r r r r r r }
%%   \multicolumn{1}{c}{$n$} & itr & \multicolumn{1}{c}{$\varepsilon_{PCG}$}  & Diag. (s) & Itr (s) & Total (s) \\
%%   \hline
%%   16,000  &    386 &  $2.88 \times 10^{-4}$    &     282 &         &         \\
%%   32,000  &    392 &  $1.65 \times 10^{-4}$    &     866 &         &         \\
%%   64,000  &    581 &  $1.89 \times 10^{-4}$    &    2013 &         &         \\
%%   128,000 &    736 &  $1.00 \times 10^{-4}$    &    6247 &         &         \\
%%   256,000 &    906 &  $1.14 \times 10^{-4}$    &   16936 &         &         \\
%%   512,000 &   1690 &  $5.95 \times 10^{-5}$    &   45892 &         &         \\
%% \end{tabular}\\
%% (a) $\btheta^2_{a} = (1/2,1/6)$, $d = 2$, $f = 12$, $\tilde{f} = 12$ \\
\corb{
(a) $\btheta_{a} = (1/2,1/6)$, $d = 2$, $f = 3$  ($p = 10$), $\tilde{f} = 12$  ($\tilde{p} = 91$)
}\\
\begin{tabular} { r r r r r r r }
  \multicolumn{1}{c}{$n$} & itr($\bC_W$) & itr($\bC$) & \multicolumn{1}{c}{$\varepsilon_{PCG}$}  & Diag. (s) & Itr (s) & Total (s) \\
  \hline
  16,000  &    330 & 3603 & $2.39 \times 10^{-3}$    &     246 &     115 &   361   \\
  32,000  &    333 & 5429 & $1.39 \times 10^{-3}$    &     750 &     251 &  1001   \\
  64,000  &    455 & 8152 & $1.32 \times 10^{-3}$    &    1947 &     589 &  2536   \\
  128,000 &    564 &    - & $7.10 \times 10^{-4}$    &    5570 &    1577 &  7147   \\
  256,000 &    619 &    - & $9.78 \times 10^{-4}$    &   15266 &    3065 & 18331   \\
  512,000 &   1230 &    - & $4.50 \times 10^{-4}$    &   42254 &   13101 & 55355   \\
\end{tabular}\\
\corb{(b) $\btheta_b = (3/4,1/6)$, $d = 2$, $f = 14$ ($p = 120$), $\tilde{f} = 14$ ($\tilde{p} = 120$)} \\
\begin{tabular} { r r r r r r r }
  \multicolumn{1}{c}{$n$}  & itr($\bC_W$) & itr($\bC$) & \multicolumn{1}{c}{$\varepsilon_{PCG}$}  & Diag. (s) & Itr (s) & Total (s) \\
  \hline
  16,000  &   965  &  26795 & $2.78 \times 10^{-3}$    &     370 &     397 &    767  \\ 
  32,000  &  1110  &  41079 & $2.04 \times 10^{-3}$    &    1125 &    1061 &   2186  \\
  64,000  &  2239  &  82166 & $1.35 \times 10^{-3}$    &    2892 &    3714 &   6606  \\
  128,000 &  3443  &      - & $1.09 \times 10^{-3}$    &    8268 &   13130 &  21398  \\
  256,000 &  4557  &      - & $7.63 \times 10^{-4}$    &   23175 &   30302 &  53477  \\
% 512,000 &        &                          &         &         &         \\
\end{tabular}\\
\corb{
(c) $\btheta_c = (1,1/6)$, $d = 2$, $f = 14$ ($p = 120$), $\tilde{f}
= 14$ ($\tilde{p} = 120$) }\\
\begin{tabular} { r r r r r r r }
  \multicolumn{1}{c}{$n$} & itr($\bC_W$) & itr($\bC$) & \multicolumn{1}{c}{$\varepsilon_{PCG}$}  & Diag. (s) & Itr (s) & Total (s) \\
  \hline
  16,000  &  2710  &  $>100,000$ &  $1.90 \times 10^{-3}$  &     553 &    1844 &    2397  \\
  32,000  &  4261  &          -  &  $1.43 \times 10^{-3}$  &    1522 &    5713 &    7235  \\
  64,000  &  8801  &          -  &    $1.00 \times 10^{-4}$  &    5022 &   23785 &   28807  \\
  128,000 & 14405  &          -  &  $7.28 \times 10^{-4}$  &   12587 &   75937 &   88524 \\
\end{tabular}\\
\end{center}
\label{SRBF:table2}
\end{table}

The diagonal preconditioner we use is one of the simplest. We plan to
extend this approach to more sophisticated preconditioners such as
block Symmetric Successive OverRelaxation (SSOR) (see
\cite{Castrillon2013}) in the future.

The residual errors are then propagated to the final estimate
$Z(\bs_0)$ around the same magnitude. However, as a final experiment
in Table \ref{SRBF:table3} we tabulate the relative $l_2$ error
between the multi-level kriging approach and the direct method for
data set \#2 with exponential covariance function $\exp(-\theta r)$,
where $\theta = 5.9915$ and $f = 3$. The PCG tolerance is set to
$10^{-5}$ and $\tilde f = 3$. Notice that the error increases with
$n$. This is expected since the unpreconditioned system error will
degrade.
\begin{table}[h!]
\caption{Tabulation of the kriging estimate relative $l_2$ error
  between the multi-level kriging approach and the direct method for
  Data Set \#1 (3D) for 1000 target points. The covariance function is
  $\exp(-\theta r)$, where $\theta = 5.9915$. The polynomial bias term in
  the Gaussian model is set to $f = 3$ ($p = 20$). The solver
  tolerance is set to $\varepsilon_{PCG} = 10^{-5}$ and accuracy
  parameter is set to $\tilde p = 20$ ($\tilde f = 3$). }
 \begin{center}
 \begin{tabular}{ r c }
 \multicolumn{1}{c}{$n$} & 
 \multicolumn{1}{c}{$l_{2}$ Relative Error} \\
 \hline 
 1,000  &   $1.53 \times 10^{-6}$ \\
 2,000  &   $6.71 \times 10^{-5}$ \\
 4,000  &   $6.42 \times 10^{-5}$ \\
 % \hline
 8,000  &   $1.01 \times 10^{-4}$ \\
 16,000 &   $9.14 \times 10^{-5}$ \\
 32,000 &   $1.05 \times 10^{-4}$ \\
 \end{tabular}
 \end{center}
 \label{SRBF:table3}
\end{table}



%% \section{Statistical Numerical Examples}
%% \label{statistical}
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% % \section{Example}

%% In this section we study the statistical performance of our method to
%% estimate the covariance parameters and the Kriging estimate. Due to
%% computational limitations we are only able to generate test data for
%% up to 32,000 observation locations. In principle we should be able to
%% estimate the covariance parameters for 100,000 size problems and the
%% Kriging estimate for up to 512,000 irregularly placed observations
%% with good accuracy ($10^{-3}$). We now describe four different tests.



%% The optimization problem from the log-likelihood function
%% \eqref{Introduction:multilevelloglikelihood} is solved using fminbnd
%% from the optimization toolbox in MATLAB with the local minimizer
%% search in the interval $[1,10]$.  We set the distance criterion to
%% $\tau_{i,j} := 2^{-i}$.

%% {\bf Test \#1} Suppose we have a single realization $Z$ of the
%% Gaussian spatial random field $Z(\bs)$ with the Mat\'{e}rn covariance. We
%% set $f = 3$ ($p = 10$) and $\tilde{f} = 7$ ($\tilde{p} = 36$) and set
%% the covariance parameters to $\btheta = (3/4,1/6)$.  For a single
%% realization $Z$ we split the observations values into $n = 1000, 2000,
%% \dots, 32000$ size problems $Z_{1} \subset \dots \subset Z_{6}$ and
%% the respective locations $\bS_1 \subset \dots \subset \bS_6$. For each
%% observation values $Z_{1}, \dots, Z_{6}$ (and locations) apply the
%% transformation $\bW$ to compute $Z_{W,1}, \dots, Z_{W,6}$ and solve
%% the six optimization problems $\corblue{\hat{\btheta}_{i} : =
%%   \argmax_{\btheta \in (0,\infty) \times (0,\infty)}
%%   \ell_{W}(Z_{W,i},\btheta)}$ is solved for $i = 1,\dots,6$.

%% The optimization problem from the log-likelihood function
%% \eqref{Introduction:multilevelloglikelihood} is solved using
%% fminsearch from the optimization toolbox in MATLAB with the local
%% minimizer search for $\nu$ in the interval $[1/2,5/4]$ and $\rho$ in
%% the interval $[1/5,1/7]$.  We set the distance criterion to $C_{\tau}
%% := 1$ for $1,000$ to $32,000$ and the fminsearch tolerance to
%% $10^{-3}$.  In Table \ref{SRBF:table4}(b) we tabulate the results for
%% the parameters estimate with respect to problem size. We notice that
%% the estimates of $(\nu,\rho)$ converge to the actual values as we
%% increase the number of observations. In particular, for $n = 32,000$
%% the estimate $(0.75,0.16)$ is very close to the parameters
%% $(0.75,1/6)$ of the covariance matrix of the data.


%% {\bf Test \#3} (Exponential covariance).  The next step is to study
%% the accuracy of the Kriging estimate. We set the covariance parameters
%% to $\theta = 5.9915$.  We also pick 1000 random locations in the unit
%% cube for the Kriging estimate locations. In Table \ref{SRBF:table4}(c) we
%% tabulate the relative $l_2$ error of the Kriging estimate locations
%% for our method and the direct computation for $\bS_1, \dots, \bS_6$
%% with $\varepsilon_{PCG} = 10^{-5}$.

%% \begin{table}[ht]
%% \caption{(a) Tabulation of the mean estimate $\bbE[ \hat{\btheta} ]$
%%   and $std(\bbE[\hat{\btheta}])$ where the covariance matrix
%%   $\bC(\btheta)$ parameter is set to $\theta = 5.9915$ and $p = 20$.
%%   The second column is the parameter estimate for test \#1.  The third
%%   column is the number of realizations for test \#2. The fourth column
%%   corresponds to the mean of the covariance parameter estimate
%%   $\hat{\btheta}$ with respect to $M$ realizations of the noise
%%   model. The fifth column corresponds to the standard deviation of the
%%   mean estimate.  $\bbE[\hat{\btheta}]$. (b) Parameter estimate of
%%   $\nu$ and $\rho$ for test \#3 (Mat\'{e}rn covariance). (c)
%%   Tabulation of the Kriging estimate relative $l_2$ error between the
%%   multi-level Kriging approach and the direct method.}
%% \begin{center}
%% \begin{tabular}{ r r r}
%% \multicolumn{1}{c}{$n$} & 
%% \multicolumn{1}{c}{$\nu                     $} &
%% \multicolumn{1}{c}{$\rho                     $} 
%%  \\ 
%% \hline 
%% 1,000  & 0.6835 & 0.1762 \\ 
%% 2,000  & 0.7398 & 0.1725 \\ 
%% 4,000  & 0.6984 & 0.1822 \\ 
%% 8,000  & 0.7354 & 0.1724 \\ 
%% 16,000 & 0.7300 & 0.1737 \\
%% 32,000 & 0.7478 & 0.1690 \\ 
%% \end{tabular}
%% \hspace{10mm}
%% \begin{tabular}{ r c }
%%  \multicolumn{1}{c}{$n$} & 
%%  \multicolumn{1}{c}{$l_{2}$ Relative Error} \\
%%  \hline 
%%  1,000  &   $1.53 \times 10^{-6}$ \\
%%  2,000  &   $6.71 \times 10^{-5}$ \\
%%  4,000  &   $6.42 \times 10^{-5}$ \\
%%  % \hline
%%  8,000  &   $1.01 \times 10^{-4}$ \\
%%  16,000 &   $9.14 \times 10^{-5}$ \\
%%  32,000 &   $1.05 \times 10^{-4}$ \\
%%  \end{tabular}
%% \\
%% (b) \hspace{50mm} (c)
%% \end{center}
%% \label{SRBF:table4}
%% \end{table}


%% \begin{table}[ht]
%% \caption{$\tilde{p} = 36$, $p = 10$}
%% \begin{center}
%% \begin{tabular}{ r r r r r r r r}
%% \multicolumn{1}{c}{$n$} & 
%% \multicolumn{1}{c}{$t$} &
%% \multicolumn{1}{c}{$i$} &
%% \multicolumn{1}{c}{$\hat{\nu}$} &
%% \multicolumn{1}{c}{$\hat{\rho}$} &
%% \multicolumn{1}{c}{Sparsity} &
%% \multicolumn{1}{c}{size} &
%% \multicolumn{1}{c}{s/itr}
%%  \\ 
%% \hline 
%%  8,000 &  4 &  4 & 0.7530   & 0.1620 & &   & \\
%%  8,000 &  4 &  3 & 0.7143   & 0.1966 & &   & \\
%%  8,000 &  4 &  2 & 0.7326   & 0.1814 & &   & \\
%% \hline
%% 16,000 &  4 &  4  &  0.7745 & 0.1509  & &   &   \\
%% 16,000 &  4 &  3  &  0.7500 & 0.1673  & &   &   \\
%% 16,000 &  4 &  2  &  0.7582 & 0.1618  & &   &   \\
%% \hline
%% 32,000 &  5 &  5  &  0.7369 & 0.1790  & &   &    7  \\
%% 32,000 &  5 &  4  &  0.7458 & 0.1715  & &   &   73  \\
%% 32,000 &  5 &  3  &  0.7516 & 0.1668  & &   &  260  \\
%% 32,000 &  5 &  2  &  0.7499 & 0.1683  & &   &  900  \\
%% 32,000 &  5 &  1  &  0.7506 & 0.1679  & &   & 1572  \\
%% \end{tabular}
%% \end{center}
%% \label{SRBF:table5}
%% \end{table}




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conclusions}
\label{multilevelconclusions}

 In this paper we developed a multi-level restricted Gaussian maximum
 likelihood method for estimating the covariance function parameters
 and the computation of the best unbiased predictor. Our
 approach produces a new set of multi-level contrasts that decouples
 the covariance parameters from the deterministic components. In
 addition, the covariance matrix exhibits fast decay independently
 from the decay rate of the covariance function.  Due to the fast
 decay of the covariance matrix only a small set of coefficients of
 the covariance matrix are computed with a level-dependent criterion.
 We showed results of our method for the Mat\'{e}rn covariance with
 highly irregular placement of the observation locations to good
 accuracy.

%% The results for the construction of the multi-level covariance
%%   matrix and the BLUP estimate are obtained from a single core on the
%%   CPU at the moment.  However, almost all the computational cost is
%%   directly attributed to the KIFMM solver and the covariance function
%%   computation. Thus we should be able to significantly increase the
%%   performance of our code by the following avenues:
%% i) Although there exists a parallel library version of the
%%   covariance independent KIFMM \citep{Lashuk2012} we have not been able
%%   to compile it. We have been in contact with the authors and expect a
%%   new version of the code that properly works for multiple processors.
%% ii) Another important avenue is to develop a fast math
%%   library for the Mat\'{e}rn covariance function that takes advantage of an
%%   optimized instruction set of the x86 CPU.
%% iii) The covariance dependent KIFMM sacrifices speed for
%%   flexibility of the covariance. Highly optimized FMM have been developed
%%   for fixed covariances \citep{Yokota2013}.
%% iv) The current implementation of our code is in 3D. We
%%   expect to re-implement our code for 2D. Since the Taylor expansion
%%   in 2D is smaller than 3D for the same polynomial order we expect to
%%   obtain the same order of decay of the multi-level covariance matrix
%%   but with significant computational improvements.

We are currently working on deriving error
estimates of the kriging estimate and determinant computation with
respect to the number of degrees of freedom $n$. We are also
contemplating extending our multi-level approach to multivariate
random fields and cokriging (e.g. \cite{Furrer2011}).

Our method also applies to
non-stationary problems if the covariance function is differentiable
to degree $\tilde{f}$. For example, if the covariance function changes
smoothly with respect to the location, Lemma
\ref{multilevelcov:lemma1} still applies and the multi-level
covariance matrix decays at the same rate as a stationary one. Now,
even if the covariance function is non differentiable everywhere with
respect to the location, Lemma \ref{multilevelcov:lemma1} still
applies, but at a lower decay rate.

\corb{Note that we have not made direct comparisons with many of the
  approaches described in Section \ref{Introduction}. These methods
  are very good at solving a particular type of problem i.e. grid-like
  geometries and/or compact covariance functions. For these situations
  we recommend using the approaches already developed in the
  literature as they work well and are easier to implement. However,
  to our knowledge we have no seen fast results where the placement of
  the observations are highly irregular and the Mat\'{e}rn covariance
  function decays slowly. We plan to include some comparison for the
  spatial covariance paper we plan to write.}

% In this approach we construct the multi-level basis with respect to
% polynomials. This choice allows a lot of flexibility in the type of
% covariance matrices we can compress. However, for any reproducing
% kernel an alternative choice is the first $p$ eigenfunctions. This
% will lead to an optimal decay of the covariance matrix and will not be
% subject to the smoothness of the covariance function. The downside is
% that the covariance matrix will not necessarily decouple.
% We are also contemplating extending our multi-level approach
%   to multivariate random fields and cokriging
%   (e.g. \cite{Furrer2011}).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% \section*{Appendix A: Multi-Level REML Covariance Estimation}
%% \label{MultiLevelREML}

%% In this section, we show how to construct the multi-level
%% log-likelihood function
%% \eqref{Introduction:multilevelloglikelihood}. The underlying approach
%% is to construct a Multi-level (or Hierarchical) Basis (MB) that
%% decomposes the space $\R^{n} $ into a series of orthogonal
%% hierarchical spaces that are adapted to the locations of the
%% observations, $\bs_1, \ldots,\bs_n$.  The MB can then be used to: (i)
%% form the multi-level REML function
%% \eqref{Introduction:multilevelloglikelihood}; (ii) sparsify the
%% covariance matrix $\bC_{W}(\btheta)$; and (iii) improve the
%% conditioning over the covariance matrix $\bC(\btheta)$.

%% First, we establish some notation and definitions:
%% \begin{itemize}
%% \item Let $\alpha:= (\alpha_{1},\dots,\alpha_{d}) \in \mathbb{Z}{^d}$,
%%   $|\alpha| := \alpha_{1}+\dots+\alpha_{d}$, $\bx : = [x_1,\dots,x_d]$
%%   and $D^\alpha_{\bx} := \frac{\partial^{\alpha_{1} + \dots +
%%       \alpha_{d}}} {\partial x_{1}^{\alpha_{1}} \dots \partial
%%     x_{d}^{\alpha_{d}}}$. For any $f \in \bbN_0$ (where
%% $\mathbb{N}_0 := \mathbb{N} \cup \{0\}$) let ${\cal Q}^d_f$ be
%%   the set of monomials $\{x_1^{\alpha_1} \dots
%%   x_d^{\alpha_d}\,\,\,|\,\,\, |\alpha| \leq f\}$. Furthermore, let
%%   $\bM_f$ be the design matrix with respect to all the monomials in
%%   ${\cal Q}^d_f$. Note that the number of monomials of degree $f$ with
%%   dimension $d$ is $\begin{pmatrix} d + f \\ f \end{pmatrix}$.
 
%% \item Let $\bC(\btheta) : = \{ \phi(r_{i,j}; \btheta) \}$ where $\phi$
%%   is the covariance function, $r_{i,j} := \|\bs_i - \bs_j\|_2$ and
%%   $\bs_i,\bs_j \in \R^{d}$ for $i,j = 1, \dots, n$. Alternatively we
%%   refer to $\phi(\bx,\by; \btheta)$ as $\phi(r; \btheta)$, where
%%   $r:=\|\bx - \by\|_2$ and $\bx,\by \in \R^d$. Suppose ${\cal
%%     P}^{p}(\mathbb{S})$ is the span of the design matrix $\bM_f$.  We
%%   now assume that $\phi(r; \btheta)$ is a positive definite function
%%   and $C^{\tilde{f}+1}(\R)$ (where the degree $\tilde{f}$ is explained
%%   later) for all $r \in \R$ except at the origin.
%% \end{itemize}



%% \begin{rem}
%% In practice instead of using the set of monomials ${\cal Q}^d_f$ we use
%% the set of Chebyshev polynomials of the first kind as these lead to a
%% more stable numerical scheme.
%% \end{rem}


%% We now show how to construct a multilevel basis for $\R^n$. However,
%% we refer the reader to \cite{Castrillon2013} for a detailed
%% description. The first step is to decompose the locations into a
%% series of multi-level cubes of dimension $d$. Without loss of
%% generality we assume that all the locations are contained in a unit
%% cube $B^{0}_{0}$ at level 0 and index 0. If the number of locations
%% inside $B^{0}_{0}$ is more than $p$ then equally subdivide $B^{0}_{0}$
%% into $2^{d}$ cubes ($B^{1}_{0}$, $B^{0}_{1}$, ...), where $d$ is the
%% number of dimensions. If the number of locations is $p$ or less then
%% stop, associate every location with the cube $B^{0}_{0}$ and denote
%% this a leaf cube.  Otherwise, for each non-empty cube $B^{q}_{k}$ at
%% level $q = 1$ and index $k$ if the number of locations is more than
%% $p$ then subdivide, otherwise associate all the locations to
%% $B^{q}_{k}$ and denote this a leaf cube.  This process is repeated for
%% all the subdivided cubes at levels $q = 2, \dots,$ until no
%% subdivisions are possible. The result is a tree structure with $0,
%% \ldots, t$ levels (See Algorithm 1 in \cite{Castrillon2013} for more
%% details).  We denote the leaf cubes as all the non-empty cubes that
%% contain at most $p$ locations, i.e., they will correspond to the leafs
%% of the tree structure.

%% Suppose that there is a one-to-one mapping between the set vectors
%% ${\cal E}:=\{\bve_{1},\dots,\bve_{n}\}$, which we denote as leaf unit
%% vectors, and the set of locations $\{ \bs_{1},\dots,\bs_{n}\}$,
%% i.e. $\bs_{i} \longleftrightarrow \bve_{i}$ for all $i = 1, \dots,
%% n$. Furthermore, for any index $i,j \in \mathbb{N}_{0}$, $1 \leq i
%% \leq n$, $1 \leq j \leq n$, we have that $\bve_{i}[j] = \delta[i-j]$,
%% where $\delta[\cdot]$ is the discrete Kronecker delta function. It is
%% clear that the span of $\bve_{1},\dots,\bve_{n}$ corresponds to the
%% space $\R^{n}$.  Now, for each leaf cube $B^{q}_{k}$ in the tree ($q$
%% and $k$ are the corresponding level and index respectively) form the
%% set of unit vectors $C^{q}_{k}$ such that $\{ \bve_i \in
%% C^{q}_{k}\,\,|\,\, \bs_i \in B^{q}_{k}\}$.  In Figure \ref{MLRLE:fig1}
%% we show an example of the multi-level domain decomposition for $d=p =
%% 2$.

%% \begin{figure}

%% \vspace{-5cm}

%% \begin{tabular}{c c}

%% \begin{tikzpicture}[scale=.75] 
%%   \begin{scope} 
%%     [place/.style={circle,draw=blue!50,fill=blue!20,thick,
%%       inner sep=0pt,minimum size=1.5mm}]
%%     \draw[step=8,gray,very thin] (0, 0) grid (8, 8);

%%     % Anchors
%%     \node[anchor=center] at (0, -4.5) {$$};
%%     \node[anchor=center] at (0,   10) {$$};
%%     % \node[anchor=center] at (16, -4.5) {$$};
%%     % \node[anchor=center] at (16,   10) {$$};

%%     \node at (0.5,7.5) [place] {};
%%     \node at (0,8.0) [] {($\bs_1$,$\bve_1)$};

%%     \node at (0.75,7.2) [place] {};
%%     \node at (0.3,6.3) [place] {};
%%     \node at (0.6,6.7) [place] {};
%%     \node at (1.3,6.3) [place] {};
%%     \node at (1.5,6.8) [place] {};

%%     \node at (2.5,4.5) [place] {};
%%     % \node at (3.5,7.5) [place] {};
%%     \node at (0.3,4.5) [place] {};
%%     \node at (1.6,5.7) [place] {};
%%     \node at (2.6,6.3) [place] {};
%%     \node at (3.2,5.2) [place] {};

%%     \node at (0.5,3.5) [place] {};
%%     \node at (1.5,2.5) [place] {};
%%     \node at (2.3,2.2) [place] {};
%%     \node at (3.3,2.5) [place] {};


%%     \node at (2.3,0.3) [place] {};
%%     \node at (2.7,0.5) [place] {};
%%     \node at (2.2,1.4) [place] {};
%%     \node at (2.6,1.4) [place] {};
%%     \node at (3.3,0.5) [place] {};
%%     \node at (3.7,0.7) [place] {};
%%     \node at (3.2,1.3) [place] {};
%%     \node at (3.6,1.5) [place] {};


%%     \node at (4.3,2.3) [place] {}; \node at (5.7,3.5) [place] {};
%%     \node at (6.2,3.4) [place] {}; \node at (7.3,2.4) [place] {};
%%     \node at (7.5,0.5) [place] {}; \node at (7.75,1) [] {($\bs_n$,$\bve_n)$};
%%     \node at (8.5,8.5) [] {$B^{0}_0$};
%%   \end{scope}


%% \end{tikzpicture}
%% \vspace{-8cm}
%% & 



%% \begin{tikzpicture}[scale=1]

%%     \node[anchor=center] at (0, -4.5) {$$};
%%     \node[anchor=center] at (0,   10) {$$};

%% \begin{scope}[xshift=5cm, yshift=4cm,
%% place/.style={circle,draw=blue!50,fill=blue!20,thick,
%%       inner sep=0pt,minimum size=1.5mm},
%% placer/.style={circle,draw=red!50,fill=red!20,thick,
%%       inner sep=0pt,minimum size=1.5mm},
%% ]
%% \Tree [.\node[placer]{$B^{0}_{0}$}; 
%%              [.\node[place]{$B^{1}_{0}$};
%%                     [.\node[place]{$B^{2}_{0}$}; 
%%                             \node[place]{$B^{3}_{0}$}; 
%%                             \node[place]{$B^{3}_{1}$}; 
%%                             \node[place]{$B^{3}_{2}$}; 
%%                     ] 
%%                     [.\node[place]{$B^{2}_{1}$}; ] 
%%                     [.\node[place]{$B^{2}_{2}$}; ] 
%%                     [.\node[place]{$B^{2}_{3}$}; ] 
%%              ]                                         
%%              [.\node[place]{$B^{1}_{1}$}; 
%%                     [      [. $\vdots$ ] ]
%%              ]   
%%              [.\node[place]{$B^{1}_{2}$};
%%                     [      [.\node[place]{$B^{2}_{8}$}; ]
%%                            [.\node[place]{$B^{2}_{10}$}; ] 
%%                            [.\node[place]{$B^{2}_{11}$}; ]
%%                                           ]]] 
%% \end{scope}
%% \end{tikzpicture} 

%% %\Tree [.$B^{0}_{0}$ {$B^{1}_{0}$} {$B^{2}_{0}$} {$B^{3}_{0}$} ]
%% \\

%% \begin{tikzpicture}[scale=.75] 
%%   \begin{scope} 
%%     [place/.style={circle,draw=blue!50,fill=blue!20,thick,
%%       inner sep=0pt,minimum size=1.5mm}]
%%     \draw[step=4,gray,very thin] (0, 0) grid (8, 8);

%%     % Anchors
%%     \node[anchor=center] at (0, -4.5) {$$};
%%     \node[anchor=center] at (0,   10) {$$};

%%     \node at (0.5,7.5) [place] {};
%%     \node at (0,8.0) [] {($\bs_1$,$\bve_1)$};

%%     \node at (0.75,7.2) [place] {};
%%     \node at (0.3,6.3) [place] {};
%%     \node at (0.6,6.7) [place] {};
%%     \node at (1.3,6.3) [place] {};
%%     \node at (1.5,6.8) [place] {};

%%     \node at (2.5,4.5) [place] {};
%%     % \node at (3.5,7.5) [place] {};
%%     \node at (0.3,4.5) [place] {};
%%     \node at (1.6,5.7) [place] {};
%%     \node at (2.6,6.3) [place] {};
%%     \node at (3.2,5.2) [place] {};

%%     \node at (0.5,3.5) [place] {};
%%     \node at (1.5,2.5) [place] {};
%%     \node at (2.3,2.2) [place] {};
%%     \node at (3.3,2.5) [place] {};


%%     \node at (2.3,0.3) [place] {};
%%     \node at (2.7,0.5) [place] {};
%%     \node at (2.2,1.4) [place] {};
%%     \node at (2.6,1.4) [place] {};
%%     \node at (3.3,0.5) [place] {};
%%     \node at (3.7,0.7) [place] {};
%%     \node at (3.2,1.3) [place] {};
%%     \node at (3.6,1.5) [place] {};


%%     \node at (4.3,2.3) [place] {}; \node at (5.7,3.5) [place] {};
%%     \node at (6.2,3.4) [place] {}; \node at (7.3,2.4) [place] {};
%%     \node at (7.5,0.5) [place] {}; \node at (7.75,1) [] {($\bs_n$,$\bve_n)$};
%%     \node at (8.5,8.5) [] {$B^{0}_0$};
%%     \node at (3.5,3.5) [scale=0.75] {$B^1_1$};
%%     \node at (7.5,3.5) [scale=0.75] {$B^1_2$};
%%     \node at (3.5,7.5) [scale=0.75] {$B^1_0$};

%%  \end{scope}
%% \end{tikzpicture} 
%% \vspace{-8cm}

%% & 

%% \begin{tikzpicture}[scale=1]

%%     \node[anchor=center] at (0, -4.5) {$$};
%%     \node[anchor=center] at (0,   10) {$$};

%% \begin{scope}[xshift=5cm, yshift=4cm,
%% place/.style={circle,draw=blue!50,fill=blue!20,thick,
%%       inner sep=0pt,minimum size=1.5mm},
%% placer/.style={circle,draw=red!50,fill=red!20,thick,
%%       inner sep=0pt,minimum size=1.5mm}
%% ]
%% \Tree [.\node[placer]{$B^{0}_{0}$}; 
%%              [.\node[placer]{$B^{1}_{0}$};
%%                     [.\node[place]{$B^{2}_{0}$}; 
%%                             \node[place]{$B^{3}_{0}$}; 
%%                             \node[place]{$B^{3}_{1}$}; 
%%                             \node[place]{$B^{3}_{2}$}; 
%%                     ] 
%%                     [.\node[place]{$B^{2}_{1}$}; ] 
%%                     [.\node[place]{$B^{2}_{2}$}; ] 
%%                     [.\node[place]{$B^{2}_{3}$}; ] 
%%              ]                                         
%%              [.\node[placer]{$B^{1}_{1}$}; 
%%                     [      [. $\vdots$ ] ]
%%              ]   
%%              [.\node[placer]{$B^{1}_{2}$};
%%                     [      [.\node[place]{$B^{2}_{8}$}; ]
%%                            [.\node[place]{$B^{2}_{10}$}; ] 
%%                            [.\node[place]{$B^{2}_{11}$}; ]
%%                                           ]]] 
%% \end{scope}
%% \end{tikzpicture} 

%% \\


%% \begin{tikzpicture}[scale=.75] 
%%   \begin{scope} 
%%     [place/.style={circle,draw=blue!50,fill=blue!20,thick,
%%       inner sep=0pt,minimum size=1.5mm}]

%%     \draw[step=4,gray,very thin] (0, 0) grid (8, 8);
%%     \draw[step=2,gray,very thin] (0, 0) grid (4, 4);
%%     \draw[step=2,gray,very thin] (0, 4) grid (4, 8);
%%     \draw[step=2,gray,very thin] (4, 0) grid (8, 4);

%%     \draw[step=1,gray,very thin] (0, 6) grid (2, 8);
%%     \draw[step=1,gray,very thin] (2, 0) grid (4, 2);
    

%%     % Anchors
%%     \node[anchor=center] at (0, -4.5) {$$};
%%     \node[anchor=center] at (0,   10) {$$};

%%     \node at (0.5,7.5) [place] {};
%%     \node at (0,8.0) [] {($\bs_1$,$\bve_1)$};

%%     \node at (0.75,7.2) [place] {};
%%     \node at (0.3,6.3) [place] {};
%%     \node at (0.6,6.7) [place] {};
%%     \node at (1.3,6.3) [place] {};
%%     \node at (1.5,6.8) [place] {};

%%     \node at (2.5,4.5) [place] {};
%%     % \node at (3.5,7.5) [place] {};
%%     \node at (0.3,4.5) [place] {};
%%     \node at (1.6,5.7) [place] {};
%%     \node at (2.6,6.3) [place] {};
%%     \node at (3.2,5.2) [place] {};

%%     \node at (0.5,3.5) [place] {};
%%     \node at (1.5,2.5) [place] {};
%%     \node at (2.3,2.2) [place] {};
%%     \node at (3.3,2.5) [place] {};


%%     \node at (2.3,0.3) [place] {};
%%     \node at (2.7,0.5) [place] {};
%%     \node at (2.2,1.4) [place] {};
%%     \node at (2.6,1.4) [place] {};
%%     \node at (3.3,0.5) [place] {};
%%     \node at (3.7,0.7) [place] {};
%%     \node at (3.2,1.3) [place] {};
%%     \node at (3.6,1.5) [place] {};


%%     \node at (4.3,2.3) [place] {}; \node at (5.7,3.5) [place] {};
%%     \node at (6.2,3.4) [place] {}; \node at (7.3,2.4) [place] {};
%%     \node at (7.5,0.5) [place] {}; \node at (7.75,1) [] {($\bs_n$,$\bve_n)$};
%%     \node at (8.5,8.5) [] {$B^{0}_0$};
%%     \node at (3.5,7.5) [scale=0.75] {$B^2_3$};
%%     \node at (7.5,3.5) [scale=0.75] {$B^2_{11}$};
%%     \node at (1.5,3.5) [scale=0.75] {$B^2_4$};

%%  \end{scope}
%% \end{tikzpicture} 
%% &

%% \begin{tikzpicture}[scale=1]

%%     \node[anchor=center] at (0, -4.5) {$$};
%%     \node[anchor=center] at (0,   10) {$$};

%% \begin{scope}[xshift=5cm, yshift=4cm,
%% placer/.style={circle,draw=red!50,fill=red!20,thick,
%%       inner sep=0pt,minimum size=1.5mm}]
%% \Tree [.\node[placer]{$B^{0}_{0}$}; 
%%              [.\node[placer]{$B^{1}_{0}$};
%%                     [.\node[placer]{$B^{2}_{0}$}; 
%%                             \node[placer]{$B^{3}_{0}$}; 
%%                             \node[placer]{$B^{3}_{1}$}; 
%%                             \node[placer]{$B^{3}_{2}$}; 
%%                     ] 
%%                     [.\node[placer]{$B^{2}_{1}$}; ] 
%%                     [.\node[placer]{$B^{2}_{2}$}; ] 
%%                     [.\node[placer]{$B^{2}_{3}$}; ] 
%%              ]                                         
%%              [.\node[placer]{$B^{1}_{1}$}; 
%%                     [      [. $\vdots$ ] ]
%%              ]   
%%              [.\node[placer]{$B^{1}_{2}$};
%%                     [      [.\node[placer]{$B^{2}_{8}$}; ]
%%                            [.\node[placer]{$B^{2}_{10}$}; ] 
%%                            [.\node[placer]{$B^{2}_{11}$}; ]
%%                                           ]]] 
%% \end{scope}
%% \end{tikzpicture} 

%% \\

%% \end{tabular}

%% \vspace{-3cm}

%% \caption{Multi-level domain decomposition of the location of
%%   observations for $d=p = 2$. All the observation locations $\bs_i$
%%   are assumed to be contained in the unit box $B^{0}_{0}$ (colored red
%%   node in tree). If the number of observations are greater than $p =
%%   2$ we subdivide into equal boxes. At the end we obtain a multilevel
%%   decomposition of the observation points.}
%% \label{MLRLE:fig1}
%% \end{figure}

% The full algorithm
% is described in detailed in \cite{Castrillon2013}.

%% The next step is to construct an orthonormal multi-level basis.
%% Without loss of generality, for each finest level leaf cube
%% $B^{t}_{l}$ let ${\cal S}$ be the number of unit vectors in $C^t_l$. A
%% new basis can be constructed from the unit vectors in $C^t_l$ by
%% taking the following steps:
%% \begin{enumerate}[(a)]
%% \item Start with $q = t$.
%% \item For each non empty cube $C^q_l$ perform
%% \begin{enumerate}[i)]
%% \item Let $\bphi^{q,l} _{j} := \sum_{\bve_i \in C^q_l} c^{q,l} _{i,j}
%%   \bve_i, \hspace{2mm} j=1, \dots, a;
%% \hspace{5mm} \bpsi^{q,l}_{j} := \sum_{\bve_i \in C^q_l}  d^{q,l}_{i,j} \bve_i, \hspace{2mm}
%% j=a+1, \dots, s$,
%% where $c^{q,l}_{i,j}$, $d^{q,l}_{i,j} \in \mathbb{R}$ and
%% for some $a \in \mathbb{N}^{+}$. We desire that the new discrete MB
%% vector $\bpsi^{q,l}_{j}$ be orthogonal to ${\cal P}^{p}(\mathbb{S})$,
%% i.e., for all $g \in {\cal P}^{p}(\mathbb{S})$,
%% \begin{equation}
%% \sum_{i=1}^{n} g[i] \bpsi^{q,l}[i] = 0.
%% \label{hbconstruction:eqn1}
%% \end{equation}
%% \item Form the matrix ${\cal M}^{q,l}_{\cal S} := \bM\T \bV$, where
%%   ${\cal M}^{q,l}_{\cal S} \in \R^{p \times {\cal S}}$, $\bV \in \R^{n
%%   \times s}$, and $\bV: = [ \dots, \bve_i, \dots ]$ for all $\bve_i
%%   \in C^q_l$. Now, suppose that the matrix ${\cal M}^{q,l}_{\cal S}$
%%   has rank $a$ and perform the Singular Value Decomposition
%%   (SVD). Denote by $\bU_{\cal S}\bD_{\cal S}\bV_{\cal S}$ the SVD
%%   of ${\cal M}^{q,l}_{\cal S}$, where $\bU_{\cal S} \in \R^{p \times \cal
%%     S}$, $\bD_{\cal S} \in \R^{{\cal S} \times {\cal S}}$, and $\bV_{\cal S} \in
%%   \R^{ {\cal S} \times {\cal S}} $. 
%% \item To satisfy the equation \eqref{hbconstruction:eqn1} we make the
%%   choice
%% \[
%%   \left[ \begin{array}{ccc|ccc}
%%       c^{q,l}_{0,1} & \dots &c^{q,l}_{a,1} & d^{q,l}_{a+1,1} & \dots &d^{q,l}_{{\cal S},1} \\
%%       c^{q,l}_{0,2} & \dots &c^{q,l}_{a,2} & d^{q,l}_{a+1,2} & \dots &d^{q,l}_{{\cal S},2} \\
%%       \vdots & \vdots & \vdots & \vdots & \vdots & \vdots   \\
%%       c^{q,l}_{0,{\cal S}} & \dots &c^{q,l}_{a,{\cal S}} & d^{q,l}_{a+1,{\cal S}} & \dots &d^{q,l}_{{\cal S},{\cal S}}
%%     \end{array}
%% \right] := \bV\T_{\cal S},
%% % \label{eqDefVspT*}
%% \]
%% \noindent where the columns $a+1$, \dots, $s$ form an orthonormal
%% basis of the nulls pace ${N_0}({\cal M}_{\cal S})$. Similarly, the
%% columns $1,\dots, a$ form an orthonormal basis of $\R^{\cal S}
%% \backslash {N_0}({\cal M}_{\cal S})$. Since the vectors in $C^q_l$ are
%% orthonormal then $\bphi^{q,l}_{1}, \dots, \bphi^{q,l}_a,
%% \bpsi^{q,l}_{a+1}, \dots, \bpsi^{q,l}_{\cal S}$ form an orthonormal
%% basis of $\mbox{span}_{\bve_i \in C^{q}_{l}} \bve_i$.  Moreover
%% $\bpsi^{q,l}_{a+1}, \dots, \bpsi^{q,l}_{\cal S}$ satisfy equation
%% \eqref{hbconstruction:eqn1}, i.e., are orthogonal to ${\cal
%%   P}^{p}(\mathbb{S})$ and are locally adapted to the locations contained
%% in the cube $B^{q}_{l}$.
%% \item Denote by $D^{q}_{l}$ the collection of all the vectors
%%   $\bpsi^{q,l}_{a+1}, \dots, \bpsi^{q,l}_{\cal S}$ for every cube
%%   $B^{q}_{l}$. Notice that the vectors $\bphi^{q,l}_{1}, \dots,
%%   \bphi^{q,l}_a$, which we denote as scaling vectors, are {\it not}
%%   orthogonal to ${\cal P}(\mathbb{S})$ and we need to process them further.
%% \item Repeat steps ii) - v) for every cube $B^{q}_{l}$ at level $q$. 
%% \item Denote by $W_{q,l}(\mathbb{S})$ the span of all the vectors in
%%   $D^{q}_{l}$ and $W_{l}(\mathbb{S}) := \cup_{q,l} W_{q,l}(\mathbb{S})$, i.e.,  for
%%   all the cubes $B^{q}_{l}$ contained at level $q$. Denote by 
%%   $V_{q,l}(\mathbb{S})$ the span of all the vectors in $C^{q}_{l}$ and
%%   $V_{l}(\mathbb{S}) := \cup_{q,l} V_{q,l}(\mathbb{S})$, i.e.,  for all the cubes
%%   $B^{q}_{l}$ contained at level $q$.
%% \end{enumerate}
%% \item The next step is to go to level $q - 1$. For each non empty cube
%%   $B^{q-1}_{k}$ at level $q-1$ identify all the cubes $B^{q}_{l}$ at
%%   level $q$ such that $B^{q}_{l} \subset B^{q-1}_{k}$. Collect all the
%%   scaling vectors associated to all cubes such that $B^{q}_{l} \subset
%%   B^{q-1}_{k}$ and denote them by $C^{q-1}_{k}$. For notational
%%   simplicity ${\cal S}$ refers to all the scaling vectors in
%%   $C^{q-1}_k$.
%% \item Let $q: = q - 1$. Repeat steps (b) to (d), but replace the leaf
%%   unit vectors with the scaling vectors contained in $C^{q}_{l}$ for
%%   all non empty cubes $B^{q}_{l}$ contained in level $q$. Stop when
%%   level $q = -1$ is reached.
%% \end{enumerate}

%% When the algorithm stops we have a series orthogonal subspaces
%% $W_{0}(\mathbb{S}),\dots, W_{t}(\mathbb{S})$ (and their corresponding
%% basis vectors) that are orthogonal to $V_{0}(\mathbb{S})$. The
%% orthonormal basis vectors of $V_{0}(\mathbb{S})$ also span the space
%% ${\cal P}^{p}(\mathbb{S})$. It is not hard to see that $\R^{n} = {\cal
%%   P}^{p}(\mathbb{S}) \oplus W_{0}(\mathbb{S}) \oplus W_{1}(\mathbb{S})
%% \oplus \dots \oplus W_{t}(\mathbb{S}).$ (See Lemma 2.2 in
%% \cite{Castrillon2013}). Now, from Theorem 2.1 in
%% \cite{Castrillon2013}, it is shown that at most ${\cal O}(nt)$
%% computational steps are needed to compute the multi-level basis for $d
%% = 3$. Note that this the multi-level construction and the theorems in
%% \cite{Castrillon2013} can be easily extended to any finite dimension
%% $d$.

%% From the basis vectors of the subspaces $({\cal P}^{p}(\mathbb{S}))^{\perp} 
%% = \cup_{i=0}^{T} W_{i}(\mathbb{S})$ we can build an orthogonal projection
%% matrix $\bW:\R^{n} \rightarrow ({\cal P}^{p}(\mathbb{S}))^{\perp}$. 
%% The dimensions of $\bW$ is $(n - p) \times n$ since the total number
%% of orthonormal vectors that span ${\cal P}^{p}(\mathbb{S})$ is
%% $p$. Conversely, the total number of orthonormal vectors that span
%% $({\cal P}^{p}(\mathbb{S}))^{\perp}$ is $n-p$.

%% Let $\bL$ be a matrix where each row is an orthonormal basis vector of
%% ${\cal P}^{p}(\mathbb{S})$. For $i = 0,\dots,t$ let $\bW_i$ be a
%% matrix where each row is a basis vector of the space $W_i(\mathbb{S})$. We
%% can now form the matrix $\bW \in \mathbb{R}^{(n - p) \times n}$ where
%% $\bW := \left[ \bW_t\T, \dots, \bW_1\T,\bW_0\T \right] \T$, therefore
%% a) The matrix 
%% $ \bP := \left[
%% \begin{array}{c}
%% \bW \\
%% \bL
%% \end{array}
%% \right
%% ]$ 
%% is orthonormal, i.e., $\bP\bP\T= \bI_n$;
%% b) Any vector $\bv \in \R^{n}$ can be written as $\bv = \bL\T\bv_{L} +
%%   \bW\T\bv_{W}$ where $\bv_{L} \in \R^{p} $ and $\bv_{W} \in \R^{n-p}$ are
%%   unique; c) The matrix $\bW$ contains at most ${\cal O}(nt)$ non-zero
%%   entries and $\bL$ contains at most ${\cal O}(np)$ non-zero
%%   entries. This implies that for any vector $\bv \in \R^{n}$ the
%%   computational cost of applying $\bW \bv$ is at most ${\cal
%%     O}(nt)$ and $\bL \bv$ is at most ${\cal O}(np)$.
%%   % (i.e. the dimension of $\bbeta$).

%% It is now clear that we have achieved property i) from Section
%% \ref{Introduction}. As for property ii) the efficiency depends on the
%% placement of the locations. Indeed for certain pathological cases $t$
%% can increase faster than $\log{n}$. 


%% \begin{rem}
%% For practical cases, $t$ increases proportionally to
%% $\log{n}$. If the inter location spacing collapses as $n^{-q}$ , where
%% $q$ is independent of $n$ , then $q\log{n}$ levels are needed, see
%% Section 4 in \cite{Beatson1997} for details.
%% \end{rem}

%% \begin{rem} There are some cases (as shown in Section 
%% \ref{numericalstudy} ) where it is desirable that equation
%% \eqref{hbconstruction:eqn1} is satisfied for all $g \in {\cal
%%   P}^{\tilde{p}}(\mathbb{S})$ for some $\tilde{p} > p$ and ${\cal
%%   P}^{p}(\bbS) \subset {\cal P}^{\tilde{p}}(\bbS)$.  If
%% $\bM_{\tilde{f}}$ is full rank and $n > 2\tilde{p}$ we can construct a
%% decomposition of $\R^{n}$ as $ {\cal P}^{p}(\mathbb{S}) \oplus
%% W_{-1}(\mathbb{S}) \oplus W_{0}(\mathbb{S}) \oplus W_{1}(\mathbb{S})
%% \oplus\dots \oplus W_{t}(\mathbb{S}).$, where
%% $W_{-1}(\mathbb{S}):=({\cal P}^{\tilde{p}}(\mathbb{S})\backslash {\cal
%%   P}^{p}(\mathbb{S}))$. Now, $\bW$ is rewritten as $\bW:=
%%             [\bW\T_{t}, \dots, \bW\T_{-1}]\T$, where $\bW_{-1}:
%%             \R^n \rightarrow ({\cal
%%               P}^{\tilde{p}}(\mathbb{S})\backslash {\cal
%%               P}^{p}(\mathbb{S}))$ corresponds to a projection matrix
%%             (see \cite{Castrillon2013} for the details of this
%%             construction). We refer to $\tilde{p}$ as the
%%             accuracy parameter of the multi-level basis.
%% \label{MultiLevel:Construction:remark}
%% \end{rem}

%% With the multilevel basis (MB) we can transform the observation data
%% vector $\bZ$ by applying the matrix $\bW$. This leads to the
%% multi-level log-likelihood function
%% \eqref{Introduction:multilevelloglikelihood}. The covariance matrix
%% $\bC(\btheta)$ is now transformed into $\bC_{W}(\btheta)$ with the
%% structure shown in Figure \ref{multilevelcov:fig1} where each of the
%% blocks $\bC^{i,j}_W = \bW_i \bC(\btheta) \bW_j \T$ for all $i,j = 0,
%% \dots, t$. This implies that %the $l^{th}$ row and $k^{th}$ column of
%% the entries of the matrix $\bC^{i,j}_W$ are formed from all the
%% interactions of the MB vectors between level $i$ and $j$. Thus for any
%% $\bpsi^{i,k}_{\tilde{l}} \in D^{i}_k$ and any $\bpsi^{j,l}_{\tilde{k}}
%% \in D^{j}_l$ there is a unique entry of $\bC^{i,j}_W$ of the form $
%% (\bpsi^{i,k}_{\tilde{k}})\T \bC(\btheta) \bpsi^{j,l}_{\tilde{l}}$.
%% \begin{figure}[t!] 
%% \psfrag{a}[l]{$\bC^{t,t}_{W}$} 
%% \psfrag{b}[c]{$\dots$}
%%   \psfrag{c}[c]{\hspace{1mm} {$\bC^{t,i}_{W}$} } 
%%   \psfrag{d}[c]{ \hspace{5mm}  $\vdots$}
%%   \psfrag{e}[l]{$\ddots$} 
%%   \psfrag{f}[c]{$\vdots$}
%%   \psfrag{g}[c]{\hspace{5mm} $\bC^{-1,t}_{W}$} 
%%   \psfrag{h}[c]{$\dots$}
%%   \psfrag{i}[c]{\small $\bC^{i,i}_{W}$}
%% \begin{center}
%% \includegraphics[width=2.7in,height=2.10in,trim=0 0 110 0,clip]{./figures/NewKmatrix}
%% \end{center}
%% \caption{Organization of the multi-level covariance matrix
%%   $\bC_W(\btheta)$. Note for this figure $i = -1$.}
%% \label{multilevelcov:fig1}
%% \end{figure}
%% Note that the blocks $\bC_W^{i,j}$, where $i=-1$ or $j=-1$, correspond
%% to the case where the accuracy term $\tilde{p} > p$ and we include the
%% interactions from the basis functions of the space $W_{-1}(\bbS)$.

%% The following lemma relates the covariance function $\phi$, the degree
%% $\tilde{f}$ (corresponding to the accuracy parameter $\tilde{p}$) of
%% the design matrix $\bM_{\tilde{f}}$ to the decay of the entries of the
%% matrix $\bC_W(\btheta)$.
%% \begin{lem} Let $B_{\ba}$ be the smallest ball in ${\mathbb R}^{d}$
%% with radii $r_{\ba}$ centered around the midpoint $\ba \in {\mathbb
%%   R}^{d}$ of the cube $B^{i}_{l}$ such that $B^{i}_{l} \subset
%% B_{\ba}$. Similarly, let $B_{\bb}$ be the smallest ball in ${\mathbb
%%   R}^{d}$ with radii $r_{\bb} \in {\mathbb R}^{d}$ centered around the
%% midpoint $\bb$ of the cube $B^{j}_{k}$ such that $B^{j}_{k} \subset
%% B_{\bb}$.  Now, since $\bpsi^{i,l}_{\tilde{l}} \in W_i(\bbS)$ and
%% $\bpsi^{j,k}_{\tilde{k}} \in W_{j}(\bbS)$ satisfy the moment
%% orthogonality condition from equation \eqref{hbconstruction:eqn1} for
%% all $g \in {\cal P}^{\tilde{p}}(\bbS)$ then the following bound holds:
%% \begin{equation}
%% |(\bpsi^{i,k}_{\tilde{k}})\T \bC(\btheta) \bpsi^{j,l}_{\tilde{l}} |
%% \leq
%% \sum_{|\alpha|=\tilde{f}+1} \sum_{|\beta|=\tilde{f}+1}
%% \frac{r_a^{\alpha}}{\alpha!} \frac{r_b^{\beta}}{\beta!}
%% \sup_{\bx \in B_{\ba}, \by \in B_{\bb}}
%%  |D^{\alpha}_{\bx} D^{\beta}_{\by} \phi(\bx,\by;\btheta)|,
%% \label{multilevelcov:eqn1}
%% \end{equation}
%% \label{multilevelcov:lemma1}
%% for $i,j = 0,\dots,t$.
%% \end{lem}
%% From Lemma \ref{multilevelcov:lemma1} we observe that the decay of the
%% entries of $\bC_W(\btheta)$ is dependent on the magnitude of the
%% derivatives of the covariance function $\phi(r;\btheta)$, the size of
%% $B_{\ba}$ and $B_{\bb}$ and the degree of ${\cal Q}^d_{\tilde{f}}$. Thus if
%% $\phi(r;\btheta)$ is smooth on $B_{\ba}$ and $B_{\bb}$ the entries of
%% $\bC_W(\btheta)$ will be small.

\section*{Appendix A: Proofs}

{\bf Lemma \ref{multilevelcov:lemma1}:}

Since $\phi(r;\btheta)$ is in $C^{\tilde{f}+1}(\R)$, then by Taylor's
theorem we have that for every $\bx \in B_{\ba}$
$\phi(\bx,\by;\btheta) = \sum_{|\alpha| \leq {\tilde f}}
\frac{D^{\alpha}_{\bx}\phi(\ba,\by;\btheta)}{\alpha
  !}(\bx-\ba)^{\alpha} + R_{\alpha}(\bx,\by;\btheta)$,
\noindent where $(\bx-\ba)^{\alpha}:=(x_{1}-a_{1})^{\alpha_{1}} \cdots
(x_{d}-a_{d})^{\alpha_{d}}$, $\alpha ! := \alpha_{1}! \cdots
\alpha_{d}!$, and $R_{\alpha}(\bx,\by;\btheta) :=
\sum_{|\alpha|=\tilde{f}+1}\frac{ (\bx-\ba)^{\alpha}}{\alpha!}
D^{\alpha}_{\bx} \phi(\ba+s(\bx-\ba),\by;\btheta)$ for some $s \in
[0,1]$.  Now, recall that $\bpsi^{i,k}_{\tilde{k}}$ is orthogonal to
${\cal P}^{p}(\mathbb{S})$ then $\sum_{h = 1}^{n}
\bpsi^{i,k}_{\tilde{k}}[h] \phi(\bs_h,\by;\btheta) = \sum_{h = 1}^{n}
\bpsi^{i,k}_{\tilde{k}}[h]$ $R_{\alpha}(\bs_h,\by;\btheta)$.  Since
$\bpsi^{j,l}_{\tilde{l}}$ is also orthogonal to ${\cal
  P}^{p}(\mathbb{S})$ then by applying Taylor's theorem centered at
$\bb \in B_{\ba}$:
\[
\begin{split}
&|(\bpsi^{i,k}_{\tilde{k}})\T \bC(\btheta)
\bpsi^{j,l}_{\tilde{l}}|
= 
\left|\sum_{h = 1}^{n} 
\sum_{e = 1}^{n} 
\bpsi^{i,k}_{\tilde{k}}[h] \bpsi^{j,l}_{\tilde{l}}[e]
\phi(\bs_h,\bs_e;\btheta) \right| 
=
\left| \sum_{h = 1}^{n} \sum_{e = 1}^{n} \left(
\sum_{|\alpha|=\tilde{f}+1}\frac{(\bs_h-\ba)^{\alpha}}{\alpha!}  \right. \right. \\
& \left. \left. 
\sum_{|\beta|=\tilde{f}+1}\frac{(\bs_e-\bb)^{\beta}}{\beta!}
D^{\alpha}_{\bx} D^{\beta}_{\by}
\phi(
a+s(\bs_h-\ba), b+t(\bs_e-\bb);\btheta) \bpsi^{i,k}_{\tilde{k}}[h] \bpsi^{j,l}_{\tilde{l}}[e]
\right) \right| \\
\leq&
\sum_{|\alpha|=\tilde{f}+1} \sum_{|\beta|=\tilde{f}+1
}
\frac{r_{\ba}^{\alpha}}{\alpha!} \frac{r_{\bb}^{\beta}}{\beta!}   
% \sup_{\bx \in B_{\ba}, \by \in B_{\bb}}
|  D^{\alpha}_{\bx} D^{\beta}_{\by} \phi(\ba+s(\bs_h-\ba), \bb+t(\bs_e-\bb);\btheta)|
\left|\sum_{r = 1}^{n} \sum_{e = 1}^{n} 
\bpsi^{i,k}_{\tilde{k}}[r] \bpsi^{j,l}_{\tilde{l}}[e]\right|\\
% \leq& 
% \sum_{|\alpha|=f+1} \sum_{|\beta|=f+1}
% \frac{r_{\ba}^{\alpha}}{\alpha!} \frac{r_{\bb}^{\beta}}{\beta!}   
% \sup_{\bx \in B_{\ba}, \by \in B_{\bb}} |D^{\alpha}_{\bx} D^{\beta}_{\by} \phi(\bx, \by;\btheta)|   
% \sum_{r = 1}^{n} |\bpsi^{i,k}_{\tilde{k}}[r]| 
% \sum_{e = 1}^{n} |\bpsi^{j,l}_{\tilde{l}}[e] \\
\leq& 
\sum_{|\alpha|=\tilde{f}+1} \sum_{|\beta|=\tilde{f}+1}
\frac{r_{\ba}^{\alpha}}{\alpha!} \frac{r_{\bb}^{\beta}}{\beta!}   
\sup_{\bx \in B_{\ba}, \by \in B_{\bb}} |D^{\alpha}_{\bx} D^{\beta}_{\by} \phi(\bx, \by;\btheta)|,
% \sup_{x \in B_a, y \in B_b}
%  |D^{\alpha}_{x} D^{\beta}_{y} \phi(x, y;\btheta)|.
\end{split} 
\]
%% \[
%% \begin{split}
%% |(\bpsi^{i,k}_{\tilde{k}})\T \bC(\btheta)
%% \bpsi^{j,l}_{\tilde{l}}|
%% =& 
%% \left|\sum_{h = 1}^{n} 
%% \sum_{e = 1}^{n} 
%% \bpsi^{i,k}_{\tilde{k}}[h] \bpsi^{j,l}_{\tilde{l}}[e]
%% \phi(\bs_h,\bs_e;\btheta) \right| 
%% \end{split}
%% \]
%% \[
%% \begin{split}
%% =&
%% \left| \sum_{h = 1}^{n} \sum_{e = 1}^{n} \left(
%% \sum_{|\alpha|=f+1}\frac{(\bs_h-\ba)^{\alpha}}{\alpha!}  
%% \sum_{|\beta|=f+1}\frac{(\bs_e-\bb)^{\beta}}{\beta!}
%% \right. \right. \\
%% & \left. \left. \int^{1}_{0} \int^{1}_{0} 
%% D^{\alpha}_{x} D^{\beta}_{y}
%% \phi(
%% a+s(x_r-a), b+t(y_e-b);\btheta) \bpsi^{i,k}_{\tilde{k}}[r] \bpsi^{j,l}_{\tilde{l}}[e]
%% \,ds\,dt \right) \right| \\
%% \leq&
%% \sum_{|\alpha|=f+1} \sum_{|\beta|=f+1}
%% \frac{r_a^{\alpha}}{\alpha!} \frac{r_b^{\beta}}{\beta!}   
%% \sup_{x \in B_a, y \in B_b}
%% \int^{1}_{0} \int^{1}_{0}
%% |  D^{\alpha}_{x} D^{\beta}_{y} \phi(a+s(x-a), b+t(y-b);\btheta)|
%% \,ds\,dt \\
%% &
%% \left|\sum_{r = 1}^{n} \sum_{e = 1}^{n} 
%% \bpsi^{i,k}_{\tilde{k}}[r] \bpsi^{j,l}_{\tilde{l}}[e]\right|\\
%% \leq& 
%% \sum_{|\alpha|=f+1} \sum_{|\beta|=f+1}
%% \frac{r_a^{\alpha}}{\alpha!} \frac{r_b^{\beta}}{\beta!}
%% \sup_{x \in B_a, y \in B_b}
%%  |D^{\alpha}_{x} D^{\beta}_{y} \phi(x, y;\btheta)|   
%% \sum_{r = 1}^{n} |\bpsi^{i,k}_{\tilde{k}}[r]| 
%% \sum_{e = 1}^{n} |\bpsi^{j,l}_{\tilde{l}}[e] \\
%% \leq& 
%% \sum_{|\alpha|=f+1} \sum_{|\beta|=f+1}
%% \frac{r_a^{\alpha}}{\alpha!} \frac{r_b^{\beta}}{\beta!}.   
%% \sup_{x \in B_a, y \in B_b}
%%  |D^{\alpha}_{x} D^{\beta}_{y} \phi(x, y;\btheta)|.
%% \end{split} 
%% \]
for some $s,t \in [0,1]$. The last inequality follows since from
Schwartz' inequality $\sum_{h = 1}^{n} |\bpsi^{i,k}_{\tilde{k}}[h]|
 \leq
\sqrt{ \sum_{h = 1}^{n} (\bpsi^{i,k}_{\tilde{k}}[h])^2} = 1$ and
$\sum_{e = 1}^{n} |\bpsi^{j,l}_{\tilde{l}}[e]| \leq
\sqrt{ \sum_{e = 1}^{n} (\bpsi^{j,l}_{\tilde{l}}[e])^2} = 1$.
\hfill $\Box$

\noindent 
{\bf Lemma \ref{multilevelkriging:lemma2}:} Since the Mat\'{e}rn
covariance function $\phi(\bx,\by; \btheta)$ is positive definite
we have that for all $\bv \neq {\bf 0}$:
\[
\sum_{i,j=1}^{n} v_i v_{j} \bC^{i,j}(\btheta) = 
\sum_{i,j=1}^{n} v_i v_{j} \phi(\bx_{i},\by_{j}; \btheta) > 0,
\]
where $\bC^{i,j}(\btheta)$ is the $(i,j)$ entry of the matrix
$\bC(\btheta)$.  The diagonal terms of $\bC_{W}$ are of the form
$(\bpsi^{i,k}_{\tilde{k}})\T \bC(\btheta) \bpsi^{i,k}_{\tilde{k}} $.
This implies that

\[
(\bpsi^{i,k}_{\tilde{k}})\T \bC(\btheta)
\bpsi^{i,k}_{\tilde{k}}  > 0.
\]
Thus, $\bD_W$ will always be positive definite.
\hfill $\Box$

\section*{Appendix B: Notation}

{\footnotesize
\corb{
\begin{tabular}{l l  l l}
\multicolumn{4}{l}{\bf Index for when the following are first defined, 
mentioned or reformulated.} \\
$d \in \mathbb{N}$               & Dimension of problem (p1)&
$n \in \mathbb{N}$               & Number of observations (p1) \\
$t \in \mathbb{N}$               & Maximum MB level (p8) &
$p \in \mathbb{N}$               & Number of columns of $\bM$ (p1) \\
$f \in \mathbb{N}$               & Polynomial degree (p7) &
$\tilde{p} \in \mathbb{N}$       & Accuracy parameter of MB (p7) \\
$\tilde{f} \in \mathbb{N}$       & Degree of multilevel basis. p(7) &
$w \in \mathbb{N}$               & Dimension of $\btheta$  (p1) \\
$\bM \in \R^{n \times p}$           & Design matrix (p1) &
$\phi(\cdot)$                    & Covariance function (p7) \\
$\mathbb{S} := \{ \bs_1,\dots,\bs_n\}$    & Locations of observations (p1) &
$\btheta:=(\nu,\rho) \in \R^{2}$  & Param. of matern kernel (p1) \\
$\bC(\btheta) \in \R^{n \times n}$  & Covariance matrix (p1) &
$\bZ \in \R^{n}$                  & Observation values (p1) \\
$\bs_{0} \in \R^{d}$               & Target points set (p1) &   
${\cal P}^{p}(\mathbb{S})$        & Span of the columns of  $\bM$ (p1) \\    
$\bbeta \in \R^{p}$               & Vector of unknowns from & 
$\hat{\bbeta} \in \R^p $,        & Estimates \\
                                 & deterministic model (p1) &   
$\hat{\bZ}(\bs_0)$               & kriging estimate at $\bs_0$ (p2) \\
$l(\bbeta,\btheta) \in \R$       & Log-likelihood function (p1) &
$B^{q}_k$                         & Cube at level $q$ and index $k$ (p8)\\
$\bL$                           &  $\R^n \rightarrow {\cal P}^{p}(\mathbb{S})$ (p4) & 
$\bW$                           &  $\R^n \rightarrow ({\cal P}^{p}(\mathbb{S}))^{\perp}$ (p4) \\
$\bP \T$                           &  
$ \left[
\begin{array}{cc}
\bW \T &
\bL \T
\end{array}
\right
]$ (p10) & $\tilde{\bZ}^i_W$ &  (p17) \\
                      & 
    &
$\tilde{\ell}^i_{W}(\btheta)$        & Multilevel log-likelihood (p17)\\
${\cal Q}^d_f$ & Set of polynomial monomials & 
$\tau \in \mathbb{N}_0$  & Level dependent \\
               & of order $f$ and dimension $d$ (p7)& & criterion constant (p13)\\
\end{tabular}
}
}
\section*{Appendix C}

Using the approach developed in this paper we can
compute the MSE at the target point $\bs_0$. Now, since $\bP\T \bP =
\bI$ we have that $\bM_f \T \bC(\btheta)^{-1} \bM_f = \bM_f \T \bP\T
\bP \bC(\btheta)^{-1} \bP \T \bP \bM_f = \bM_f \T \bP\T
\bK_{W}(\btheta)^{-1}\bP \bM_f$, $\bK_W(\btheta) := \bP \bK(\btheta)
\bP \T$, $\bK_{W}(\btheta):= \left[
\begin{array}{c c}
 \bC_W(\btheta) & \ba_W \\
\ba_W \T & \bS_W(\btheta)
\end{array}
\right ]$ 
where $\bS_W(\btheta) := \bL \bC(\btheta) \bL \T \in \R^{p
  \times p}$ and $\ba_W(\btheta) := \bW \bC(\btheta) \bL \T \in \R^{
  (n - p) \times p}$. Let
$\tilde{\bS}_W(\btheta):=(\bS_W(\btheta) - \ba_W \T \bC_W(\btheta)^{-1} \ba_W)^{-1}$,
then
\[
\bK_{W}(\btheta)^{-1}:= \left[
\begin{array}{c c}
\bC_W(\btheta)^{-1}  +
\bC_W(\btheta)^{-1} \ba_W \tilde{\bS}_W(\btheta) \ba_W \T \bC_W(\btheta)^{-1}
& -\bC_W(\btheta)^{-1} \ba_W \tilde{\bS}_W(\btheta) \\
-\tilde{\bS}_W(\btheta) \ba_W\T \bC_W(\btheta)^{-1}
& \tilde{\bS}_W(\btheta) 
\end{array}
\right].
\]
Given that $\bW \bM_f  = {\bf 0}$ then $ \bM_f \T \bC(\btheta)^{-1} \bM_f
= (\bL \bM_f)^T \tilde{\bS}_W(\btheta) (\bL \bM_f) $.  Following a
similar argument we have that equation \eqref{Kriging:MSE} becomes
\begin{equation}
1 + 
\tilde{\bu}
\T
\bM_f \T 
(\bL \bM_f) \T \tilde{\bS}_W(\btheta) (\bL \bM_f)
\tilde{\bu}
-\bc_{W}\T\bK_{W}^{-1}(\btheta)\bc_{W},
\label{Kriging:MSE2}
\end{equation}
where $\bc_W := \bP \bc$.  By using matrix-vector products with the
PCG method
% method we have that $\tilde{\bu}$ can be compute in ${\cal O}(kp
% + p^3)$ and $\bK_{W}^{-1}(\btheta)\bc_{W}$.
 $\tilde{\bS}_W(\btheta)$ can be computed in ${\cal O}((k(t+2)+1)np +
p^d)$. Thus the term $(\bL \bM_f) \T \tilde{\bS}_W(\btheta) (\bL
\bM_f) $ can be computed in ${\cal O}((k(t+2)+1)np + p^d + (t+2)pn)$.  Now, by
also using matrix-vector products the term
$\bK_{W}^{-1}(\btheta)\bc_{W}$ can be computed in ${\cal O}(kn(t+2) + p^d
)$. Thus the total cost for computing equation \eqref{Kriging:MSE2} is
${\cal O}((k(t+2)+1)np + p^d + (t+2)pn)$.\\

\noindent {\large {\bf Acknowledgements:}} We appreciate the help and
advice from Jun Li and Lisandro Dalcin in getting the C++ code working
properly, and to Stefano Castruccio for giving us feedback on our
manuscript.


\baselineskip 22pt
\pagestyle{empty}


\input{Multi-Level-REML-Avignon-Arxiv.bbl}

\clearpage
\pagenumbering{gobble} 
\end{document}





Multi-Level-REML-Avignon-Arxiv.bbl
\begin{thebibliography}{23}
\newcommand{\enquote}[1]{``#1''}
\expandafter\ifx\csname natexlab\endcsname\relax\def\natexlab#1{#1}\fi

\bibitem[{Anitescu et~al.(2012)Anitescu, Chen, and Wang}]{Anitescu2012}
Anitescu, M., Chen, J., and Wang, L. (2012), \enquote{A Matrix-Free Approach
  for Solving the Parametric Gaussian Process Maximum Likelihood Problem,}
  \textit{SIAM Journal on Scientific Computing}, 34, 240--262.

\bibitem[{Balay et~al.(2013{\natexlab{a}})Balay, Brown, Buschelman, Eijkhout,
  Gropp, Kaushik, Knepley, McInnes, Smith, and Zhang}]{petsc-user-ref}
Balay, S., Brown, J., Buschelman, K., Eijkhout, V., Gropp, W.~D., Kaushik, D.,
  Knepley, M.~G., McInnes, L.~C., Smith, B.~F., and Zhang, H.
  (2013{\natexlab{a}}), \enquote{{PETS}c Users Manual,} Tech. Rep. ANL-95/11 -
  Revision 3.4, Argonne National Laboratory.

\bibitem[{Balay et~al.(2013{\natexlab{b}})Balay, Brown, Buschelman, Gropp,
  Kaushik, Knepley, McInnes, Smith, and Zhang}]{petsc-web-page}
Balay, S., Brown, J., Buschelman, K., Gropp, W.~D., Kaushik, D., Knepley,
  M.~G., McInnes, L.~C., Smith, B.~F., and Zhang, H. (2013{\natexlab{b}}),
  \enquote{{PETSc} {W}eb {P}age,} Http://www.mcs.anl.gov/petsc.

\bibitem[{Balay et~al.(1997)Balay, Gropp, McInnes, and Smith}]{petsc-efficient}
Balay, S., Gropp, W.~D., McInnes, L.~C., and Smith, B.~F. (1997),
  \enquote{Efficient Management of Parallelism in Object Oriented Numerical
  Software Libraries,} in \textit{Modern Software Tools in Scientific
  Computing}, eds. Arge, E., Bruaset, A.~M., and Langtangen, H.~P.,
  Birkh{\"{a}}user Press, pp. 163--202.

\bibitem[{Beatson and Greengard(1997)}]{Beatson1997}
Beatson, R. and Greengard, L. (1997), \enquote{A Short Course on Fast Multipole
  Methods,} in \textit{Wavelets, Multilevel Methods and Elliptic PDEs}, Oxford
  University Press, pp. 1--37.

\bibitem[{Castrill\'{o}n-Cand\'{a}s et~al.(2013)Castrill\'{o}n-Cand\'{a}s, Li,
  and Eijkhout}]{Castrillon2013}
Castrill\'{o}n-Cand\'{a}s, J., Li, J., and Eijkhout, V. (2013), \enquote{A
  Discrete Adapted Hierarchical Basis Solver for Radial Basis Function
  Interpolation,} \textit{BIT Numerical Mathematics}, 53, 57--86.

\bibitem[{Chen et~al.(2008)Chen, Davis, Hager, and Rajamanickam}]{Chen2008}
Chen, Y., Davis, T.~A., Hager, W.~W., and Rajamanickam, S. (2008),
  \enquote{Algorithm 887: CHOLMOD, Supernodal Sparse Cholesky Factorization and
  Update/Downdate,} \textit{ACM Trans. Math. Softw.}, 35, 22:1--22:14.

\bibitem[{Davis and Hager(1999)}]{Davis1999}
Davis, T. and Hager, W. (1999), \enquote{Modifying a Sparse Cholesky
  Factorization,} \textit{SIAM Journal on Matrix Analysis and Applications},
  20, 606--627.

\bibitem[{Davis and Hager(2001)}]{Davis2001}
--- (2001), \enquote{Multiple-Rank Modifications of a Sparse Cholesky
  Factorization,} \textit{SIAM Journal on Matrix Analysis and Applications},
  22, 997--1013.

\bibitem[{Davis and Hager(2005)}]{Davis2005}
--- (2005), \enquote{Row Modifications of a Sparse Cholesky Factorization,}
  \textit{SIAM Journal on Matrix Analysis and Applications}, 26, 621--639.

\bibitem[{Davis and Hager(2009)}]{Davis2009}
Davis, T.~A. and Hager, W.~W. (2009), \enquote{Dynamic Supernodes in Sparse
  Cholesky Update/Downdate and Triangular Solves,} \textit{ACM Trans. Math.
  Softw.}, 35, 27:1--27:23.

\bibitem[{Elden et~al.(2004)Elden, Wittmeyer-Koch, and Nielsen}]{Elden2004}
Elden, L., Wittmeyer-Koch, L., and Nielsen, H. (2004), \textit{Introduction to
  Numerical Computation - analysis and Matlab illustrations},
  Studentlitteratur.

\bibitem[{Furrer and Genton(2011)}]{Furrer2011}
Furrer, R. and Genton, M.~G. (2011), \enquote{Aggregation-Cokriging for
  Highly-Multivariate Spatial Data,} \textit{Biometrika}, 98, 615--631.

\bibitem[{Furrer et~al.(2006)Furrer, Genton, and Nychka}]{Furrer2006}
Furrer, R., Genton, M.~G., and Nychka, D. (2006), \enquote{Covariance Tapering
  for Interpolation of Large Spatial Datasets,} \textit{Journal of
  Computational and Graphical Statistics}, 15, 502--523.

\bibitem[{George(1973)}]{George1973}
George, A. (1973), \enquote{Nested Dissection of a Regular Finite Element
  Mesh,} \textit{SIAM Journal on Numerical Analysis}, 10, 345--363.

\bibitem[{Gilbert and Tarjan(1987)}]{Gilbert1987}
Gilbert, J.~R. and Tarjan, R.~E. (1987), \enquote{The Analysis of a Nested
  Dissection Algorithm,} \textit{Numerische Mathematik}, 50, 377--404.

\bibitem[{Kaufman et~al.(2008)Kaufman, Schervish, and Nychka}]{Kaufman2008}
Kaufman, C.~G., Schervish, M.~J., and Nychka, D.~W. (2008), \enquote{Covariance
  Tapering for Likelihood-Based Estimation in Large Spatial Datasets,}
  \textit{Journal of the American Statistical Association}, 103, 1545--1555.

\bibitem[{Stein et~al.(2012)Stein, Chen, and Anitescu}]{Stein2012}
Stein, M.~L., Chen, J., and Anitescu, M. (2012), \enquote{Difference Filter
  Preconditioning for Large Covariance Matrices,} \textit{SIAM Journal on
  Matrix Analysis and Applications}, 33, 52--72.

\bibitem[{Stein et~al.(2013)Stein, Chen, and Anitescu}]{Stein2013}
--- (2013), \enquote{Stochastic Approximation of Score Functions for Gaussian
  Processes,} \textit{Annals of Applied Statistics}, 7, 1162--1191.

\bibitem[{Stein et~al.(2004)Stein, Chi, and Welty}]{Stein2004}
Stein, M.~L., Chi, Z., and Welty, L.~J. (2004), \enquote{Approximating
  Likelihoods for Large Spatial Data Sets,} \textit{Journal of the Royal
  Statistical Society, Series B}, 66, 275--296.

\bibitem[{Sun et~al.(2012)Sun, Li, and Genton}]{Sun2012}
Sun, Y., Li, B., and Genton, M.~G. (2012), \enquote{Geostatistics for Large
  Datasets,} in \textit{Space-Time Processes and Challenges Related to
  Environmental Problems}, eds. Porcu, M., Montero, J.~M., and Schlather, M.,
  Springer, pp. 55--77.

\bibitem[{Sun and Stein(2015)}]{Sun2015}
Sun, Y. and Stein, M.~L. (2015), \enquote{Statistically and Computationally
  Efficient Estimating Equations for Large Spatial Datasets,} \textit{Journal
  of Computational and Graphical Statistics}, Posted online 07 Nov 2014.

\bibitem[{Ying et~al.(2004)Ying, Biros, and Zorin}]{ying2004}
Ying, L., Biros, G., and Zorin, D. (2004), \enquote{A Kernel-Independent
  Adaptive Fast Multipole Method in Two and Three Dimensions,} \textit{Journal
  of Computational Physics}, 196, 591--626.

\end{thebibliography}
