Merging Multigrid Optimization with SESOP
Tao Hong, Irad Yavneh, Michael Zibulevsky
TL;DR
This paper introduces a novel optimization method that combines SESOP and multigrid techniques, enhancing convergence efficiency for quadratic function optimization through a hybrid approach and optimal parameter tuning.
Contribution
It presents a new combined framework of SESOP and multigrid optimization, including convergence analysis and parameter optimization for quadratic problems.
Findings
Numerical experiments show improved optimization performance.
Derived approximately optimal fixed parameters for SESOP-TG.
Effective reduction in computational overhead for quadratic functions.
Abstract
A merger of two optimization frameworks is introduced: SEquential Subspace OPtimization (SESOP) with MultiGrid (MG) optimization. At each iteration of the algorithm, the search direction implied by the coarse-grid correction process of MG is added to the low dimensional search-space of SESOP, which includes the preconditioned gradient and search directions involving the previous iterates, called {\em history}. Numerical experiments demonstrate the effectiveness of this approach. We then study the asymptotic convergence factor of the two-level version of SESOP-MG (dubbed SESOP-TG) for optimization of quadratic functions, and derive approximately optimal fixed parameters, which may reduce the computational overhead for such problems significantly.
| Bilinear | Bicubic | Idealized | ||||||
|---|---|---|---|---|---|---|---|---|
| TG-1 | SESOP | Opt | TG-1 | SESOP | Opt | |||
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsSparse and Compressive Sensing Techniques · Matrix Theory and Algorithms · Advanced Optimization Algorithms Research
\newsiamremark
remarkRemark \newsiamremarkhypothesisHypothesis
\newsiamthmclaimClaim \headersMerging Multigrid Optimization with SESOPT. Hong & I. Yavneh & M. Zibulevsky
\externaldocumentex_supplement
Merging Multigrid Optimization with SESOP††thanks: Submitted to the editors DATE.
\fundingThis work was funded by ****
Tao Hong
Irad Yavneh
and Michael Zibulevsky
Abstract
A merger of two optimization frameworks is introduced: SEquential Subspace OPtimization (SESOP) with MultiGrid (MG) optimization. At each iteration of the algorithm, the search direction implied by the coarse-grid correction process of MG is added to the low dimensional search-space of SESOP, which includes the preconditioned gradient and search directions involving the previous iterates, called history. Numerical experiments demonstrate the effectiveness of this approach. We then study the asymptotic convergence factor of the two-level version of SESOP-MG (dubbed SESOP-TG) for optimization of quadratic functions, and derive approximately optimal fixed parameters, which may reduce the computational overhead for such problems significantly.
keywords:
SESOP, acceleration, multigrid, optimization, linear and nonlinear large-scale problems
{AMS}
90C06, 90C25, 65N55, 65K05, 65N22, 65H10
1 Introduction
Multigrid (MG) methods are widely considered to be an efficient approach for solving elliptic partial differential equations (PDEs) and systems, as well as other problems which can be effectively represented on a hierarchy of grids or levels [2, 5, 20, 25]. However, it is often challenging to design efficient stand-alone MG methods for difficult problems, and therefore MG methods are often used in combination with acceleration techniques (e.g., [18]).
In this paper we study MG in an optimization framework and seek robust solution methods by merging this approach with so-called SEquential Subspace OPtimization (SESOP) [28]. SESOP is a general framework for iteratively solving large-scale optimization problems, as described in the next section. The combined framework is called SESOP-MG, and its two-grid (TG) version is called SESOP-TG. We then analyze the asymptotic convergence factor (ACF) of a fixed-stepsize version of SESOP-TG for quadratic optimization problems, and estimate the expected acceleration due to SESOP by means of the so-called -ellipticity measure [3, 22]. Resorting to local Fourier analysis (LFA) [4, 24], we propose two methods to estimate optimal fixed stepsizes of SESOP-TG for quadratic optimization problems cheaply in cases where LFA is applicable. Numerical tests demonstrate the relevance of the theoretical analysis in practice.
The paper is organized as follows. The standard TG and SESOP algorithms are briefly described in the remainder of this section. The merged SESOP-TG/MG algorithm is proposed and tested in Section 2. An analysis of the ACF for a fixed-stepsize version of the SESOP-TG method for quadratic problems is presented in Section 3. There, we also show how to estimate the optimal fixed stepsizes cheaply for some specific problems. Numerical tests validate our analysis and the effectiveness of the proposed method for estimating the optimal fixed stepsizes. Conclusions are drawn in Section 4.
We adopt the following notation. denotes the unknown solution vector, and a convex function we aim to minimize. In the two-grid case, we use superscripts and to denote the fine and coarse grid, e.g., and denote the fine and coarse functions, respectively. In general, we use boldface font to denote vectors and matrices, and T denotes the transpose operator.
1.1 Multigrid (MG)
We consider MG as a method for convex optimization. For an extensive review of MG for PDE optimization, see Borzì and Schulz [1] and references therein. Several authors developed MG optimization for specific problems in the 1990’s, and at the end of the decade Nash [15] formulated MG as a general optimization framework called MG/OPT based directly on the well-known full approximation storage (FAS) scheme of Brandt [2]. Following this, Lewis and Nash applied this framework to systems governed by differential equations [13]. In similar vein, Wen and Goldfard proposed a line search MG method to solve unconstrained convex and nonconvex problems [23]. Toint et al. merged MG optimization with the trust region approach, applying MG to the series of linear subproblems arising in each step of trust region methods. This was applied to nonlinear convex or nonconvex problems, including bound constraints [10, 11, 21, 9]. Recently, Calandra et al. applied MG optimization to reduce the cost of step computation in high-order optimization [6].
Consider the following unconstrained problem defined on the fine-grid:
[TABLE]
where is a convex and differentiable function and the solution set of (1) is not empty. We describe a single iteration of the two-grid algorithm next, with MG obtained by straightforward recursion. Denote by the approximation to the fine-grid solution after the th iteration. Assume that we have defined two full-rank operators, a restriction , and a prolongation , where is the size of the coarse-grid. Furthermore, let denote a coarse-grid approximation to (for example, we may use , but other choices may be used as well). The coarse-grid problem is then defined as follows:
[TABLE]
where is a coarse approximation to , and . The correction term , adapted from FAS, is used to enforce the same first-order optimality condition [15] on the fine and coarse levels. After solving (2) (exactly in TG, and approximately and recursively in MG), the coarse-grid correction CGC direction is defined by
[TABLE]
Finally, the CGC is added to the current fine-grid approximation:
[TABLE]
This may be followed by additional relaxation steps, yielding the st approximation . The two-grid algorithm is written in Algorithm 1.
1.2 SEquential Subspace OPtimization (SESOP)
SESOP [14, 7, 29, 28] is a framework for solving smooth large-scale unconstrained minimization problems such as (1), by sequential optimization over affine subspaces , spanned by the current descent direction (typically preconditioned gradient) and previous propagation directions of the method. If is convex, SESOP yields the optimal worst-case convergence factor of [14], while achieving efficiency of the quadratic Conjugate Gradients (CG) method when the problem is close to quadratic or in the vicinity of the solution.
The affine subspace at the th iteration is defined by
[TABLE]
where is the th iterate, the matrix contains the spanning directions in its columns, the preconditioned gradient and last steps ,
[TABLE]
The new iterate is obtained via optimization of over the current subspace ,
[TABLE]
By keeping the dimension of low, (6) can be solved efficiently with a Newton-type method. SESOP is relatively efficient if the objective function can be represented as ; the evaluation of is cheap, but calculating is expensive. In this case, we can avoid re-computation of during the subspace minimization by storing and . A typical class of problems of this type is the well-known optimization [29].
SESOP can yield faster convergence if more efficient descent directions are added to the subspace . This may include a cumulative or parallel coordinate descent step, a separable surrogate function or expectation-maximization step, Newton-type steps and other methods [7, 29, 28]. In this work we suggest adding the CGC direction provided by the MG framework.
2 Merging multigrid with SESOP
2.1 SESOP-TG
We begin this section by introducing a two-grid version of our scheme, SESOP-TG, and later extend it to the multilevel version, SESOP-MG. In its basic form, the idea is to add the CGC of (3) into the affine subspace by replacing with an augmented subspace and computing the locally optimal to obtain the next iterate, . Our intention is to combine the efficiency of MG that results from fast reduction of smooth error by the CGC, with the robustness of SESOP that results from the local optimization over the subspace. That is, even in difficult cases where the CGC is inadequate, the algorithm should still converge at least as fast as the standard SESOP, because an inefficient direction simply results in small (or conceivably even negative) coefficient. SESOP-TG is presented in Algorithm 2. Note that we add the option of two additional steps to the usual SESOP algorithm, the pre- and post-relaxation at steps 2 and 13 in Algorithm 2, commonly applied in MG algorithms. This allows us to advance the solution with low computational expense whenever the coarse-grid direction is inefficient (due the the fact that previous use of the CGC direction greatly reduced the smooth error). For the relaxation we typically use nonlinear Jacobi or Gauss-Seidel applied to the gradient equation. The flowchart of Algorithm 2 is presented in Figure 1.
Remark 2.1**.**
Evidently, if we select to contain only , with , then SESOP-TG- reduces to TG. Note that in Algorithm 1, we need to valuate to formulate but we do not use it further. Here, however, we put together with the CGC direction in the augmented subspace for further use. The main drawback of Algorithm 2 is of course the need to solve the subspace minimization problem for at each iteration. However, in certain cases of problems with special structure or large size, the cost of the subspace minimization may become negligible; see numerical tests in Section 2.2.
Remark 2.2**.**
*Similarly to standard multigrid methods, the multilevel algorithm called SESOP-MG- is derived by recursively treating (2) at step 10 of Algorithm 2. Compared with the subspace in (2), the subspaces in the coarse problems only contain two directions: the preconditioned gradient and the coarse grid correction from the coarser level. *
2.2 Numerical tests
We demonstrate the SESOP-MG- algorithm performance for two nonlinear problems and compare it with steepest descent (SD), Nesterov acceleration [16], and limited-memory BFGS (LBFGS) [17]. The MG method proposed in [23]—denoted “MG-Line”—is adapted for this comparison, to illustrate the effectiveness of using history. We use SD as the relaxation for “MG-Line”. For SESOP-MG-, we employ Newton’s method for the subspace minimization at each iteration. All the tests are performed on a laptop with 2.3GHz Intel Core i.
Example I
Consider the following variational problem on a uniform grid [23],
[TABLE]
where is the gradient, , , and
[TABLE]
Through the Euler-Lagrange equation, the PDE formulation of (7) reads
[TABLE]
where denotes the Laplacian.
Equation (7) is discretized by finite differences, using first order forward differences for in (7), resulting in a five-point stencil for in (8). The finest grid size is , and we employ levels, coming down to grids of size . The coarse problems are defined by rediscretization, and full-weighted residual transfers and bilinear interpolation are used as the restriction and prolongation. The relaxation sweep parameters are and for MG-Line. BFGS with up to ten iterations is used for solving the problem on the coarsest level. The minFunc toolbox [19] is used for BFGS and Newton’s method. The maximal number of iterations of SESOP-MG- and MG-Line are set to be and , respectively. The other methods are allowed to use iterations. For clarity of display, we denote by the minimal objective value over all methods after running the maximal allowed number of iterations.
From Figure 2, we clearly observe that SESOP-MG- is the fastest algorithm in terms of both iteration count and CPU time. Moreover, we note that SD, Nesterov, and LBFGS converge fast initially, but then slow down. This well-known phenomenon is a result of the fact that these methods cannot efficiently eliminate low-frequency error, unlike MG methods which use CGC. Note also that SESOP-MG- is significantly faster than MG-Line, demonstrating the effectiveness of introducing history for acceleration.
Example II
Our second example is the -Laplacian,
[TABLE]
where . The corresponding PDE of (9) is:
[TABLE]
The parameter is introduced to maintain differentiability (hence a positive denominator). The function is defined by substituting into (10). Note that solving (9) becomes especially challenging when is close to . We choose to study the performance of our approach. The experimental setting is the same as that of Example I, except the maximal number of iterations allowed. As seen in Figure Figure 3, SESOP-MG- is the fastest of the five methods, followed by MG-Line, as in Example I.
Dependence of SESOP-MG- on
Next we show how the choice of the search-space size affects the performance of SESOP-MG-. To this end, we apply SESOP-MG- with different to (9) with . In Figure 4(a), we see that using a larger yields faster convergence rates. However, a larger also increases the complexity of the subspace minimization results in higher CPU times per iteration. From Figure 4(b), we find that is a good compromise for achieving low CPU times for this test. We also see that results in higher times than because, in practice, the larger size of the subspace introduces some numerical difficulties.
The numerical tests demonstrate the potential advantage of SESOP-TG/MG- for such types of optimization problems. However, SESOP-TG/MG- comes with the cost of solving a subspace minimization problem at each iteration. From Figure 4, we observe that the performance of SESOP-TG/MG- deteriorates for large . In the case of quadratic optimization problems, we may be able to avoid the subspace minimization, by using fixed nearly-optimal stepsizes for SESOP-TG-. Indeed, we derive such fixed stepsizes, and show that they yield ACF’s that are comparable to those obtained by subspace minimization. In such cases, we get the acceleration nearly for free, provided that we can estimate the fixed parameters efficiently. To this end, we propose two heuristic methods, based on local Fourier analysis (LFA) and smoothing analysis, to estimate the optimal fixed stepsizes cheaply.
3 Convergence factor analysis of SESOP-TG- for quadratic problems
To gain insight, we analyze SESOP-TG- for quadratic optimization problems, which are equivalent to the solution of linear systems. We first derive a fairly general formulation for SESOP-TG-. Then we explore, under certain simplifying assumptions, a fixed-stepsize variant of SESOP-TG-. In this analysis we assume for simplicity no pre- or post-relaxation steps.
Consider the linear system
[TABLE]
where is a symmetric positive-definite (SPD) matrix, and we omit superscripts for notational simplicity. Evidently, solving (11) is equivalent to the following quadratic minimization problem:
[TABLE]
Given iterates and , the next iterate produced by SESOP-TG- is given by
[TABLE]
where are the optimal weights associated with the three directions comprising : multiplies the so-called history, that is, the difference between the last two iterates; multiplies the preconditioned gradient; multiplies the CGC direction . Here, represents the coarse-grid matrix approximating , which is most commonly defined by the Galerkin formula, , or simply by rediscretization on the coarse-grid in the case where is the discretization of an elliptic PDE on the fine-grid. Subtracting from both sides of (13), and denoting the error by , we get
[TABLE]
Rearranging (14) yields
[TABLE]
where
[TABLE]
and denotes the identity matrix. Define the vector . Then, (15) implies the following relation:
[TABLE]
and the asymptotic convergence factor (ACF) of SESOP-TG- is given by the spectral radius of .
To analyze the ACF of SESOP-TG-, we continue under the assumption that the coefficients , , are fixed, and compute the optimal coefficients, yielding the smallest ACF. Denote by an eigenvalue of with eigenvector ,
[TABLE]
Hence, and . This yields
[TABLE]
Thus, is an eigenvector of with eigenvalue leading to
[TABLE]
with solutions
[TABLE]
The spectral radius of is therefore
[TABLE]
where runs over the eigenvalues of and is a sign function which is defined to be 1 for non-negative arguments and -1 for negative arguments. For the remainder of our analysis we focus on the case where the eigenvalues of are all real.
3.1 The case of real b
In many practical cases, the eigenvalues of are all real, which simplifies the analysis and yields insight. We focus next on a common situation where this is indeed the case.
Lemma 3.1**.**
Assume:
The prolongation has full column-rank and the restriction is the adjoint of the prolongation: . 2. 2.
The coarse-grid operator is SPD. 3. 3.
The preconditioner is SPD.
*Then the eignevalues of are all real. *
Proof 3.2**.**
Because is SPD, there exists a SPD matrix such that . The matrix
[TABLE]
*is similar to so they have the same eigenvalues. Moreover, is evidently symmetric, so its eigenvalues are all real. *
The first two assumptions of this lemma are satisfied commonly, including of course both the case where is defined by rediscretization on the coarse grid and the case of Galerkin coarsening. The preconditioner is SPD for some commonly used MG relaxation methods, including Richardson (where is the identity matrix), Jacobi (where is the inverse of the diagonal of ), and symmetric Gauss-Seidel.
We next adopt the change of variables , . This yields
[TABLE]
with
[TABLE]
Lemma 3.3**.**
*The eigenvalues of are real and positive for any . *
Proof 3.4**.**
*This follows from the fact that , and are SPD. *
We henceforth denote the eigenvalues of by , and assume (that is, and are of the same sign), so the ’s are all real and positive. We then proceed by fixing and optimizing and so as to minimize the spectral radius of . In light of Lemma 3.1, we have , and the ACF as a function of is given by:
[TABLE]
where . For convenience, we use (respectively, ) to denote considered as a single-variable function with a fixed (respectively, ). Our parameter optimization problem is now defined by
[TABLE]
The following three lemmas provide us with a closed-form solution of this minimization problem, if the smallest and largest eigenvalues of are given.
Lemma 3.5**.**
*For any fixed and , is maximized at either or , the smallest and largest eigenvalues of , respectively. *
Proof 3.6**.**
*See Appendix A. *
Lemma 3.7**.**
*For any fixed , . *
Proof 3.8**.**
*See Appendix B. *
Plugging into (20) yields
[TABLE]
where with the condition number of . Lemma 3.9 provides the optimal , which minimizes in (21).
Lemma 3.9**.**
*. *
Proof 3.10**.**
*See Appendix C. *
Substituting the expression of into (21), we get the worst-case convergence factor with optimal and
[TABLE]
Note that indicates that the history term becomes significant if the problem is ill-conditioned.
Remark 3.11**.**
The optimal ACF of the fixed-coefficient variant of SESOP-TG- is found to be , where is the condition number of , optimized over . (The optimal choice of will be discussed further below.) By using the definition of , the optimal coefficient for the preconditioned gradient term as specified in Lemma 3.7 can be written as
[TABLE]
By setting , we retain SESOP-TG-[math] and then the optimal becomes
[TABLE]
Evidently, the asymptotic convergence factor of SESOP-TG-[math] is
[TABLE]
*Comparing (22) with (25), we see the significant improvement provided by the use of a single history direction, with the condition number replaced by the square root. This result is reminiscent of the conjugate gradients (CG) method compared to steepest descent (SD) for quadratic problems, and indeed there is an equivalence between SESOP and CG (respectively, SD) for the case of single history (respectively, no history) with . We note that the condition number in our scheme is that of , not , which differs from the case where we do not use the CGC direction. The advantage of SESOP is in allowing the addition of various search directions, at the cost of optimizing the coefficients. This study is partly aimed at reducing this cost. *
3.2 Towards optimizing the condition number of
The previous analysis indicates that we should aim to minimize , the condition number of . In certain cases, particularly when is a circulant matrix (typically the discretization of an elliptic PDE with constant coefficients on a rectangular or infinite domain), this can be done by means of Fourier analysis [24], as discussed in Section 3.3. Here, we begin with a more general discussion to gain insight into this matter. We consider the case of Galerkin coarsening, , and . Guided by [8], we consider the case where the columns of the prolongation matrix are comprised of a subset of the eigenvectors of .
Lemma 3.12**.**
Let , denote the eigenvalues and eigenvectors of . Assume that the columns of the prolongation matrix are comprised of a subset of the eigenectors of , and denote by the range of the prolongation, that is, the subspace spanned by the columns of . Denote
[TABLE]
Then, the condition number of in (19) with is given by
[TABLE]
Proof 3.13**.**
*See Appendix D. *
We see that the second term in , which corresponds to the direction provided by the CGC, increases the eigenvalues associated with the columns of the prolongation by . It thus follows from (26) that, to obtain any advantage at all from the CGC direction in reducing , the eigenvector associated with the smallest must be included amongst the columns of the prolongation, and therefore . Similarly, considering the numerator in (26), it is clearly advantageous that the eigenvector corresponding to the largest not be included in the range of the prolongation. With these assumptions, we obtain the following result for the optimal which minimizes .
Theorem 3.14**.**
Assume and . Then, the condition number is minimized by choosing
[TABLE]
resulting in the optimal ,
[TABLE]
Proof 3.15**.**
*See Appendix E. *
Remark 3.16**.**
Notice that is either equal to or at most “slightly larger” because
[TABLE]
*That is, even in the regime , the optimal condition number is increased by less than 1. Note that yields a convergence factor (with no history) of , which matches that of the classical TG algorithm with optimally weighted Richardson relaxation followed by CGC. In Section 3.4, we will show the connection of presented here with the so-called -ellipticity measure, which is a qualitative criterion for the existence of local smoothers for a given elliptic PDE [22]. Finally, we note that the optimal prolongation is obtained by choosing the columns of to be the eigenvectors associated with the smallest eigenvalues of (similarly to the classical TG case [8]). This clearly minimizes the ratio . Furthermore, this choice yields , so if (so ), then we obtain . *
3.3 Optimizing the condition number of in practice
As our analysis indicates, optimizing so as to minimize is a crucial step towards obtaining an optimal convergence factor. In general, the complexity of minimizing is very high and the explicit results of Section 3.2 are only valid when the prolongation is comprised of eigenvectors of , which is impractical, because the prolongation needs to be very sparse for efficiency. In certain cases, such as when results from discretizing elliptic PDEs with constant or slowly varying coefficients, Local Fourier Analysis (LFA) [2, 24] can be used in conjunction with the analysis of Section 3.2 to yield effective approximate results.
LFA, also called Local Mode Analysis, is a useful quantitative tool for estimating the asymptotic convergence factor ACF of MG methods for elliptic PDEs with constant coefficients. We describe LFA here briefly, and for further details refer the reader to standard multigrid textbooks or to [24], a comprehensive book on LFA for MG. We consider the two-dimensional case for simplicity. For that is a discretization of a constant-coefficient elliptic PDE on an infinite or doubly periodic domain of uniform mesh-size , grid-based functions of the form with and , are eigenfunctions of . Furthermore, if standard shift-invariant prolongation is used, such as bilinear or bicubic interpolation, then appropriate subsets of dimension four of the eigenfunctions form subspaces that are invariant under multiplication by (as in standard two-level analysis). The upshot is that we can compute the eigenvalues of at a fairly moderate cost, and use linesearch to find which optimizes the condition number of . This approach is demonstrated later in numerical examples, including the option of saving computational cost by optimizing on relatively coarse grids.
We can furthermore obtain a very cheap approximation to the optimal and by making the simplifying assumptions that are commonly used in computing the so-called smoothing factor of relaxation (known as smoothing analysis) as follows. Denote by the eigenvalue of associated with the grid-function , and partition the domain into low- and high-frequencies as in Definition 3.17.
Definition 3.17** (Low- and High-frequency components).**
[TABLE]
Smoothing analysis simplifies by assuming that the CGC acts as a projection onto the high-frequency subspace, that is, it has no effect on high-frequency error components, while it eliminates exactly low-frequency error components. Under these simplifying assumptions, we obtain
[TABLE]
Moreover, the eigenvalues of the second term in are given by with the eigenvalues of for [26]. Now, we can write as:
[TABLE]
where and denote the maximum and minimum value of , respectively. Similarly to Theorem 3.14, the optimal is given by
[TABLE]
We note that, for elliptic PDEs, the accuracy of using the approximations (28) and depend on the order of intergrid transfers [12, 26]. In the following numerical tests, we evaluate experimentally the accuracy of using (30) to estimate the fixed stepsizes in practice.
3.4 A connection with the -ellipticity measure
Using Definition 3.17, we define an “idealized” MG method as a two-grid method that affects high-frequency error components only on the fine grid and eliminates all low-frequency error components via the CGC, corresponding to (28). Then, the ACF of an idealized SESOP-TG- is with , corresponding to the discussion of Section 3.2 for ill-conditioned problems. Alternatively, we can express this idealized ACF in terms of , which is known in the literature as the -ellipticity measure, obtaining . For SESOP-TG-[math], in contrast, the ACF becomes , which is the well-known smoothing factor of optimally damped Jacobi relaxation for symmetric problems [22].
3.5 Numerical tests—continued
In this section, we first test the accuracy of using (17) to predict the ACF of SESOP-TG- with fixed stepsizes. Then, compare the ACF with minimization over the subspace to that obtained with optimized fixed stepsizes. Finally, we examine the performance of the two proposed heuristic methods (cf. Section 3.3) for estimating the fixed stepsizes. The rotated anisotropic diffusion problem is chosen as the test problem,
[TABLE]
where and denote the second partial derivatives of in the coordinate system. Denote by the angle between and the grid-aligned coordinate system, . We rewrite (31) as
[TABLE]
where and . Using a nine-point stencil,
[TABLE]
to discretize (32) on a uniform grid with mesh-size and prescribed boundary conditions, we get a linear system of the form
[TABLE]
ACF Prediction
In our first test, we simply choose , and , and then use Lemmas 3.7 and 3.9 to compute and . The resulting algorithm is referred to as “SESOP-TG--Fixed”. We discretize (32) on a grid, imposing Dirichlet boundary conditions, and we employ bilinear prolongation and full-weighted residual transfers. The coarse-grid operators are defined by direct rediscretization. Finally, we set , i.e., no preconditioning. The ACF achieved in practice by SESOP-TG--Fixed is evaluated as the geometric mean of the convergence factor per iteration in the last iterations, which are terminated at iterations or when the residual norm is smaller than , whichever comes first. The convergence factor at the th iteration is defined by the ratio of the successive residual norms , where and is the approximate solution at the th iteration. In Figure 5, we see that the ACF predicted by the spectral radius of in (17) matches the practical results well.
Comparison between subspace minimization and optimal fixed stepsizes
Next, we compare the ACFs of three different approaches to determining the fixed stepsizes: 1) SESOP-TG--Fixed defined above; 2) subspace minimization (classical SESOP); 3) optimized fixed stepsizes. The optimized stepsizes are obtained by minimizing the true condition number of in (19) by employing linesearch over . The test problem remains unchanged, except that we test both bilinear and bicubic prolongations and use grids. The ACF achieved in practice is estimated as above by the geometric mean of the last iterations when the algorithm reaches iterations or the residual norm is smaller than . Note also that of Section 3.4, i.e., one over the -ellipticity measure, is used here to represent the idealized ACF for comparison.
From Table 1, we clearly see that optimized fixed stepsizes yield a lower ACF than SESOP-TG--Fixed for the rotated anisotropic diffusion problem, and in fact its ACF is at least as good as that obtained by subspace minimization. This suggests that, if we can estimate the optimal stepsizes efficiently, then we can significantly reduce computation time because the subspace minimization, in general, is expensive. Also, we see that bicubic prolongation yields a lower ACF than bilinear prolongation, consistent with [12].
The relative disadvantage of SESOP-TG--Fixed for rotated anisotropic diffusion stems mainly from fixing , whereas the optimal is higher for this problem, consistent with the analysis of [26].
Optimizing of in practice
We next study the efficacy of the two heuristic methods proposed in Section 3.3—estimating on a coarse grid or using (30)—to select the fixed stepsizes.
In the following tests we use grids to discretize (32). Denote by the ACF obtained by optimizing of on a grid. Then, we define the “Deterioration Factor” (DF),
[TABLE]
to measure the deterioration of the ACF incurred by optimizing on a coarser grid. For example, means that, asymptotically, it takes twice as many iterations of the algorithm which uses fixed stepsizes optimized on coarser grids to achieve the same error reduction as a single iteration with the true optimal stepsize. Additionally, we compare between bilinear and bicubic interpolation [12]. In Figure 6, we see for , which means that determining the stepsizes on grids results in just a asymptotic increase in the required iterations. Moreover, we observe that the use of higher order interpolation tends to reduce the DF.
Now, we study the accuracy of using (30) for selecting the fixed stepsizes. From [26, Equation (14)], we know that the accuracy of (29) for approximating of depends on the order of intergrid transfers. In Figure 7, we show a comparison of ACFs with different orders of intergrid transfers. Note that the intergrid transfers used here have the same order for low and high-frequency modes. From Figure 7, we find that (30) becomes more accurate and very close to the idealized estimate when the order of restriction increases. This is due to the fact that higher order intergrid transfer operators filter out the high-frequency modes, and therefore the idealized assumptions are more closely satisfied. Note that the optimized algorithm is even better than the the idealized one, as also observed in Table 1, because the optimized algorithm takes into account all the modes globally when minimizing . Moreover, for the optimized version we see that increasing the order of the restriction begins increase the ACF slightly, bringing it closer to the idealized value. We finally note that the purpose of this test is to academically study the relationship between the accuracy of (30) and the use of intergrid transfers and, in practice, it is not cost-effective to use very high order intergrid transfers.
Remark 3.18**.**
*In this part we have examined the accuracy of two heuristic methods for determining the fixed stepsizes. In practice, the method shown in Figure 6 may be more attractive because the increase of iterations due to computing the stepsizes on a coarse grid is modest. However, for some practical problems, (30) is also attractive because its computation is cheaper than working on a coarse grid. Moreover, for many practical problems, we need to solve (33) with multiple , and then the additional computation for selecting the stepsizes is negligible. *
Practical tests
Now we study the performance of applying the multilevel version (V-cycle) of the proposed scheme to (32). The multilevel version of Algorithm 1 is denoted by “MG”. Using preconditioned conjugate gradients (PCG) with MG as the preconditioner is denoted by “PCG-MG”. The fixed-stepsize version of SESOP-MG- is denoted by SESOP-MG--Opt when is computed on a coarse grid and by SESOP-MG--(30) when (30) is used. We use grids to discretize (32) and obtain for SESOP-MG--Opt on grids. The coarse problems are obtained by direct rediscretization, and the bilinear prolongation and full-weighted residual transfers are employed. For coarse problems, we use Jacobi relaxation with optimal damping factor as estimated in [27] and and . Note that on the finest level, we set and and only one evaluation of the gradient is allowed.
In Figure 8(a) (), we see SESOP-MG- and PCG-MG perform identically. However, in Figure 8(c) (), we find that SESOP-MG- becomes faster than PCG-MG. Indeed, in this case, it helps to scale the residual before restricting to the coarse for MG (also see the deterioration of MG in Figure 8(c)), but SESOP-MG- can automatically scale the residual via the subspace minimization [4]. Moreover, from Figures 8(b) and 8(d), we observe that SESOP-MG--Opt and SESOP-MG--(30) are the fastest methods in terms of CPU time because these two methods inherit the effectiveness of SESOP-MG- but need much less CPU time. It is interesting to note that SESOP-MG--Opt and SESOP-MG--(30) work almost the same in these two tests, demonstrating the effectivenss of the proposed strategies for determining the stepsizes.
4 Conclusion
In this paper, we merge multigrid (MG) optimization with SESOP optimization. The numerical experiments on linear and nonlinear problems illustrate the effectiveness and robustness of our scheme. Moreover, for linear problems, if only one history is used, we derive optimal fixed stepsizes that allow us to avoid the expensive subspace minimization and save computation. Specifically, for elliptic PDEs with constant coefficients, we propose two heuristic methods to estimate the fixed stepsizes based on LFA and smoothing analysis and the efficiency of these two methods is demonstrated in numerical tests.
Appendix A Proof of Lemma 3.5
Consider in (20) as a continuous function of a positive variable , with fixed and :
[TABLE]
We distinguish between two regimes: (I): , the square-root term is imaginary, and we simply get ; (II): ,the square-root term is real and the derivative of with respect to is
[TABLE]
We ignore the irrelevant choice , for which the method is obviously not convergent. Using , we get
[TABLE]
Notice that is a symmetric function of and it is furthermore convex because its derivative is strictly negative for and positive for . It follows that, regardless of the sign of throughout the regime , there exists no local maximum of . Hence,
Appendix B Proof of Lemma 3.7
Consider as a continuous function of a real variable with and fixed. For , the derivative of with respect to is
[TABLE]
As in Lemma 3.5, the meeting point of and , which lies between and , is where is minimized with respect to and then the optimal enforces resulting in leading to For , we have which is irrelevant to .
Appendix C Proof of Lemma 3.9
Consider two regimes: (I) ; (II) . For (I), we have and where are the two solutions of . Evidently, is the solution to minimize .
For (II), we have . We first note that
[TABLE]
and then . Evidently, the derivative of with respect to is
[TABLE]
Note that is positive for and negative for , hence, (34) is positive for , whereas for we have
[TABLE]
The last inequality is due to the fact that that and . Evidently, the optimal to minimize is either or . However, for that is the solution. Substituting the expression of into , we get the desired result.
Appendix D Proof of Lemma 3.12
Any eigenvector (respectively, ) is an eigenvector of the CGC matrix , with eigenvalue 1 (respectively, 0), because if then it can be written as (that is, it is the th column for some ), so
[TABLE]
whereas if then it is orthogonal to the columns of , so
[TABLE]
It follows that the eigenvectors of are , with eigenvalues given by
[TABLE]
By using the definition of , , , , and , the desired result is derived.
Appendix E Proof of Theorem 3.14
To minimize , the optimal should minimize the numerator and maximize the denominator. Note that and are the ones to minimize the numerator and maximize the denominator, respectively. For , we have and then:
2. 2.
3. 3.
,
Evidently, the optimal can be any one between and yielding . By applying the same reasoning to , the optimal achieves at . Summarizing, is the optimal solution for both cases. Substituting into (26), we get the desired result.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] A. Borzì and V. Schulz , Multigrid methods for PDE optimization , SIAM Review, 51 (2009), pp. 361–395.
- 2[2] A. Brandt , Multi-level adaptive solutions to boundary-value problems , Mathematics of Computation, 31 (1977), pp. 333–390.
- 3[3] A. Brandt , 1984 multigrid guide with applications to fluid dynamics . Monograph, GMD-Studie 85, GMD-FIT, Postfach 1240, D-5205, St. Augustin 1, West Germany, 1985.
- 4[4] A. Brandt and O. E. Livne , Multigrid Techniques: 1984 Guide with Applications to Fluid Dynamics, Revised Edition , vol. 67, SIAM, 2011.
- 5[5] W. L. Briggs, V. E. Henson, and S. F. Mc Cormick , A multigrid tutorial , SIAM, second ed., 2000.
- 6[6] H. Calandra, S. Gratton, E. Riccietti, and X. Vasseur , On high-order multilevel optimization strategies , SIAM Journal on Optimization, 31 (2021), pp. 307–330.
- 7[7] M. Elad, B. Matalon, and M. Zibulevsky , Coordinate and subspace optimization methods for linear least squares with non-quadratic regularization , Applied and Computational Harmonic Analysis, 23 (2007), pp. 346–367.
- 8[8] R. D. Falgout, P. S. Vassilevski, and L. T. Zikatanov , On two-grid convergence estimates , Numerical Linear Algebra with Applications, 12 (2005), pp. 471–494.