1 Introduction↩︎

Nonlinear GMRES (NGMRES) was first proposed in 1997 by Washio and Oosterlee as a Krylov subspace based acceleration technique for solving difficult nonlinear equations [1]. Shortly after, they extended the approach to nonlinear multigrid and to solving recirculating incompressible flow problems [2]. There have been many works since then using NGMRES to accelerate or enable convergence in nonlinear solvers for many types of application problems, including image restoration [3], [4], tensor decomposition [5]–[7], and CUTEst problems [8]. NGMRES takes the following form: for a given fixed point operator \(q(x)\) used to solve \(g(x)=0\), depth \(m\) NGMRES can be written at iteration \(k+1\) as \[x_{k+1} = q(x_k) + \sum_{i=0}^{m_k}\beta_i^{k+1} \left(q(x_k)- x_{k-i} \right), \label{eq:min-NG1}\tag{1}\] where \(m_k=\min\{k,m\}\), and \(\beta^{k+1}=\big(\beta_0^{k+1},\beta_1^{k+1},\cdots, \beta_{m_k}^{k+1}\big)\) are determined by solving the least-squares optimization problem \[\label{eq:min-NG2} \min_{\beta^{k+1}} \| g(q(x_k))+\sum_{i=0}^{m_k} \beta_i^{k+1} \left(g(q(x_k))-g(x_{k-i}) \right)\|^2.\tag{2}\] NGMRES can equivalently be written in a constrained optimization form, which is more natural for our analysis herein: \[x_{k+1} = \alpha^{k+1}_{k+1} q(x_k) + \sum_{j=k-m_k}^k \alpha^{k+1}_i x_i,\] \[ where\alpha^{k+1} = \underset{\sum_{i=k-m_k}^{k+1} \alpha_i^{k+1} = 1}{\operatorname{arg}\,\operatorname{min}}\; \| \alpha^{k+1}_{k+1} g(q(x_k))+ \sum_{j=k-m_k}^k \alpha^{k+1}_i g(x_i) \|^2.\] While the norm historically used in the NGMRES optimization problems above is most often \(\ell^2\), the convergence theory in this paper requires it to be the range norm of \(g\). Numerical tests in this paper and [9] show that while sometimes similar convergence behavior is observed with these two norms, other times using the range norm of \(g\) produces significantly better results.

A convergence theory and proof of acceleration for NGMRES was recently established in [9] for NGMRES applied to the specific problem of solving the steady Navier-Stokes equations (NSE) with the Picard iteration, and under a small data assumption to enforce that Picard is contractive with rate \(\kappa_{nsepic}<1\) for both the fixed point residuals and the nonlinear residuals [9]. To our knowledge, this was the first work to mathematically prove that NGMRES accelerates convergence, albeit for a specific contractive iteration and nonlinear system.

The purpose of this paper is to generalize and extend the NSE/Picard-specific results to general nonlinear systems \(g=0\) where \(g:H\rightarrow H'\) (where \(H\) is a Hilbert space and \(H'\) is the dual space of \(H\)) with associated fixed point iteration function \(q:H\rightarrow H\). We provide herein a complete single step convergence analysis for NGMRES applied to FPIs for nonlinear systems that reveals how it accelerates and even enables convergence. Our theory covers both the contractive and noncontractive cases, as well as general depths \(m\ge 0\).

The results we find, under some reasonable smoothness assumptions on \(g\) and \(q\) (given in section 2) and denoting convergence ratio of the nonlinear residuals for the unaccelerated FPI to be \(\kappa\) (which the theory allows to be greater than 1), take the form: \[\begin{align} \| g(x_{k+1}) \|_{H'} & \le \min\{ \gamma_{k+1} \kappa,\theta_{k+1} \} \| g(x_{k}) \|_{H'} +higher order terms, \end{align}\] where \(\theta_{k+1}\) and \(\gamma_{k+1}\) are defined by \[\begin{align} \theta_{k+1} & := \frac{ \| \alpha^{k+1}_{k+1} g( q(x_k)) + \alpha^{k+1}_k g(x_k) + ... + \alpha^{k+1}_{k-m}g(x_{k-m}) \|_{H'} }{ \| g(x_k) \|_{H'} }, \\ \gamma_{k+1}& := \frac{ \| \alpha^{k+1}_{k+1} g(q(x_k) ) + \alpha^{k+1}_k g(x_k) + ... + \alpha^{k+1}_{k-m}g(x_{k-m}) \|_{H'} }{ \| g(q(x_k)) \|_{H'} }. \end{align}\] Note that both \(\theta_{k+1}\le 1\) and \(\gamma_{k+1}\le 1\) due to the structure of the optimization problem, and are only equal to one in the special case that the (constrained form) optimization problem returns \((0,1,0,...,0)\) or \((1,0,...,0)\), respectively. The only difference between \(\theta_{k+1}\) and \(\gamma_{k+1}\) is the denominator in their respective definitions. The term \(\gamma_{k+1}\) is observed to be the ‘gain of the optimization problem,’ as it is the ratio of what the optimization objective function gives to what a FPI with no acceleration gives. A key point is that \(\gamma_{k+1}\) scales the linear convergence rate \(\kappa\) of the nonlinear residuals to \(\gamma_{k+1}\kappa\), which identifies the gain of the optimization problem as the mechanism responsible for acceleration by NGMRES. Hence up to higher order terms, this result shows that NGMRES accelerates convergence precisely through the gain of the optimization problem. The term \(\theta_{k+1}\) in the convergence results directly approximates the linear convergence rate and moreover, since \(\theta_{k+1}\le 1\) (and only equals 1 in very special cases), it establishes convergence of NGMRES, up to the effect of higher order terms. Numerical tests in Section 4 show that \(\theta_{k+1}\) is a remarkably accurate predictor of the convergence rate when the higher order terms are negligible, and thus it can be used to determine when the nonlinear contributions are negligible and thus determine an adaptive depth strategy to speed up convergence even more.

