text:
\begin{document}

\markboth{A. Brandt, J. Brannick, K. Kahl, and I. Livshits}{Bootstrap Algebraic Multigrid}
\title{Bootstrap Algebraic Multigrid: status report, open problems, and outlook}
\author[Achi Brandt, James Brannick, Karsten Kahl, and Ira Livshits]{Achi Brandt\affil{1}, James Brannick\affil{2}\comma\corrauth, Karsten Kahl\affil{3}, and Ira Livshits\affil{4}}
\address{\affilnum{1}\ Microsoft Corporation, Mountainview, CA, 94043, USA.
 \\
 \affilnum{2}\ Department of Mathematics, The
       Pennsylvania State University, University Park, PA 16802, USA.
\\
 \affilnum{3}\ Department of Mathematics, University of Wuppertal, Wuppertal,
 Germany.
 \\
  \affilnum{4}\ Department of Mathematical Sciences, Ball State University, Muncie IN, 47306, USA.
       }

\emails{
 {\tt acbrandt@microsoft.com} (Achi Brandt),
  {\tt brannick@psu.edu} (James Brannick),
 {\tt kkahl@math.uni-wuppertal.de} (Karsten Kahl),
 {\tt ilivshits@bsu.edu} (lra Livshits)
 }

\begin{abstract}
This paper provides an overview of the main ideas driving the bootstrap algebraic multigrid methodology, including compatible relaxation and algebraic distances for defining effective coarsening strategies, the least squares method for computing accurate prolongation operators and the bootstrap cycles for computing the test vectors that are used in the least squares process.
We review some recent research in the development, analysis and application of bootstrap algebraic multigrid and point to open problems in these areas.  Results from our previous research as well as some new results for some model diffusion problems with highly oscillatory diffusion coefficient are presented
to illustrate the basic components of the BAMG algorithm.
\end{abstract}

\keywords{Bootstrap algebraic multigrid, compatible relaxation, algebraic distance, least squares interpolation}

 \ams{65F10, 65N22, 65N55}

\maketitle
\section{Introduction}
{Numerous scientific and engineering problems find their formulation in terms of (systems)
of partial differential equations, which in turn require the solution of
large-scale finite element, finite difference, or finite volume equations.
Modern applications involve large-scale parallel processing of
(linear) systems with millions or even billions of unknowns, for which multigrid (MG)
methods often provide solvers that are optimal with respect to their computational complexity
and, hence, their parallel scalability.

Multigrid solvers for sparse systems of linear equations
$Ax=b$ are all based on two complimentary processes:
a local relaxation scheme (smoother) that eliminates certain
components of the error by working on the fine level and a coarse-level
correction that treats the remaining global error.}

Generally, the design and analysis of these two MG components are
based on the following smoothing property of relaxation.
For any $0 < \rho < 1$, an error vector $e$ is called $\rho$-smooth if
all its {\it normalized} residuals\footnote{The vector of {\it
normalized residuals} is $\ta e$, where $\ta$, the {\it normalized
matrix}, is scaled so that the $l_1$ norm of each of its rows is
$1$. Correspondingly, we define the {\it normalized eigenvalues}
of $A$ to be the eigenvalues of $\ta$. The magnitude of any
normalized eigenvalue is at most $1$.} are smaller than $\rho
\|e\|$. The basic observation \cite{B83} is that the convergence
of a proper relaxation process\footnote{Kaczmarz relaxation can
always be used (but better schemes are often available),
supplemented when needed by local relaxation steps around
exceptionally large residuals (as suggested in \cite[\S
A.9]{BAMG1}).} slows down only when the current error is
$\rho$-smooth with $\rho \ll 1$, the smaller the $\rho$ the slower
the convergence. This basic observation implies that when
relaxation slows down, the error vector $e$ can be approximated in
a much lower-dimensional subspace. Very efficient {\it ``geometric
multigrid''} solvers have been developed for the case that the
lower-dimensional subspace corresponds to functions on a
well-structured grid (the {\it coarse} level), on which the smooth
errors can be approximated by easy-to-derive equations, based for
example on discretizing the same continuum operator that has given
rise to the fine-level equations $Ax=b$ {to define the coarse-level operator
$A_c$.} The coarse-level
equations are solved using recursively a similar combination of
relaxation sweeps and still-coarser-level approximations to the
resulting smooth errors.

The basic two-grid method for solving $Ax=b$, from which a multigrid method is defined by recursion,
involves a stationary linear iterative
method (the smoother) applied to the fine-grid system,
and a coarse-grid correction: given an approximation $w\in \mathbb
C^n$, compute an update $v\in \mathbb C^n$ by \medskip
\begin{enumerate}
\item Pre-smoothing: $y = w+M(f-Aw)$,
\item Correction: $v = y + PA_c^{-1}R(f-Ay)$.

\end{enumerate}\medskip
Here, $M$ is the approximate inverse of $A$ that defines the multigrid smoother and
$R: \mathbb{C}^{n} \mapsto \mathbb{C}^{n_c}$ and $P: \mathbb{C}^{n_c} \mapsto \mathbb{C}^n$ with $n_c < n$ are the restriction and
interpolation operators that map information between the coarse grid of size $n_c$ and the
fine grid of size $n$.

To deal with much more general situations, where the
fine-level system may not be defined on a well-structured grid
nor perhaps arise from any continuum problem, {\it ``algebraic
multigrid''} (AMG) methods were developed to derive the
set of coarse-level {\it variables} and coarse-level {\it
equations} directly from the given matrix $A$. The basic
approach (developed in \cite{BMR83, BMR84, RS87} and called today
``classical AMG'' or RS-AMG) involves the following two steps. \bi
\item[(1)] Choosing the coarse-grid variables, e.g., as a subset $C$ of
the set of fine-grid variables, $\Omega$, such that each variable in $F=\Omega
\setminus C$ is
``strongly connected'' to variables in $C$. \item[(2)]
Approximating the fine-level residual problem $Ae=r$ by the
coarse-level equations $A_c e_c=r_c$ using the Galerkin
definition $A_c = RAP$ with $r_c=Rr$, so that $e \approx P e_c$.\footnote{{
A common choice when $A$ is HPD is to set $R=P^{H}$, which leads to a variational
Galerkin scheme.}}
\ei
In this classical approach, the ``strength of connection'' measure as well as the
interpolation matrix $P$, and the {restriction matrix} $R$
are all defined directly in terms of the elements of the matrix
$A$.
The overall algorithm gives satisfactory results only when the (properly scaled) matrix $A$ has a dominant diagonal and (with small possible exceptions) all its off-diagonal elements have the same sign. Even then, the produced interpolation in many cases has limited accuracy, insufficient for full multigrid efficiency.

{The Bootstrap AMG (BAMG) algorithm, introduced in \cite[\S
17.2]{B00}, and the closely related Adaptive AMG algorithms introduced in~\cite{MBrezina_RFalgout_SMacLachlan_TManteuffel_SMcCormick_JRuge_2003,MBrezina_2005}, were developed
to extend AMG to an even wider class of problems for which the assumptions on smooth error
made in classical AMG are violated.

The main idea in Bootstrap and Adaptive AMG is to
apply the current solver to appropriately formulated
test problems to explicitly compute approximations of the
smooth error components.
These components are then incorporated
into the coarse space by the AMG setup algorithm
as needed.

The adaptive algebraic multigrid methods developed originally in \cite{MBrezina_RFalgout_SMacLachlan_TManteuffel_SMcCormick_JRuge_2003,MBrezina_2005},
use a sequential setup procedure to adaptively identify the components of the space that need to be represented in the range of the coarse grid operator.  While such an approach insures that test vectors are incorporated into the coarse space only as they are needed, this filtering procedure requires reconstructing the entire multilevel hierarchy at each iteration, i.e., recomputing interpolation and the coarse grid matrices on all levels.  Thus, these approaches are efficient when one (or few) components are needed to define an effective interpolation operator, but can be costly in cases where many such components are needed to obtain full MG efficiency.

The two main ideas in bootstrap AMG that distinguish the approach from other Adaptive AMG methods are: (1) the idea of defining restriction and interpolation to fit a collection of smooth vectors in a least squares sense, and (2) the idea of computing the collection of test vectors simultaneously in a multigrid scheme, where the bulk of the work is done on coarse, i.e., on cheap levels.
The BAMG algorithm, described below, combines the general devices of (A) compatible relaxation and algebraic distances to choose the coarse variables
with (B) a least squares procedure and (C) a bootstrap cycling scheme for computing the test vectors
used in the LS process to construct high-quality restriction and interpolation operators
and thereby a multigrid hierarchy.
Before giving an overview of those Bootstrap techniques, we remark here that they should be used as needed, e.g., in cases where geometric MG and AMG solvers are not available and/or the additional costs they introduce can be amortized, such as problems where the same system has to be solver for multiple right hand sides.}

\bi
\item[(A)]  The coarse variables are selected using compatible relaxation
and algebraic distances, as tools for gauging the quality of the coarse variable set
and guiding the coarsening process, respectively.  \ei

A general criterion for measuring the quality of the set of coarse
variables $C$ is the fast convergence of compatible relaxation
(CR), as introduced in \cite{etna2000} (being a special case of
choosing coarse variables for much more general types of systems
\cite{SU}, introduced first for problems in statistical mechanics
\cite{RMG}).  Generally, CR refers to a relaxation scheme that keeps
the coarse-grid variables invariant.  {A general approach introduced by
Brandt and Livne (see also~\cite{HCR}), that can be applied
to any relaxation scheme is
``habituated CR'' (HCR), which uses the full fine-grid relaxation scheme and
then imposes certain constraints on the coarse-grid variables.
The HCR criterion allows choosing a set of coarse variables which
is not necessarily a subset of the set of fine-level variables. In
particular, coarse variables can represent some {\it averages} of
fine variables, as needed for example for equation homogenization
or for coarsening nonlinear problems. Another important type of a
set of coarse variables whose choice can be guided by CR arises in
{\it wave equations}, where it should include several different
subsets, each for example corresponding to a different wave
direction.
Moreover, given a suitable definition of the averages used in defining
coarse variables,
the HCR convergence rate yields an
accurate prediction of the convergence rate that can be achieved
by a two-grid solver for a given relaxation scheme and the
proposed set of coarse variables~\cite{JBrannick_2005a,JBrannick_RFalgout}. In this way, the HCR quantitative prediction
is also very useful in designing and debugging the actual solver.}

Algebraic distance gives a general approach
for measuring the connectivity among nodes in a graph.
In AMG, the algebraic distance measure
 is based on a notion of strength of connectivity among variables that is derived by minimizing a local least squares (LS) problem which is also used to derive interpolation (described in Section~\ref{sec:LS}). The approach first constructs direct (caliber-one) LS interpolation for a given set of test vectors (representatives of algebraically smooth error obtained by the bootstrap process described
in Section~\ref{sec:BS}) and then defines the algebraic distance between a fine point and its neighboring points in terms of the values of the local least squares functionals resulting from the so defined interpolation. The algebraic distance measure thus aims to address the issue of strength of connections in a general context -- the goal being to determine explicitly those degrees of freedom from which high quality least squares interpolation for some given set of smooth test vectors can be constructed.

An algorithm using algebraic distance was first developed for combinatorial optimization problems on graphs~\cite{SRB09,RonSB11}.   In these works, algebraic distance is developed as a measure of the closeness between nodes that is then used in an aggregation algorithm to build multilevel representations of graphs.
In~\cite{BAMG_Aniso}, algebraic distance was developed as a measure of strength of connectivity
among variables for coarsening non-grid-aligned anisotropic diffusion problems.  Further, it was shown that combing the measure with compatible relaxation results in a robust algorithm for choosing
interpolation operators for such problems.  And an affinity measure based on algebraic distance
together with an energy measure
are used to form aggregates in the Lean AMG solver developed in~\cite{LAMG_Report} for the graph Laplacian problem.
Related work is found in~\cite{brannick_local_stab_2011}.

\bi
\item[(B)] {The interpolation operators (as well as restriction operators for non-Hermitian problems) are derived by a least-squares
fit to a set of $\rho$-smooth Test Vectors (TVs) with
sufficiently small $\rho$}. \ei

