The idea of using polynomial methods to improve simple smoother iterations within a multigrid method for a symmetric positive definite (SPD) system is revisited. When the single-step smoother itself corresponds to an SPD operator, there is in particular a very simple iteration, a close cousin of the Chebyshev semi-iterative method, based on the Chebyshev polynomials of the fourth instead of first kind, that optimizes a two-level bound going back to Hackbusch. A full V-cycle bound for general polynomial smoothers is derived using the V-cycle theory of McCormick. The fourth-kind Chebyshev iteration is quasi-optimal for the V-cycle bound. The optimal polynomials for the V-cycle bound can be found numerically, achieving an about 18% lower error contraction factor bound than the fourth-kind Chebyshev iteration, asymptotically as the number of smoothing steps $k \to \infty$. Implementation of the optimized iteration is discussed, and the performance of the polynomial smoothers are illustrated with a simple numerical example.

§ INTRODUCTION

Given a multigrid method for solving a linear system that uses simple linear iterations for smoothers, a natural question is whether Krylov (i.e., polynomial) acceleration techniques can be applied to improve the smoother iterations. General-purpose Krylov methods don't necessarily have good smoothing properties, since they are typically designed to optimize some norm of the total error and not just its “un-smooth” components. In particular, a method like conjugate gradients can dramatically worsen the performance of a smoother iteration on the “high-frequency” error components in return for a marginal improvement in the performance on the “smooth” error components, which should of course be left for the coarse grid correction to handle. There are, however, polynomial methods specifically designed to be used as smoothers. One popular strategy for symmetric positive definite (SPD) systems, suggested by Adams, Brezina, Hu, and Tamaro <cit.>, is to use the standard Chebyshev semi-iteration, but adjusted to target only an upper portion $[\lambda^*, \lambda_{\max}]$ of the range of eigenvalues $[\lambda_{\min}, \lambda_{\max}]$, with $\lambda^*$ an ad-hoc parameter. While this method works well in practice, and it is reported that the performance is not very sensitive to the particular choice of $\lambda^*$, the presence of the ad-hoc parameter is still unaesthetic.
We will see that we can design methods that avoid such ad-hoc parameters, by picking polynomials that optimize multigrid convergence bounds. We will restrict our attention to SPD systems and single-step SPD smoothers (since polynomial methods tend to be most effective in the SPD setting), and examine both a two-level bound for general polynomial smoothers going back to Hackbusch <cit.>, as well as develop a full V-cycle convergence bound for general polynomial smoothers starting from the V-cycle theory developed by McCormick <cit.>. We will see that a very simple Chebyshev-type iteration, which avoids the need for any ad-hoc parameter, optimizes the two-level bound and is quasi-optimal for the V-cycle bound. We will also see that a small modification to the iteration, though somewhat complicated to implement, optimizes the V-cycle bound. The effectiveness of these smoother acceleration techniques are demonstrated with a numerical example.

To fix notation, let $A$ be a symmetric positive definite matrix, and consider the “coarse-to-fine” half V-cycle multigrid iteration $x_{n+1} = x_n + M (b - A x_n)$ determined by the iteration matrix

\[ E = I - M A = G K \]

where $G$ is an iteration matrix for the smoother, e.g., $G = (I - B A)^k$ for $k$ steps of the iteration $x_{i+1} = x_i + B(b - A x_i)$, and $K$ is the coarse grid correction

\[ K \isdef I - P M_c P^* A,\]

where the prolongation operator $P$ is an injection from a space of smaller dimension, and
where $M_c$ is defined recursively in the same way as $M$, but for the system governed by the Galerkin coarse operator

\[ A_c \isdef P^* A P, \]