To our knowledge, these are the first theoretical results for NGMRES that prove acceleration and identify the mechanism responsible for it, for both contractive and noncontractive iterations, for any depth, and for general nonlinear systems. Existing theory for NGMRES is all quite recent, and include important results from [10] where it is established the NGMRES converges \(q\)-linearly (and \(r\)-linearly when \(m=0\)), that full depth NGMRES applied to the Richardson iteration is equivalent to classical GMRES [11], an extension of this result to preconditioned Richardson and preconditioned GMRES in [12], and convergence analysis for alternating NGMRES is given for linear systems in [13]. Hence we believe the results herein establishing NGMRES acceleration for general nonlinear systems and associated contractive and noncontractive iterations helps to fill a gap in NGMRES theory.

NGMRES is related to the recently popularized Anderson acceleration (AA) [14], [15] method in that it has a similar extrapolation that defines the next iterate, but it significantly differs in that the optimization problem of NGMRES uses nonlinear residuals while AA uses fixed point residuals. Although for many problems AA and NGMRES provide similar improvement to convergence behavior [7], [16], NGMRES can outperform AA in highly anisotropic problems, steep optimization landscapes, and systems with sharp discontinuities [5], [7], [17]. Still, comparisons in the literature are scarce and it is not yet possible to draw any strong conclusions. It is an open question to classify what types of applications and iterations will have one of AA or NGMRES significantly outperform the other. We partially address this question herein, by showing that NGMRES is better suited for use with superlinearly converging iterations.

We note that convergence results for AA applied to general FPIs were established in [15], [18] that are similar to NGMRES results herein as they show the gain of the AA optimization problem is the mechanism responsible for the acceleration provided by AA, and AA adds higher order terms to the fixed point residual expansion. A key difference compared to NGMRES results is that AA theory is fixed-point residual based, while NMGRES theory is based on the true nonlinear residuals. Another key difference in our NGMRES theory is that for equal depths, higher order terms for AA for lower order than for NGMRES; in particular for depth 1, NGMRES has quadratic higher order terms while AA’s higher order terms are of order 1.62 [19]. Hence NGMRES appears better suited for use with superlinearly converging methods, and we give a numerical test that illustrates this. Yet another difference in AA and NGMRES theories is the simplicity of the NGMRES proofs: while the proofs that establish AA theory are rather long and technical [15], [18]–[22], the NGMRES convergence and acceleration theory herein is much shorter and straight-forward.

For the choice of the optimization norm, we find analogous results as found for AA in [23]: using the correct norm in the optimization problem matters. For many years, \(\ell^2\) was the most common choice for AA and NMGRES optimization norms since it was efficient and simple to use, and there was no theory yet that suggested to use a different norm. However, this is no longer true. The NMGRES theory herein requires the range norm of \(g\) be used for the optimization problem, and tests herein and in [9] show that using \(\ell^2\) instead can give worse (or no) convergence. AA theory relies on the Hilbert space norm associated with \(q\) being used for the AA optimization problem, and in [23] it is shown that using this norm instead of \(\ell^2\) can provide better convergence results.

This paper is arranged as follows. In section 2 we give some mathematical preliminaries to set the notation, present assumptions, and recall known results to allow for a smooth analysis to follow. Section 3 gives a convergence analysis for NGMRES for contractive and noncontractive fixed point iterations, and also shows how NGMRES provides acceleration. Section 4 gives results of several numerical tests which illustrate the new theory and show the effectiveness of the method.

2 Mathematical Preliminaries↩︎