Denote by $C_i$ the set of coarse variables used in
interpolating to a fine grid variable $i\in\finevar$.
Then the main idea in the least squares interpolation approach is to compute for each
$i\in \finevar$ a row of interpolation that approximates
the given set of test vectors at the points in $C_i$, minimizing the interpolation error for these
vectors in a least squares sense.
It follows from the
satisfaction of the CR criterion that, with a proper choice of
$C_i$, there exists an interpolation that will have only $O(\rho)$
errors in reproducing $\xk_i$ for {\it any} normalized vector $x^{(k)}$ which
is $\rho$-smooth. The size $|C_i|$ of this set should in principle
increase as $\rho$ decreases, but in practice a preset
interpolation {\it ''caliber''} $\caliber$ yields small enough
errors. For example, for a discretized scalar PDE in
$d$-dimensions, a proper set with $|C_i|=d+1$ points would give
$O(h^2+\rho)$ accuracy, which may suffice if $\rho$ is small
enough (see below). Typically, the number of TVs need not be larger than $1.5
\caliber$ for this accuracy to be attained for almost every
$\xk_i$.

The set $C_i$ can often be adequately chosen by geometric
considerations, such as the set of geometrical neighbors.
If the chosen set is inadequate, the
least-squares procedure will show {\it bad fitness} (interpolation
errors large compared with $\rho$), and the set can then be
improved. The least-squares procedure can also detect variables in
$C_i$ that can be discarded without a significant loss in accuracy.
Thus, this approach allows creating interpolation with whatever
needed accuracy which is {\it as sparse as possible}. Moreover, the LS
derivation of interpolation is applicable to {\it general types}
of coarse variables (e.g., averages). \bi \item[(C)] {The TVs
are constructed in a bootstrap manner}, in which several tentative
AMG levels generated by fitting interpolation to
moderately smooth TVs are used to produce {\it improved} (much
smoother) TVs, repeating this if needed until fully efficient AMG
levels have been generated. \ei

Quantitatively, the algorithm is driven by the following
relations.
The first TVs are each produced by relaxing the homogeneous system
$Ax=0$ with a different starting vector. This quickly leads to a
set of $\rho$-smooth TVs with $\rho \ll 1$.  A combination of random
vectors and additional geometrically derived smooth vectors can generally be
used as starting vectors.  In the case of discretized isotropic
PDEs, if geometrically smooth vectors that satisfy the homogeneous
boundary conditions are used, relaxation may not be needed at all.
In many other cases relaxation can be confined to the neighborhood
of boundaries and discontinuities. As mentioned, the
least-squares fitting to such a $\rho$-smooth set produces an
interpolation scheme whose errors in interpolating {\it any} normalized
$\rho$-smooth vector are at most $O(\rho)$. The coarse-level
Galerkin system based on such an interpolation will then yield an
$O(\rho^2)$-accurate approximation to every normalized eigenvalue
of a $\rho$-smooth eigenvector. Increasingly coarser AMG levels
can recursively be built in this way, all typically having such
$O(\rho^2)$ accuracy for eigenvectors with small normalized eigenvalues
and their linear combinations.

The AMG cycles based on these levels can further be applied to the
homogeneous system $Ax=0$.

In many cases, the convergence (to
$x=0$) will exhibit the desired multigrid efficiency (close to the
HCR prediction), indicating that the produced AMG structure is
ready for production. If the given method needs further improvement, these cycles will quickly produce
$O(\rho^2)$-smooth TVs, which can then be used for a second round
of constructing the AMG hierarchy, potentially (with suitable
sets $C_i)$ providing $O(\rho^4)$-accuracy. More such rounds of TV creation, to obtain
still higher accuracy, is seldom needed.

Some important
applications do result in nearly singular systems, e.g., the discretized Dirac equation arising in
lattice quantum chromodynamics.
For such problems the system matrix has several almost
zero modes (AZMs), which are a few special, nearly-singular eigenvectors
with normalized eigenvalues much smaller than usual. For example, in discretizing a $q$-order PDE on a (structured or unstructured) grid with $O(h)$ meshsizes, the lowest normalized eigenvalues will all usually have size $O(h^q)$; but by a special adjustment of a parameter in the PDE, one normalized eigenvalue can be made arbitrarily small, having size $\eta \ll h^q$.

Usually, the potential existence of one, or sometimes several,
AZMs is known in advance, as a property of the problem at hand. To
obtain $O(\eta)$ approximations to normalized eigenvectors would
generally require substantially increasing the number of interpolation points, or caliber
denoted by $\caliber$. However, since the number of AZMs is typically small, the
interpolation can be tailored to fit them individually
sufficiently well, without increasing its caliber. This will in
fact automatically happen as soon as the AZMs emerge (in the
cycles applied to the homogeneous system), provided the
least-squares fitting weights each test vector
proportionally to the inverse of its normalized residual norm.
We note that a multilevel eigensolver can be integrated into
the bootstrap setup process to further accelerate the convergence
of the setup for the AZMs, as introduced in~\cite{BAMG2010}.

In the remainder of the paper, an overview of the main components of the BAMG approach are given.
Section \ref{sec:CR} contains a review of compatible relaxation and algebraic distances.
The least squares approach for computing interpolation is discussed in Section \ref{sec:LS}
and the bootstrap procedure aimed at computing increasingly accurate TVs is
discussed in Section \ref{sec:BS}.  A main focus in the presentation in these
sections is to point to open research questions in the development, analysis and
application of BAMG and to the most pressing research directions on these topics.
Promising future research areas in BAMG are outlined in Section \ref{sec:FD}.

\section{Compatible relaxation and algebraic distance}\label{sec:CR}
The first task in coarsening a given system is the choice of coarse-level variables, which in the BAMG
framework is accomplished using compatible relaxation and algebraic distances.
Compatible relaxation coarsening algorithms have been developed and analyzed for both classical
AMG~\cite{HCR,PanayotRob_2003,JBrannick_RFalgout,Brannick_Trace_06} and aggregation-based coarsening~\cite{Jacob}.  These works have demonstrated the utility of
CR as a tool for guiding the coarsening process and assessing the quality of a given coarse variable set.
The works in~\cite{HCR,JBrannick_RFalgout} focused on adding robustness to classical AMG using CR coarsening approaches, in particular for anisotropic diffusion problems with highly oscillatory and strongly anisotropic diffusion coefficient.  The method developed in~\cite{Jacob} does not use CR directly in the setup algorithm, however, CR was used as the main tool for the analysis and design of the resulting algorithm, which
was developed specifically for general non-grid-aligned anisotropic diffusion problems.

An outline of a generic coarsening process based on CR for the linear system $Ax=b$ is given in Algorithm
\ref{alg:CR}.  Here, $\rho_{cr}$ is an approximation to the asymptotic convergence rate of CR computed after the $\nu$ iterations are applied and $\delta$ is a tolerance that can be set by the user or more generally it can be optimized within the algorithm.  In~\cite{HCR}, the CR coarsening algorithm is designed to approximately optimize $\beta = \rho_{cr}^{\frac{1}{W}}$, where $W$ denotes the total work of a single AMG cycle that is estimated using the  coarsening ratios in different stages of a CR coarsening algorithm.
Starting with an initial set $C_0$ of coarse variables, whose values are held fixed at zero, a few sweeps of compatible relaxation (with vanishing right-hand side so that the solution should converge to zero, but starting of course from a suitable non-zero first approximation) are applied to detect convergence slowness if the coarse-variable set is inadequate. In this case, a set of variables that should be added to the coarse set needs to be
determined.   The typical approach is to select a subset of the set of variables which are slow to converge (to zero) and add these to the set of coarse variables.  As in AMG, the additional variables to be added to $C_0$ can be chosen using a notion of connectivity of a variable to its neighbors, which is defined via algebraic distance in BAMG.

\begin{algorithm}
\caption{: Compatible relaxation coarsening, \textit{Input:} $\coarsevar_{0}$ \textit{Output:} $\coarsevar$\label{alg:CR}}
\begin{algorithmic}
\STATE Initialize $\coarsevar = \coarsevar_{0}$
\STATE Perform $\nu$ CR iterations
\WHILE{$\rho_{cr} > \delta$}
\STATE Update $\coarsevar$
\STATE Perform $\nu$ CR iterations
\ENDWHILE
\end{algorithmic}
\end{algorithm}