and with the recursion stopping by taking $M_c = A_c^{-1}$ at the coarsest level.
We will restrict our attention primarily to bounds of $\norm{E}_A^2$, the square of the contraction factor in the energy norm of this half V-cycle. As is well-known (see, e.g., <cit.>), this suffices for the analysis of general V-cycles.
If we denote the energy inner product adjoint for an operator $T$ by $T^{\dagger}$, so that
$  A T^{\dagger} = T^* A, $
i.e., $T^{\dagger} \isdef A^{-1} T^* A$, then $E^{\dagger}$ is the iteration matrix for a “fine-to-coarse” half V-cycle iteration, with equal norm $\norm{E^{\dagger}}_A = \norm{E}_A$. A general V-cycle $E_V$ decomposes as some “fine-to-coarse” half V-cycle $(E')^\dagger$ followed by some “coarse-to-fine” half V-cycle $E$,
$ E_V = E (E')^{\dagger}, $
so that

\[ \norm{E_V} \le \norm{E}_A \norm{E'}_A . \]

For a symmetric V-cycle with $E' = E$ this inequality strengthens to the equality $\norm{E_V}_A = \norm{E}_A^2$ .

Instead of using the simple smoothing iteration $G = (I - BA)^k$
we will consider polynomial iterations

\[ G = p_k(BA) \]

for some polynomial $p_k(x)$ of degree $k$ with $p_k(0) = 1$. Picking $p_k$ to optimize $\norm{G}_A$ leads to Krylov acceleration methods like the Chebyshev semi-iterative method and the method of conjugate gradients. Our central concern will be instead to investigate how to choose $p_k$ to optimize the multigrid convergence rate, through achieving optimal or nearly optimal bounds on $\norm{E}_A^2$. As polynomial methods tend to work best for symmetric iterations, we will confine our attention to the case of symmetric positive definite $B$. For example, Jacobi corresponds to $B = D^{-1}$ where $D$ is the diagonal part of $A$. Furthermore, so as to remove some redundancy, we will assume $B$ is scaled such that

\[ \rho(B A) = 1 \]

where $\rho$ denotes the spectral radius. In particular, these assumptions imply that the single step ($k=1$) smoother

\[ G = I - \omega B A \]

has contraction factor $\norm{G}_A = \rho(G) < 1$ whenever the damping parameter $\omega$ is in the range $0 < \omega < 2$.

Our results will all be expressed in terms of the approximation property

\begin{equation} \label{eq:C-def} C \isdef \norm{A^{-1} - P A_c^{-1} P^*}_{A,B}^2 \isdef \sup_{\norm{f}_B \le 1} \norm{(A^{-1} - P A_c^{-1} P^*)f}_{A}^2 . \end{equation}

Hackbusch <cit.> proved the bound

\begin{equation}\label{eq:hackbusch-v-cycle-bound} \norm{E}_A^2 \le \frac{C}{C+2k} \end{equation}

for the case $G = (1 - BA)^k$, i.e., with $\omega = 1$. Using as our main tool the general V-cycle theory developed by McCormick <cit.>, we will generalize this to bounds depending on $C$ and $p_k$. We will show that using a simple Chebyshev-type iteration for the smoother, based on a three-term recurrence, achieves the bound

\begin{equation} \label{eq:cheb-smoother-simple-bound} \norm{E}_A^2 \le \frac{C}{C+\frac{4}{3}k(k+1)} \; . \end{equation}

The polynomials associated with this Chebyshev-type iteration are the lesser-known Chebyshev polynomials of the third and fourth kinds, also known as “airfoil” polynomials within the aerodynamics community, where they have been extensively studied <cit.>. Within the multigrid community, they have found extensive use in smoothed aggregation methods, originally as the optimal polynomials for constructing the smoothing operator to apply to some tentative prolongator, but also used for the smoothing iteration itself. The first appearance traces back to the thesis of Brezina <cit.>, with continued development in papers by Vanêk, Brezina, Tezaur, and  Krızková <cit.>, particularly <cit.>. The use of these polynomials for smoothing iterations outside of the context of smoothed aggregation was suggested by Adams, Brezina, Hu, and Tamaro <cit.>, though at the time, the authors were apparently only aware of the existence of these polynomials for specific degrees (1, 4, and 13 being the first three). As mentioned above, Adams et al. also suggested the use of the usual Chebyshev semi-iteration based on the polynomials of the first kind, adjusted to target the range of
eigenvalues $[\lambda^*, \lambda_{\max}]$ instead of $[\lambda_{\min}, \lambda_{\max}]$, with $\lambda^*$ an ad-hoc parameter.

Our main contribution is to point out that there is a very simple iteration based on the three-term recurrence satisfied by the fourth-kind Chebyshev polyonimals, namely

\[ \begin{aligned}
x_k &= x_{k-1} + z_k, \\
z_0 &= 0, \\
z_k &= \frac{2k-3}{2k+1} z_{k-1} + \frac{8k-4}{2k+1} \frac{1}{\rho(BA)} B (b - A x_{k-1}),
\end{aligned}\]

and to give full multi-level V-cycle convergence bounds, particularly (<ref>), when this iteration is used as the smoother. Note the iteration is just as simple, if not simpler, than the standard (first-kind) Chebyshev semi-iterative method, and for use as a smoother avoids the ad-hoc choice of parameter $\lambda^*$. The fourth-kind Chebyshev polynomials and associated iteration deserves more wide-spread recognition in the broader multigrid community, as a simple method for accelerating an arbitrary multigrid V-cycle for SPD systems.

§ FOURTH-KIND CHEBYSHEV ITERATION

Let us first consider the bound due to Hackbusch <cit.> for the two-level iteration contraction factor (taking $M_c = A_c^{-1}$),

\begin{equation} \label{eq:weak-bound}
\begin{split}
\norm{E}_A &= \norm{E^{\dagger}}_A = \norm{(I - P A_c^{-1} P^* A) G}_A \\
&\le \norm{A^{-1} - P A_c^{-1} P^*}_{A,B} \norm{AG}_{B,A} \\
&= C^{1/2} \sup_{\lambda \in \sigma(BA)} \lambda^{1/2} |p_k(\lambda)| \\
&\le C^{1/2} \sup_{0 < \lambda \le 1} \lambda^{1/2} |p_k(\lambda)|,
\end{split}
\end{equation}

where, recall, we have taken the smoothing iteration matrix to be $G = p_k(B A)$, with $p_k$ a polynomial of degree $k$ and $p_k(0) = 1$. This bound has its shortcomings: it is only for a two-level method, and it need not be smaller than one, whereas we always have $\norm{E}_A < 1$ provided $\norm{G}_A < 1$. We will be considering other bounds without these particular shortcomings later. The reason for considering the bound (<ref>) is that the problem of picking optimal polynomials for it has an elegant solution.

Optimizing the bound (<ref>) involves solving a weighted minimax polynomial approximation problem, with weight $w(x) = x^{1/2}$. Luckily, the solution is known. The unique monic polynomials minimizing the maximum of $(1+x)^{1/2} |p(x)|$ on $[-1, 1]$ are $2^{-n}V_n(x)$, where $V_n(x)$ are the Chebyshev polynomials of the third kind, see Mason <cit.>. Up to scaling, these are Jacobi polynomials with $\alpha=-\tfrac{1}{2}$, $\beta=\tfrac{1}{2}$, but, like the Chebyshev polynomials of the first and second kinds, enjoy particularly simple forms of the various special properties of Jacobi polynomials (see <cit.>). For our purposes, it will be slightly more convenient to use the Chebyshev polynomials of the fourth kind $W_n(x)$, which are essentially the same: $W_n(x) = (-1)^n V_n(-x)$. They are given by

\[
W_n(x) = \frac{\sin \, (n+\tfrac{1}{2})\theta}{\sin \tfrac{1}{2}\theta}
\]

where $x = \cos \theta$. Let us map $x \in [-1, 1]$ to $\lambda \in [0, 1]$ via $x = 1-2\lambda$. The half angle formula for sine gives $\sin \tfrac{1}{2}\theta = \lambda^{1/2}$, and we see that the weighted polynomials

\begin{equation} \label{eq:eqos}
\lambda^{1/2} W_n(1 - 2\lambda) = \sin \, (n+\tfrac{1}{2})\theta
\end{equation}

equioscillate between $\pm 1$. To enforce the constraint $p_k(0) = 1$ we need to scale by the inverse of
$ W_n(1) = 2n+1, $

\begin{equation}\label{eq:cheb-pk} p_k(\lambda) = \frac{1}{2k+1} W_k(1 - 2\lambda). \end{equation}

Using (<ref>) we see that

\begin{equation}\label{eq:cheb-weak-bound} \sup_{0 \le \lambda \le 1} \lambda^{1/2} |p_k(\lambda)| = \frac{1}{2k+1}. \end{equation}

The first few $p_k(\lambda)$ are shown in <ref> and the corresponding weighted polynomials $\lambda^{1/2} p_k(\lambda)$, exhibiting the equioscillation property, in <ref>. <Ref> anticipates the multilevel bound of the next section. Explicitly, the first four polynomials are

\[\begin{aligned}
p_0(\lambda) &= 1, \\
p_1(\lambda) &= 1 - \tfrac{4}{3} \lambda, \\
p_2(\lambda) &= 1 - 4\lambda + \tfrac{16}{5} \lambda^2, \\
p_3(\lambda) &= 1 - 8\lambda + 16 \lambda^2 - \tfrac{64}{7}\lambda^3 .
\end{aligned}
\]

ymin=-.4, ymax=1.1,
xmin=0, xmax=1,
grid = both,

cycle multi list=Dark2,
trig format plots=rad,
+[ domain=0:pi, samples=201,] (lambda(x), pcheb(#1, x));

[][$\lambda^{1/2} p_k(\lambda)$]

ymin=-.4, ymax=.4,
xmin=0, xmax=1,
grid = both,
cycle multi list=Dark2,
trig format plots=rad,
+[ domain=0:pi, samples=201,] (lambda(x), sqrt(lambda(x))*pcheb(#1, x));

[][$\lambda^{1/2} p_k(\lambda) / \sqrt{1 - p_k(\lambda)^2}$]

ymin=-0.4, ymax=0.65, ytick distance=0.2,
xmin=0, xmax=1,
grid = both,

cycle multi list=Dark2,
trig format plots=rad,
cycle list shift=1
+[ domain=0:pi, samples=201,] (lambda(x), fcheb(#1, x));

The first few scaled and shifted Chebyshev polynomials of the fourth kind, (a) unweighted, (b) weighted by $\lambda^{1/2}$, (c) weighted by $\sqrt{\lambda/(1-p_k(\lambda)^2)}$

The polynomials $W_n$ satisfy the same recurrence as the Chebyshev polynomials $T_n$, $U_n$ of the first and second kinds, just with a different initial condition <cit.>.

\[ W_n(x) = 2xW_{n-1}(x) - W_{n-2}(x), \qquad W_0(x) = 1, W_1(x) = 2x + 1 . \]

It follows that the scaled and shifted polynomials $p_k$ satisfy the recurrence

\[ p_k(\lambda) = \frac{4k-2}{2k+1} (1 - 2 \lambda) p_{k-1}(\lambda) - \frac{2k-3}{2k+1} p_{k-2}(\lambda). \]

We can derive an iteration based on this recurrence, exactly as is done to derive the Chebyshev semi-iterative method. The result is the iteration

\begin{equation}\label{eq:cheb-iter} \begin{aligned}
x_k &= x_{k-1} + z_k, \\
z_0 &= 0, \\
z_k &= \frac{2k-3}{2k+1} z_{k-1} + \frac{8k-4}{2k+1} \frac{1}{\rho(BA)} B (b - A x_{k-1}),
\end{aligned}\end{equation}

satisfying $A^{-1} b - x_k = p_k(B A) (A^{-1} b - x_0)$. For the sake of easy reference, we have included the scaling of $B$ in (<ref>) so that it works without our assumption $\rho(B A) = 1$. The bound

\begin{equation} \label{eq:cheb-weak-bound-alt} \norm{A p_k(B A)}_{B,A} \le \frac{1}{2k+1} \end{equation}

follows from (<ref>), leading to the two-level contraction factor estimate

\begin{equation} \label{eq:cheb-two-level-bound} \norm{E}_A \le \frac{C^{1/2}}{2k+1} . \end{equation}

§ CONJUGATE GRADIENTS FOR SMOOTHERS

The polynomials $p_k$ and corresponding iteration (<ref>) of the previous section were chosen to minimize $\norm{AG}_{B,A}$ for any potential $B$ with $\sigma(B A) \in (0, 1]$. If we allow $G$ to be nonlinear, we can instead consider Krylov methods that directly minimize $\norm{A G e_0}_B$, i.e., the $B$-norm of the residual.

One possibility is to apply preconditioned conjugate gradients (CG) to the transformed system

\begin{equation}\label{eq:ABAx=ABb} A B A x = A B b \end{equation}

with preconditioner $A^{-1}$. Since CG minimizes the energy norm of the error $x - x_k$ where $x = A^{-1} b$ is the exact solution, applied to the transformed system CG minimizes $\norm{x - x_k}_{ABA} = \norm{b - A x_k}_{B}$, i.e., the $B$-norm of the residual $r_k = b - A x_k$. By taking $A^{-1}$ as preconditioner, the CG iterates $x_k$ will be in the Krylov space spanned by $(BA) e_0, (BA)^2 e_0, \dotsc, (BA)^k e_0$ where $e_0$ is the initial error. As the Chebyshev iterates (<ref>) are also in this space, the bound on the $B$-norm of their residuals

\[ \norm{r_k}_B \le \frac{\norm{r_0}_{A^{-1}}}{2k+1}, \]

a restatement of (<ref>), must also hold for the residuals of the CG iteration.

After some minor manipulation, applying the standard CG algorithm to (<ref>) using preconditioner $A^{-1}$ leads to the equivalent iteration

\begin{equation}\label{eq:CG} \begin{aligned}
p_0 &= z_0 = B r_0, \quad g_0 = A z_0, \\
\alpha_k &= \frac{\norm{z_k}_A^2}{\norm{g_k}_B^2}, \\
x_{k+1} &= x_k + \alpha_k p_k, \\
z_{k+1} &= z_k - \alpha_k B g_k, \\
\beta_{k+1} &= \frac{\norm{z_{k+1}}_A^2}{\norm{z_{k}}_A^2}, \\
p_{k+1} &= z_{k+1} + \beta_k p_k, \\
g_{k+1} &= A p_{k+1} = A z_{k+1} + \beta_k g_k .
\end{aligned}\end{equation}

Note that the action of the preconditioner $A^{-1}$ is not actually needed. The cost of (<ref>) is similar to the usual preconditioned CG method: each iteration requires one application of $A$ to form $A z_k$, one application of $B$ to form $B g_k$, and two inner products.

§ MULTILEVEL ANALYSIS

Let us recap the core ideas of the general V-cycle theory developed by McCormick <cit.>. The central idea is to decompose the solution space $\mathsf{V} = \mathbb{R}^n$ into the direct sum of the coarse space (the range of $P$) and its $A$-orthogonal complement,

\[ \mathsf{V} = \ran(P) \oplus \ran(P)^{\perp_A} . \]

Let us denote the pair of projections corresponding to this decomposition by

\[ \pi_c \isdef P A_c^{-1} P^* A, \quad\text{and}\quad \pi_f \isdef I - \pi_c. \]

Being complementary projections, these satisfy the relations $\pi_c^2 = \pi_c$, $\pi_f^2 = \pi_f$, $\pi_c \pi_f = \pi_f \pi_c = 0$, and $\pi_c + \pi_f = I$. Being $A$-orthogonal, they additionally are self adjoint in the energy inner product, $A \pi_c = \pi_c^* A$, $A \pi_f = \pi_f^* A$, and satisfy the Pythagorean theorem

\[ \norm{u}_A^2 = \norm{\pi_f u}_A^2 + \norm{\pi_c u}_A^2 . \]

Note that $\pi_f = I - P A_c^{-1} P^* A$ is the coarse grid correction with an exact coarse grid solve, and that the approximation property constant $C$ defined by (<ref>) can be written as

\[ C = \norm{\pi_f A^{-1}}_{A, B}^2. \]

We will soon need an alternate expression for $C$, which we derive by first noting

\[\begin{split}
C &= \norm{\pi_f A^{-1}}_{A, B}^2 = \rho(B^{-1} A^{-1} \pi_f^* A \pi_f A^{-1}) \\
&= \rho(B^{-1} \pi_f A^{-1} \pi_f^*)
= \rho(A^{-1} \pi_f^* B^{-1} \pi_f) \\
&= \sup_{u \ne 0} \frac{\norm{\pi_f u}_{B^{-1}}^2}{\norm{u}_A^2}
= \sup_{u \ne 0} \frac{\norm{\pi_f u}_{B^{-1}}^2}{\norm{\pi_f u}_A^2 + \norm{\pi_c u}_A^2}
\end{split}\]

where we have used that $\pi_f A^{-1} = A^{-1} \pi_f^*$ and that the spectral radius is preserved by cyclic permutations of factors,
$\rho(TU) = \rho(UT)$. Next, observe that the maximum is attained when $\pi_c u = 0$, so that we may restrict $u$ to $\nul(\pi_c) = \ran(\pi_f) = \ran(P)^{\perp_A}$, giving

\begin{equation} \label{eq:C-def-alt}
C = \norm{\pi_f A^{-1}}_{A, B}^2 = \sup_{u \in \ran(\pi_f) \setminus \{ 0 \}} \frac{\norm{u}_{B^{-1}}^2}{\norm{u}_{A}^2} \ge 1,
\end{equation}

where the lower bound of one follows from our assumption $\rho(BA) = 1$. Note that, without the restriction of $u$ to $\ran(\pi_f)$, the supremum of $ \norm{u}_{B^{-1}}^2/\norm{u}_{A}^2$ is the reciprocal of the smallest eigenvalue of $BA$, which is the condition number of $BA$ as $\rho(BA) = 1$. So, an (imprecise) way to think of $C$ is as the condition number of $BA$ “restricted” to the $A$-orthogonal complement of the coarse space.

The decomposition $e = \pi_f e + \pi_c e$ is useful for the multilevel analysis because the multilevel coarse grid correction

\[ K \isdef I - P M_c P^* A \]

leaves the space $\ran(\pi_f) = \ran(P)^{\perp_A}$ invariant, while acting on the coarse space by $E_c \isdef I - M_c A_c$,

\[ K \pi_f = \pi_f, \quad K P = P E_c . \]

Hence, for any $e$,

\[ K e = \pi_f e + P E_c e_c, \quad\text{where}\quad P e_c = \pi_c e, \]

i.e., $e_c = A_c^{-1} P^* A e$. If we now assume the smoother satisfies

\begin{equation} \label{eq:smoothing-factor}
\norm{G e}_A^2 \le \alpha \norm{\pi_f e}_A^2 + \norm{\pi_c e}_A^2
\end{equation}

for some $\alpha < 1$, which McCormick calls the “smoothing factor,” we then have

\[ \norm{E e}_A^2 = \norm{G K e}_A^2 = \norm{G (\pi_f e + P E_c e_c)}_A^2 \le \alpha \norm{\pi_f e}_A^2 + \norm{P E_c e_c}_A^2, \]

but since

\[ \norm{P E_c e_c}_A = \norm{E_c e_c}_{A_c} \le \norm{E_c}_{A_c} \norm{e_c}_{A_c} = \norm{E_c}_{A_c} \norm{\pi_c e}_{A}, \]

it follows that

\[ \norm{E e}_A^2 \le \max(\alpha, \norm{E_c}_A^2) (\norm{\pi_f e}_A^2 + \norm{\pi_c e}_A^2)
= \max(\alpha, \norm{E_c}_A^2) \norm{e}_A^2. \]

If the smoother on every level satisfies the appropriate version of (<ref>), we see by induction that $\norm{E}_A^2$ is bounded by the largest smoothing factor constant $\alpha$ across all levels.

Let us assume the smoother iteration matrix $G$ is a contraction in the energy norm,

\[ \norm{G}_A < 1 . \]

Then the assumption (<ref>) is satisfied if we take

\[ \alpha = \sup_{\pi_f e \ne 0} \frac{\norm{G e}_A^2 - \norm{\pi_c e}_A^2}{\norm{\pi_f e}_A^2}
= 1 - \inf_{\pi_f e \ne 0} \frac{\norm{e}_A^2 - \norm{G e}_A^2}{\norm{\pi_f e}_A^2} . \]

Define the symmetrized smoother $N$ so that the corresponding iteration matrix is

\begin{equation}\label{eq:1-NA} I - N A = G^{\dagger} G = A^{-1} G^* A G , \end{equation}

i.e., the smoother followed by its adjoint. Explicitly, $N \isdef A^{-1} - A^{-1} G^* A G A^{-1}$ and we see that

\begin{equation}\label{eq:N<Ai} 0 < N \le A^{-1} . \end{equation}

That is, $N$ is SPD and $A^{-1} - N$ is symmetric positive semidefinite.

\[ \norm{e}_A^2 - \norm{G e}_A^2 = ([A - G^* A G]e, e) = \norm{A e}_N^2, \]

we can take

\[ \alpha = 1 - C_N^{-1} \]

where we define

\[ C_N \isdef \sup_{e \ne 0} \frac{\norm{\pi_f e}_A^2}{\norm{A e}_N^2} =
\sup_{r \ne 0} \frac{\norm{\pi_f A^{-1} r}_A^2}{\norm{r}_N^2} = \norm{\pi_f A^{-1}}_{A, N}^2, \]

i.e., the approximation property constant $C$ but with the symmetrized smoother $N$ in place of $B$. Note that (<ref>) holds for $C_N$ with $N$ in place of $B$. In particular we have the lower bound $C_N \ge 1$ since (<ref>) implies $A \le N^{-1}$.

Summarizing, we have the following lemma. This is essentially McCormick's Lemma 2.3 <cit.>, though with an alternate expression for $C_N$ (corresponding to the reciprocal of McCormick's $\beta$) anticipating some of the ideas from McCormick's Theorem 3.4. The result in the same form can also be found in Vassilevski <cit.> as Theorem 5.16 and Lemma 5.17.
Note that the lemma holds for general smoother iterations $G$, as we did not make use of any additional assumptions (e.g., symmetry).

If $\norm{G}_A < 1$ holds on each level, then the V-cycle contraction factor is bounded by

\[ \norm{E}_A^2 \le 1 - C_N^{-1} \]

where $C_N$ is the largest value across all levels of the approximation property constant

\[ \norm{\pi_f A^{-1}}_{A, N}^2 = \sup_{u \in \ran(\pi_f) \setminus \{ 0 \}} \frac{\norm{u}_{N^{-1}}^2}{\norm{u}_{A}^2} \ge 1 \]

defined in terms of the symmetrized smoother $N \isdef A^{-1} - A^{-1} G^* A G A^{-1}$.

It is now fairly straight-forward to prove the following V-cycle bound for our main case of interest, polynomial smoothers $G = p(BA)$. Like the two-level bound (<ref>), this bound only depends on $C = \norm{\pi_f A^{-1}}_{A,B}^2$ and $p$. The proof is quite similar to that of Theorem 3.4 of McCormick <cit.>.

Let the smoother iteration (on each level) be given by

\[ G = p(BA) \]

where $B$ is symmetric positive definite and scaled such that $\rho(BA) = 1$, and $p(x)$ is a polynomial satisfying $p(0) = 1$ and $|p(x)| < 1$ for $0 < x \le 1$, possibly different on each level. Then the V-cycle contraction factor $\norm{E}_A^2$ is bounded by the largest value across all levels of

\[ \frac{C}{C+\gamma^{-1}} \]

where $C$ is the approximation property constant $\norm{\pi_f A^{-1}}_{A, B}^2$ and

\begin{equation} \label{eq:gamma-defn} \gamma = \sup_{0 < \lambda \le 1} \frac{\lambda \, p(\lambda)^2}{1 - p(\lambda)^2} . \end{equation}

We have

\[\begin{split}
C_N &= \sup_{u \in \ran(\pi_f) \setminus \{ 0 \}} \frac{\norm{u}_{N^{-1}}^2}{\norm{u}_{A}^2}
= 1 + \sup_{u \in \ran(\pi_f) \setminus \{ 0 \}} \frac{\norm{u}_{N^{-1}}^2 - \norm{u}_{A}^2}{\norm{u}_{A}^2} \\
&\le 1 + \left(\sup_{u \ne 0} \frac{\norm{u}_{N^{-1}}^2 - \norm{u}_{A}^2}{\norm{u}_{B^{-1}}^2} \right)
\left(\sup_{u \in \ran(\pi_f) \setminus \{ 0 \}} \frac{\norm{u}_{B^{-1}}^2}{\norm{u}_{A}^2} \right) \\
&\le 1 + \gamma C
\end{split}\]

\[\begin{split}
\sup_{u \ne 0} \frac{\norm{u}_{N^{-1}}^2 - \norm{u}_{A}^2}{\norm{u}_{B^{-1}}^2}
&= \rho(B[N^{-1} - A]) \\
&= \rho(BA [(I - p(BA)^2)^{-1} - I]) \\
&= \sup_{\lambda \in \sigma(BA)} \lambda \, [(1 - p(\lambda)^2)^{-1} - 1] \le \gamma
%= \sup_{\lambda \in \sigma(BA)} \frac{\lambda \, p(\lambda)^2}{1 - p(\lambda)^2} \\
%&\le \gamma
\end{split}\]

where we have used that $N^{-1} = A (I - G^2)^{-1}$, which follows from (<ref>) as $G^{\dagger} = G$. The result now follows by applying the previous lemma.

The reason this is a useful bound is that, for reasonable polynomials, the ratio $\lambda p(\lambda)^2 / (1 - p(\lambda)^2)$ in the supremum defining $\gamma$ is significantly less than one throughout the entire interval $(0, 1]$. The presence of $p(\lambda)^2$ in the numerator makes the ratio small wherever $p(\lambda)$ is small, while for $\lambda \to 0$ where $p(\lambda) \to 1$,

\begin{equation}\label{eq:lim-0}
\lim_{\lambda \to 0} \frac{\lambda \, p(\lambda)^2}{1 - p(\lambda)^2}
= \frac{1}{-2p'(0)}
\end{equation}

by an application of L'Hôpital's rule.

The key step in the above proof, which comes from McCormick's proof of Theroem 3.4 in <cit.>, is to bound the ratio $(\norm{u}_{N^{-1}}^2 - \norm{u}_{A}^2)/\norm{u}_{A}^2$ determining $C_N-1$ instead of $\norm{u}_{N^{-1}}^2/\norm{u}_{A}^2$ determining $C_N$. Without this step, the obvious approach might have been to bound the factors of

\[ \frac{\norm{u}_{N^{-1}}^2}{\norm{u}_{A}^2} =  \frac{\norm{u}_{B^{-1}}^2}{\norm{u}_{A}^2} \frac{\norm{u}_{N^{-1}}^2}{\norm{u}_{B^{-1}}^2}, \]

leading to

\[
C_N \le C \sup_{0 < \lambda \le 1} \frac{\lambda}{1 - p(\lambda)^2} .
\]

But without $p(\lambda)^2$ in the numerator, the supremum is greater than one, rendering the bound useless.

As a quick sanity check, we note that, as a direct consequence, this lemma recovers Hackbusch's bound (<ref>) on the contraction factor for the V-cycle with a simple iteration.

If the smoother is $k$ steps of a simple damped iteration,

\[ G = (I - \omega BA)^k, \]

for possibly different $0 < \omega < 2$ and $k$ on each level,
then $\norm{E}_A^2$ is bounded by the largest value of

\[ \frac{C}{C + 2 \omega k} \]

\begin{equation}\label{eq:omega-condition} (1 - \omega)^{2k} \le \frac{1}{1+2\omega k} . \end{equation}

In particular, $\omega \le \tfrac{3}{2}$ is necessary for $k=1$ and sufficient for all $k$. Asymptotically, the largest permissible $\omega$ for the bound to hold is

\[ \omega = 2 - \frac{\log 4k}{2k} + o(k^{-1}) \quad\text{as } k \to \infty . \]

Here $p(\lambda) = (1-\omega \lambda)^k$ and $f(\lambda) \isdef \lambda p(\lambda)^2 / (1 - p(\lambda)^2)$ is convex on $[0, 1]$, so that the supremum of $f$ is attained at either $\lambda = 1$ where $f(1)^{-1} = (1-\omega)^{-2k}-1$ or (in the limit) at $\lambda = 0$, where $\lim_{\lambda \to 0} f(\lambda)^{-1} = 2 \omega k$. The latter is smaller when (<ref>) holds.

§.§ Fourth-kind Chebyshev Iteration

The bound for the V-cycle using the fourth-kind Chebyshev iteration as smoother also follows directly from <ref>.

If (<ref>) is used for the smoother iteration, so that

\[ G = p_k(BA), \]

where $p_k$ are the scaled and shifted Chebyshev polynomials of the fourth kind given by (<ref>), then $\norm{E}_A^2$ is bounded by the largest value of

\begin{equation}\label{eq:cheb-simp-bound} \frac{C}{C + \tfrac{4}{3} k (k+1) }. \end{equation}

This follows by applying (<ref>), since the supremum determining $\gamma$ for these polynomials is attained in the limit $\lambda \to 0$, as illustrated in <ref>.

It is possible to somewhat improve the bound of <ref> and thereby <ref>.
As the derivation is quite technical, and the improvement is only marginal unless the approximation property $C$ is quite small, e.g., $C \approx 2$, we leave the derivation to <ref> and <ref>. Here we simply state the improved bound for the fourth-kind Chebyshev iteration.

If (<ref>) is used for the smoother iteration, so that

\[ G = p_k(BA), \]

where $p_k$ are the scaled and shifted Chebyshev polynomials of the fourth kind given by (<ref>), then $\norm{E}_A^2$ is bounded by the largest value of

\begin{equation}\label{eq:cheb-bound-sharp} \frac{C'}{C' + \tfrac{4}{3} k (k+1) }, \qquad C' = C - \beta f, \end{equation}

with $C = \norm{\pi_f A^{-1}}_{A,B}^2$ the approximation property constant, $\beta$ the constant

\[ \beta \isdef 1 - \tfrac{1}{3} \left[ \inf_{4<\phi\le\tfrac{3}{2}\pi} \csc^2 \phi - \phi^{-2} \right]^{-1}
= 0.650914713503148\ldots \ge \frac{41}{63}
\]

and $0 \le f \le \tfrac{2}{3} \beta^{-1} < 1.03$ given by

\[
f = \begin{cases}
\left(1+\frac{65}{15n^2-33}\right)\left(1-\frac{40 C}{10n^2-n+19}\right), & C < \frac{n^2}{7.97} + 0.512 \\
\frac{1+2C}{32C^2 - 10} (n^2 +1.79 +n^{-2}), & \text{otherwise,}
\end{cases}
\]

where $n = 2k+1$. In particular, $f \to 1$ as $k \to \infty$.

The sharpened bound (<ref>) of <ref> is of exactly the same form as (<ref>) of <ref>, but with the approximation constant $C$ “discounted” by a little less than $2/3$ as the number of steps $k \to \infty$.

§.§ Optimized Polynomials

The scaled and shifted Chebyshev polynomials of the fourth kind (<ref>) are not optimal for the multilevel bound of <ref>, as they minimize the supremum of $\lambda p(\lambda)^2$ and not $\lambda p(\lambda)^2/(1-p(\lambda)^2)$. The optimal polynomials for the multilevel bound can be found numerically (see <ref> for details). Explicitly, the first two are

\[\begin{aligned}
p_1(\lambda) &= 1 - \tfrac{3}{2} \lambda, \\
p_2(\lambda) &= (1 - \tfrac{\sqrt{5}}{2} \lambda)(1-\tfrac{5+\sqrt{5}}{2} \lambda) = 1 - (\tfrac{5}{2} + \sqrt{5})\lambda + \tfrac{5}{4}(1+\sqrt{5}) \lambda^2 .
\end{aligned}\]

The first few optimal polynomials are shown in <ref>, and they can be seen to exhibit the equioscillation property in <ref>. Numerical evidence (see <ref>) supports the conjecture that, asymptotically, these polynomials achieve

\[
\gamma_{\text{opt}}^{-1} \sim 4\left(\frac{2k+1}{\pi}\right)^2 - \frac{2}{3} + \frac{1}{60} \left(\frac{2k+1}{\pi}\right)^{-2} + O(k^{-4}) \quad\text{as } k \to \infty,
\]

with, moreover, the first two terms of the series giving a fairly good lower bound of $\gamma_{\text{opt}}^{-1}$. This would imply the error contraction factor for the V-cycle using these optimal polynomials as smoothers is bounded by the largest value of

\begin{equation}\label{eq:opt-bound}
\frac{C}{C+\tfrac{4}{\pi^2}(2k+1)^2-\tfrac{2}{3}}
\end{equation}

across all levels. In the limit $k \to \infty$, this is an improvement over the bound (<ref>) for the fourth-kind Chebyshev iteration by the factor $\frac{\pi^2}{12} \approx 0.822$, i.e., a reduction in the contraction factor bound by about 18%.

On the other hand, the improved bound (<ref>) for the fourth-kind Chebyshev iteration is asymptotically better than (<ref>) when
$C < \beta/(1-\pi^2/12) \approx 3.67$. Note that (<ref>) is the bound for the polynomials optimizing the bound of <ref>, whereas (<ref>) was based on the sharper bound of <ref>. Finding the polynomials to optimize the sharper bound of <ref>, even numerically, appears to be a much more difficult problem, and it is not obvious the benefit would be more than marginal.

ymin=-0.6, ymax=1.1,
xmin=0, xmax=1,
grid = both,

cycle multi list=Dark2,
trig format plots=rad,
cycle list shift=1
table[x=x,y=p#1] optp.dat;

[][$\lambda^{1/2} p_k(\lambda) / \sqrt{1 - p_k(\lambda)^2}$]

ymin=-0.65, ymax=0.65, ytick distance=0.2,
xmin=0, xmax=1,
grid = both,

cycle multi list=Dark2,
trig format plots=rad,
cycle list shift=1
table[x=x,y=p#1] optpw.dat;

The optimal polynomials $p_k$ for the bound of <ref>, (a) unweighted, (b) weighted by $\sqrt{\lambda/(1-p_k(\lambda)^2)}$

2c$\frac{4}{\pi^2}(2k+1)^2 - \frac{2}{3}$
1     3             2   981        1   $91 \times 10^{-2}$     1   $83 \times 10^{-2}$

2     9   4721      9   465        6   $68 \times 10^{-3}$     6   $58 \times 10^{-3}$

3     19   1957     19   1923      3   $38\times 10^{-3}$     3   $36 \times 10^{-3}$

4     32   1634     32   1614      2   $04\times 10^{-3}$     2   $03 \times 10^{-3}$

5     48   3742     48   3728      1   $36\times 10^{-3}$     1   $36 \times 10^{-3}$

10     178   0643     178   0639      3   $73\times 10^{-4}$     3   $73 \times 10^{-4}$

100     16373   241899     16373   241895      4   $07\times 10^{-6}$     4   $07 \times 10^{-6}$

1000     1   $623\times 10^{6}$     1   $623\times 10^{6}$      4   $11\times 10^{-8}$     4   $11 \times 10^{-8}$

Optimal $\gamma^{-1}$ for <ref>

§ IMPLEMENTATION ISSUES

The iteration (<ref>) is the natural implementation for a smoother using the Chebyshev polynomials of the fourth kind, derived from the three-term recurrence they satisfy. The numerically determined optimal polynomials for the V-cycle bound of <ref> do not satisfy a three-term recurrence, and so require a different implementation.

One obvious approach is to translate the factorization into monomials involving the roots $r_i$,

\[ p_k(\lambda) = \prod_{i=1}^k \left(1 - \frac{\lambda}{r_i}\right) \, ,\]

directly into the iteration

\[ x_{i+1} = x_i + \frac{1}{r_i} B (b - A x_i) . \]

While probably fine for small $k$, a concern of this approach is that the polynomials corresponding to the intermediate steps are poorly behaved. Depending on the chosen ordering of roots, the iteration may first dampen some error modes by several orders of magnitudes before then amplifying them by several orders of magnitude, or vice versa.

Instead, we propose using an iteration based on the expansion in the basis of fourth-kind Chebyshev polynomials

\[p_k(\lambda) = \sum_{i=0}^k \alpha_i W_i(1 - 2\lambda), \]

exploiting the fact that we expect the orthogonal (w.r.t. a weight function) Chebyshev polynomials of the fourth kind $W_i$ to be a good basis for representing the optimal $p_k$ and the fact that we have a nice iteration, namely (<ref>), for computing the basis $W_i$. The Gaussian quadrature rule associated with the polynomials $W_i$ can be used to determine the coefficients $\alpha_i$, requiring only that $p_k$ be evaluated at the quadrature nodes (see Appendix A for details). We propose then the iteration

\begin{equation}\label{eq:opt-iter} \begin{aligned}
z_0 &= 0, \quad r_0 = b - A x_0, \\
z_i &= \frac{2i-3}{2i+1} z_{i-1} + \frac{8i-4}{2i+1} \frac{1}{\rho(BA)} B r_{i-1}, \\
x_i &= x_{i-1} + \beta_i z_i, \quad r_i = r_{i-1} - A z_i.
\end{aligned}\end{equation}

This is equivalent to (<ref>) when we take $\beta_i = 1$. Note that $\beta_i$ does not appear in the “residual” update, so that it is not the case that $r_i = b - A x_i$. Instead, $r_i$ tracks the residual of the unmodified fourth-kind Chebyshev iteration, and the vectors $z_i$ are identical in both iterations. From the properties of (<ref>), we know that $z_i = [\tfrac{1}{2i-1}W_{i-1}(BA) - \tfrac{1}{2i+1}W_i(BA)] e_0$ where $e_0 = A^{-1} b - x_0$. It follows that the polynomial $p_k(\lambda)$ corresponding to the iteration (<ref>), i.e., so that $A^{-1} b - x_k = p_k(BA) e_0$, is given by

\[p_k(\lambda) = \sum_{i=0}^k \frac{\beta_i-\beta_{i+1}}{2i+1} W_i(1-2\lambda) \]

provided we let $\beta_0 = 1$ and $\beta_{k+1} = 0$. Hence, the coefficients $\beta_i$ can be found from the coefficients $\alpha_i$ via the recursion

\[ \beta_{i+1} = \beta_i - (2i+1)\alpha_i . \]

For the polynomials that optimize the multilevel bound (<ref>), the resulting $\beta_i$ (which depend on $k$) are in the interval $[1, 1.6)$, for $k$ up to several hundred. We can thus think of the optimized iteration (<ref>) as an “over-relaxed” version of the fourth-kind Chebyshev iteration (<ref>). Explicit Matlab code for computing the coefficients $\beta_i$ is given in <ref>.

A natural question that arises is how to pick the optimal number of smoothing steps $k$. For the same total work, one must choose between more V-cycles with fewer smoothing steps  per cycle versus fewer cycles with more steps per cycle. If we use total number of smoothing steps as a proxy for total work, then the obvious thing to compare is two V-cycles using $k$ smoothing steps versus a single V-cycle with $2k$ smoothing steps. Using the bound for the simple iteration of <ref>, since

\[
\left(\frac{C}{C+2\omega k}\right)^2 \le \frac{C}{C+2\omega(2k)}\quad\text{for $k \ge 0$},
\]

it would seem to always be preferable to take fewer smoothing steps in the simple iteration, suggesting $k=1$ is optimal. In contrast, using the bound of <ref> for the fourth-kind Chebyshev iteration, we find that

\[
\frac{C}{C+\tfrac{4}{3}(2k)(4k+1)} \le \left(\frac{C}{C+\tfrac{4}{3}k(k+1)}\right)^2\quad\text{when $\tfrac{2}{3}(k+1)^2 \le C$,}
\]

suggesting $k$ should be chosen large enough so that $\tfrac{2}{3}(k/2+1)^2 \ge C$. Now, this simple analysis is flawed, especially because we are ignoring the fixed cost of the grid transfer operators (restriction and prolongation) in each cycle, which would tend to push the optimal $k$ higher. In particular, even for the simple iteration, the optimal $k$ may be greater than $1$. But it is striking that, for the fourth-kind Chebyshev iteration, the optimal $k$ grows with $C$, even ignoring the grid transfer costs.

§ NUMERICAL RESULTS

[][$C \approx 2$]

log ticks with fixed point,
legend style=font=,
legend pos=south west,
+[mark=+] table[x=k,y=w43] data-1.dat; $\omega=4/3$
+[mark=o] table[x=k,y=w32] data-1.dat; $\omega=3/2$
+[mark=*] table[x=k,y=cheb] data-1.dat; Chebyshev
+[mark=triangle*] table[x=k,y=opt] data-1.dat; optimized

cycle list shift=-4

+[ style=dashed, domain=1:10, samples=101,] 2/(2 + 2*(4/3)*x);
+[ style=dashed, domain=1:10, samples=101,] 2/(2 + 2*(3/2)*x);
+[ style=dashed, domain=1:10, samples=101,] 2/(2 + (4/3)*x*(x+1));
+[ style=dashed, domain=1:10, samples=101,] 2/(2 + ((2/3.1415926535897932384626)*(2*x+1))^2 - 2/3 );

[][$C \approx 8$]

log ticks with fixed point,
legend style=font=,
legend pos=south west,
+[mark=+] table[x=k,y=w43] data-2.dat; $\omega=4/3$
+[mark=o] table[x=k,y=w32] data-2.dat; $\omega=3/2$
+[mark=*] table[x=k,y=cheb] data-2.dat; Chebyshev
+[mark=triangle*] table[x=k,y=opt] data-2.dat; optimized

cycle list shift=-4

+[ style=dashed, domain=1:10, samples=101,] 8/(8 + 2*(4/3)*x);
+[ style=dashed, domain=1:10, samples=101,] 8/(8 + 2*(3/2)*x);
+[ style=dashed, domain=1:10, samples=101,] 8/(8 + (4/3)*x*(x+1));
+[ style=dashed, domain=1:10, samples=101,] 8/(8 + ((2/3.1415926535897932384626)*(2*x+1))^2 - 2/3 );

[][$C \approx 32$]

log ticks with fixed point,
legend style=font=,
legend pos=south west,
+[mark=+] table[x=k,y=w43] data-4.dat; $\omega=4/3$
+[mark=o] table[x=k,y=w32] data-4.dat; $\omega=3/2$
+[mark=*] table[x=k,y=cheb] data-4.dat; Chebyshev
+[mark=triangle*] table[x=k,y=opt] data-4.dat; optimized

cycle list shift=-4

+[ style=dashed, domain=1:10, samples=101,] 32/(32 + 2*(4/3)*x);
+[ style=dashed, domain=1:10, samples=101,] 32/(32 + 2*(3/2)*x);
+[ style=dashed, domain=1:10, samples=101,] 32/(32 + (4/3)*x*(x+1));
+[ style=dashed, domain=1:10, samples=101,] 32/(32 + ((2/3.1415926535897932384626)*(2*x+1))^2 - 2/3 );

[][$C \approx 128$]

ymin=0.1, ymax=1,
log ticks with fixed point,
legend style=font=,
legend pos=south west,
+[mark=+] table[x=k,y=w43] data-8.dat; $\omega=4/3$
+[mark=o] table[x=k,y=w32] data-8.dat; $\omega=3/2$
+[mark=*] table[x=k,y=cheb] data-8.dat; Chebyshev
+[mark=triangle*] table[x=k,y=opt] data-8.dat; optimized

cycle list shift=-4

+[ style=dashed, domain=1:10, samples=101,] 128/(128 + 2*(4/3)*x);
+[ style=dashed, domain=1:10, samples=101,] 128/(128 + 2*(3/2)*x);
+[ style=dashed, domain=1:10, samples=101,] 128/(128 + (4/3)*x*(x+1));
+[ style=dashed, domain=1:10, samples=101,] 128/(128 + ((2/3.1415926535897932384626)*(2*x+1))^2 - 2/3 );

V-Cycle error contraction factors

<Ref> illustrates the performance of the fourth-kind Chebyshev iteration (<ref>), the optimized iteration, and the simple iteration with damping parameters $\omega=4/3$ and $\omega=3/2$ as smoothers in a multigrid V-cycle. The damping parameters for the simple iterations were chosen so that they coincide with the two polynomial smoothers for $k=1$ smoothing step. The test problem is the 2D Poisson equation with Dirichlet boundary conditions discretized by the finite element method (FEM) on a $1023 \times 1023$ grid. For $B$ we chose the (rescaled) Jacobi iteration, which for this case is the same as the Richardson iteration as the diagonal of $A$ is constant. The approximation constant $C$ for this problem is approximately twice the element aspect ratio squared (as can be seen from a local Fourier analysis), and the four subfigures correspond to aspect ratios of 1, 2, 4, and 8, and thus to $C\approx2,8,32,128$. The solid lines with markers are the measured V-cycle contraction factors for different numbers of smoothing steps $k$. For each solid line, the dashed line of the same color is the bound derived from <ref>, i.e., <ref> for the simple iteration, <ref> for the Chebyshev iteration, and (<ref>) for the optimized iteration.

In all cases, the fourth-kind Chebyshev iteration significantly out-performs the simple iterations. While there is a significant gap between the measured contraction factors and the corresponding bounds, apart from the optimized iteration the measured contraction factors and bounds qualitatively follow the same trends, for the most part. The optimized iteration performs worse than even the simple iteration for the ideal case of aspect ratio $1$ for $k<5$. However, it does deliver the expected improvement over the fourth-kind Chebyshev iteration as suggested by the improved bound, for $C$ large enough and $k$ small enough.  While the performance of the optimized iteration can be relatively poor for small $C$ and/or large $k$, we observe that it approaches the bound derived from <ref> quite closely at some points, showing that the bound can be rather sharp in some cases.

Recall that the sharpened bound of <ref> for the fourth-kind Chebyshev iteration was asymptotically better than the bound for the optimized polynomials for $C < 3.66$. This helps to explain why the fourth-kind Chebyshev iteration outperforms the optimized iteration in the $C \approx 2$ case for all numbers of smoothing steps $k$. Since the polynomial iterations coincide with the simple iterations at $k=1$ smoothing step, this also explains why in the $C \approx 2, k=1$ case the lower damping factor $\omega=4/3$ outperforms the higher damping factor $\omega=3/2$ contrary to the bound of <ref>. The curves crossover near $k=2$, with the higher damping factor performing better for $k>2$, in line with <ref>.

[][$C \approx 2$]

log ticks with fixed point,
legend style=font=,
legend pos=south west,
+[mark=*] table[x=k,y=cheb] data-1.dat;
cycle list shift=-1

+[ style=dashed, domain=1:10, samples=101,] 2/(2 + (4/3)*x*(x+1));
+[ domain=1:10, samples=201,] 1/(1 + (4/3)*x*(x+1)/cprime2(2,x));
+[mark=+] table[x=k,y expr=1-1/a1] CN-cheb.dat;
+[ style=dashed, domain=1:10, samples=101,] 2/(2*x+1)^2;

[][$C \approx 8$]

log ticks with fixed point,
legend style=font=,
legend pos=south west,
+[mark=*] table[x=k,y=cheb] data-2.dat;
cycle list shift=-1

+[ style=dashed, domain=1:10, samples=101,] 8/(8 + (4/3)*x*(x+1));
+[ domain=1:10, samples=201,] 1/(1 + (4/3)*x*(x+1)/cprime2(8,x));
+[mark=+] table[x=k,y expr=1-1/a2] CN-cheb.dat;
+[ style=dashed, domain=1:10, samples=101,] 8/(2*x+1)^2;

[][$C \approx 32$]

log ticks with fixed point,
legend style=font=,
legend pos=south west,
+[mark=*] table[x=k,y=cheb] data-4.dat;
cycle list shift=-1

+[ style=dashed, domain=1:10, samples=101,] 32/(32 + (4/3)*x*(x+1));
+[ domain=1:10, samples=201,] 1/(1 + (4/3)*x*(x+1)/cprime2(32,x));
+[mark=+] table[x=k,y expr=1-1/a4] CN-cheb.dat;
+[ style=dashed, domain=1:10, samples=101,] 32/(2*x+1)^2;

[][$C \approx 128$]

ymax=1, ymin=0.1,
log ticks with fixed point,
legend style=font=,
legend pos=south west,
+[mark=*] table[x=k,y=cheb] data-8.dat; Chebyshev
cycle list shift=-1

+[ style=dashed, domain=1:10, samples=101,] 128/(128 + (4/3)*x*(x+1)); $C/(C+\tfrac{4}{3}k(k+1))$
+[ domain=1:10, samples=201,] 1/(1 + (4/3)*x*(x+1)/cprime2(128,x));
+[mark=+] table[x=k,y expr=1-1/a8] CN-cheb.dat; $1-C_N^{-1}$
+[ style=dashed, domain=1:10, samples=101,] 128/(2*x+1)^2; $C/(2k+1)^2$

V-cycle error contraction factor for Chebyshev iteration and theoretical bounds

The improvement of the sharpened bound of <ref> for the fourth-kind Chebyshev iteration over the bound of <ref> is illustrated in <ref>, which contains the same data as <ref>, but restricted to the Chebyshev iteration, and includes more bounds. The constant $C_N$ of the McCormick bound $1-C_N^{-1}$ (<ref>) was measured empirically. It can be seen that the improvement of the bound $C'/(C'+\tfrac{4}{3}k(k+1))$ of <ref> over $C/(C+\tfrac{4}{3}k(k+1))$ of <ref> is quite marginal, except in the case $C \approx 2$ where it is fairly significant. Indeed, in the $C \approx 2$ case, the improved bound is actually a fairly sharp estimate of the McCormick bound (<ref>) on which both <ref> and <ref> are based. Also illustrated is the two-level bound (<ref>), which, while not a bound of the multilevel V-cycle, nonetheless seems to give a good prediction of the V-cycle behavior for large enough $k$, at least on this test problem.

§ CONCLUSION

The fourth-kind Cheybshev iteration (<ref>) is a very simple way to accelerate the simple smoother iteration $(I-\omega BA)^k$ in a multigrid V-cycle for an SPD system, provided the smoother $B$ is also SPD. All that is required is an estimate of the largest eigenvalue of $BA$, $\rho(BA)$. No further parameters are involved. In particular, compared to the usual  (first-kind) Chebyshev semi-iteration, the fourth-kind iteration avoids the need to pick a minimum eigenvalue determining the range of eigenvalues to target. The iteration improves the full multilevel V-cycle error contraction factor bound from (the maximum over all levels of)
\[\frac{C}{C+2\omega k}\quad\text{to}\quad\frac{C}{C+\tfrac{4}{3}k(k+1)}\]

where $k$ is the number of smoothing steps and where $C$ is the approximation property involving the smoother $B$,

\[C = \norm{A^{-1} - P A_c^{-1} P^*}_{A, B}^2 = \sup_{u \in \ran(P)^{\perp_A} \setminus \{0\}} \frac{\norm{u}_{B^{-1}}^2}{\norm{u}_A^2} \ge 1, \]

assuming $B$ has been scaled so that $\rho(BA) = 1$. Note that if the largest eigenvalue of $BA$ (before scaling) is overestimated by some fraction, this effectively increases $C$ by the same fraction. E.g., overestimating $\rho(BA)$ by 1% effectively increases $C$ by 1%, thus also increasing the contraction factor bound by 1% in the limit as $k \to \infty$. The fourth-kind Chebyshev iteration is quasi-optimal for the multilevel bound. The optimal polynomials for the bound improve it by about 18% in the limit as $k \to \infty$. These can be implemented with a small modification to the fourth-kind Chebyshev iteration with the same run-time cost. The implementation is more complicated as the optimal polynomials need to be found numerically. The extra complication may be worthwhile, especially for large $C$, although the optimized iteration can actually perform worse, particularly when $C$ is small or $k$ is large.

Missing 'biblatex' package
The bibliography requires the 'biblatex' package.

abstractGauss-Seidel method is often the smoother of choice within multigrid applications. In the context of unstructured meshes, however, maintaining good parallel eﬃciency is diﬃcult with multiplicative iterative methods such as Gauss-Seidel. This leads us to consider alternative smoothers. We discuss the computational advantages of polynomial smoothers within parallel multigrid algorithms for positive deﬁnite symmetric systems. Two particular polynomials are considered: Chebyshev and a multilevel speciﬁc polynomial. The advantages of polynomial smoothing over traditional smoothers such as Gauss-Seidel are illustrated on several applications: Poisson’s equation, thin-body elasticity, and eddy current approximations to Maxwell’s equations. While parallelizing the Gauss-Seidel method typically involves a compromise between a scalable convergence rate and maintaining high ﬂop rates, polynomial smoothers achieve parallel scalable multigrid convergence rates without sacriﬁcing ﬂop rates. We show that, although parallel computers are the main motivation, polynomial smoothers are often surprisingly competitive with Gauss-Seidel smoothers on serial machines.
journaltitleJournal of Computational Physics
shortjournalJournal of Computational Physics
shorttitleParallel multigrid smoothing
titleParallel multigrid smoothing: polynomial versus Gauss-Seidel
dams et al. - 2003 - Parallel multigrid smoothing polynomial versus Ga.pdf:application/pdf

University of Colorado at Denver

abstractWe propose and analyze three multilevel iterative solvers of both domain decomposition and multigrid type. All of these methods are algebraic, allowing almost or fully black-box implementation. Their development was motivated by the need to solve large algebraic systems of equations resulting from finite element discretizations of self-adjoint, second order uniformly elliptic problems on unstructured three-dimensional meshes. Two of the methods discussed perform a simple, but effective domain decomposition as a part of the solving process. This allows for a remarkable adaptivity, where the decomposition is generated depending on the difficulty of the problem without requiring an input of a different decomposition. We focus on achieving robustness features that allow using the new methods as a replacement of direct solvers for solving these systems. The new methods are superior in terms of computational complexity and storage requirements. On serial architectures, the asymptotic computa...
titleRobust Iterative Methods on Unstructured Meshes

titleTables of properties of airfoil polynomials
typeNASA Reference Publication

booktitleMultigrid Methods
noteJournal Abbreviation: Multigrid Methods
titleMulti-grid convergence theory
ackbusch - 1982 - Multi-grid convergence theory.pdf:application/pdf

abstractMason, J.C., Chebyshev polynomials of the second, third and fourth kinds in approximation, indefinite integration, and integral transforms, Journal of Computational and Applied Mathematics 49 (1993) 169-178.
journaltitleJournal of Computational and Applied Mathematics
shortjournalJournal of Computational and Applied Mathematics
titleChebyshev polynomials of the second, third and fourth kinds in approximation, indefinite integration, and integral transforms
AF5LSED/Mason - 1993 - Chebyshev polynomials of the second, third and fou.pdf:application/pdf

CRC press

titleChebyshev polynomials

abstractIn two earlier papers [SIAM J. Numer. Anal., 19 (1982), pp. 924–929; 21 (1984), pp. 255–262], we developed an algebraic convergence theory for a class of multigrid methods applied to positive definite self-adjoint linear operator equations. The purpose of the present paper is to extend these results by eliminating an earlier approximation order restriction, developing additional rate estimates and allowing for very general relaxation schemes, including those that are nonstationary and nonsymmetric. These results apply to most well-known iterative methods and preconditioners.
journaltitleSIAM Journal on Numerical Analysis
notePublisher: Society for Industrial and Applied Mathematics
shortjournalSIAM J. Numer. Anal.
shorttitleMultigrid Methods for Variational Problems
titleMultigrid methods for variational problems: general theory for the V-cycle
CKBN/0722039.html:text/html;McCormick - 1985 - Multigrid Methods for Variational Problems Genera.pdf:/home/jlottes/Zotero/storage/FS6BYVE6/McCormick - 1985 - Multigrid Methods for Variational Problems Genera.pdf:application/pdf

abstractAn improved convergence bound for the polynomially accelerated two-level method of Brousek et al. [Electron. Trans. Numer. Anal., 44 (2015), pp. 401-442, Section 5] is proven. This method is a reinterpretation of the smoothed aggregation method with an aggressive coarsening and massive polynomial smoothing of Vaněk, Brezina, and Tezaur [SIAM J. Sci. Comput., 21 (1999), pp. 900-923], and its convergence rate estimate is improved here quantitatively. Next, since the symmetrization of the method requires two solutions of the coarse problem, a modification of the method is proposed that does not have this disadvantage, and a qualitatively better convergence result for the modification is established. In particular, it is shown that a bound of the convergence rate of the method with a multiply (k-times) smoothed prolongator is asymptotically inversely proportional to d2k, where d is the degree of the smoothing polynomial. In earlier works, this acceleration effect is only quadratic. Finally, for another modified multiply smoothed method, it is proved that this convergence improvement is not limited only to an asymptotic regime but holds true everywhere.
journaltitleETNA - Electronic Transactions on Numerical Analysis
shortjournalETNA - Electronic Transactions on Numerical Analysis
titleImproved convergence bounds for two-level methods with an aggressive coarsening and massive polynomial smoothing
ull Text:/home/jlottes/Zotero/storage/NGQW4YJF/Vanek et al. - 1996 - Two-level method with coarse space size independen.pdf:application/pdf

abstractThis monograph is the first to provide a comprehensive, self-contained and rigorous presentation of some of the most powerful preconditioning methods for solving finite element equations in a common block-matrix factorization framework. The book covers both algorithms and analysis using a common block-matrix factorization approach which emphasizes its unique feature. Topics covered include the classical incomplete block-factorization preconditioners, the most efficient methods such as the multigrid, algebraic multigrid, and domain decomposition. This text can serve as an indispensable reference for researchers, graduate students, and practitioners. It can also be used as a supplementary text for a topics course in preconditioning and/or multigrid methods at the graduate level.
noteGoogle-Books-ID: KuV0_LNndj8C
shorttitleMultilevel Block Factorization Preconditioners
titleMultilevel Block Factorization Preconditioners: Matrix-based Analysis and Algorithms for Solving Finite Element Equations
Mathematics / Counting & Numeration,Mathematics / Differential Equations / General,Mathematics / Numerical Analysis,Mathematics / Probability & Statistics / Stochastic Processes

§ COMPUTING OPTIMAL POLYNOMIALS

The optimal polynomials for <ref> can be found numerically by applying Newton's method to a system of nonlinear equations expressing the equioscillation property. Let us represent the polynomial $p$ by its $k$ roots $r_1, \dotsc, r_k$,

\[ p = \prod_{i=1}^k \left(1-\frac{x}{r_i}\right). \]

Note that $p$ satisfies the required condition $p(0) = 1$. We seek to minimize $\gamma = \sup_{0 < x \le 1} f^2$ where

\[ f \isdef w^{1/2} p, \quad w \isdef \frac{x}{1-p^2}, \quad x > 0, \]

and where we take $w(0) \isdef [-2p'(0)]^{-1}$ so that $w$ and $f$ are continuous at $x=0$.
For the optimal $p$, $f$ will equioscillate,

\begin{equation} \label{eq:f-equioscillate} f(0) = |f(x_1)| = \dotsb = |f(x_{k-1})| = |f(1)| \end{equation}

where $x_i$ are the $k-1$ extrema of $f$, $f'(x_i) = 0$, located between the $k$ roots $r_i$ of $p$. For convenience,

\[ x_k \isdef 1 \]

and let $F:\mathbb{R}^k \to \mathbb{R}^k$ be the function of the vector of roots $r = \begin{bmatrix} r_1 & \hdots & r_k \end{bmatrix}^T$ with components

\[ F_i = f(0) - |f(x_i)| . \]

Then (<ref>) is equivalent to $F(r)=0$ and the optimal $r$ can be found with Newton's method. A natural initial guess for $r$ is the vector of roots of the shifted Chebyshev polynomial of the fourth kind,

\[
r_i^{(0)} = \tfrac{1}{2}-\tfrac{1}{2}\cos \frac{i}{k+\tfrac{1}{2}} \pi, \quad i=1,\dotsc,k .
\]

A second, nested, Newton iteration can be used to compute the extrema $x_1, \dotsc, x_{k-1}$ of $f$. For the first outer iteration, a good initial guess for the extrema are those of the initial $p$,

\[
x_i^{(0)} = \tfrac{1}{2}-\tfrac{1}{2}\cos \frac{i+\tfrac{1}{2}}{k+\tfrac{1}{2}} \pi, \quad i=1,\dotsc,k-1 .
\]

Subsequent iterations can use the extrema from the previous iteration for the initial guess.

The total logarithmic derivative of $p$ is

\[ \frac{dp}{p} = \sum_{i=1}^k \frac{dx - \tfrac{x}{r_i} dr_i}{x - r_i} = q\,dx - \sum_{i=1}^k \frac{x}{x - r_i} \frac{ dr_i}{r_i}, \]

\[ q(x) \isdef \frac{p'(x)}{p(x)} = \sum_{i=1}^k \frac{1}{x-r_i}, \]

and the total logarithmic derivative of $f$ is

\begin{equation}\label{eq:f-total-log-deriv} \frac{df}{f} = \frac{dx}{2x} + \frac{w}{x} \frac{dp}{p} = (\tfrac{1}{2} + w q) \frac{dx}{x} - \sum_{i=1}^k \frac{w}{x-r_i} \frac{dr_i}{r_i} . \end{equation}

It follows that the extrema $x_i$ of $f$ are determined by the condition $\tfrac{1}{2} + w q = 0$, or equivalently by $g(x_i) = 0$ where

\[ g(x) \isdef \frac{x}{2 w(x)} + x q(x) = \frac{1-p(x)^2}{2} + x q(x) . \]

The negative derivative of $g(x)$ is

\[ -g'(x) = \sum_{i=1}^k \frac{1}{x-r_i} \left[ p(x)^2 + \frac{r_i}{x-r_i} \right], \]

giving the Newton step to update the extrema, $x_i \gets x_i - g(x_i) / g'(x_i)$.

Suppose $x$ is constrained to equal $x_i$, i.e., to be an extremum of $f$ when $i<k$ or else fixed at $1$ when $i=k$. In the first case $\tfrac{1}{2}+wq=0$ and in the second case $dx = 0$, so that in both cases (<ref>) reduces to

\[ \frac{df}{f} = - \sum_{j=1}^k \frac{w}{x-r_j} \frac{dr_j}{r_j}, \]

from which we can read off $\partial f(x_i)/\partial r_j$. Combining with $\partial f(0) / \partial r_j = f(0)^3 / r_j^2$, which can be derived from $f(0) = [-2q(0)]^{-1/2}$, we find that the components of the Jacobian $J$ of $F$ are

\[ J_{ij} = \frac{f(0)^3}{r_j^2} + \frac{w(x_i) |f(x_i)|}{r_j \, (x_i - r_j)}. \]

The Newton step for updating the vector of roots is then just $r \gets r + \Delta r$ where $J \Delta r = -F$. A full Matlab implementation of the algorithm is shown in <ref>. Empirically, the initial guesses suggested above are good enough for convergence of the pair of Newton iterations, at least up to $k = 200$.

style              = Matlab-editor,
basicstyle         = ,
escapechar         = ",
mlshowsectionrules = true,

A Matlab implementation of Newton's method for finding the roots of the optimal polynomials for <ref>

Once the roots of the optimal $p_k(\lambda)$ have been determined, we can compute the coefficients $\alpha_i$ of the expansion

\begin{equation}\label{eq:p-exp} p_k(\lambda) = \sum_{i=0}^k \alpha_i W_i(1-2\lambda) \end{equation}

of $p_k$ in the fourth-kind Chebyshev polynomial basis, so that

\[ \beta_0 = 1, \quad \beta_{i+1} = \beta_i - (2i+1) \alpha_i \]

then gives the coefficients $\beta_i$ for implementing $p_k$ in the iteration (<ref>). By combining the orthogonality relation for the fourth-kind Chebyshev polynomials <cit.>

\[ \frac{1}{\pi} \int_{-1}^1 \sqrt{\frac{1-x}{1+x}} W_i(x)W_j(x) \; dx = \delta_{ij} \]

with the associated Gaussian quadrature rule <cit.>,

\[ \frac{1}{\pi} \int_{-1}^1 \sqrt{\frac{1-x}{1+x}} p(x) \; dx = \sum_{i=1}^k \frac{1-x_i}{k+\tfrac{1}{2}} p(x_i), \quad x_i = \cos \frac{i}{k+\tfrac{1}{2}} \pi, \]

which holds for any polynomial $p(x)$ of degree at most $2k-1$, it follows that if (<ref>) holds then

\[ \alpha_j = \sum_{i=1}^k \frac{1-x_i}{k+\tfrac{1}{2}} W_j(x_i) p\left(\frac{1-x_i}{2}\right) \quad\text{for $j \le k-1$.}\]

Note that $\alpha_k$ is not needed to compute $\beta_k$. A Matlab implementation of the computation of the coefficients $\beta_i$ is shown in <ref>.

Matlab code for computing the coefficients $\beta_i$ for the optimized iteration (<ref>)

§ SHARPER ESTIMATES

It is possible to derive a slightly sharper bound than <ref>, particularly for the case of small $C$, e.g., $C = 2$. While providing marginally little improvement for modest to large $C$, the improved bound does predict some of the features illustrated in <ref>. In particular, at $C=2$, it predicts a lower error contraction factor for the fourth-kind Chebyshev iteration than the optimized iteration.

One of the steps of the proof of <ref> was to bound

\[
\sup_{u \in \ran(\pi_f) \setminus \{0\}} \frac{\norm{u}_{N^{-1}}^2 - \norm{u}_{A}^2}{\norm{u}_{B^{-1}}^2}
\]

by removing the restriction of $u$ to $\ran(\pi_f)$, the complement of the coarse space, and taking the supremum over arbitrary nonzero $u$ instead. This latter supremum was then bounded, effectively, by considering potential eigenvectors $u$ of $BA$, $BA u = \lambda u$ for some eigenvalue $\lambda \in (0, 1]$. In applying the bound to both the simple iteration and the Chebyshev iteration (<ref>), the supremum over potential eigenvalues was attained in the limit $\lambda \to 0$. But the eigenvectors corresponding to eigenvalues near $0$ are exactly those we expect to be “most excluded” from $\ran(\pi_f)$ in some sense, suggesting we can sharpen these bounds.

Let us make this intuition precise, as follows. Let $0 \le \mu < C^{-1}$ be an arbitrary parameter and let

\begin{align*}
\sV_l &= \spn \{ u \mid BA u = \lambda u, \lambda \le \mu \}, \\
\sV_h &= \spn \{ u \mid BA u = \lambda u, \lambda > \mu \},
\end{align*}

i.e., the spans of the eigenvectors with corresponding eigenvalues smaller than and larger than $\mu$. These spaces are $BA$-invariant, e.g., $BA \sV_l = \sV_l$, and give the decomposition

\[ \sV = \mathbb{R}^n = \sV_l \oplus \sV_h . \]

Let $\pi_l$ and $\pi_h = I - \pi_l$ be the projections corresponding to this decomposition. Because $BA$ is self adjoint in both the $A$- and $B^{-1}$-inner products, the $BA$-invariant spaces $\sV_l$ and $\sV_h$ are both $A$- and $B^{-1}$-orthogonal, so that we have the two Pythagorean theorems

\begin{equation}\label{eq:pythag-AB} \begin{aligned}
\norm{u}_A^2 &= \norm{\pi_l u}_A^2 + \norm{\pi_h u}_A^2, \\
\norm{u}_{B^{-1}}^2 &= \norm{\pi_l u}_{B^{-1}}^2 + \norm{\pi_h u}_{B^{-1}}^2 .
\end{aligned}\end{equation}

By construction, the Rayleigh quotients for the subspaces lie in the intervals $(0, \mu]$ and $(\mu, 1]$,

\begin{align*}
\norm{\pi_l u}_A^2 &\le \mu \norm{\pi_l u}_{B^{-1}}^2, \\
\mu \norm{\pi_h u}_{B^{-1}}^2 < \norm{\pi_h u}_A^2 &\le \norm{\pi_h u}_{B^{-1}}^2 .
\end{align*}

Since, by (<ref>), $\norm{u}_{B^{-1}}^2 \le C \norm{u}_A^2$ for any $u \in \ran(\pi_f)$, we see that

\[
\norm{\pi_l u}_{B^{-1}}^2 + \norm{\pi_h u}_{B^{-1}}^2 \le C (\mu \norm{\pi_l u}_{B^{-1}}^2 + \norm{\pi_h u}_{B^{-1}}^2),
\quad\text{for $u \in \ran(\pi_f)$},
\]

and thus

\begin{equation}\label{eq:ratio-bound} \frac{ \norm{\pi_l u}_{B^{-1}}^2 }{ \norm{\pi_h u}_{B^{-1}}^2 } \le \frac{C - 1}{1 - \mu C}, \quad\text{for $u \in \ran(\pi_f),$} \end{equation}

where we used the assumption $\mu < C^{-1}$.

The bound (<ref>) on the ratio of $B^{-1}$-norms of the components of an arbitrary $u \in \ran(\pi_f)$ in the subspaces $\sV_l$ and $\sV_h$ is a precise statement of how $C$ characterizes the extent to which the eigenspaces associated with smaller eigenvalues are “excluded” from  $\ran(\pi_f)$. We can use the bound to sharpen <ref>, by replacing the supremum of $\lambda p(\lambda)^2 / (1 - p(\lambda)^2)$ over $\lambda \in (0, 1]$ by a weighted average of this supremum with that taken over $\lambda \in (\mu, 1]$.

Let the smoother iteration (on each level) be given by

\[ G = p(BA) \]

where $B$ is symmetric positive definite and scaled such that $\rho(BA) = 1$, and $p(x)$ is a polynomial satisfying $p(0) = 1$ and $|p(x)| < 1$ for $0 < x \le 1$, possibly different on each level. Define

\[ \gamma_{\mu} \isdef \sup_{\mu < \lambda \le 1} \frac{\lambda \, p(\lambda)^2}{1 - p(\lambda)^2} . \]

Then the V-cycle contraction factor $\norm{E}_A^2$ is bounded by the largest value across all levels of

\[ \frac{C}{C+\gamma^{-1}} \]

where $C$ is the approximation property constant $\norm{\pi_f A^{-1}}_{A, B}^2$ and

\begin{equation}\label{eq:gamma-weighted-avg} \gamma = \frac{1 - C^{-1}}{1 - \mu} \gamma_0 + \frac{C^{-1} - \mu}{1 - \mu} \gamma_{\mu} \end{equation}

for any $0 \le \mu \le C^{-1}$, provided (<ref>) is taken to mean $\gamma = \gamma_{\mu}$ when $C=1$ particularly at the singularity at $\mu = 1$.

Let $N$ be the “symmetrized” smoother defined as in the proof of <ref>, so that $I - NA = G^2$ and

\[ N^{-1} = A (I - G^2)^{-1} . \]

We need to show

\[ \sup_{u \in \ran(\pi_f) \setminus \{0\}} \frac{\norm{u}_{N^{-1}}^2 - \norm{u}_{A}^2}{\norm{u}_{B^{-1}}^2} \le \gamma , \]

for then the result follows by a small modification to the proof of <ref>.

Note that

\begin{equation} \label{eq:gamma-bounds} \frac{\norm{\pi_l u}_{N^{-1}}^2 - \norm{\pi_l u}_A^2}{\norm{\pi_l u}_{B^{-1}}^2} \le \gamma_0,
\qquad \frac{\norm{\pi_h u}_{N^{-1}}^2 - \norm{\pi_h u}_A^2}{\norm{\pi_h u}_{B^{-1}}^2} \le \gamma_{\mu}
\end{equation}

follows similarly as in the proof of <ref>, here making use of the fact that the spectrum of $BA$ restricted to $\sV_h$ is a subset of $(\mu, 1]$ by construction.

When $C=1$, (<ref>) implies $\ran(\pi_f) \subseteq \sV_h$, and since (<ref>) reduces to $\gamma = \gamma_{\mu}$, the result is immediate. So we can assume $C > 1$. When $\mu = 0$ or $\mu = C^{-1}$, (<ref>) reduces to $\gamma = \gamma_0$, recovering <ref>, so we can assume $0 < \mu < C^{-1}$.

Because $(I - G^2)^{-1}$ leaves $\sV_l$ and $\sV_h$ invariant, and thus commutes with the projections $\pi_l$ and $\pi_h$, it follows that
$\sV_l$ and $\sV_h$ are orthogonal also in the $N^{-1}$-inner product, as, e.g.,

\[ \pi_l^* N^{-1} \pi_h = \pi_l^* A (I-G^2)^{-1} \pi_h = \pi_l^* A \pi_h (I-G^2)^{-1} = 0. \]

Hence, we also have a Pythagorean theorem for the $N^{-1}$ norm,

\begin{equation}\label{eq:pythag-N} \norm{u}_{N^{-1}}^2 = \norm{\pi_l u}_{N^{-1}}^2 + \norm{\pi_h u}_{N^{-1}}^2 . \end{equation}

Combining (<ref>), (<ref>), and (<ref>) yields

\[\begin{split}\frac{ \norm{u}_{N^{-1}}^2 - \norm{u}_A^2 }{ \norm{u}_{B^{-1}}^2}
\le  \frac{\gamma_0 \norm{\pi_l u}_{B^{-1}}^2 + \gamma_{\mu} \norm{\pi_h u}_{B^{-1}}^2}{ \norm{\pi_l u}_{B^{-1}}^2 + \norm{\pi_h u}_{B^{-1}}^2 }
= \gamma_{\mu} + \delta (\gamma_0 - \gamma_{\mu})
\end{split}\]

\[ \delta = \frac{\norm{\pi_l u}_{B^{-1}}^2}{\norm{\pi_l u}_{B^{-1}}^2 + \norm{\pi_h u}_{B^{-1}}^2} . \]

From (<ref>) it follows that

\[
\delta \le \frac{1 - C^{-1}}{1 - \mu}\quad\text{for $u \in \ran(\pi_f)$} .
\]

Hence, since $\gamma_0 \ge \gamma_{\mu}$,

\[
\sup_{u \in \ran(\pi_f) \setminus \{0\}} \frac{\norm{u}_{N^{-1}}^2 - \norm{u}_{A}^2}{\norm{u}_{B^{-1}}^2} \le
\gamma_u + \frac{1 - C^{-1}}{1 - \mu} (\gamma_0 - \gamma_{\mu}) = \gamma .
\]

§ SHARPER CHEBYSHEV ITERATION BOUND

Let us apply <ref> to the fourth-kind Chebyshev iteration by finding estimates of the optimal $\mu$ for the bound together with an estimate of

\[ \gamma_{\mu} \isdef \sup_{\mu < \lambda \le 1} f^2, \qquad f^2 = \frac{\lambda p_k^2}{1-p_k^2}, \]

\[ p_k = \frac{1}{2k+1} \frac{\sin (k+\tfrac{1}{2}) \theta}{\sin \tfrac{1}{2}\theta}, \quad 1-2\lambda = \cos \theta, \]

is the scaled and shifted Chebyshev polynomial of the fourth kind.

ymin=0, ymax=1,
xmin=0, xmax=1,

cycle multi list=Dark2,
trig format plots=rad,

+[ domain=0:pi, samples=201,] (lambda(x), ((2*3+1)^2-1)/3*fcheb(3, x)^2);
+[ domain=0:0.07914261619345, samples=2,] 1-8.2805671833*x;
cycle list shift=-1
+[ domain=0.07914261619345:1, samples=2,] 0.344655085885;

$\tfrac{n^2-1}{3} f(\lambda)^2$ and bound (<ref>) for $k=3$

It will be convenient to let $n \isdef 2k+1$, $\phi \isdef (k+\tfrac{1}{2})\theta$. Then, $\lambda = \sin^2 \tfrac{\phi}{n}$ and

\[ \frac{1}{n^2 f^2} = \csc^2 \phi - \frac{1}{n^2} \csc^2 \frac{\phi}{n} . \]

As illustrated in <ref>, $f^2$ attains its maximum of $3/(n^2-1)$ at $\lambda=0$ and has local maxima of decreasing magnitude, all again larger than $f^2$ at the endpoint $\lambda=1$. This follows from noticing that, for $\phi > \pi$, $f^2$ is not larger than the strictly decreasing envelope $(n^2 - \csc^2 \frac{\phi}{n})^{-1}$ which it contacts near each local maximum. The first (and largest) local maximum of $f^2$ occurs near $\phi = \tfrac{3}{2}\pi$, with the value $1/(n^2 w)$ where

\[ w \isdef \inf_{4<\phi\le\tfrac{3}{2}\pi} \csc^2 \phi - \frac{1}{n^2} \csc^2 \frac{\phi}{n} . \]

That is, we have $f^2 \le 1/(n^2 w)$ for $\phi > \pi$. In the region $\phi<\pi$, $f^2$ is again equal to $1/(n^2 w)$ near $\phi=2$. Let the specific value of $\phi$ where this occurs be denoted $\phi^*$. That is, define $\phi^*$ as the solution of

\[ \csc^2 \phi^* - \frac{1}{n^2} \csc^2 \frac{\phi^*}{n} = w \]

closest to $2$. Also, let $\lambda^* \isdef \sin^2 \tfrac{\phi^*}{n}$ be the corresponding $\lambda$. Because $f^2$ is convex as a function of $\lambda$ from $\lambda = 0$ up to the first zero at $\phi = \pi$, it follows that we can bound $\gamma_{\mu}$ as

\begin{equation}\label{eq:gamma-mu-bound} \gamma_{\mu} \le \tfrac{3}{n^2-1} \max \{ 1 - (1-y) \frac{\lambda}{\lambda^*}, y \}, \quad y \isdef \frac{n^2 -1}{3 n^2 w}. \end{equation}

The case $k=3$ ($n=7$) is illustrated in <ref>.

It is possible to find a series expansion for $w$ in powers of $1/n^2$, by using the series

\[ \csc^2 \phi - \frac{1}{n^2} \csc^2 \frac{\phi}{n}
= \csc^2 \phi - \frac{1}{\phi^2} - \frac{1}{3 n^2} - \frac{\phi^2}{15 n^4} + \dotsb, \]

expanding $\phi$ as a series, then minimizing order by order. We find

\[ w \approx 0.954876 - \frac{1}{3n^2} - \frac{1.47441}{n^4} - \frac{5.17589}{n^6} - \frac{16.1247}{n^8} + \dotsb . \]

The leading order term $w_0 = 0.95487649\ldots$ is the minimum of $\csc^2 \phi - \phi^{-2}$ near $\phi=\tfrac{3}{2}\pi$. Similarly, we can find an expansion for $\phi^*$,

\[ \phi^* \approx 1.99661 - \frac{0.898596}{n^4} - \frac{3.72312}{n^6} - \frac{13.1537}{n^8} + \dotsb . \]

These in turn lead to series expansions for $y$ and $\lambda^*$,

\begin{align*}
y &\approx 0.349085 - \frac{0.227225}{n^2} + \frac{0.459695}{n^4} + \frac{1.70183}{n^6} + \frac{5.96713}{n^8} + \dotsb, \\
\lambda^* &\approx \frac{3.98647}{n^2} - \frac{5.29731}{n^4} - \frac{0.772624}{n^6} - \frac{6.1326}{n^8} - \frac{19.6675}{n^{10}} + \dotsb .
\end{align*}

The leading term of $y$ here is $y_0 = 1/(3w_0)$.

Recall that <ref> bounds the V-cycle contraction factor by

\begin{equation}\label{eq:app-ml-poly-bound-sharp} \frac{C}{C+\gamma^{-1}}, \qquad \gamma = \frac{1 - C^{-1}}{1 - \mu} \gamma_0 + \frac{C^{-1} - \mu}{1 - \mu} \gamma_{\mu}. \end{equation}

If we let
\[ C' \isdef \frac{n^2-1}{3} \gamma C \]

then we can write the V-cycle contraction factor bound as

\[ \frac{C}{C+\gamma^{-1}} = \frac{C'}{C'+\tfrac{1}{3}(n^2-1)} = \frac{C'}{C'+\tfrac{4}{3}k(k+1)}. \]

Applying (<ref>) to (<ref>), we find

\[ C - C' \ge \frac{1 - C \mu}{1 - \mu} \frac{1-y}{\lambda^*} \min \{ \mu, \lambda^* \} . \]

We are free to pick any positive $\mu < C^{-1}$, and since the bound on $C - C'$ is maximized at

\[ \mu^* = \begin{cases} 1-\sqrt{1-C^{-1}}, & C^{-1} \le \lambda^* (2-\lambda^*), \\
\lambda^*, & \text{otherwise},
\end{cases} \]

which is always smaller than $C^{-1}$, we may pick $\mu = \mu^*$, giving the bound

\begin{equation}\label{eq:C'-exact}
C - C' \ge \begin{cases} \frac{1-y}{\lambda^*} g(C), & C^{-1} \le \lambda^* (2-\lambda^*), \\
\frac{1 - \lambda^* C}{1-\lambda^*}(1-y), & \text{otherwise,}
\end{cases}
\end{equation}

\[ g(C) \isdef 1+2(C-1)\left(1-\frac{1}{\sqrt{1-C^{-1}}}\right) \ge \frac{2 + 4 C}{16 C^2 - 5}. \]

Since $\lambda^* \to 0$ as $n \to \infty$, we have

\[ \lim_{n \to \infty} C-C' \ge \beta \isdef 1-y_0 = 1-(3w_0)^{-1} = 0.6509147135\ldots. \]

We can rewrite (<ref>) as

\begin{equation}\label{eq:C'-exact2}
C - C' \ge \beta \begin{cases} b - a C, & C < C^*, \\
c \; g(C), & \text{otherwise,}
\end{cases}
\end{equation}

\[ b = \frac{1}{\beta} \frac{1 - y}{1-\lambda^*} \approx 1 + \frac{4.33555}{n^2} + \frac{11.28}{n^4} +\dotsb
\ge 1 + \frac{65}{15n^2-33}, \]

\[ a = b \lambda^* \approx  \frac{3.98647}{n^2} + \frac{11.9862}{n^4} +\dotsb
\le \left(1 + \frac{65}{15n^2-33}\right)\frac{40}{10n^2-n+19}, \]

\[ c=\frac{1-y}{\lambda^*} \approx 0.250849 n^2 + 0.420901 + \dotsb
\ge \tfrac{1}{4}(n^2+1.79+n^{-2}), \]

\[ C^* = \frac{1}{\lambda^*(2-\lambda^*)} \approx 0.125424 n^2 + 0.416667 + \dotsb
\le \frac{n^2}{7.97} + 0.512 . \]

<Ref> now follows from (<ref>) by applying the estimates above. Note the reason we can replace $C^*$ with an upper bound is that $b - a C \le c \; g(C)$ holds for all $C \ge 1$ (with equality holding at $C = C^*$). These estimates were chosen in an ad-hoc way to balance simplicity with sharpness for all $n$. They can be verified to hold at leading order, implying they hold for $n$ large enough, as well as checked explicitly for small values of $n$.