We now present notation and mathematical preliminaries to be used in the sections to follow. Let \((H,\|\cdot \|_H)\) be a Hilbert space with inner product \((\cdot,\cdot)\), and \((H', \|\cdot \|_{H'})\) be the dual space of \(H\) with \(\displaystyle \|y\|_{H'} = \sup\limits_{0\neq x\in H} \frac{\langle y,x\rangle}{\|x\|}\), where \(\langle \cdot,\cdot \rangle\) represents the duality pairing.

This paper is concerned with solving nonlinear systems of the form \(g(x)=0\) for \(g:H\rightarrow H'\) and \(f\in H'\), where \(x^*\in H\) denotes a root. The associated fixed point function is denoted by \(q:H\rightarrow H\), and the goal of NGMRES is to accelerate or enable convergence of the FPI \(x_{k+1}=q(x_k)\). We make the following regularity assumptions on \(g\) and \(q\), and set notation for the associated constants. Together, these assumptions amount to requiring \(g\) be Lipschitz continuously differentiable on \(H\), and \(q\) to be continuous and satisfy that fixed point residuals cannot grow beyond the nonlinear residuals without bound. Despite the weakened assumptions on \(q\) relative to \(g\), many fixed point functions \(q\) are derived from \(g\) (e.g. \(q(x)=g(x)+x\)) so the assumed regularity on \(g\) is inherited by \(q\). A key feature of our analysis is that it covers the contractive and noncontractive cases, and note that we make no assumptions on \(\kappa_0\) and \(\kappa_1\) other than they are finite.

Note that since \(g\) is Fréchet differentiable, \(g'\) is linear with respect to the second argument. Consequently we have \[\begin{align} \label{dglinear} g'(\cdot; t_1h_1+ t_2h_2) = t_1g'(\cdot ; h_1) + t_2 g'(\cdot, h_2),\;\; \forall t_1,t_2\in \mathbb{R}, h_1,h_2\in H. \end{align}\tag{3}\]

We now recall the depth \(m\) NGMRES algorithm for accelerating a fixed-point iteration defined by \(q(x)\) that is used to solve \(g(x)=0\). Here we use the constrained form of the optimization problem, which is equivalent to the more commonly used formulation 1 2 .

The depth \(m\) can vary at each step, and \(m\) can be no larger than \(k\). While it is common to use a constant depth \(m\) after the initial \(m\) iterations, we discuss an adaptive strategy for dynamically picking \(m\) in Section 4.

While historically the \(\ell^2\) norm is commonly used for the optimization norm, our results herein rely on using the \(H'\) norm. Moreover, in [9] it is shown that for \(g\) representing the NSE and \(q\) the associated Picard iteration, using the dual norm for the optimization problem gives much better results compared to using \(\ell^2\) for certain problems, in particular in 3D.

We will assume that the \(\{ \alpha_j^k \}\) coefficients are uniformly bounded, which is implied by the optimization problems being (uniformly numerically) well-posed. Note that this can be checked on the fly, and \(m\) can be reduced until a well-posed optimization problem is obtained.

We define the following important quantities, which arise in our convergence analysis. \[\begin{align} \tag{4} \theta^m_{k+1}\mathrel{\vcenter{:}}= \frac{\left\|\alpha_{k+1}^{k+1} g(q(x_k)) + \sum\limits_{j=k-m}^k \alpha^{k+1}_j g(x_j) \right\|_{H'} }{\|g(x_k)\|_{H'}}, \\ \tag{5} \gamma^m_{k+1} \mathrel{\vcenter{:}}= \frac{\left\|\alpha_{k+1}^{k+1} g(q(x_k)) + \sum\limits_{j=k-m}^k \alpha^{k+1}_j g(x_j) \right\|_{H'}}{\|g(q(x_k))\|_{H'}} \end{align}\] From ?? , it holds that \(0\le \gamma^m_{k+1},\theta_{k+1} \le 1\). Observe that the numerator of \(\gamma^m_{k+1}\) is the optimization objective function evaluated at the minimizer \(\alpha^{k+1}=\{ \alpha^{k+1}_i \}_{i=k-m}^{k+1}\), and the denominator is the optimization objective function evaluated at \((1,0,...,0)\). Hence \(\gamma^m_{k+1}=1\) only when the optimization problem returns \(\alpha^{k+1}=(1,0,...,0)\) as optimal, and so except in these very rare instances it holds that \(0\le \gamma^m_{k+1}<1\) (and our numerical tests illustrate this). Similarly, \(\theta^m_{k+1} = 1\) only if \(\alpha^{k+1}=(0,1,0,...,0)\). Thus, we can generally expect \(\gamma^m_{k+1},\theta_{k+1} < 1\) except in diabolical cases. Note that if one does not perform NGMRES and simply does a fixed point iteration, this is equivalent to running NGMRES with optimization coefficient \(\alpha^{k+1}=(1,0,...,0)\). Hence we refer to \(\gamma^m_{k+1}\) as the ratio ‘gain of the optimization problem’ at Step \(k+1\).

3 Convergence↩︎

We now present a new convergence theory for NGMRES. This is a single step analysis, and we set \(m=\min \{m_{max},k \}\). We first consider the base case of \(m=0\), and prove a convergence result that holds for contractive and noncontractive fixed point iterations. We then prove an analogous result for general \(m\ge 1\) in the contractive case, and finally a result for \(m\ge 1\) in the noncontractive case. In each case the main result is that NGMRES scales the linear convergence rate of the nonlinear residuals (i.e. \(\kappa_1\)) by the gain of the optimization problem \(\gamma^m_{k+1}\), and \(\theta^m_{k+1} \le 1\) directly predicts the NGMRES nonlinear residuals’ linear convergence rate (up to higher order terms).

3.1 NGMRES convergence theory for \(m=0\)↩︎

Consider Algorithm 2 with \(m=0\), which reduces to the following algorithm.

3.2 NGMRES convergence theory for general \(m\): the contractive case↩︎

We now extend the NGMRES convergence and acceleration results to the case of general \(m\). Here, we split the analysis into two cases: contractive and noncontractive. While for \(m=0\) our results hold for both the contractive and noncontractive cases, for \(m\ge 1\) the structure of NGMRES leads us to different analyses for these cases. We now give results for the contractive case, after a preliminary lemma.

3.3 NGMRES convergence theory for general \(m\): the noncontractive case↩︎

Now we consider the general depth case in the noncontractive case. Here, we use a different preliminary lemma.

4 Numerical Experiments↩︎

We now give results for four numerical tests that illustrate the theory above for NGMRES. These tests show the effectiveness of NGMRES in accelerating and enabling convergence, show that \(\theta_k^m\) is a very accurate approximation of the nonlinear residuals’ linear convergence rate and how this can be used to develop an adaptive depth strategy, how restarts can significantly improve convergence in the noncontractive case, and that better convergence can be obtained if the optimization norm matches that of the theory. We also show how NGMRES is well suited to use with known deflation techniques to directly find distinct solutions for multi-solution PDEs, and that NGMRES is effective for superlinearly convergent iterations. Unless otherwise mentioned, \(H'\) is used as the optimization norm in our tests.

4.1 Nonlinear Helmholtz equation↩︎

The nonlinear Helmholtz (NLH) equation arises in nonlinear optics as a model for propagation of linearly-polarized, time-harmonic electromagnetic waves in Kerr-type dielectrics. i.e. “lasers.” Following [24], [25], the system can be written using two-way boundary conditions as: Find \(u:[0,1]\rightarrow \mathbb{C}\) satisfying \[\begin{align} \frac{d^2 u}{dx^2} + k_0^2 \left( 1 + \epsilon(x) |u|^2 \right) u & = 0, \;\;\;0<x<1, \\ \frac{du}{dx} + i k_0 u & = 2 i k_0,\;\;\;x=0,\\ \frac{du}{dx} - i k_0 u & = 0,\;\;\;x=1. \end{align}\] \(k_0\) represents the linear wave number in the surrounding medium, and we take \(\epsilon(x)\) on \([0,1]\) to be a piecewise constant function tho approximate a grated Kerr medium by \[\epsilon(x) = \left\{ \begin{array}{lc} 0 & 0\le x<\frac{1}{3} \\ 0.5 & \frac{1}{3}\le x < \frac{2}{3} \\ 1 & \frac{2}{3} \le x\le 1 \end{array} \right. .\] Solutions for \(k_0\)=20, 40 and 60 (found using our methods below) are shown in Figure 1.

Even in 1D, NLH is a very challenging problem, especially for larger values of \(k_0\) since this increases the strength of the cubic nonlinearity [24]–[26]. NLH has recently become a benchmark for testing of nonlinear solvers and acceleration methods, e.g. [15], [18], [27].

We consider the Picard iteration for solving NLH, which is given by \[\begin{align} \frac{d^2 u_{j+1}}{dx^2} + k_0^2 \left( 1 + \epsilon(x) |u_j|^2 \right) u_{j+1} & = 0, \;\;\;0<x<1, \tag{8} \\ \frac{du_{j+1}}{dx} + i k_0 u_{j+1} & = 2 i k_0,\;\;\;x=0, \tag{9} \\ \frac{du_{j+1}}{dx} - i k_0 u_{j+1} & = 0,\;\;\;x=1. \tag{10} \end{align}\] The Picard iteration 8 10 is a fixed point iteration, since we can write \(u_{j+1}=q(u_j)\) where \(q\) is the solution operator of 8 10 . Following [24], we take \(u_0=e^{ik_0 x}\). Note that due to the second derivative in NLH (and the Picard iteration), the natural function space is \(H = H^1(\Omega)\) with two way boundary conditions. The Picard system is discretized using second order finite differences (which is equivalent to using \(P_1\) finite elements in this setting) and uniform point spacing of \(h=0.002\) using to the iteration. We approximate the discrete dual norm for the NGMRES optimization problem using \[\| \phi \|_{H'} \approx (\phi, S^{-1} \phi)^{1/2},\] where \(S\) is the finite difference stiffness matrix with endpoint degrees of freedom removed. We note with this discretization and problem setup, the Picard iteration converges for \(k_0 \le 17\), and diverges for \(k_0>18\) (tests omitted). Linear solves are performed with a direct solver for this problem.

4.1.1 Illustration of NGMRES convergence theory↩︎

For our first test of NGMRES using NLH and its associated Picard iteration, we consider \(k_0=20\) and run NGMRES with \(m\)=2, 5, 10 and 20, as well as the Picard iteration with no acceleration. Results are shown in Figure 2 for convergence, \(\gamma_k\) and \(\theta_k\). For convergence, we observe that unaccelerated Picard fails to converge, while all NMGRES-Picard tests converge. Convergence improves as \(m\) is increased, except \(m=10\) and \(m=20\) perform about the same. Acceleration coefficients are shown in Figure 2 at center, and as expected we see that overall, the acceleration coefficients get smaller as \(m\) gets larger. Shown at right in Figure 2 are plots of the actual convergence rate \(\frac{\| g(u_{k}) \|_{H'}}{ \| g(u_{k-1}) \|_{H'}}\) and predicted linear convergence rate \(\theta_k\) for \(m=5\). We observe that after about the tenth iteration, \(\theta_k\) is a very good predictor of the actual convergence rate; once the nonlinear residual is sufficiently small, the higher order terms in the convergence estimate become negligible and \(\theta_k\) becomes a very good predictor of the convergence rate at iteration \(k\).

4.1.2 Improving NGMRES with restarts↩︎

Theorem 10 shows that in the noncontractive case, the nonlinear residual bound higher order term sum goes back to the initial iterate. If sharp, this suggests that gains made on the linear terms can be dominated by the higher order terms as more iterations are performed. Hence restarts can be a way to mitigate this effect. To test this, we pick \(k_0=40\) and repeat the same NLH test as above with varying \(m\). Results are shown in Figure 3, and we see that while each of \(m\)=5, 10, 20 and 50 converge (note \(m<5\) failed; plots omitted), convergence is quite erratic when no restarts are used. Using restarts every 50 iterations improves convergence for each choice of \(m\), and also yields somewhat smoother convergence behavior. Even more improvement is found when restarts are used every 25 iterations.

4.1.3 Dual norm versus \(\ell^2\) for the optimization norm↩︎

Our analysis requires the use of the dual norm \(H'\) for the optimization problem, however historically NGMRES users choose the \(\ell^2\) norm for convenience and simplicity (but in their defense, until now there was no general convergence theory that suggested using a different norm than \(\ell^2\)).

We test NGMRES-Picard for NLH using the \(H'\) and \(\ell^2\) optimization norms, under three parameter sets: \([k_0\)=40, \(m\)=10, restart=25]; \([k_0\)=60, \(m\)=10, restart=25]; \([k_0\)=60, \(m\)=50, restart=none]. Results are shown in Figure 4, and we observe that for the easier case of \(k_0=40\), results are similar for the \(H'\) and \(\ell^2\) optimization norm choices. However, for the harder case of \(k_0\)=60, the results using the \(H'\) norm are much better than those using \(\ell^2\); for \(m\)=10 and restart=25, convergence is achieved in 475 iterations using \(H'\) optimization norm while convergence is never reached using \(\ell^2\) (in 2000 iterations, plot only shows up to 700 for comparison purposes). For \(m\)=50 and no restarts, using optimization norm \(H'\) provides convergence in about 450 iteration while it takes 625 for convergence using \(\ell^2\).

These results are similar to what was found in [9] for NGMRES applied to the Picard iteration for NSE. That, for easier problems, using the dual norm for the optimization problem or using \(\ell^2\) gives similar convergence results. However for harder problems, especially in 3D, using the dual norm gives much better convergence results.

4.2 3D Buoyancy driven flow (Boussinesq equations)↩︎

We next consider tests on flows driven by natural convection or buoyancy. These flows occur in applications such as ventilation, solar collectors, window insulation, cooling in electronics, and many others [28], and are modeled by the Boussinesq system: in \(\Omega\subset\mathbb{R}^d\) (\(d\)=2 or 3) by \[\begin{align} \label{eqn:bouss} u_t + (u\cdot \nabla ) u - \nu \Delta u + \nabla p &= Ri [ 0; \theta ] + f, \nonumber \\ \nabla \cdot u &= 0, \nonumber \\ T_t + (u\cdot \nabla) T - \kappa \Delta T &= F, \end{align}\tag{11}\] with \(u\) representing the velocity field, \(T\) the temperature (or density), \(p\) the pressure, and \(f\) and \(F\) the external forces. The kinematic viscosity \(\nu>0\) is defined to be the inverse of the Reynolds number (\(Re=\nu^{-1}\)), and the thermal conductivity \(\kappa\) is given by \(\kappa=Re^{-1}Pr^{-1}\), with \(Pr\) representing the Prandtl number and \(Ri\) the Richardson number accounting for the gravitational force. The associated Rayleigh number is defined by \[Ra=Ri\cdot Re^2 \cdot Pr.\] Larger \(Ra\) leads to more complex physical behavior as well as more difficulties in numerically solving the equations.

For a test example, we use the 3D differentially heated cavity benchmark problem, where \(\Omega=(0,1)^3\), there is no external forcing (\(f=0\), \(F=0\)), no slip velocity boundary conditions are enforced on all sides, and for temperature we strongly enforce \(T(0,y,z)=1\), \(T(1,y,z)=0\) and perfect insulation \((\nabla T\cdot n=0)\) is weakly enforced on the other four sides. We set \(Ri=1\), \(Pr=0.71\) to model air, and then vary \(Re\) in order to vary \(Ra\).

We solve directly for steady solutions (i.e. \(u_t=0\) and \(T_t=0\)) using the Picard iteration: Given \(u_0\), define \(u_{k},T_{k}\), \(k\ge 1\), to be solutions of the (decoupled) system \[\begin{align} (u_{k-1}\cdot \nabla ) u_k - \nu \Delta u_k + \nabla p_k & = Ri [ 0;T_{k} ], \tag{12} \\ \nabla \cdot u_k &= 0, \tag{13} \\ (u_{k-1}\cdot \nabla) T_k - \kappa \Delta T_k & = 0, \tag{14} \end{align}\] and also satisfy the boundary conditions discussed above.

Note that 14 can be decoupled and solved to obtain \(T_k\), and then an Oseen solve of 12 13 can be performed for \(u_k\) and \(p_k\). We note that for Boussinesq, one iteration of Picard is much cheaper than one Newton iteration, which would include the terms \(+(u_k-u_{k-1})\cdot\nabla u_{k-1}\) on the left hand side of the momentum equation and \(+(u_k-u_{k-1})\cdot\nabla T_{k-1}\) on the left hand side of the energy equation. Hence, the Newton linear systems would be fully coupled and much more difficult to efficiently and reliably solve. Thus, accelerating the Picard iteration for Boussinesq is highly desirable.

We discretize velocity-pressure using \((P_3,P_2^{disc})\) Scott-Vogelius (SV) elements on a tetrahedral grid constructed with a rectangular box grid generated by \(12^3\) Chebyshev points that is then refined by splitting each box into 6 tetrahedra and finally barycenter refining each of those tetrahedra. This SV pair is known to be inf-sup stable on such a grid [29]. For the temperature, we use \(P_3\) on this same mesh, and in total the discretization yields 954K degrees of freedom (dof). The temperature linear solves are done directly, and for the Oseen solves we use an augmented Lagrangian type preconditioner for GMRES following [30]–[32], with direct inner solves.

The natural function spaces of velocity and temperature are (after equivalently switching to homogeneous Dirichlet boundary conditions for temperature on the walls and adjusting the temperature equation forcing appropriately), \[V = \{ v\in H^1_0(\Omega)^3, \;\nabla \cdot v=0\},\;\;\; X = \{ v\in H^1(\Omega), v(x,y,z)=0forx=0,1\}.\] We denote \(q_B: V\times X \rightarrow V\times X\) as the solution operator of 12 14 and thus also the fixed point function: \([u_{k+1},T_{k+1}]=q_B([u_k,T_k])\). Note the pressure and pressure space is removed since the SV pair is inf-sup stable and thus the divergence-free velocity subspace can be used. We denote \(g_B: V\times X \rightarrow V'\times X'\) to be the nonlinear residual operator.

We approximate the velocity and temperature dual space norms for the NGMRES optimization using \[\| \phi \|_{V'} \approx (\phi,A_h^{-1} \phi)^{1/2},\;\;\| \chi \|_{X'} \approx (\chi,S_h^{-1} \chi)^{1/2},\] where \(A_h\) is the discrete Stokes operator and \(S_h\) is the stiffness matrix, each with its associated Dirichlet boundary dof removed. Hence the optimization dual norm used is \[\| (\phi,\chi) \|_{B'} := \left( \| \nabla \phi \|_{V'}^2 + \| \nabla \chi \|_{X'}^2 \right)^{1/2}.\] This norm is also used for measuring convergence in our tests. We note that the Picard iteration converges for \(Ra\)=10000 in 85 iterations, and for \(Ra\ge 25000\) it does not converge in 200 iterations (and shows no indication that it will ever converge).

4.2.1 Illustration of NGMRES convergence theory↩︎

For our first test with the Boussinesq system, we compute NGMRES-Picard with varying \(m\) and \(Ra\)=10000 and 100000. Results are shown in Figure 6. For \(Ra\)=10000, we observe that the unaccelerated Picard iteration converges in about 85 iterations, while NGMRES-Picard converges in about 15 iterations for each of \(m=1,5,10\). Acceleration coefficients \(\gamma_k\) for \(Ra\)=10000 are as expected in that we see an overall decrease as \(m\) increases. We also observe that \(\theta_k\) becomes an excellent predictor for the convergence rate at each step, after about the 7th iteration.

For \(Ra\)=100000, unaccelerated Picard fails to converge. However, NGMRES-Picard converges for each of \(m=3,5,10\), and their respective convergence behaviors are similar. Worse converge is observed for \(m=20\), with convergence in 90 iterations, while for \(m=50\) there is no convergence at all. The plot of \(\gamma_k\)’s shows an overall decrease as \(m\) increases, but as our theory shows, this only affects the linear convergence rate and for large \(m\) the higher order terms can play a significant role in slowing/preventing convergence. We also observe that for \(m=10\), \(\theta_k\) agrees well with the actual convergence ratio by about iteration 50.

4.2.2 Using \(\theta_k\) for an adaptive depth strategy↩︎

As observed in our theory and in the numerical tests above, it is often best to use smaller \(m\) early in the iteration when the nonlinear residual is large and then larger \(m\) once the residual is small. The key idea is that increasing \(m\) improves the linear convergence rate, but at the expense of higher order terms having a greater effect. Hence once the nonlinear residuals are moderately small, then the higher order terms become negligible and so increasing \(m\) will only help.

Based on our theory above in the contractive regime, we can use \(\theta_k\)’s sharp prediction of the linear convergence rate to better choose when to increase \(m\). Following along the lines we proposed in [16], we increase \(m\) in the following way. First, a small \(m\) is chosen and NGMRES is run as usual. At each step, we check the inequality \[\bigg| \theta_k - \frac{ \| g_B(u_k,T_k)\|_{B'}}{ \| g_B(u_{k-1},T_{k-1})\|_{B'}}\bigg| <0.001.\] If this is satisfied, it strongly suggests that the higher order terms are not affecting the convergence rate, and so we increase \(m\) by one at the next iteration.

Convergence results for this strategy are shown in Figure 7 for \(Ra\)=100000. We observe that by starting with \(m=3\), the adaptive strategy leads to convergence in 53 iterations, and starting with \(m=5\) gives convergence in 58 iterations. Convergence with constant \(m\)=3,5,10 are also shown for comparison (these are the same as in Figure 6 at bottom left), and we observe that the simple adaptive strategy gives a significant speedup in convergence.

4.3 Incompressible Navier-Stokes equations↩︎

For our third test problem, we consider the incompressible NSE and the channel expansion problem studied by Farrell et al [33] which as domain \(\Omega=(0,2.5)\times (-1,1) \cup (2.5,150)\times(-6,6)\) representing a narrow channel of width 2 entering into a channel of width 12. The NSE take the form: \[\label{NS1} \left\{\begin{align} -\nu \Delta u+u\cdot\nabla u+ \nabla p&=0 \quad \text{in}~\Omega,\\ \nabla\cdot {u}&=0\quad \text{in}~\Omega, \end{align}\right.\tag{15}\] where \(u\) and \(p\) represent the unknown velocity and pressure, respectively, and \(\nu>0\) is the kinematic viscosity which for this test is the inverse of the Reynolds number \(Re = \nu^{-1}.\). Dirichlet boundary conditions are enforced as parabolic inflow with max inlet velocity of 1 at \(x=0\) and as no-slip on the walls, and zero traction is weakly enforced at the outflow. We choose \(Re\)=50 for this test, for which the NSE is known to have 5 distinct solutions [33], although two of these are symmetric reflections of others (with one solution being self-reflectionally symmetric). The 5 distinct solutions are shown in Figure 8.

State of the art techniques for finding multiple solutions typically use variations of (deflated) Newton-type continuation methods [33], and so one must start with very low \(Re\) and find all solutions there before ‘climbing the ladder’ to slightly higher \(Re\), repeating the process until all solutions are found at the desired \(Re\). While generally effective, this process can be slow since many solves need done at lower \(Re\). Hence we aim to find solutions directly at \(Re\)=50 using NGMRES-Picard, and avoid continuation methods altogether.

The Picard iteration for the NSE is given by \[\label{Pic1} \left\{\begin{align} -\nu \Delta u_{k+1}+u_k\cdot\nabla u_{k+1}+ \nabla p_{k+1}&=0, \\ \nabla\cdot {u}_{k+1}&=0, \end{align}\right.\tag{16}\] together with the boundary conditions given above. We now use NGMRES-Picard together with deflation and to directly find all five solutions (and note it is enough to find 3 solutions not reflectionally symmetric of each other).

To use NGMRES-Picard to directly find still unknown solutions, we simply rescale the NGMRES nonlinear residual at each step by \[\| g(w) \|_{V'} \rightarrow \frac{ \| g(w) \|_{V'}}{\sum_{j=1}^n \| w - u_i \|_{V}^3 },\] where \(V\) is the divergence-free velocity space and \(V'\) is its dual, and \(n\) is the number of velocity solutions \(\{ u_i \}_{i=1}^n\) already found. Such a deflation rescaling is suggested in [33], and we find for this problem the cubic power works well. The general idea is that if the iteration gets close to \(u_i\), the residual will grow and push the iteration away through the optimization problem.

Hence the algorithm is to run NGMRES-Picard until a solution is found, and then rerun using the adjusted nonlinear residual until all solutions are found. Initial guesses are random perturbations of the Stokes solution. We set a NGMRES-Picard nonlinear residual convergence tolerance to \(10^{-5}\) in the \(V'\) norm, and this solution is then passed to Newton to quickly zoom in on the solution. We used \(m=10\) with restarts every 30 iterations, with 100 max iterations per run.

We discretize using \((P_2,P_1^{disc})\) SV elements on a barycenter refinement of a Delaunay mesh refined more heavily near the inlet. A plot of the mesh is shown in Figure 9, and this discretization provides for 97K total dof. Since the system is small, we use direct solves for the linear systems.

The procedure works well overall, and was able to find all 3 solutions that are not reflectionally symmetric of each other in 4 total runs. The first run found the symmetric solution quickly, converging in 45 iterations, the second run found the solution at the bottom in Figure 8 in 69 iterations, the third run did not converge in 100 iterations, and the fourth run found the solution second from top in Figure 8 in 67 iterations. Undoubtedly, our algorithm could be improved e.g by optimizing NGMRES parameters and restarts, allowing more than 100 iterations at each run, or by using smarter Newton algorithms (we used classical, unstabilized Newton), but still NGMRES was able to easily incorporate deflation and was quite effective.

4.4 NGMRES for superlinearly converging iterations↩︎

For a final test, we illustrate the potential for \(m=0\) NGMRES to be used with superlinearly converging iterations. NGMRES (and AA) is likely to help a superlinear solver at early iterations, before the iterates are close enough for superlinear convergence to kick in. Once the iterates get close to a root, how NGMRES affects convergence is not completely understood. However, our theory above in Theorem 4 shows that the higher order terms are quadratic, which suggests that \(m=0\) NGMRES will not slow down an order \(1<r\le 2\) iteration. While this is not a proof, NGMRES theory does compare well with depth 1 AA theory, whose fixed point residual bound has higher order terms that are less than quadratic, and furthermore AA is known to reduce convergence order \(r>1\) to \(\frac{r+1}{2}\) [19], [22].

We now test convergence of \(m=0\) NGMRES and depth 1 AA (using the standard AA from e.g. [15]) on Newton and Secant iterations to solve \(g(x) = x^2 -x-2 =0\): \[q_N(x) = x - \frac{g(x)}{g'(x)}, \;\; q_S([x,y]) = [x - \frac{g(x)}{ \frac{g(x)-g(y)}{x-y}},x].\] Using a 1D test allows us to focus on the higher order term effects.

Results are shown in Figure 10 for initial guess \([0,100]\) and stopping tolerance of \(10^{-10}\) for the nonlinear residual. For both tests, we observe that NGMRES improves convergence, while AA slows it down. The improvement from NGMRES comes early in the iteration, as expected, but NGMRES does not appear to slow down the superlinear asymptotic convergence rates of the Newton and Secant iterations. While this is a simple 1D test, it suggests that further study of NGMRES applied to superlinearly converging methods is warranted.

5 Conclusions↩︎

This paper has provided the first proof of NGMRES acceleration for contractive and noncontractive fixed point iterations used to solve general nonlinear systems. We have established the gain of the optimization problem as the mechanism responsible for the acceleration, and have also identified an additional quantity arising from the optimization problem that accurately predicts the linear convergence rate of the nonlinear residuals at each step. The newly identified quantity, since it is at most 1, also shows the NGMRES only reduces the nonlinear residual (up to higher order terms). Several numerical tests illustrate the theory and also show how NGMRES can be improved for certain problems by using restarts and/or an adaptive depth strategy based on our newly discovered linear convergence rate predictor . Additionally, we have shown that NGMRES has the potential to improve superlinearly converging iterations (an advantage over AA).

For future work, we believe there are several important directions. First, further improvements can be made to the adaptive depth strategy given above, e.g. by considering what is the best depth between \(1\) and \(m\) to use at any given iteration, based on weighing higher order term effects against improvement in the linear convergence rate. Second, NGMRES appears to help with superlinearly converging iterations (as opposed to AA, which makes them worse); this should be explored both analytically and numerically. Third, more documented comparisons of AA and NGMRES need performed to better classify when one is better than the other; the analysis herein may give some hints in this direction. Finally, it may be possible that our result for the noncontractive case with \(m\ge 1\) may be improvable in the higher order terms, since our NLH tests found convergence sometimes after 500 iterations.

References↩︎