{In order to illustrate the flow of the above CR coarsening process, we report results from~\cite{JBrannick_RFalgout} for a specific implementation of the algorithm applied to the anisotropic diffusion problem
\begin{equation}\label{eq:jumps}
 - \nabla \cdot \mathcal{K}(x,y) \nabla u(x,y) = f(x,y),
 \end{equation}
where ${\mathcal{K}}$ is a tensor of the form
\begin{equation} \label{eqn-2D-tensor}
{\mathcal K} = \left[
\begin{array}{cc}
a & c \\
c & b
\end{array}
\right] ,~~~
\left\{
\begin{array}{l}
a = \cos^2(\theta) + \epsilon \sin^2(\theta) \\
b = \epsilon \cos^2(\theta) + \sin^2(\theta) \\
c = (1-\epsilon) \cos(\theta) \sin(\theta)
\end{array}
\right. .
\end{equation}
The parameter $0 < \epsilon \leq 1$ specifies the strength of anisotropy in the
problem, while parameter $\theta$ specifies the direction of anistropy.
We use 2D domain $\Omega = [0,1] \times [0,1]$ with Dirichlet boundary conditions
with the following choice of problem parameters
\begin{equation*}
\begin{array}{ll}
\epsilon=1,~ d=10^4                &~~~ (x,y) \in {[0,1/2]} \times {[0,1/2]}, \\
\epsilon=1,~ d=0                   &~~~ (x,y) \in {(1/2,1]} \times {[0,1/2]}, \\
\epsilon=0.01,~ \theta=0,~ d=0     &~~~ (x,y) \in {[0,1/2]} \times {(1/2,1]}, \\
\epsilon=0.01,~ \theta=\pi/2,~ d=0 &~~~ (x,y) \in {(1/2,1]} \times {(1/2,1]}.
\end{array}
\end{equation*}
The discrete system is obtained by a finite element discretization of this problem with bilinear basis functions on uniform quadrilateral mesh.  The results of the algorithm for this test problem are provided in Figure~\ref{fig-cr-4reg};
we refer to~\cite{JBrannick_RFalgout} for additional
numerical experiments and results obtained using CR in a classical AMG
setup algorithm.
In the reported results, an $F$-point Gauss Seidel
variant of CR is used, the initial set $\coarsevar_0=\emptyset$ and the updating of $\coarsevar$ in the subsequent stages of the algorithm is based on adding an independent set of coarse-point candidates that are slow to converge in the CR process.
The precise implementation of the algorithm is given in~\cite{JBrannick_RFalgout}, see Algorithm 3.1
in Section 3.

Figure~\ref{fig-cr-4reg} illustrates the
coarsening stages and the final coarse grid chosen by the algorithm.  Because $\rho_{cr}
> 0.7$, the algorithm enters the first coarsening stage.  Figure~\ref{fig-cr-4reg}(a) contains a plot of the
values of the candidate measures, scaled pointwise error values, that are used in selecting
coarse-grid candidates.  The candidate measures are large
in three regions of the domain, signaling that coarsening is needed there.  The
algorithm also correctly detects that coarsening is not needed in the bottom
left region of the domain where pointwise Gauss-Seidel is effective. The first
$C$-point set chosen is depicted by dotted squares in
Figure~\ref{fig-cr-4reg}(b).  For these updated sets, we
obtain $\rho_{cr} = 0.91$, which leads to a second stage.

The candidate measures at stage two of the algorithm are plotted in
Figure~\ref{fig-cr-4reg}(b) and the
algorithm signals that still further coarsening is needed only in the top two regions
of the domain.  The
result of the second stage is the coarse grid depicted in
Figure~\ref{fig-cr-4reg}(c).  The new CR convergence factor is $\rho_{cr}  = 0.49$,
so there are no additional stages in the algorithm.  Notice that the algorithm gives
semi-coarsening in $x$ in the top-left region where the problem is anisotropic
in $x$, semi-coarsening in $y$ in the top-right where the problem is anisotropic
in $y$, full-coarsening in the bottom-right where the problem is isotropic, and
no coarsening in the bottom-left where the mass term dominates.  Further,
we mention that the algorithm does not make explicit use of an strength of connectivity
meausre.

\begin{figure}[ht!]
\begin{center}
\begin{tabular}{@{}c@{}}
\includegraphics[width= 0.48\linewidth]{cr-4reg-dir-qua-033-stage_1} \\ (a)
\end{tabular}
\hfill
\begin{tabular}{@{}c@{}}
\includegraphics[width= 0.48\linewidth]{cr-4reg-dir-qua-033-stage_2} \\ (b)
\end{tabular}
\vspace{1em} \\
\begin{tabular}{@{}c@{}}
\includegraphics[width= 0.48\linewidth]{cr-4reg-dir-qua-033-grid}
\end{tabular}
\end{center}
\caption{
Candidate measures computed during
stage one (a) and two (b) for the 2D anisotropic diffusion problem.  Each square depicts the value
of the scaled pointwise error at grid point $i$ according to the color bar.  Dotted squares
indicate points currently in $\coarsevar$.  The final coarse grid is
given by the larger red circles in (c).}
\label{fig-cr-4reg}
\end{figure}
}

Generally, while CR algorithms have been successfully applied and developed for various model problems, there remain several design issues that must be addressed and remain open in developing a robust implementation
of the above CR algorithm.
One issue that arises is formulating an appropriate relaxation scheme for the problem
being solved.  If specific features of the problem are known,
e.g., it comes from a system of PDEs, then an appropriate smoother
can be designed using knowledge of the problem.  More generally,
the smoother can be adapted within the CR coarsening process.  For example,
in cases where $|\coarsevar|$ or $\rho_{cr}$ are too large, a block relaxation scheme can be used and the
blocks can be adapted using CR to test the efficiency of the smoother.
This is research topic which has not been explored much in the AMG setting,
in particular, together with CR and algebraic distance.   {We mention the recent work
on adaptive smoothers in~\cite{Adap_Smooth} that shows the potential of such
an approach.}

The two main research questions in CR that require further investigation and with which we conclude
our discussion of CR are:

\begin{itemize}
\item  How to define an appropriate form of CR and related convergence rate estimate from $\rho_{cr}$
that accurately predicts the performance of the resulting two-grid method for a given interpolation strategy?
\item  How to choose $C_0$ and how to identify suitable local modifications to a given set of coarse variables
whenever the CR convergence rate is slow?
\end{itemize}

The fact that fast CR rates imply the suitability of a given coarsening scheme has been established
theoretically~\cite{PanayotRob_2003} and by the extensive numerical experiments found in~\cite{JBrannick_RFalgout,Jacob}.  However, a precise relation between the convergence rate
of a given CR iteration and the resulting two-grid convergence is not known in general.
Specific versions of CR have been proposed that do provide
accurate predictions of the convergence rate of the associated two grid method, {\color{black}for example,
$F$-relaxation and} Habituated CR.  A single HCR iteration is given by relaxation steps applied
to the entire fine-grid system, followed
by imposing constraints on the resulting error restricted to the coarse variables:
$$e_C := Xe = 0,$$
where $X$ is a matrix that defines
the averaging process for deriving coarse variables.
The solution to this problem can be calculated using the Kaczmarz iteration, starting with the error $e$ that results from the full relaxation steps as an initial guess or it
can be computed using the $\ell_2$ projection of $e$ on $\text{range} (X)$.  The latter approach
gives the overall HCR error propagation operator
\begin{equation}\label{eq:HCR}
 E = (I-\pi)(I - MA)^{\nu}, \quad \pi = X^H(XX^H)^{-1}X.
 \end{equation}
Unlike the $\finevar$-relaxation version of CR which can only be applied with the special coarsening used in classical AMG, and then only for simple relaxation schemes, HCR is applicable in general setting since it only involves
the fine-grid relaxation scheme followed by a simple projection based on the averaging matrix $X$.   Hence, in addition to giving a better qualitative prediction of the performance of the two-grid solver, HCR is more flexible and easier to process.

{\color{black}The sharpness of HCR in predicting the convergence rate
  of a two-grid method depends critically on the averaging matrix $X$
  used in enforcing the constraints and, as
shown in~\cite{JBrannick_RFalgout}, a suitable choice can be defined
by $X=P^H$, where $P$ denotes a sufficiently accurate and stable
interpolation operator. We note that this $P$ does not need to be the
same interpolation operator that is used in the AMG solver.
Generally, the accuracy of the iteration in~\eqref{eq:HCR} in
predicting the convergence rate of the two-grid method depends on the
energy norm of $\pi$, i.e., $\|\pi\|_A\geq 1$
(see~\cite{JBrannick_RFalgout}).  We mention, in addition, the work
in~\cite{SMacLYS07}, where it is shown that if the $C/F$ partitioning
is such that the weighted Jacobi $F$-relaxation converges quickly,
then a bound on the two-grid convergence rate follows.}

The issue of identifying suitable local modifications to the coarse variable set
whenever the CR convergence rate is too slow is a more subtle issue.  Assume
that the relaxation scheme is fixed, $C$ is nonempty,
and the convergence rate of CR is slow, indicating that additional variables need to be added
to $C$.
The usual approach is to find a candidate set, $\bar{F}$, of
$F$-points from which an independent set is added to the existing set
$C$.
The primary issue is to determine a set of candidates $\bar{F}$.
The basic observation from~\cite{etna2000} that provides insight on this topic is that
a set of variables that should be added to the coarse variable set is exposed by the CR process; it is a subset of the set of variables which are slow to converge to zero when CR is applied.
This observation motivated the CR algorithms developed and analyzed in~\cite{HCR,JBrannick_RFalgout}.

In~\cite{JBrannick_RFalgout}, a candidate set measure based on the values of the normalized pointwise
errors that result from the CR process was studied.  The set $\bar F$ is formed by adding to it $i\notin \coarsevar$ for which this error is large.   Intuitively, the approach is based on choosing additional coarse grid points in regions where the CR process is inefficient.
In~\cite{HCR}, a similar algorithm was developed for the matrix $\tilde A$ that is first normalized such that each of its rows has $\ell_1$-norm equal to $1$.
While the effectiveness of these strategies has been demonstrated by the numerous numerical experiments reported in~\cite{HCR,JBrannick_RFalgout},
there is still no rigorous understanding of or justification for using pointwise errors of CR to determine the variables that will most effectively reduce the convergence rate of the CR process when added to $C$.

Very recent work~\cite{BAMG_Aniso} on this topic has focused on incorporating the use of algebraic distance as a measure of strength of connectivity into the CR coarsening process.
The main idea in this context is to use algebraic distance to form a
strength of connectivity graph that identifies the coupling among
variables which can then be used in an independent set algorithm to
update $C$.
In the simplest form, the notion of algebraic distances is
{\color{black}straightforward: For any given pair of fine-grid variables $i,j \in \Omega$ compute
\begin{equation}\label{eq:AD}
 \mu_{ij}  =  \frac{
\sum_{\kappa = 1}^k\left(v_{i}^{(\kappa)} -  p_{ij} v_{j}^{(\kappa)}\right)^{2}}
{\sum_{\kappa = 1}^k \bigg(v^{(\kappa)}_i\bigg)^2},
\end{equation}
where $p_{ij}$ is the coefficient that minimizes this least squares
problem for caliber-one interpolation of the form $v_i = p_{ij} v_j$.
Here, the set $\tvV = \{v^{(1)}, \ldots,
v^{(k)}\}$ consists of properly scaled smooth TVs, and the choice of
the denominator ensures that the expression is symmetric with respect
to the indices $i,j$.
Statistically, the measure $\mu_{ij}$ is based on the idea that
the nodes $i$ and $j$ in a graph are strongly connected only if $v_i$ and $v_j$ are highly correlated for all TVs.

Another, more direct, way to write the algebraic distance measure is given by solving the LS problem \eqref{eq:AD} for the coefficient
$p_{ij}$ and substituting the minimizer into the definition of $\mu_{ij}$.
The resulting
$p_{ij} = \tfrac{\langle V_{(i)},V_{(j)}\rangle}{\langle V_{(j)},V_{(j)}\rangle}$ gives a rise to  an equation for the minimum of \eqref{eq:AD}:
\begin{equation}\label{eq:AD2}
  \mu_{ij} = 1- \frac{\big|\big\langle V_{(i)},V_{(j)}\big\rangle\big|^2}{\|V_{(i)}\|^2\|V_{(j)}\|^2},
  \end{equation}
  where
with  $V = [v^{(1)} \: ... \:v^{(k)}] \in \mathbb{C}^{n\times k}$,
where    $\langle V_{(i)},V_{(j)}\rangle := \sum_{\kappa = 1}^k
v_i^{(\kappa)}v_j^{(\kappa)}$ measures the angle between $ V_{(i)}$
and $ V_{(j)}$.  From this explicit form of the minimum
  it is now easy to see that the measure $\mu_{ij}$ is symmetric with respect to $i,j$, as noted above.
}
{The algebraic distance measure is a specific case of the least squares interpolation process described in more detail in the next section.}
More generally, the algebraic distance measure can be computed for sets of neighboring coarse points $\coarsevar_i$, using the more general
LS problem defined in~\eqref{eq:LSfuncrowi}.  In this way, $\mu_{ij}$ can be used as an a posteriori measure of the quality of the interpolatory set $\coarsevar_i$ and the resulting row of interpolation.
For cases where an algebraic notion of dependence is not sufficiently well defined, e.g., in the early stages of the BAMG setup process when the TVs may not be smooth enough to yield an accurate algebraic distance measure, one can assume any non-vanishing dependence to be strong in the CR coarsening process; this may result in slow coarsening, but then the
scheme can be applied again in an iterative manner, to this new set instead of $C_0$.  The latter CR approach that does not rely on strength of connectivity to choose $\coarsevar$ was studied in the classical AMG setting in~\cite{JBrannick_RFalgout} and its extension that incorporates the algebraic distance measure in \eqref{eq:AD} was developed for anisotropic problems in~\cite{BAMG_Aniso}.

Generally, such algorithms that use iterative strategies for choosing $C$ together with CR often lead to irregular alignment of coarse variables, even for structured grid problems, which in turn can result in an increase in the number of nonzeros in the coarse-grid matrix $A_c$.   This issue can be treated in latter stages of the BAMG process by using more aggressive coarsening, once the TVs become more smooth and the algebraic distance measure becomes more precise.  This is the case, for example, for strongly anisotropic problems where neighbors in the characteristic direction are much closer than in the other directions.  For such problems the appropriate long range connections in these characteristic directions are revealed only by very low energy TVs, which normally need to be the product of subsequent bootstrap setup cycles.
An algorithm that  both removes and adds variables to $C$ in each CR coarsening stage can be used to further control complexities, as needed.  We note that in our experience, for many problems
only the first coarsening stage in which $C_0$ is formed is needed, given sufficient accuracy of the algebraic distance measure.

\section{Least squares interpolation}\label{sec:LS}
After the coarse grid $\coarsevar$, has been selected by compatible relaxation and algebraic distance,
the interpolation operator is computed using a least
squares process. The least squares interpolation operator $P$ is defined to fit
collectively a set test vectors that should characterize the near kernel of the
system matrix.
Assuming the sets of interpolatory variables, $\coarsevar_i$, for each $i\in \finevar$, and a set of test vectors, $\tvV = \{v^{(1)}, \ldots,
v^{(k)}\}$, have been determined, the $i$th
row of $P$, denoted by $p_i$,
is defined as the minimizer of the local least squares problem:
\begin{equation}\label{eq:LSfuncrowi}
\mathcal{L}(p_i) =
\sum_{\kappa=1}^k\omega_k\left(v_{\{i\}}^{(\kappa)} - \sum_{j\in \coarsevar_{i}} \left(p_{i}\right)_{j} v_{\{j\}}^{(\kappa)}\right)^{2} \rightarrow \min.
\end{equation}
Here, the notation $v_{\widetilde{\Omega}}$
denotes the canonical restriction of the vector $\widetilde{v}$ to
the set $\widetilde{\Omega} \subset \Omega$, e.g.,
$v_{\{i\}}$ is simply the $i$th entry of $v$.
Conditions on the uniqueness of the solution to minimization problem
\eqref{eq:LSfuncrowi} and an explicit form of the minimizer have been
derived in \cite{BAMG2010}.

The weights $\omega_{\kappa}>0$ are defined using the $A$-norms of the test vectors $\|v\|^2_A = \langle
Av,v\rangle$ when the system matrix $A$ is Hermitian and
positive definite (HPD)).

 A main issue with this choice of weighting is that
the LS process uses the entries of the test vectors restricted to
{\em local} coarse neighborhoods to define the LS fit while the weights $\omega_{\kappa}$ are defined in terms
of the {\em global} test vectors.  Hence, it may happen that a given test vector is globally smooth (i.e.,
produces a large weight and is emphasized in the LS process), but it does not provide an adequate local representation of smooth error.  This leads to a loss in accuracy of interpolation and slows down the convergence
of the resulting two-grid method.

\begin{figure}
  \centering
  \begin{picture}(0,0)
\includegraphics{interprelbw.pdf}
\end{picture}
\setlength{\unitlength}{3947sp}

\begingroup\makeatletter\ifx\SetFigFont\undefined
\gdef\SetFigFont#1#2#3#4#5{
  \reset@font\fontsize{#1}{#2pt}
  \fontfamily{#3}\fontseries{#4}\fontshape{#5}
  \selectfont}
\fi\endgroup
\begin{picture}(6777,3624)(2389,-3973)
\put(9151,-1636){\makebox(0,0)[lb]{\smash{{\SetFigFont{14}{16.8}{\rmdefault}{\mddefault}{\updefault}{\color[rgb]{0,0,0}$i \in \finevar$}
}}}}
\put(9151,-1036){\makebox(0,0)[lb]{\smash{{\SetFigFont{14}{16.8}{\rmdefault}{\mddefault}{\updefault}{\color[rgb]{0,0,0}$i \in \coarsevar$}
}}}}
\end{picture}

\caption[Full coarsening and interpolation relations]{Full coarsening and interpolation relations for $i\in \finevar \setminus \coarsevar$.}
\label{fig:nr:scalar:interprel}
\end{figure}
{In~\cite{iBAMG} (see also~\cite{Brezina:2012:RCB})}, a residual-based LS form of BAMG interpolation (LSR) was introduced as a technique for
addressing this issue of using global weights.  The basic idea in the LSR process
is to implicitly smooth the test vectors to produce an improved local representation
of smooth error.  In this respect, the approach is related to the adaptive relaxation scheme described by Brandt in~\cite{BAMG1} (and suggested in~\cite{etna2000}) that applies relaxation selectively to points exhibiting especially large residuals, as they occur for example in problems with singularities.
Specifically, in LSR interpolation the test vectors are (temporarily) updated as follows:
\begin{equation}\label{eq:LSfsmooth}  v_{\{i\}}^{(\kappa)} \leftarrow  v_{\{i\}}^{(\kappa)} - \frac{1}{d_{ii}}  r_{\{i\}}^{(\kappa)}, \end{equation}
and then used in computing $p_i$, again as the solution to \eqref{eq:LSfuncrowi}.

Table~\ref{tab:nr:lsinterpketaN01} contains results from~\cite{BAMG2010} for a variational two-grid method resulting from using standard (full) coarsening and interpolation with $|\coarsevar_i| \leq 4$ in \eqref{eq:LSfuncrowi} (see Figure~\ref{fig:nr:scalar:interprel}) for the original LS formulation and the LSR version
applied to the central finite difference discretization of the Poisson problem with homogenous Dirichlet boundary conditions on a uniform quadrilateral grid.

The reported results are the approximate asymptotic convergence rates for different values of the number relaxed test vectors, $k_r$, and the number of Gauss Seidel iterations, $\eta$, used to compute them.  {Note that here the total number of TVs is set to $k=k_r$, whereas, when a multilevel eigensolver is combined with relaxation
{\color{black}(an approach that we consider in the next section)}
the total number of TVs is given by $k=k_r+k_e$, with $k_e$ the number of TVs computed by the multilevel eigensolver described in the following section.}  The test vectors are computed by applying Gauss Seidel iterations with lexicographic ordering to the Laplace problem starting with $k$ distinct initial guesses, generated randomly with a normal distribution with expectation zero and variance one, $N(0,1)$, and the update in \eqref{eq:LSfsmooth} used in the LSR formulation is applied only to $20\%$ of the entries of the TVs for which the associated values of the residual $r_{\{i\}}^{(\kappa)}$ are largest in absolute value.

\begin{table}
  \begin{center}
    \begin{tabular}{|c||c|c|c|c|c|c|c|}
      \hline
     $k$ & $6$ & $7$ & $8$ & $10$ & $12$\\
             \hline
        \hline
       $\eta =2$ & $.934$    $(.866)$  &  $.916$    $(.827)$  &  $.904$ $(.806)$ &   $.878$    $(.786)$    & $.861$    $(.759)$  \\
      $\eta =4$ & $.735$    $(.524)$ &   $.681$    $(.453)$  &  $.648$    $(.403)$  &  $.608$    $(.269)$ & .575 (.316) \\
      $\eta =6$ & $.522$    $(.294)$   & $.451$    $(.242)$  &  $.425$    $(.205)$   & $.381$    $(.154)$ &.305 (.136) \\
      $\eta =8$ & $.372$    $(.169)$  &  $.308$    $(.140)$  &  $.287$    $(.123)$   & $.238$   $(.054)$  & .216 (.061)\\
      \hline
    \end{tabular}
    \caption[Two-grid LS interpolation forPoisson probelm ($k$ vs.\ $\eta$)]{Asymptotic convergence rate estimates, $\rho$, of the two-grid methods with two pre- and two post-smoothing steps constructed using LS (residual-based LS) interpolation for the standard finite difference discretization of the Poisson problem with homogenous Dirichlet boundary conditions on a uniform rectangular grid with $h=1/64$.
     }\label{tab:nr:lsinterpketaN01}
  \end{center}
\end{table}

The improved performance of the method that results from increasing $\eta$
is perhaps expected, since applying more relaxation sweeps to the test vectors should
give a better approximation of smooth error.  The situation gets more interesting and more  complicated when we fix the number of Gauss Seidel iterations and then increase the sampling space (i.e, vary the number of test vectors), which results in an ad hoc (randomized) approach for sampling the dominant part of the space that is not significantly reduced by the smoother.
Recently, randomized methods for range approximation
have been given particular attention in data mining applications and further studies
of the LS process along these lines are still needed.  {The paper~\cite{antrg} summarizes the latest
results obtained on this topic and highlights the gains which are possible with judicious modifications of the
ad hoc approach. }

\begin{table}[!ht]
  \begin{center}
    \begin{tabular}{|c|c|c|c|c|c|c|}
      \hline
      $h$ & $1/32$ & $1/64$ & $1/128$ & $1/256$ & $1/512$ \\
      \hline
    $\rho$ & $.267$ $(.075)$ & $.648$ $(.403)$ & $.886$
      $(.758)$ & $.966$ $(.917)$ & $.990$ $(.977)$\\
            \hline
    \end{tabular}
    \caption[Two-grid LS interpolation ($k_{r}=8, \eta=4$) for FD Laplace]{Asymptotic convergence rate estimates, $\rho$, of the LS (LSR) based two-grid methods with two pre- and two post-smoothing steps applied to the same FD Laplace problem as in Table 1.  The solver is constructed using the LS (LSR) schemes with $\eta=4$ and $k_{r}=8$.\label{tab:nr:lsinterp84N01}}
     \end{center}
\end{table}

In Table~\ref{tab:nr:lsinterp84N01}, we report results from~\cite{BAMG2010} for a two-grid method applied to the same central finite difference Poisson problem for different choices of the grid spacing $h$.  The solver is constructed using the LS and LSR formulations for defining $P$ with a fixed number of relaxed test vectors $k_{r} = 8$ and smoothing iterations
$\eta = 4$ independent of the choice of $h$.  Here, we see that the convergence rates of both the LS and LSR based two-grid methods increase
as the size of the problem is increased, suggesting that applying a constant number of relaxation steps to compute a fixed number of test vectors does not produce a sufficiently accurate local representation of  the algebraically smooth error for defining least squares interpolation.  {The numerical results reported in~\cite{iBAMG,Brezina:2012:RCB}
further support this claim.}

In Table~\ref{tab:nr:lsinterp714N01}, we present results of the two-grid solvers obtained by adding the constant vector, ${\bf 1}$, to $\mathcal{V}$ and again applying the LS and LSR schemes to the FD Poisson problem.  Here a fixed number of
test vectors $k_{r} = 7+1$ are used and $\eta =4$ relaxation steps are used to compute them.
\begin{table}[!ht]
  \begin{center}
    \begin{tabular}{|c|c|c|c|c|c|c|}
      \hline
      $h$ & $1/32$ & $1/64$ & $1/128$ & $1/256$ & $1/512$ \\
      \hline
       $\rho$ & $.089$ ($.039$) & $.121$ ($.040$) & $.130$
      ($.042$) & $.148$ ($.042$) & $.150$ ($.043$)\\
      \hline
    \end{tabular}
    \caption[Two-grid LS interpolation ($k_{r}=7+\mathbf{1}, \eta=4$) for FD Laplace]{Asymptotic convergence rate estimates, $\rho$, of the LS (LSR) based two-grid methods with two pre- and two post-smoothing steps applied to the same FD Laplace problem as in Table 1.  The solver is constructed using the LS (LSR) schemes using $\eta=4$ Gauss Seidel iterations to compute the seven initially random test vectors; the additional test vector is the constant vector chosen a priori.
      \label{tab:nr:lsinterp714N01}}
    \end{center}
\end{table} With this modification the performance of the LS and LSR based
two-grid solvers improve significantly.  In fact, both the LS and LSR methods perform nearly optimally for the problem sizes
considered.  This clearly demonstrates the ability of the LS process to incorporate geometric and spectral information from the problem when it is available.  We note, in addition, that with our choice of full coarsening and interpolation relations (see Figure~\ref{fig:nr:scalar:interprel}), the
resulting interpolation is not derived from immediate neighbors, which
further demonstrates the flexibility of the least squares formulation in that it naturally
accommodates long-range interpolation.
We refer to the paper~\cite{BAMG_Aniso} in which algebraic distance and long-range LS interpolation were studied for
non-grid-aligned anisotropic problems for additional details and results on this topic.

The LSR process provides a simple technique for improving the accuracy of interpolation in cases where the global weights are large but the TVs are not sufficiently smooth in some local regions.  As the numerical experiments show, the main point is to obtain extra smoothness (i.e., smaller residuals) at the points
where interpolation is defined.  And such large residuals in general are the result of special local configurations in the random vectors used as initial guesses.  This issue can also be inexpensively treated using adaptive relaxation together with several short $V$-cycles based on the
interpolation that is constructed for the first few levels.

The related problem of determining appropriate local weights in the LS definition of interpolation is more
complicated and has not been resolved in a general algebraic setting, except for finite element discretizations of elliptic
boundary value problems.  For FE systems, the local weights can be computed using the local stiffness matrices.
The problem of determining local weights is also related to the question
of deriving local equations for computing the test vectors, as in spectral AMG~\cite{rhoAMGe}.

\section{A bootstrap algorithm for computing the test vectors}\label{sec:BS}
The BAMG process provides a general approach for computing the TVs used in
algebraic distance and LS interpolation.  For definiteness, we provide an overview
of one specific version of bootstrap AMG designed and analyzed in~\cite{BAMG2010}.
The algorithm begins with
relaxation applied to the homogenous system,
\begin{equation}\label{homog}
A_lx_l =0,
\end{equation}
on each grid, $l = 0,..., L-1$; assuming that a priori knowledge of the
algebraically smooth error is not available, these vectors are initialized randomly on the finest grid, whereas on all coarser grids
they are defined by restricting the existing test vectors computed on the previous finer grid.
Given interpolation constructed using CR, algebraic distances and the LS process,
the coarse-grid operators are computed, e.g., using the variational definition when $A$ is HPD.
Once an initial MG hierarchy has been computed,
the current sets of TVs are further enhanced on all grids using the
existing multigrid structure.  Specifically, the given hierarchy is used to formulate a multigrid eigensolver
which is then applied to an appropriately chosen generalized eigenproblem to compute additional test vectors.
{We note that the use of a coarse-grid eigensolver to compute test vectors has been studied theoretically in the context
of reduction-based AMG in~\cite{MacRed} and in the context of adaptive smoothed aggregation
in~\cite{GESA}.  The approach outlined here differs in its use of LS interpolation together with
a multilevel eigensolver as considered in~\cite{BAMG2010}.}
This overall process is then repeated with the current AMG method applied in addition to (or replacing)
relaxation as the solver for the homogenous systems in \eqref{homog}.
Figure \ref{fig:boot:setupcycle} provides an schematic outline of the bootstrap $V$- and $W$-cycle setup algorithms.  {In general, $V^m$ and $W^m$ denote setup algorithms that use $m$ iterations of the $V$- and $W$-cycles, respectively.}

The rationale behind the multilevel generalized eigensolver (MGE) is as follows.  Assume an initial multigrid hierarchy has been constructed.  Given the initial Galerkin operators $A_{0}, A_{1}, \ldots, A_{L}$ on each level
and the corresponding interpolation operators $P_{l+1}^{l}, l =
0,\ldots,L-1$, define the composite interpolation operators as
$P_{l} = P_{1}^{0}\cdot \ldots \cdot P_{l}^{l-1},\ l = 1, \ldots, L$.  Then,
for any given vector $x_{l} \in \mathbb{C}^{n_l}$ we have
$
  \innerprod[A_{l}]{x_{l}}{x_{l}} = \innerprod[A]{P_{l}x_{l}}{P_{l}x_{l}}.
$
Furthermore, defining $T_{l} = P_{l}^{H}P_{l}$ we obtain
\begin{equation*}
  \frac{\innerprod[A_{l}]{x_{l}}{x_{l}}}{\innerprod[T_{l}]{x_{l}}{x_{l}}} = \frac{\innerprod[A]{P_{l}x_{l}}{P_{l}x_{l}}}{\innerprod{P_{l}x_{l}}{P_{l}x_{l}}}.
\end{equation*}
This observation in turn implies that
on any level $l$, given a vector $x^{(l)} \in \mathbb{C}^{n_l}$ and $\lambda^{(l)}\in
  \mathbb{C}$ such that
$  A_l x^{(l)} = \lambda^{(l)}T_l x^{(l)},
$
we have Rayleigh quotient (RQ)
\begin{equation}\label{eq:evalapprox}
  \mbox{rq}(P_l x^{(l)}) :=
  \frac{\innerprod[A]{P_{l}x^{(l)}}{P_{l}x^{(l)}}}{\innerprod{P_{l}x^{(l)}}{P_{l}x^{(l)}}}
  = \lambda^{(l)}.
\end{equation}
This provides a relation among the eigenvectors and eigenvalues
of the operators in the multigrid hierarchy on all levels with
the eigenvectors and eigenvalues of the finest-grid operator $A$
that can be used to measure the accuracy of interpolation for smooth
modes.
Note that the eigenvalue approximations in~\eqref{eq:evalapprox} are continuously updated within the algorithm so that the overall approach resembles an inverse Rayleigh-Quotient
iteration found in eigenvalue computations (cf.~\cite{JWilkinson_1965}).
The overall MGE process is motivated by the following observation.

Consider the generalized eigenvalue problem
  \begin{equation}\label{eq:eva}
    A_{j}x_{j} = {\color{black}\lambda^{(j)}}T_{j}x_{j},
  \end{equation} on two subsequent grids, $j=l, l-1$, with
  $\ {\color{black}\lambda^{(l)}}$ denoting the RQ of a given approximation of an eigenvector on the coarser grid and $ {\color{black}\lambda^{(l-1)}}$ denoting the RQ of the vector obtained by applying relaxation to the system~\eqref{eq:eva} with this vector interpolated to the next finer grid as the initial guess.  Define the \emph{eigenvalue approximation measure} $\tau^{(l,l-1)}_{\lambda}$ by
    \begin{equation}\label{eq:EWapproxmeasure}
  {\color{black}   \tau^{(l,l-1)}_{\lambda} =
    \frac{|\lambda^{(l)}-\lambda^{(l-1)}|}{|\lambda^{(l-1)}|}}.
  \end{equation}
Then, at any stage of the MGE iteration, a large value of $\tau^{(l,l-1)}_{\lambda}$ indicates that the hierarchy should be recomputed to incorporate this relaxed eigenvector approximation.  Otherwise, relaxation has not significantly changed this vector and it is accurately represented by the existing interpolation operator $P_l^{l-1}$.  In this way, the MGE serves as an approach for efficiently identifying components that must  be interpolated accurately (e.g., the low modes of $A$), and with the use of the eigenvalue approximation measure a technique for determining if the current $P$ approximates them sufficiently well.
For additional details of the algorithm and its implementation we refer to the paper~\cite{BAMG2010}.

To illustrate the effect of the MGE in the bootstrap process, we provide results from~\cite{BAMG2010}, in which a $V^2$-cycle bootstrap setup using four pre- and post-smoothing steps is applied to compute the set of relaxed vectors $\mathcal{V}^r$ and set of bootstrap vectors $\mathcal{V}^e$ coming from the MGE process, with $k_r=|\mathcal{V}^r| = 8$ and $k_e= |\mathcal{V}^e|=8$, respectively.  The sets $\mathcal{V}^r$ and $\mathcal{V}^e$ are then combined to form the set of TVs $\mathcal{V}$ that is used to compute the least squares interpolation operator on each level.
In these tests, only relaxation is applied to the homogenous systems in both setup cycles to update the sets $ \mathcal{V}^r$.
The coarse grids and sparsity structure of interpolation are defined as in the previous tests on all levels
and the problem is coarsened to a coarsest level with $h=1/16$.
\begin{figure}
\begin{center}
     \tikzstyle{greenpoint}=[circle,inner sep=0pt,minimum size=2mm,draw=black!100,fill=black!30]
    \tikzstyle{whitepoint}=[rectangle,inner sep=0pt,minimum size=2mm,draw=black!100,fill=black!80]
    \tikzstyle{redpoint}=[diamond,inner sep=0pt,minimum size=2.55mm,draw=black!100,fill=black!80]
    \tikzstyle{blackpoint}=[circle,inner sep=0pt,minimum size=2mm,draw=black!100,fill=black!100]
    \tikzstyle{bluepoint}=[circle,inner sep=0pt,minimum size=2mm,draw=black!100,fill=black!0]
    \resizebox{\textwidth}{!}{\begin{tikzpicture}
      \draw [sharp corners] (0,0) node[blackpoint] {} --
      ++(300:1cm) node[blackpoint] {} --
      ++(300:1cm) node[blackpoint] {} --
      ++(300:1cm) node[blackpoint] {} --
      ++(300:1cm) node[redpoint] (L4) {} --
      ++(60:1cm) node[bluepoint] (L3) {} --
      ++(60:1cm) node[bluepoint] (L2) {} --
      ++(60:1cm) node[bluepoint] (L1) {} --
      ++(60:1cm) node[whitepoint] (L0) {} --
      ++(0:1cm) node[greenpoint] {} --
      ++(300:1cm) node[greenpoint] (end) {};
      \draw[dashed] (end) --  ++(300:1cm);

      \draw (L0) + (3cm,0) node[blackpoint,label=right:{\small Relax on $Av=0, v \in \mathcal{V}^r$, compute $P$}] {};
      \draw (L1) + (3.5cm,0) node[redpoint,label=right:{\small Compute $v$, s.t., $Av=\lambda Tv$, update $\mathcal{V}^e$}] {};
      \draw (L2) + (4.0cm,0) node[bluepoint,label=right:{\small Relax on $Av = \lambda T v, v \in \mathcal{V}^e$}] {};
      \draw (L3) + (4.5cm,0) node[greenpoint,label=right:{\small Relax on and / or solve using MG $Av=0, v \in \mathcal{V}^r$ and $Av = \lambda T v, v \in \mathcal{V}^e$, recompute $P$}] {};
      \draw (L4) + (5cm,0) node[whitepoint,label=right:{\small Test
        MG method, update $\mathcal{V}$}] {};

      \draw (L3) +(0,4cm) node (start) {};
      \draw[sharp corners] (start) ++(300:-3cm) node[blackpoint] {} --
      ++(300:1cm) node[blackpoint] {} --
      ++(300:1cm) node[blackpoint] {} --
      ++(300:1cm) node[redpoint] {} --
      ++(60:1cm) node[greenpoint] {} --
      ++(300:1cm) node[redpoint] {} --
      ++(60:1cm) node[bluepoint] {} --
      ++(60:1cm) node[greenpoint] {} --
      ++(300:1cm) node[greenpoint] {} --
      ++(300:1cm) node[redpoint] {} --
      ++(60:1cm) node[greenpoint] {} --
      ++(300:1cm) node[redpoint] {} --
      ++(60:1cm) node[bluepoint] {} --
      ++(60:1cm) node[greenpoint] {} --
      ++(300:1cm) node[greenpoint] {} --
      ++(300:1cm) node[redpoint] {} --
      ++(60:1cm) node[bluepoint] {} --
      ++(60:1cm) node[bluepoint] {} --
      ++(60:1cm) node[whitepoint] {} --
      ++(0:1cm) node[greenpoint] {} --
      ++(300:1cm) node[greenpoint] {} --
      ++(300:1cm) node[greenpoint] {} --
      ++(300:1cm) node[redpoint] {} --
      ++(60:1cm) node[greenpoint] (end) {};
      \draw[dashed] (end) --  ++(300:1cm);
    \end{tikzpicture}}
  \caption{Galerkin Bootstrap AMG $W$-cycle and $V$-cycle setup schemes.\label{fig:boot:setupcycle}}
\end{center}
\end{figure}

 As the results reported in Table~\ref{tab:nr:lsinterpMG84N01V2} show, using the MGE to enhance the
test vectors consistently improves the performance of the resulting solvers
when compared to the results reported in Table~\ref{tab:nr:lsinterp84N01},
especially for the LSR scheme.
\begin{table}[!ht]
  \begin{center}
    \begin{tabular}{|c|c|c|c|c|c|c|}
      \hline
      $h$ & $1/32$ & $1/64$ & $1/128$ & $1/256$ & $1/512$ \\
      \hline
     $\rho$ & $.041$ ($.038$) & $.062$ ($.041$) & $.075$
      ($.043$) & $.125$ ($.043$) & $.971$ ($.043$)\\
      \hline
    \end{tabular}
    \caption[$V^{2}$-BootAMG ($k_{r}=8, \eta=4$) for FD Laplace]{
   Asymptotic convergence rate estimates, $\rho$, of the LS (LSR) based $V(2,2)$ BAMG solvers applied to the FD Laplace problem.  The BAMG solvers are constructed using a $V^2$-cycle setup with $\eta=4$ and $k_{r}=k_e=8$.
\label{tab:nr:lsinterpMG84N01V2}}
  \end{center}
\end{table}
Here, the LSR approach yields an optimal method, whereas for the LS scheme, as the problem size increases the convergence rate of the solver also increases substantially -- from $\rho \approx.125$ to $\rho \approx.971.$

Next, we consider results of the $W$-cycle setup, illustrated at the top of Figure~\ref{fig:boot:setupcycle}.
In Table~\ref{tab:nr:lsinterpMG84N01W}, we present results
for this scheme again with $k_{r}=k_e=8$ and using only relaxation to solve the homogenous systems in both setup cycles to update the sets $ \mathcal{V}^r$.  We notice a marked improvement in the performance of the LS-based solver constructed using a $W$-cycle setup over that of the solver
constructed using a $V$-cycle setup and note that the LSR scheme again scales optimally for the given problem sizes.  The improved performance of the $W$-cycle setup indicates the effectiveness of using additional  coarse level processing in computing accurate TVs.

\begin{table}[ht!]
  \begin{center}
    \begin{tabular}{|c|c|c|c|c|c|c|}
      \hline
      $h$ & $1/32$ & $1/64$ & $1/128$ & $1/256$ & $1/512$ \\
      \hline
    $\rho$ & $.042$ ($.038$) & $.062$ ($.041$) & $.073$
      ($.043$) & $.130$ ($.043$) & $.161$ ($.044$)\\
      \hline
    \end{tabular}
    \caption[$W$-BootAMG ($k=j=8, \eta=4$) for FD Laplace]{Asymptotic convergence rate estimates, $\rho$, of the LS (LSR) based $V(2,2)$ BAMG solvers applied to the FD Laplace problem.  The BAMG solvers are constructed using a $W$-cycle setup (see~Figure~\ref{fig:boot:setupcycle}) with $ \eta=4$ and $k_{r}=k_e =8$.\label{tab:nr:lsinterpMG84N01W}}
  \end{center}
\end{table}

We conclude the numerical experiments with results of the BAMG setup algorithm applied to the
same bilinear finite element discretization of the 2D diffusion problem considered earlier
\begin{equation}
 - \nabla \cdot \mathcal{K}(x,y) \nabla u(x,y) = f(x,y),
 \end{equation}
 with homogenous Dirichlet boundary conditions.
However, now we assume that $\mathcal{K}(x,y)$ denotes the two-scale permeability:
$$\mathcal{K} = \bigg\{ \begin{array}{lr}
1 & (x,y) \in \Omega_1  \\
10^\nu  & (x,y) \in \Omega_2
\end{array},$$
with $\Omega_1$ given by the union of the black subdomains in Figure~\ref{fig:jumps}
and $\Omega_2 = \Omega \setminus \Omega_1$.
The problem is discretized on a uniform quadrilateral grid with various choices of the grid spacing
$h$ so that the grid aligns with the interfaces
for the given problem sizes (see the plot on the right of Figure~\ref{fig:jumps}).  We consider the values $\nu=-2,-8$
and various tilings, denoted by $\tau=1,4,8,16$ (see Figure~\ref{fig:jumps}) in testing the BAMG algorithm.

The BAMG method uses standard (full) coarsening with nearest-neighbor interpolation on all levels of the hierarchy, and the problem is coarsened until the coarsest systems has grid spacing $h=1/16$.  We note that on coarse levels, the grids and the interpolation patterns do not align with the jumps at the interfaces of $\Omega_1$ and $\Omega_2$.
We consider tests of the $W^2$-cycle setup with $\eta =4$ relaxation steps to compute the $k_{r}=k_e=8$ TVs.   For the results reported in the bottom two tables of Table~\ref{tab:aBAMG} we use the solver resulting from the first $W$-cycle setup in an adaptive step in which we solve the homogenous system $Ax=0$ using five $V(2,2)$ cycles on the finest level only.
Here, the initial guess is the TV from $ \mathcal{V}^r$ with the smallest energy, which is then replaced
by the resulting approximation as the new TV in $ \mathcal{V}^r$.
\begin{figure}[!t]
\centering{
\includegraphics[width=.65\textwidth]{jumps}   }
\caption{Tiling of the domain for the interface problem \eqref{eq:jumps}.  The figure on the left corresponds to $\tau = 4$ and the one on the right to $\tau =1$. \label{fig:jumps}}
\end{figure}

\begin{table}[!ht]
    \begin{center}\small
      \begin{tabular}{lll}
        \begin{tabular}{|c|c|c|c|c|c|c|c|}
          \hline
          $\tau$\\
          \hline
          $1$\\
          \hline
          $4$\\
          \hline
          $8$ \\
          \hline
          $16$\\
          \hline
           $h$\\
          \hline
        \end{tabular}

        \begin{tabular}{|c|c|c|c|c|c|c|c|}
          \hline
          \multicolumn{4}{|c|}{$\nu=-2$} \\
          \hline
          .25 & .33 & .57 & .88\\
          \hline
          .24 &  .32  & .55 & .63\\
          \hline
          .24 & .31 & .50 & .59 \\
          \hline
          * & .29 & .47& .57 \\
          \hline
          1/32 & 1/64 & 1/128 & 1/256 \\
          \hline
        \end{tabular}
        \begin{tabular}{|c|c|c|c|c|c|c|c|}
          \hline
          \hline
        \end{tabular}
        \begin{tabular}{|c|c|c|c|c|c|c|c|}
          \hline
          \multicolumn{4}{|c|}{$\nu=-8$} \\
          \hline
          .29 & .47 & .64 & .91\\
          \hline
          .26 & .57 & .79 & .88\\
          \hline
          .25 & .41 & .71 & .83\\
          \hline
          * &.37 & .68 & .70 \\
          \hline
          1/32 & 1/64 & 1/128 & 1/256\\
          \hline
        \end{tabular}
      \end{tabular}
    \end{center}

    \begin{center}\small

       \begin{tabular}{lll}
        \begin{tabular}{|c|c|c|c|c|c|c|c|}
          \hline
          $\tau$\\
          \hline
          $1$\\
          \hline
          $4$\\
          \hline
          $8$ \\
          \hline
          $16$\\
          \hline
           $N$\\
          \hline
        \end{tabular}

        \begin{tabular}{|c|c|c|c|c|c|c|c|}
          \hline
          \multicolumn{4}{|c|}{$\nu=-2$} \\
          \hline
          .20 & .21 & .23 & .23\\
          \hline
          .20 &  .22  & .22 & .22\\
          \hline
          .20 & .21 & .22 & .22 \\
          \hline
          * & .18 & .20& .20 \\
          \hline
          1/32 & 1/64 & 1/128 & 1/256 \\
          \hline
        \end{tabular}
        \begin{tabular}{|c|c|c|c|c|c|c|c|}
          \hline
          \hline
        \end{tabular}
        \begin{tabular}{|c|c|c|c|c|c|c|c|}
          \hline
          \multicolumn{4}{|c|}{$\nu=-8$} \\
          \hline
          .20 & .21 & .23 & .23\\
          \hline
          .21 & .22 & .22 & .22\\
          \hline
          .21 & .22 & .23 & .23\\
          \hline
          * &.18 & .20 & .20 \\
          \hline
          1/32 & 1/64 & 1/128 & 1/256\\
          \hline
        \end{tabular}
      \end{tabular}
 \caption[$V^{2}$-BootAMG ($k_{r}=8, \eta=4$) for FD Laplace]{
   Asymptotic convergence rate estimates, $\rho$, of the LSR based $V(1,1)$ BAMG solver applied to the FE discretization of the diffusion problem \eqref{eq:jumps}.  The solver is constructed using a $W^2$-cycle BAMG setup with $\eta=4$ and $k_{r}=k_e=8$.  For the results in the bottom tables, an additional intermediate adaptive step is applied to the TV with lowest energy on the finest level only before the second $W$-cycle setup.
\label{tab:aBAMG}}
    \end{center}
  \end{table}

 As the results illustrate, the $W^2$-cycle setup does not yield a scalable solver
 for this problem.\footnote{{In Table 6, * indicates that the fine grid
does not resolve the variations in the permeability field.}}  However, adding an adaptive step to the algorithm before the second $W$-cycle
 does result in an optimal order method.  We note that using additional $W$-cycles in the setup, i.e., a $W^m$-cycle
 with $m > 2$, yields similar results as obtained with this additional adaptive step, but that these
 additional cycles require significantly more work than does applying the adaptive step a few times only on the finest level, as implemented in the reported experiments.  Generally, using this simple modification reduces the number of bootstrap setup cycles needed to obtain an efficient solver by at least one and in most
 cases two or more for this problem, assuming the number of relaxation steps in the setup cycles is kept fixed.   {We note, in addition, that in~\cite{MBrezina_2005} adaptive AMG was  successfully applied to the same test problem and that an adaptive approach may be preferable for this problem since only a single smooth-error component is needed to computed effective AMG interpolation.}

The use of adaptive cycles in the bootstrap setup process has also been considered in various
 contexts, including the work in~\cite{BAMG_Schwinger} where a BAMG solver was designed for the non-Hermitian Dirac equation of lattice QCD and the work in~\cite{brannick_markov_2011,BAMG_Markov_2}, where a BAMG method for computing state vectors of Markov chains was developed.  These works have shown that
combining the least squares interpolation with the bootstrap setup for computing
the TVs leads to an efficient and integrated self learning process for defining algebraic multigrid interpolation.
While this research has resulted in an increased understanding of these techniques, a systematic study of how to combine them in an optimal and robust way is still needed.

\section{Future directions and open problems}\label{sec:FD}

Research in the design and analysis of the individual components of the BAMG solver
has led to efficient and robust multilevel algorithms for a variety of problems,
including indefinite Helmholtz and Schr\"odinger equations \cite{IL08, IL11}, non-Hermitian PDE systems such as the Dirac equation in lattice QCD~\cite{BAMG_Schwinger}, anisotropic diffusion problems with highly oscillatory and strongly anisotropic
diffusion coefficient~\cite{BAMG_Aniso}, and the Graph Laplacian problem~\cite{LAMG_Report}.  More generally, this work has led to an increased understanding of the techniques of
compatible relaxation, algebraic distances, least squares interpolation,
and the bootstrap process for algebraically computing representatives of
smooth error, the test vectors.   These works have demonstrated the potential
of the BAMG methodology to deliver an efficient and robust algebraic multilevel solver
for a large class of sparse matrix equations.

The basic components of the BAMG process have been studied and are (mostly) understood
for HPD problems .
The main issue to consider in developing BAMG as a fast solver for general HPD problems concerns the integration of the BAMG components in order to obtain an overall algorithm for use in practice.  Specifically, a systematic study of the BAMG processes for a carefully chosen and diverse set of test problems (that includes
systems of PDEs) is needed.  As noted in the previous sections,
the most pressing research directions as the authors currently see them are
\begin{itemize}
\item An efficient  mechanism for  updating coarse variable set $C$ based on a local  slowness of CR relaxation, see Section \ref{sec:CR};
\item A local procedure for finding weights $\omega_\kappa$ in (\ref{eq:LSfuncrowi}) -- to prevent an emphasize on a globally smooth test vector, that will otherwise be  heavily weighted everywhere,  in the regions  where it has large residuals; such approach will be a generalization of the LSR scheme which is now used concurrently with  global weights, see Section \ref{sec:LS};
\item A balance between work and accuracy when computing  TVs, see Section \ref{sec:BS};
\item A procedure  for incorporating all available information about a given problem, in particular, its  geometric information;
\item An integration of all BAMG components to obtain an easy to tune and robust multigrid solver.
\end{itemize}

Although the basic components of BAMG can in principle also
be applied to non-Hermitian problems, there are still several basic research questions in this
setting that require further investigation.  One main unresolved issue that arises in a Galerkin-based algebraic
multigrid concerns characterizing the algebraically smooth errors that need to be approximated
by restriction and interpolation.  The typical approach has been to compute approximations to singular values with small singular vectors, {e.g., in~\cite{GESA} such a strategy was developed for conviction diffusion problems.}  These approaches are based on the observation that the smallest
singular value is less than the smallest in magnitude eigenvalue.
However, for many problems, the smallest eigenvalues of a non-Hermitian matrix are much smaller than the second smallest singular value.  Thus, in general the near kernel of the system matrix can not be characterized
by the singular vectors (or eigenvectors) with small singular values, unless the matrix is normal.

The basic approach in BAMG is to iterate on the homogenous system using the AMG smoother and
the emerging multilevel structure in order to expose the smooth errors that need
to be treated on coarse levels.  Given an appropriate relaxation scheme,
e.g., the Kaczmarz iteration, which provides a general approach for
non-Hermitian problems, the subspace of smooth error is in
the most general sense related to the field of values
\begin{equation*}
  \mathcal F =\left\{ \frac{\langle Ax,x\rangle}{\langle x,x \rangle} \bigg|  x \in \mathbb{C}^n, x \neq 0\right\}.
\end{equation*}
We note that the field of values may contain the origin even if all eigenvalues of the matrix
have a positive real part.  Thus, it is possible that the Galerkin coarse
level operator can have a zero eigenvalue, even when the fine-level matrix
does not.
We have observed in practice that in such
cases, the resulting Galerkin method can
diverge.
 Here, an underestimate of the small singular values
or eigenvalues on some grid can lead to an over-correction
from the coarse grid and divergence of the
resulting AMG solver.
Thus, although the BAMG process provides techniques for constructing
restriction $R$ and interpolation $P$, each in their
own setup scheme, these processes must be carefully tied together locally
in order to obtain a suitable Galerkin coarse-grid matrix.

The challenges presented by indefinite problems, e.g., the Helmholtz and Schr\"odinger operators, are again related to the local character of their low eigenspace, which consists of many eigenfunctions that are locally different from one another.  Generally, we refer to such problems as non-elliptic systems.
The key issues is that such eigenfunctions cannot be accurately interpolated by a single prolongation operator (\ref{eq:LSfuncrowi}) with a small caliber, starting with some sufficiently fine scales. There and on all coarser scales several prolongation operators have to be employed, resulting in multiple-Galerkin approach with many coarse-grid operators. To guarantee optimal computational costs, all prolongation operators have to be locally orthogonal, and  achieving that  generates  additional computational difficulties. To this date, the BAMG-like approach has been implemented only for the one-dimensional Helmholtz \cite{IL08} and Schr\"odinger \cite{IL11} problems, which require only {\sl two} prolongation operators (instead of one).   In higher dimensions, the number of prolongation operators and coarse spaces increases, but the ideas carry over in a similar way from 1d.
Two important future research topics are thus to continue to investigate developing BAMG solvers for these and other non-elliptic problems and to determine how far these techniques can be extended as fast solvers for general sparse matrix equations.

\bibliographystyle{plain}
\bibliography{bamg_overview}

\end{document}

config: arxiv


text:
\begin{document}

\title{Reasonable non--Radon--Nikodym idealss}
\author{Vladimir Kanovei\thanks
{IITP, Moscow.
Partially supported by RFFI grants 06-01-00608 and 07-01-00445.}
\and
Vassily Lyubetsky\thanks
{IITP, Moscow.
Partially supported by RFFI grant 07-01-00445. }
}
\date{\today}
\maketitle

\begin{abstract}
For any abelian Polish \dd\sg compact group $\dH$
there exist an $\fs$ ideal $\cZ\sq\cP(\dN)$ and a Borel
\dd\cZ approximate homomorphism $f:\dH\to\dH^\dN$ which
is not \dd\cZ approximable by a continuous true
homomorphism $g:\dH\to\dH^\dN$.
\end{abstract}

\subsubsection*{Introduction}

Let $G,H$ be abelian Polish groups,
and $\cZ$ be an ideal over a countable set $A$.
We consider $H^A$ as a product group.
For $s,t\in H^A$ put
\dm
\raz{s}{t}=\ens{a\in A}{s(a)\ne t(a)}\,.
\dm
Suppose that $\cZ$ is an ideal over $A$.
A map $f:G\to H^A$ is a {\it\dd\cZ approximate \hm\/} iff
$\raz{f(x)+f(y)}{f(x+y)}\in\cZ$ for all $x,y\in G$.
Thus it is required that the set of all $a\in A$ such that
$f_a(x)+f_a(y)\ne f_a(x+y)$ belongs to $\cZ$.
Here $f_a:G\to H$ is the \dd ath co-ordinate map of the map
$f:G\to H^A$.

And $\cZ$ is a {\it Radon--Nikodym\/} ideal
(for this pair of groups)
iff for any measurable \dd\cZ approximate \hm\ $f:G\to H^\dN$
there is a continuous exact \hm\ $g:G\to H^\dN$ which
\dd\cZ approximates $f$ in the sense that
$\raz{f(x)}{g(x)}\in\cZ$ for all $x\in G$.
Here the measurability condition
can be understood as Baire measurability,
or, if $G$ is equipped with a \dd\sg additive Borel measure, as
measurability with respect to that measure.

The idea of this (somewhat loose) concept is quite clear: the
Radon--Nikodym ideals are those which allow us to approximate
non-exact \hm s by true ones.
This type of problems appears in different domains of mathematics.
Closer to the context of this note,
Velickovic \cite{v} proved that any Baire-measurable
\kon approximate Boolean-algebra automorphism $f$ of $\cP(\dN)$
(so that the symmetric differences between $f(x)\cup f(y)$
and $f(x\cup y)$ and between $f(\dN\bez x)$ and $\dN\bez f(x)$
are finite for all $x,y\sq\dN$) is \kon approximable by a true
automorphism $g$ induced by a bijection betveen two cofinite
subsets of $\dN$.
Kanovei and Reeken proved that any Baire measurable
\dd\dQ approximate homomorphism
$f:\dR\to\dR$ is \dd\dQ approximable by a homomorphism of the
form $f(x)=cx$, $c$ being a real constant.
See also some results in \cite{f,kr1,kr2}.

The term ``Radon--Nikodym ideal'' was introduced by
Farah~\cite{f,aq} in the context of Baire measurable
Boolean algebra \hm s of $\cP(\dN)$.
Many known Borel ideals were demonstrated to be Radon--Nikodym,
see \cite{f,aq,kr1,kr2}.
Suitable counterexamples, again in the context of Boolean
algebra \hm s, were defined by Farah on the base of so called
pathological submeasures.
A different and, perhaps, more transparent counterexample,
related to \hm s $\dT\to\dT^\dN$
(where $\dT=\dR/\dN$),
is defined in \cite{kr2} as a modification of an ideal
introduced in \cite{so}.
The next theorem generalizes this result.

\bte
\label m
Suppose that\/ $\dH$ is an uncountable abelian Polish group.
Then there is an analytic ideal\/ $\cZ$ over\/ $\dN$ that is
not a Radon -- Nikodym ideal for maps\/ $\dH\to\dH^\dN$ in the
sense that there is a Borel and\/ \dd\cZ approximate \hm\/
$f:\dH\to\dH^\dN$ not\/ \dd\cZ approximable by a
continuous \hm\/ $g:\dH\to\dH^\dN$.
If moreover\/ $\dH$ is\/ \dd\sg compact then\/ $\cZ$ can be
chosen to be\/ $\fs$.
\ete

Note that the theorem will not become stronger if we require
$g$ to be only Baire-measurable, or just measurable with respect
to a certain Borel measure on $\dH$ --- because by the Pettis
theorem any such a measurable group \hm\ must be continuous.

The remainder of the note contains the proof of Theorem~\ref m.
It would be interesting to prove the theorem for non-abelian
Polish groups.
(The assumption that $\dH$ is abelian is used in the proof
of Lemma~\ref{l1}.)
And it will be interesting to find non--Radon--Nikodym ideals
for \hm s $G\to H^\dN$ in the case when the Polish groups
$G$ and $H$ are not necessarily equal.

\subsubsection{Countable subgroup}
\label{pr}

Let us fix a group $\dH$ as in the theorem, that is, an
uncountable abelian Polish group.
By $\bbo$ we denote the neutral element,
by $\opl$ the group operation,
by $d$ a compatible complete separable distance
(and we do not assume it to be invariant).

The first step is to choose a certain countable subgroup
$D\sq\dH$ of ``rational elements''.

It is quite clear that there exists
a countable dense subgroup $D\sq\dH$
satisfying the following requirement of
elementary equivalence type.

\bit
\item[$(\ast)$]
Suppose that $n\ge 1$, $c_1,\dots,c_n\in D$,
$\ve$ is a positive rational, $U_i=\ens{x\in\dH}{d(x,c_i)\le\ve}$,
and $P(x_1,\dots,x_n)$ is a finite system of linear equations
with integer coefficients, unknowns $x_1,\dots,x_n$, and constants
in $D$, of the form:
\dm
b_1x_1\opl\dots\opl b_nx_n=r\,,\quad \text{where}\quad
b_i\in\dZ\;\text{ and }\;r\in D\,.
\dm
Suppose also that this system $P$ has a solution
$\ang{x_1,\dots,x_n}$ in $\dH$ such that $x_i\in U_i$ for all $i$.
Then $P$ has a solution in $D$ as well.
(That is, all $x_i$ belong to $D\cap U_i$.)
\eit

Let us fix such a subgroup $D$.

\subsubsection{The index set}
\label{is}

Let

{\it rational ball\/} mean any subset of $\dH$ of the form
$\ens{x\in\dH}{d(c,x)<\ve}$, where $c\in D$ (the center),
and $\ve$ is a positive rational number.

\bdf
\label{3a}
Let $A$, the index set, consist of all objects $a$ of the
following kind.
Each $a\in A$ consists of:
\bit
\itsep
\item[$-$]
an open non-empty set $U^a\sneq\dH$,

\item[$-$]
a partition $U^a=U^a_1\cup\dots\cup U^a_n$ of $U^a$ onto
a finite number $n=n^a$ of pairwise disjoint non-empty

rational balls $U^a_i\sq\dH$, and

\item[$-$]
a set of points $r^a_i\in U^a_i\cap D$
such that, for all $i,\,j=1,2,\dots,n$:

\ben
\aenu
\itla{31}
{\it either\/} $r^a_i\opl r^a_j$ is $r^a_{k}$ for some $k,$
and $(U^a_i\opl U^a_j)\cap U^a\sq U^a_k$,

\itla{32}
{\it or\/} $(U^a_i\opl U^a_j)\cap U^a=\pu$.\qed

\vyk{\item[$-$]
{\it or\/} $i=j$ and $s_i=0$.  }
\een
\eit

\eDf

Under the conditions of Definition~\ref{3a},
if $\bbo\in U^a_i$ then $s_i=\bbo$: for take $j=i$.

\ble
\label{ai}
$A$ is an infinite (countable) set.
\ele
\bpf
For any $\ve>0$ there is $a\in A$ such that
$U^a$ a set of diameter $\le\ve$:
just take $n^a=1$, $r^a_1=\bbo$, and let
$U^a=U^a_1$ be the \dd{\frac\ve2}nbhd of $\bbo$ in $\dH$.
\epf

The next lemma will be used below.

\ble
\label{l*}
If\/ $y_1,\dots,y_n\in\dH$ are pairwise distinct then there
exists\/ $a\in A$ such that\/ $n^a=n$ and\/
$y_i\in U^a_i$ for all\/ $i=1,\dots,n$.
\ele
\bpf
As the operation is continuous, we can pick pairwise disjoint
rational balls $B_1,\dots,B_n$ such that $y_i\in B_i$ for all
$i$ and the following holds:
If $1\le i,j\le n$ then either there exists $k$ such that
$(B_i\opl B_j)\cap B\sq B_k$, where $B=B_1\cup\dots\cup B_n$,
or just $(B_i\opl B_j)\cap B=\pu$.
Put $U^a_i=B_i$.

To obtain a system of points $r^a_i$ required,

\vyk{
Consider the finite set
$$
X=\ens{y_i}{1\le i\le n}
\cup\ens{y_i\opl y_j}{1\le i,j\le n}\,.
$$
Let $\ve$ be the smallest positive real of the form
$d(x,y)$, where $x,y$ are two distinct elements in $X$.
Put $n^a=n$ and let each $U^a_i$ ($1\le i\le n$)
be a rational ball in $\dH$ of diameter $\frac\ve4$
whose center $c_i\in D$ is closer to $x_i$ than $\frac\ve{16}$.
To define elements $r^a_i$,
}

let $P(x_1,\dots,x_n)$ be the system of all equations of
the form $x_i+x_j=x_k$ with unknowns $x_i,x_j,x_k$, where
$1\le i,j,k\le n$ and in reality $y_i+y_j=y_k$.
It follows from the choice of $D$ that this system has a
solution $\ang{r_1,\dots,r_n}$ such that
$r_i\in U^a_i\cap D$ for all $i$.
In other words we have: $r_i+r_j=r_k$ whenever $y_i+y_j=y_k$.
Let $r^a_i=r_i$.
This ends the definition of $a\in A$ as required.

(An extra care to guarantee that $U^a=\bigcup_{1\le i\le n}U^a_i$
is a proper subset of $\dH$ is left to the reader.)

\epf

\subsubsection{The ideal}
\label{z}

Let $\cZ$ be the set of all sets $X\sq A$ such that
there is a finite set $u\sq\dH$ satisfying the
following:
for any $a\in X$ we have $u\not\sq U^a$.

The idea of this ideal goes back to Solecki~\cite{so},
where a certain ideal over the set
$\Om$ of all clopen sets $U\sq2^\dN$ of measure $\frac12$
(also a countable set) is considered.
In our case the index set $A$ is somewhat more complicated.

\ble
\label p
$\cZ$ is an ideal containing all finite sets\/ $X\sq A$,
but\/ $A\nin\cZ$.
\ele
\bpf
If $a\in A$ then the singleton $\ans a$ belongs to $\cZ$.
Indeed by definition $U^a$ is a non-empty subset of $\dH$.
Therefore there is a point $x\in \dH\bez U^a$.
Then $u=\ans x$ witnesses $A\in \cZ$.

To see that $\cZ$ is closed under finite unions, suppose
that finite sets $u,v\sq\dH$ witness that resp.\ $X,Y$
belong to $\cZ$.
Then $w=u\cup v$ obviously witnesses that $Z=X\cup Y\in\cZ$.
Finally by Lemma~\ref{l*} for any finite
$u=\ans{x_1,...,x_n}\sq\dH$ there is an element $a\in A$
such that $u\sq U^a$.
This implies that $A$ itself does not belong to $\cZ$.
\epf

\bpro
\label{p1}
$\cZ$ is an analytic ideal.
If\/ $\dH$ is\/ \dd\sg compact then\/ $\cZ$ is\/ $\fs$.
\epro
\bpf
We claim that
$X\in\cZ$ iff there are a natural $n$ and a partition
$X=\bigcup_{1\le k\le n}X_k$ such that for any $k$ the set
$X_k\sq A$ satisfies $\bigcup_{a\in X_k}U^a\ne\dH$.
Indeed suppose that $X\in\cZ$ and this is witnessed by a
finite set $u=\ans{x_1,\dots,x_n}\sq \dH$, that is,
$u\not\sq U^a$ for all $a\in X$.
It follows that $X=\bigcup_{1\le k\le n}X_k$, where
$X_k=\ens{a\in X}{x_k\nin U^a}$.
Clearly $x_k\nin \bigcup_{a\in X_k}U^a$.
To prove the converse suppose that
$X=\bigcup_{1\le k\le n}X_k\sq A$ and
$\bigcup_{a\in X_k}U^a\ne\dH$ for all $k$.
Let us pick arbitrary points
$x_k\in\dH\bez\bigcup_{a\in X_k}U^a$ for all $k$.
Then $u=\ans{x_1,\dots,x_n}$ witnesses $X\in\cZ$, as required.

It easily follows that $\cZ$ is analytic.

Now suppose that $\dH=\bigcup_{\ell\in\dN} H_\ell$,
where all sets $H_\ell$ are compact.
Then the inequality $\bigcup_{a\in X_k}U^a\ne\dH$ is equivalent
to $\sus\ell\:(H_\ell\not\sq\bigcup_{a\in X_k}U^a)$.
And by the compactness, the non-inclusion
$H_\ell\not\sq\bigcup_{a\in X_k}U^a$ is equivalent to
the following statement: $H_\ell\not\sq\bigcup_{a\in X'}U^a$
for every finite $X'\sq X_k$.

Fix an enumeration $A=\sis{a_n}{n\in\dN}$.
Put $A\res m=\ens{a_j}{j<m}$.
Using K\"onig's lemma, we conclude that
$X\in\cZ$ iff there exist natural $\ell,n$ such that
for any $m$ there exists a partition
$X\cap({A\res m})=\bigcup_{k<n}X_k$,
where for every $k$ we have $H_\ell\not\sq\bigcup_{a\in X_k}U^a$.
And this is a $\fs$ definition for $\cZ$.

\epf

\subsubsection
{The main result}

Here we prove Theorem~\ref{m}.
Define a Borel map $f:\dH\to\dH^A$ as follows.
Suppose that $x\in\dH$ and $a\in A$, $n^a=n$.
If\/ $x\in U^a_i$, $1\le i\le n$, then put $f_a(x)=x\omi r^a_i$.
($\omi$ in the sense of the group $\dH$.)
If\/ $x\nin U^a$ then put simply $f_a(x)=\bbo$.

Finally define $f(x)=\sis{f_a(x)}{a\in A}$.
Clearly $f$ is a Borel map.

The maps $f_a$ do not look like \hm s $\dH\to\dH$.
Nevertheless their combination surprisingly turns out to be an
approximate \hm!

\ble
\label{l1}
$f:\dH\to\dH^A$ is a Borel and\/ \ddz approximate \hm.
\ele
\bpf
Let $x,\,y\in\dH$ and $z=x\opl y$.
Prove that the set
$$
C_{xy}=\ans{a:f_a(x)\opl f_a(y)\ne f_a(z)}
$$
belongs to $\cZ$.
We assert that this is witnessed by the set $u=\ans{x,y,z}$,
that is, if $a\in C_{xy}$ then at least one of the
points $x,y,z$ is not a point in $U^a$.
Or, equivalently, if $a\in A$ and $x,y,z$ belong to $U^a$
then $f_a(x)\opl f_a(y)=f_a(z)$.

To prove this fact suppose that $a\in A$ and $x,y,z\in U^a$.
By definition, $U^a=U^a_1\cup\dots\cup U^a_n$, where $n=n^a$
and $U^a_i$ are disjoint
rational balls in $\dH$.
We have $x\in U^a_i$, $y\in U^a_j$, $z\in U^a_k$, where
$1\le i,j,k\le n$.
Then by definition
$$
f_a(x)=x\omi r^a_i\,,\quad f_a(y)=y\omi r^a_j\,,\quad
f_a(z)=z\omi r^a_k\,.
$$
Therefore $f_a(x)\opl f_a(y)=x\opl y\omi(s_i\opl s_j)$.
(Here we clearly use the assumption that the group is
abelian.)
We assert that $r^a_i\opl r^a_j=r^a_k$ --- then obviously
$f_a(x)\opl f_a(y)=f_a(z)$ by the above, and we are done.

Note that $z=x\opl y\in U^a$, hence
$(U^a_i\opl U^a_j)\cap U^a\ne\pu$.
We conclude that \ref{32} of Definition \ref{3a} fails.
Therefore \ref{31} holds,
$r^a_i\opl r^a_j=r^a_{k'}$ for
some $k'$ and $(U^a_i\opl U^a_j)\cap U^a\sq U^a_{k'}$.
But the set $(U^a_i\opl U^a_j)\cap U^a$ obviously contains $z$, and
$z\in U^a_k$.
It follows that $k'=k$, $r^a_{k'}=r^a_k$, $r^a_i\opl r^a_j=r^a_k$,
as required.
\epF{Lemma}

\ble
\label{l2}
The approximate \hm\/
$f$ is not\/ \ddz approximable by a continuous \hm\/
$g:\dH\to\dH^A$.
\ele
\bpf
Assume towards the contrary  that $g:\dH\to\dH^A$ is
a continuous \hm\ which \ddz approximates $f$.
Thus if $x\in\dH$ then the set
$\Da_x=\ans{a:f_a(x)\ne g_a(x)}$ belongs to $\cZ,$
where, as usual, $g_a(x)=g(x)(a)$.
Note that all of these projection maps $g_a:\dH\to\dH$ are
continuous group \hm s since such is $g$ itself.

Thus if $x\in\dH$ then $\Da_x\in\dZ$, and hence there is a
finite set $u_x\sq D$ satisfying the following:
if $a\in A$ and $u_x\sq U^a$ then $a\nin \Da_x$, that is,
$f_a(x)=g_a(x)$.
Put
\dm
X_u=\ens{x\in\dH}
{\kaz a\in A\:(u\sq U^a\imp f_a(x)=g_a(x))}
\dm
for every finite $u\sq D$.
These sets are Borel since so are maps $f,g$
(and $g$ even continuous).
Moreover $\dH=\bigcup_{u\sq D\;\text{finite}}X_u$ since every
$x\in\dH$ belongs to $X_{u_x}$.
Thus at least one of the sets $X_u$ is not meager, therefore,
is comeager on a certain rational ball $B\sq\dH$.
Fix $u$ and $B$.
By definition for comeager-many $x\in B$ and all
$a\in A$ satisfying $u\sq U^a$ we have $f_a(x)=g_a(x)$.

Arguing as in the proof of Lemma~\ref{l*}, we obtain an
element $a\in A$ satisfying the following properties:
$u\sq U^a$, $U^a\cap B\ne\pu$, but the set $B\bez U^a$
is non-empty and moreover is not dense in $B$.
Fix such $a$.
Thus there exists a non-empty rational ball $B'\sq B$ that
does not intersect $U^a$.
By definition $f_a(x)=\bbo$ for all $x\in B'$, and hence
$g_a(x)=\bbo$ for comeager-many $x\in B'$ by the choice of $B$.
We conclude that $g_a(x)=\bbo$ for all $x\in B$ in general,
because $g$ is continuous.

Now, let $n^a=n$.
Then $U^a=U^a_1\cup\dots\cup U^a_n$.
Recall that the intersection $B\cap U^a$ of two open sets is
non-empty by the choice of $a$.
It follows that there exists an index $i$, $1\le i\le n$, and
a non-empty rational ball $B''\sq B\cap U^a_i$.
Then by definition $f_a(x)=x\omi r$ for all $x\in B''$,
where $r=r^a_i$.
Therefore $g_a(x)= x\omi r$ for comeager-many $x\in B''$,
and then $g_a(x)= x\omi r$ for all $x\in B''$ since $g$ is
continuous.

To conclude, $g_a$, a continuous group \hm, is constant $\bbo$
on a non-empty open set $B'$, and is bijective on another
non-empty open set $B''$.
But this cannot be the case.
\epF{Lemma}

Lemmas \ref{l1} and \ref{l2} complete the proof of
Theorem~\ref{m}.

\small


TITLE: Relation between subgroups of a cyclic group of order $p^n.$
QUESTION [1 upvotes]: Let $G$ be a finite cyclic group of order $p^n$, where $p$ is a prime and $n \geq 0 .$ If $H \& K$ are subgroups of $G$, then show that either $H \subset K$ or $K \subset H.$
My attempt
Let $G$ be of order $p^n$ and let $a$ be one of it's generator i.e $G=\langle a\rangle.$ Since $G$ is cyclic, so are $H$ and $K$. Let $|H|=p^h$ and $|K|=p^k$ (We know that order should be of this form due to Langrange's theorem). Without loss of generality assume $h \le k$. We wish to show that $H \subset K$.
When $h=k$ this is true since there is only one subgroup of order $p^h$.
Can you complete the proof when $h < k$. If you have another approach, please post.

REPLY [1 votes]: The lattice of subgroups of a cyclic group of order $n$ is isomorphic to the lattice of divisors of $n$. When $n$ is a power of a prime, this lattice is a chain and "divides" then becomes the same as "less than or equal to".


TITLE: Riemann surface for Gaussian integral $\int \exp(z^2)dz$
QUESTION [2 upvotes]: I am having difficulty understanding the Riemann surface for the inverse of the Gaussian integral $\int_0^z e^{\zeta^2}d\zeta$ in the complex plane.
Let $g(z)=\int_0^z e^{\zeta^2}d\zeta$, so $g'$ is never zero, and so is locally 1-1. Also, it seems that $g:\mathbb{C}\rightarrow\mathbb{C}$ must be surjective, by the following argument. Since $g$ is odd, if its range omits $w$ it must also omit $-w$, so by the little Picard theorem it can omit only 0, but it doesn't omit 0.
So $g$ is a nonconstant locally 1-1 entire function mapping $\mathbb{C}$ onto $\mathbb{C}$. Now, it seems not hard to prove that the Riemann surface of the inverse of $g$ has a single sheet and covers $\mathbb{C}$, and so gives a global inverse. (I will give more details below.) But then $g$ is a holomorphic homeomorphism, and the only such maps from $\mathbb{C}$ to itself are the affine maps $z\mapsto az+b$. (Note that $\lim_{z\rightarrow\infty}g(z)=\infty$ by a compactness argument, so setting $g(\infty)=\infty$ gives a meromorphic function by the Riemann removable singularities theorem.)
Here is the argument about the existence of the global inverse. Points on the Riemann surface correspond to function elements (power series), so say $g(z_0)=w_0$ and say $h_1$ and $h_2$ are two local inverses to $g$ both centered at $w_0$ and converging in a small disk about it. Since $g(h_1(w))=w=g(h_2(w))$ for all $w$ in the disk and $g$ is locally 1-1, it follows that $h_1\equiv h_2$ in the disk and so define the same function element. So every $w$ has at most one point above it in the Riemann surface. However, since $g$ is onto and an open map, there is a holomorphic local inverse around each $w$, and so every point has a point above it in the surface. (If you prefer germs to function elements, you can make pretty much the same argument.)
Where is the flaw? Can someone provide a reference discussing this surface?
Edit [Slaps forehead] I see the problem. Let $R$ be the Riemann surface. It is true that as a covering of the $w$-plane (domain of $g$), $R$ is single-sheeted. But as a surface mapping onto the $z$-plane (range of $g$), it is not, although it's a locally 1-1 holomorphic map onto the $z$-plane. So it can't be a covering of the $z$-plane.
It's an interesting surface, and if someone has a reference giving more info about it, I'd appreciate it.

REPLY [1 votes]: It is a very nice surface but the map is not 1-1. At infinity the map has an essential singularity so every value except possibly one is taken infinitely many times (Picard's Great Theorem). The symmetry argument you give shows every value is taken.
The Riemann surface argument (for the 1-1 correspondance) only works if any bounded path in the image, starting say at the origin, $w=\gamma(t)$, $0\leq t \leq 1$, $\gamma(0)=0$  may be lifted to a bounded path in the source. The lifted path is the solution of the ode:
  $$ \frac{dz}{dt}=\frac{d\gamma(t)}{dt} \frac{1}{g'(z(t))}, \ \ z(0)=0.$$
Any simply connected neighborhood $U$ of the origin in the image then lifts to a simply connected nghb $V$ of the origin in the source and $g:V\rightarrow U$ is indeed 1-1. The problem here is that the path from
e.g. $0$ to $i\sqrt{\pi}/2$ lifts to a path going to infinity and there is no unique lift beyond, so a 1-1 correspondence  is no longer guaranteed (and fails here).


TITLE: Show that for each integer $n ≥ 2$, $F_n$ is a finite index subgroup of $F_2$, where $F_n$ is the free group on $n$ generators.
QUESTION [4 upvotes]: I'm trying to answer a problem which asks me to show (using covering spaces) that for each integer $n ≥ 2$, $F_n$ is a finite index
subgroup of $F_2$, where $F_n$ is the free group on $n$ generators.
I can see that if I can find some covering space with fundamental group $F_n$ of the wedge $S^1 \vee S^1$ then this will give me the existence of the subgroup, and I know that this will be of finite index if this covering map sends finitely many points to the vertex joining the two circles in $S^1 \vee S^1$.
The problem I'm having is I'm struggling to find such a covering space - the only space I know of with fundamental group $F_n$ is the wedge of $n$ circles, but I can't see any covering map from that to $S^1 \vee S^1$ so I assume I must need another.
If anyone could suggest such a covering space (and if so how I can show it has fundamental group $F_n$) I'd really appreciate the help!

REPLY [1 votes]: You might be interested in reading the first few pages of Hatcher's section on covering spaces (starts page 56, 65 in pdf). I'll be taking the pictures from there.
let the wedge $S^1 \vee S^1$ be labeled as follows:

And consider the space of $n$ circles connected together $\bigcirc\!\!\bigcirc\cdots \bigcirc\!\!\bigcirc$ as in the post linked in the comments. This will be our covering space.
We're going to label this space as follows:

(This is a particular case, but you get the idea)
Hopefully it's clear that if I send each labeled edges (preserving orientation) to their corresponding edges in $S^1 \vee S^1$, we get a covering map. It remains to show that this space has fundamental group of $F_n$.
To do this, imagine taking the bottom left $b$ edge and gradually contracting it until until it disappears (the two endpoint are now one). The upper left $b$ is now a circle, so we have two circles (left $a$, upper left $b$) wedged with a figure eight. This process doesn't change the fundamental group (fill in the details).
Do the same thing with the bottom right $a$. We now have a $4$-rose! Because the process doesn't change the fundamental group, this means the space we started with has fundamental group $F_4$.
It seems like you know how to complete the rest of the proof about finite index subgroups.