Multi-Trajectory Parallel-in-Iteration Optimization

Introduction

Consider the unconstrained optimization problem of minimizing a convex function $F:\mathbb{R}^{N}\to\mathbb{R}$, i.e., solving $\min_{\mathbf{u}\in\mathbb{R}^{N}} F(\mathbf{u})$. Suppose that for solving this problem we employ a 1-step iterative optimization algorithm of the form

where $\Phi(\cdot)$ is a linear or nonlinear operator and $\mathbf{w}_{i}$ is a generic forcing term that does not depend on $\mathbf{u}_{i}$. For example, if $F:=f+g$, where $f$ is $L$-smooth and $g$ is nonsmooth, then we may use proximal gradient [5, 4] as the optimization method, which takes the form of (1) with

where $\partial g$ denotes the subgradient of $g$, and $\mathbf{w}_{i}=\mathbf{0}$ for all $i>0$. If $F:=f$, i.e., the objective function is $L$-smooth, then (2) reduces to the gradient descent operator [18, 19], another prototypical method of the form (1). These algorithms may converge slowly; for general convex problems, gradient descent and proximal gradient have a convergence rate of $O(1/i)$, such that tens of thousands or more iterations [7, 21, 5] may be required to get an accurate approximation to the optimum. Since the algorithms are sequential in $i$, a long wall-clock time may then be required to solve the problem.

Assume $N_t$ iterations of the method (1) starting from an initial guess $\mathbf{u}_{0}\in\mathbb{R}^{N}$ are required in order to approximate $\mathbf{u}^{*}$, a minimizer of $F$, to the desired accuracy. Computing this sequence of iterates is equivalent to solving the “all-at-once” system,

where $\mathbf{w}_0 \equiv\mathbf{u}_0$ is the initial guess. Parallel-in-time algorithms, such as multigrid reduction-in-time (MGRIT) [11], are designed to solve systems of the form (3) for time-dependent evolution equations [13, 20]. Further inspired by the fact that certain optimization algorithms can be interpreted as discretizations of time-dependent ODEs [26, 3, 23], in [2] we re-purposed parallel-in-time methods, namely MGRIT, to parallelize over iterations of optimization algorithms by writing those iterations in the form of (3). We showed that this approach can substantially reduce the predicted wall-clock time compared with executing the optimization algorithm sequentially. We note that related ideas have also since been used to parallelize gradient flows with Parareal [27].

While parallel-in-time methods approximate entire exact trajectories with high accuracy, the goal for optimization problems is to sufficiently reduce the gradient of the initial guess in order to obtain an approximate minimizer of the objective function. In our parallel-in-iteration approach [2], MGRIT is used to approximate the exact sequence of iterates produced by a fixed number of steps of the sequential optimization method (1). Since the intermediate iterates generated by (1) are irrelevant to the final solution, expending computational resources to accurately approximate the entire trajectory is inefficient. Instead, we can consider multiple trajectories, and approximate them inexactly with just one or a few MGRIT iterations, where each new trajectory is started from the final point of the previous trajectory.

We call this restarted MGRIT algorithm R-MGRIT. Given an optimization problem with minimizer $\mathbf{u}^*$ and optimization method (1), we consider applying MGRIT to the associated all-at-once system $\mathcal{A}(\mathbf{u}) = \mathbf{w}$ (3) with $N_t + 1$ time points. After one MGRIT iteration, we obtain approximate solutions $\mathbf{v}_i \approx\mathbf{u}_i$, $i = 0, 1, \ldots, N_t$. If the approximation at the final point, $\mathbf{v}_{N_t}$, is not sufficiently accurate, we formulate a new all-at-once system,

where the new initial guess is updated to $\widetilde{\mathbf{w}}_0 = \widetilde{\mathbf{u}}_0 := \mathbf{v}_{N_t}$. This updated system determines a new trajectory yielded by (1) starting from the improved initial guess $\widetilde{\mathbf{u}}_0$. Rather than accurately approximating a single long trajectory, R-MGRIT approximates a sequence of potentially shorter trajectories, using the final point of each to initialize the next. We therefore refer to R-MGRIT as a multi-trajectory parallel-in-iteration algorithm.

Our motivation is to investigate whether doing a one MGRIT iteration per trajectory leads to a more efficient parallel-in-iteration method than running MGRIT to convergence on a single trajectory, assuming the trajectory has enough optimization steps to reach the desired optimization accuracy. We prove, based on existing MGRIT error propagation theory, that standard MGRIT and R-MGRIT yield approximately the same gradient reduction after a fixed number of iterations. However, we find that we can shorten the consecutive trajectories without increasing the number of trajectories needed for convergence, which increases the parallel speedup while at the same time reducing the number of processors required.

The rest of this paper is organized as follows. In § 2, we introduce our model problems and the optimization methods considered in this paper. In § 3, we present the MGRIT algorithm in an optimization context and define our standard and multi-trajectory parallel-in-iteration schemes. In § 4, we compare the standard and restarting MGRIT schemes theoretically, which justifies our approach of truncating the time domain for multi-trajectory solves. § 5 provides numerical results, in which we observe that both standard and restarted schemes converge in approximately the same number of iterations; moreover, the number of R-MGRIT trajectories required for convergence is either preserved or reduced by truncating the time domain to a fixed fraction of the original number of time points. In § 6, we provide theoretical speedup estimates for the truncated restarting scheme, showcasing its superiority over the standard MGRIT scheme. In § 7, we apply Anderson acceleration to accelerate R-MGRIT, and concluding remarks are given in § 8.

Model problems

In the following, we briefly summarize our model problems and present the optimization algorithms of our choice for each problem; for a more complete discussion, see [2]. The optimization algorithms are employed in the parallel-in-iteration framework, such that in what follows we speak of a “fine-grid” optimization method, and a complementary “coarse-grid” optimization method; see discussion in subsection 3.1. We let $\mathbf{x} \in\mathbb{R}^{d}$ denote a coordinate on $\Omega\subset\mathbb{R}^{d}$, such that $\mathbf{x}=x$ when $d=1$ and $\mathbf{x}=(x,y)$ when $d=2$. The matrices which appear in the definition of our model problems arise from finite-difference discretizations of (negative) Laplacian operators, and are defined as

for $d=1$ and $d=2$, respectively, with $h=1/(n+1)$. Here, $T\in\mathbb{R}^{n\times n}$ is a tridiagonal matrix with elements $-1$, $4$ and $-1$ on the sub-, main and superdiagonals, respectively, and $I\in\mathbb{R}^{n\times n}$ is the identity matrix. The 2-norm of $A_{1}$ is $\|A_{1}\|_{2}=\frac{1}{h^{2}}4\sin^{2}\left(\frac{n\pi}{2(n+1)}\right)\approx\frac{4}{h^{2}}$ for $h\ll1$, while the 2-norm of $A_{2}$ is $\|A_{2}\|_{2}=\frac{1}{h^{2}}8\sin^{2}\left(\frac{n\pi}{2(n+1)}\right)\approx\frac{8}{h^{2}}$ for $h\ll1$. Throughout the paper, we refer to $A_{1}$ and $A_{2}$ as 1D and 2D Laplacian matrices, respectively.

Quadratic minimization problem (MP1)

Our first model problem, MP1 (7), is the following straightforward quadratic minimization,

where $A_{1}$ is the 1D Laplacian matrix (5), such that $N=n$, and $b\in\mathbb{R}^{N}$ is some prescribed vector. We consider solving this problem with a fine-grid gradient descent iteration, which uses the update

with stepsize $s=1/L$ for $L=\|A_{1}\|_{2}$. On the coarse-grid we are required to step over the $2/L$ stability limit of gradient descent, and as such we use the proximal point method, with updates given by

Elastic obstacle problem (MP2)

Our second model problem, MP2 (10), is an approximate linearized discretization of the elastic obstacle problem (EOP) [6, 22] with an $\ell^1$ penalty [28]. Given some discretized obstacle $\phi\in\mathbb{R}^N$ and dimension $d \in\{1,2\}$, letting $\mathbf{b}=-A_d\phi$, MP2 is given by

where $A_d$ is the $d$-dimensional Laplacian matrix from (2.1), $\lambda>0$ is a predefined penalty parameter, and $\sum_i(-(\mathbf{u})_i)_+=\lVert(-\mathbf{u})_+\rVert_1$, where $(a)_+=\max\{0,a\}$ and $\lVert\cdot\rVert_1$ denotes the $\ell^1$ norm. For sufficiently large $\lambda$, the solution of (10) is an approximate solution of the continuous EOP up to a change of variables [16, 12].

As in [2], we consider proximal gradient descent as the fine-grid method for (10), with updates taking the form

where $\nabla f$ denotes the gradient of $f:=\frac{1}{2}\langle A_d\mathbf{u},\mathbf{u}\rangle-\langle\mathbf{b},\mathbf{u}\rangle$ and $\partial g$ denotes the subgradient of $g:=\lambda\lVert(-\cdot)_+\rVert_1$; the proximal operator $P_{sg}=(I+s\partial g)^{-1}$ is given component-wise by

We take $s=1/L$ with $L=\lVert A_d\rVert_2$, $d=1,2$. Analogously to MP1, on coarse grids we must step over the $2/L$ stability limit of the proximal gradient descent method, so that on coarse grids we use the alternating proximals method:

Parallel-in-iteration optimization methods

In [2], we presented our parallel-in-iteration framework for optimization problems based on re-purposing parallel time integration algorithms for time-dependent differential equations, motivated by the fact that optimization algorithms often have interpretations as discretizations of time-dependent differential equations (such as gradient flow). In this section, we recall the MGRIT algorithm [11], which we present in the context of optimization. Then, we define our multi-trajectory R-MGRIT parallel-in-iteration scheme. While we use MGRIT in our exposition of the multi-trajectory framework, we note that other parallel-in-time methods, such as PFASST [17, 10], for example, could be employed instead.

The MGRIT algorithm

Given an initial guess $\mathbf{u}_0\in\mathbb{R}^N$ and 1-step iterative method of the form (1), we consider the associated all-at-once system of the form (3), where $\mathbf{w}_0=\mathbf{u}_0$ is the initial guess and $\mathbf{u}$ is the exact solution. The associated residual vector, given an approximation $\mathbf{v}\approx\mathbf{u}$, is defined as

We assume that $N_t$ is sufficiently large such that $\mathbf{u}_{N_t}$ is a sufficiently accurate approximation of the minimizer $\mathbf{u}^*$. Instead of solving (3) sequentially with time-stepping or forward substitution, MGRIT aims to solve it iteratively and in parallel via the introduction of a coarser approximation of the problem. Let $t_i = i\Delta t$, $i = 0, 1, \ldots, N_t$ be a fine temporal mesh with constant spacing $\Delta t$; in an optimization setting, we take it to be the fixed stepsize $s > 0$ of the underlying optimization method (1). Then, given an integer coarsening factor $m > 1$, we introduce a coarse temporal mesh $T_j = j\Delta T$, $j = 0, 1, \ldots, N_T$, where $N_T = N_t/m$, with constant spacing $\Delta T = m\Delta t$. We denote the $N_T + 1$ points appearing exclusively on the coarse mesh by $C$-points, while all remaining points are $F$-points. Figure 1 illustrates the fine- and coarse-grid temporal meshes.

MGRIT has both linear [11] and Full Approximation Scheme (FAS) variants [14], suitable for linear and nonlinear $\Phi$, respectively, and these variants differ in the way that coarse-grid correction is performed. First, we consider the linear case, in which we assume $\Phi$ is a linear operator, such that $\Phi(\mathbf{u}_{i-1}) = \Phi\mathbf{u}_{i-1}$ is a matrix-vector product with matrix $\Phi$. We also write the action of the space-time operator $\mathcal{A}(\cdot)$, meaning it operates simultaneously on all spatial coordinates and time points of the solution $\mathbf{u}$, as a matrix-vector product $\mathcal{A}(\mathbf{u}) = A\mathbf{u}$. From (3), we obtain an equivalent system written only in terms of the $C$-points, which we write in block-matrix form as

where $\widetilde{\mathbf{w}}_{jm} = \mathbf{w}_{jm} + \Phi\mathbf{w}_{jm-1} + \cdots+ \Phi^{m-1}\mathbf{w}_{(j-1)m+1}$; note that $\Phi^m$ plays the role of time-stepping operator for this coarse system. While system (15) has a factor of $m$ fewer equations than (3), it is no less expensive to solve because (naively) computing the action of $\Phi^m$ is $m$ times as expensive as computing that of $\Phi$. So, instead, MGRIT solves an approximation to the ideal coarse-grid system (15) arising from replacing the ideal coarse-grid time-stepping operator $\Phi^m$ with an approximation:

Then, MGRIT solves the approximate system $\mathcal{B}_{\Delta}\mathbf{u}_{\Delta} = \mathbf{w}_{\Delta}$, where the approximation $\mathcal{B}_{\Delta} \approx\mathcal{A}_{\Delta}$ replaces the $\Phi^m$ terms in (15) by $\Psi$.

Throughout the paper, we assume MGRIT is initialized with an initial approximate solution vector $\mathbf{v}^{(0)} = [\mathbf{u}_0\ \ \mathbf{u}_0\ \ \cdots\ \ \mathbf{u}_0]^T$ which simply extends the initial guess $\mathbf{u}_0$ to all other time points, which will be important to maintain progress in the restarted method. Letting $\mathbf{r}^{(0)}$ denote the initial MGRIT residual (14), we define the linear MGRIT algorithm in Algorithm 3.1 with the standard MGRIT pre- and post-relaxation schemes [11, 9], where $R_I(\cdot)$ stands for the injection operator. Note that we omit the post F-relaxation that follows the correction at C-points, since it is made redundant by the FCF-relaxation in the following iteration. We note that while the MGRIT convergence criterion is given in terms of the residual

Algorithm 3.1 2-level MGRIT($\Phi$, $\Psi$, $\mathbf{u}_0$, $\mathbf{w}$, $\varepsilon$).

1. Initialize $\mathbf{u}$ given initial guess $\mathbf{u}_0 \in\mathbb{R}^{N}$: $\mathbf{u} = [\mathbf{u}_0 \quad\cdots\quad\mathbf{u}_0]^T \in\mathbb{R}^{N(N_t+1)}$.

2. repeat

3. Apply FCF-relaxation to $A\mathbf{u} = \mathbf{w}$ using $\Phi$.

4. Inject residual to the coarse grid: $\mathbf{r}_{\Delta} = R_I(\mathbf{w} - A\mathbf{u})$.

5. Solve coarse-grid correction problem using $\Psi$: $\mathcal{B}_{\Delta}\mathbf{e}_{\Delta} = \mathbf{r}_{\Delta}$.

6. Correct approximate solution at C-points: $\mathbf{u}_{\Delta} \leftarrow\mathbf{u}_{\Delta} + \mathbf{e}_{\Delta}$.

7. until $\lVert\mathbf{r}\rVert/\lVert\mathbf{r}^{(0)}\rVert\leq\varepsilon$.

(14), in [2 Lemma 5.1] we proved that reducing the norm of the initial residual by a factor $\varepsilon$ implies that the norm of the gradient is also reduced to a similar factor; i.e., $\lVert\nabla f(\mathbf{v}_{N_t})\rVert/\lVert\nabla f(\mathbf{u}_0)\rVert\approx\varepsilon$, where $\mathbf{v}_{N_t}$ denotes the MGRIT approximate solution at the final time point $N_t$.

Now let us consider the MGRIT coarse-grid correction in the more general FAS case, where $\Phi(\cdot)$ in (3) is not (necessarily) a linear operator. Let $\mathbf{v} \approx\mathbf{u}$ be an approximate solution of $A(\mathbf{u}) = \mathbf{w}$ in (3). After performing the FCF-relaxation step, both the fine-grid approximation $\mathbf{v}$ and nonlinear residual $\mathbf{w} - A(\mathbf{v})$ are injected to the coarse-grid:

The ideal FAS coarse-grid problem is then to solve the system

for the coarse-grid solution $\mathbf{u}_{\Delta}$, where similarly to (15),

The coarse-grid system (18) results from the elimination of F-points variables from the residual equation $A(\mathbf{u}) - A(\mathbf{v}) = \mathbf{r}(\mathbf{v})$. Similarly to the linear case, the MGRIT iteration solves an approximation of (18),

where $\mathcal{B}_{\Delta}(\cdot) \approx A_{\Delta}(\cdot)$ replaces the $\Phi^{m}(\cdot)$ terms in $A_{\Delta}(\cdot)$ by a computationally cheaper $\Psi(\cdot) \approx\Phi^{m}(\cdot)$.

We define the 2-level FAS-MGRIT iteration in Algorithm 3.2, once again omitting the redundant post-F relaxation. Similarly to the linear case, in [2 Lemma 5.1] we proved that $\lVert\mathbf{r}\rVert/\lVert\mathbf{r}^{(0)}\rVert\leq\varepsilon$ implies that $\lVert G_{sF}(\mathbf{v}_{N_t})\rVert/\lVert G_{sF}(\mathbf{u}_0)\rVert\approx\varepsilon$, where $G_{sF}(\cdot)$ denotes the generalized gradient of $F$, given by

Algorithm 3.2 2-level FAS-MGRIT($\Phi$, $\Psi$, $\mathbf{u}_0$, $\mathbf{w}$,$\varepsilon$).

1. Initialize $\mathbf{u}$ given initial guess $\mathbf{u}_0 \in\mathbb{R}^{N}$: $\mathbf{u} = [\mathbf{u}_0 \quad\cdots\quad\mathbf{u}_0]^T \in\mathbb{R}^{N(N_t+1)}$.

2. repeat

3. Apply FCF-relaxation to $\mathcal{A}(\mathbf{u}) = \mathbf{w}$ using $\Phi$.

4. Inject approximate solution and residual to coarse grid: $\mathbf{u}_{\Delta} = R_I(\mathbf{u})$, $\mathbf{r}_{\Delta} = R_I(\mathbf{w} - \mathcal{A}(\mathbf{u}))$.

5. Solve coarse-grid correction problem: $\mathcal{B}_{\Delta}(\hat{\mathbf{u}}_{\Delta}) = \mathcal{B}_{\Delta}(\mathbf{u}_{\Delta}) + \mathbf{r}_{\Delta}$.

6. Compute coarse-grid error approximation: $\mathbf{e}_{\Delta} = \hat{\mathbf{u}}_{\Delta} - \mathbf{u}_{\Delta}$.

7. Correct approximate solution at C-points: $\mathbf{u}_{\Delta} \leftarrow\mathbf{u}_{\Delta} + \mathbf{e}_{\Delta}$.

8. until $\lVert\mathbf{r}\rVert/\lVert\mathbf{r}^{(0)}\rVert\leq\varepsilon$.

for some $\mathbf{v}_i \in\mathbb{R}^{N}$. Note that $P_{sg} \circ G_{sf}$ is the same as the proximal gradient operator (11); in particular, $P_{sg} \circ G_{sf} = I - sG_{sF}$.

While we have described 2-level MGRIT methods, multilevel MGRIT schemes are obtained by recursively applying the 2-level schemes. In our numerical tests, we apply the standard MGRIT scheme for MP1 (7) and the FAS-MGRIT scheme for MP2 (10). For both MP1 and MP2, MGRIT is iterated until the 2-norm of the residual (14) is reduced by 8 orders of magnitude from its initial value.

The restarted MGRIT algorithm for optimization problems (R-MGRIT)

We now describe the R-MGRIT algorithm, Algorithm 3.3. As introduced in § 1, instead of approximating a single exact trajectory yielded by the underlying optimization method (1), we treat the initial guess $\mathbf{u}_0$ of the all-at-once system $\mathcal{A}(\mathbf{u}) = \mathbf{w}$ (3) as a parameter to be reinitialized. Each new trajectory is approximated using one MGRIT iteration.¹ Thus, R-MGRIT can be interpreted as a restarting scheme that approximates multiple trajectories, each hopefully starting from a point closer to the true optimum. It is important to propagate $\mathbf{u}_0$ across the C-points in the new trajectory because a random initialization of the C-points may reintroduce large gradients. We note that one R-MGRIT iteration costs the same as one MGRIT iteration from Algorithms 3.1 and 3.2.

Algorithm 3.3 R-MGRIT($\Phi$, $\Psi$, $\mathbf{u}_0$, $\mathbf{w}$,$\varepsilon$,$N_t$).

1. repeat

2. Initialize $\mathbf{u}$ given initial guess $\mathbf{u}_0 \in\mathbb{R}^{N}$: $\mathbf{u} = [\mathbf{u}_0 \quad\cdots\quad\mathbf{u}_0]^T \in\mathbb{R}^{N(N_t+1)}$.

3. Run 1 iteration of MGRIT on all-at-once system $\mathcal{A}(\mathbf{u}) = \mathbf{w}$.

4. Set $\mathbf{u}_0 \leftarrow\mathbf{u}_{N_t}$, $\mathbf{w}_0 \leftarrow\mathbf{u}_0$.

5. until $\lVert G_{sF}(\mathbf{u}_{N_t})\rVert/\lVert G_{sF}(\mathbf{u}_0)\rVert\leq\varepsilon$.

Asymptotic final-point analysis of MGRIT and R-MGRIT

In the following, we compare standard linear MGRIT and R-MGRIT theoretically. We use the index $k$ to denote the number of iterations for standard MGRIT and the number of restarts (i.e., trajectories) for R-MGRIT. Based on linear MGRIT error propagation theory [9, 25], we demonstrate that if the approximate solutions at the final time point $N_t = N_T m$ of the $k$th MGRIT iteration and $k$th R-MGRIT trajectory are approximately the same, converging to the same asymptotic limit as $N_T \to\infty$.

Consider an iterative method of the form (1) with a linear operator $\Phi$. Furthermore, we assume $\mathbf{w}_i=\mathbf{c}$ for $i>0$, such that, given an initial guess $\mathbf{u}_0$, we have

Writing the $i$th iterate yielded by (1) in terms of the initial guess $\mathbf{u}_0$ yields

Since $\mathbf{u}^*$ is a fixed point of the iteration, i.e., $\mathbf{u}^*=\Phi\mathbf{u}^*+\mathbf{c}$, we may shift coordinates by defining $\mathbf{v}_j=\mathbf{u}_j-\mathbf{u}^*$, with the resulting iteration being $\mathbf{v}_i=\Phi\mathbf{v}_{i-1}$. Thus, without loss of generality, we simplify the above iteration to

For concreteness, we may consider $\Phi:=I-sA_1$ in the case of gradient descent (8) for MP1, (7).

In the following, we first explicitly determine the MGRIT error propagation matrix for $k$ MGRIT iterations. Assuming a constant initialization strategy, we then characterize the MGRIT and R-MGRIT approximate solutions, for the case of $k$ MGRIT iterations and $k$ R-MGRIT trajectories with 1 MGRIT iteration per trajectory, at the final time point in terms of the MGRIT parameters. Finally, we prove that the MGRIT and R-MGRIT end-point approximate solutions converge to the same asymptotic limit for fixed $k$ as $N_T\to\infty$.

MGRIT error propagation after $k$ iterations

Let $\Phi$ and $\Psi$ denote the fine- and coarse-grid operators. After $k$ iterations, the C-point error propagation block matrix $\mathcal{E}$ propagates the error at C-points via

where $\mathbf{e}_{\mathrm{C}}^{(0)}$ and $\mathbf{e}_{\mathrm{C}}^{(k)}$ denote the C-point initial error and the C-point error after $k$ iterations, respectively. For MGRIT with FCF-relaxation, $\mathcal{E}$ is given by [9, 25]

or block-wise, by

By taking the $k$th power of $\mathcal{E}$ in (26), we explicitly determine the error propagation matrix $\mathcal{E}^k$ for $k$ MGRIT iterations. We now introduce an assumption to hold uniformly throughout this section; while not strictly necessary for all of the results, applying the assumption uniformly aids in the presentation.

Proof. By induction on $k$. The base case $k=1$ trivially matches (26) since $\binom{i-j-2}{0}=1$. Now, assume the formula holds for $k-1$. We compute the $(i,j)$-th block of $\mathcal{E}^k=\mathcal{E}\mathcal{E}^{k-1}$ via the block matrix product

By (28) holding for $\mathcal{E}$ and $\mathcal{E}^{k-1}$, the addends in this sum are non-zero only if $(\mathcal{E})_{i,p}\ne0 \Longleftrightarrow i-p\ge2$ and $(\mathcal{E}^{k-1})_{p,j}\ne0 \Longleftrightarrow p-j\ge2(k-1)$. Adding these two inequalities, we have that $(\mathcal{E}^k)_{i,j}=0$ if $i-j<2k$. Otherwise, for $i\ge j+2k$, substituting (28) into (29) yields

where we used the commutativity of $\Phi$ and $\Psi$. Recall the hockey-stick (or Fermat’s) identity for non-negative integers $n\ge r$,

Applying the change of variables $q=p-j-k$, we get that

where the final equality comes from applying (30). ∎

From Proposition 4.2, we generalize [9 Lemma 3.2] to determine the C-point errors after $k$ iterations.

Proof. The error at the $i$-th C-point, $e_{im}^{(k)}$, is determined by the product of the $i$-th block row of $\mathcal{E}^{k}$ with the initial C-point error vector $\left[e_{0}^{(0)}\ e_{m}^{(0)}\ \cdots\ e_{N_{T}m}^{(0)}\right]^{T}$:

The statement then follows from substituting (28) into (34). $\square$

Characterizing (R-)MGRIT final-point approximate solutions

We now characterize the approximate solutions at the final time point $N_{T}$ for standard MGRIT and R-MGRIT. For both methods, we use the constant initialization strategy described in Algorithms 3.1 and 3.3.

First, consider standard MGRIT applied to a trajectory with initial guess $\mathbf{u}_{0}$. Since the exact solution at the $j$th C-point is $\mathbf{u}_{jm}=\Phi^{jm}\mathbf{u}_{0}$, the constant initialization gives an initial error $\mathbf{e}_{jm}^{(0)}=(\Phi^{jm}-I)\mathbf{u}_{0}$. Substituting this into Lemma 4.3, the error after $k$ MGRIT iterations at the final time point is $\mathbf{e}_{N_{T}}^{(k)}=E_{N_{T}}^{(k)}\mathbf{u}_{0}$, where

Since the exact solution at the final time point is $\mathbf{u}_{N_{t}}=\Phi^{N_{t}}\mathbf{u}_{0}$, the corresponding MGRIT approximation is

where $\underbrace{\Phi^{N_{t}}-E_{N_{T}}^{(k)}}_{=:\widetilde{V}_{N_{t}}^{(k)}}$.

Next, consider R-MGRIT with one MGRIT iteration per trajectory. Since the same constant initialization strategy is used for each trajectory, the operator mapping its initial guess to its approximate final point is always the one-iteration MGRIT operator

Let $\widetilde{\mathbf{u}}_{0}^{(k)}$ denote the initial guess of the $k$th R-MGRIT trajectory, with $\widetilde{\mathbf{u}}_{0}^{(1)}=\mathbf{u}_{0}$. The approximate final point of this trajectory is $\widetilde{\mathbf{v}}_{N_{t}}^{(k)}=\widetilde{V}_{N_{t}}\widetilde{\mathbf{u}}_{0}^{(k)}$. By construction, this becomes the initial guess for the next trajectory, i.e., $\widetilde{\mathbf{u}}_{0}^{(k+1)}=\widetilde{\mathbf{v}}_{N_{t}}^{(k)}$. Applying this relation recursively, we obtain

Thus, standard MGRIT applies $k$ iterations to a single trajectory, yielding the final-point approximation $\mathbf{v}_{N_{t}}^{(k)}=V_{N_{t}}^{(k)}\mathbf{u}_{0}$, whereas R-MGRIT applies the same single-iteration final-point operator $k$ times, yielding $\widetilde{\mathbf{v}}_{N_{t}}^{(k)}=\left(\widetilde{V}_{N_{t}}\right)^{k}\mathbf{u}_{0}$. In the following sections, we compare these two approximations as $N_{T}\to\infty$.

Error operators

Let $Q=\left[\mathbf{q}_{1}\ \mathbf{q}_{2}\ \cdots\ \mathbf{q}_{N}\right]$ be the unitary matrix that diagonalizes both $\Phi$ and $\Psi$ under Assumption 4.1, such that

is satisfied with $\phi_n$ and $\psi_n$ denoting the $n$th eigenvalues of $\Phi$ and $\Psi$, respectively. We denote by $u_{0,n}=\mathbf{q}_n^*u_0$ the $n$th modal component of the initial guess $u_0$, where $\mathbf{q}_n$ is the $n$th column of $Q$. Then, after an eigenvector expansion of $e_c^{(1)}=\mathcal{E}e_c^{(0)}$ on the basis of eigenvectors given by the columns of $Q$, it follows that the action of $\mathcal{E}$ on the $n$th mode is given by $\mathcal{E}_n\in\mathbb{C}^{N_T+1\times N_T+1}$, given by [9, 25]

Following similar steps, it also follows from Proposition 4.2 that the action of $E_{N_T}^{(k)}$ from (35) on the $n$th mode (i.e., the $n$th modal component of $u_0\in\mathbb{R}^N$) is

such that the modal equivalents of (36) and (38) are

and

respectively, where $v_{N_T,n}^{(k)}$ is the $n$th modal component of $v_{N_T}^{(k)}$, and $\widetilde{v}_{N_T,n}^{(k)}$ is the $n$th modal component of $\widetilde{v}_{N_T}^{(k)}$. We now establish the asymptotic modal error operator as $N_T\to+\infty$.

Proof. Applying the change of variables $p=N_T-2k-j$ reverses the summation order, rewriting the modal error operator from (41) as:

Let $M=N_T-2k$. We can split the summation into two distinct parts and evaluate their limits as $M\to\infty$. For the first part involving the $\phi_n$ term, let

To bound this sum, let $\rho= \max(\lvert\psi_n\rvert,\lvert\phi_n^m\rvert)$. By hypothesis, both $\lvert\phi_n\rvert< 1$ and $\lvert\psi_n\rvert< 1$, which implies $\rho< 1$. Taking the absolute value yields:

Using the hockey-stick identity (30), the sum of the binomial coefficients evaluates to $\binom{M+k}{k}$, giving the bound $\lvert S_M\rvert\le\binom{M+k}{k}\rho^M$. The expression $\binom{M+k}{k}$ is a polynomial in $M$ of degree $k$, while $\rho^M$ is an exponentially decaying term. Since exponential decay dominates polynomial growth as $M \to\infty$, we have $\lim_{M\to+\infty}\binom{M+k}{k}\rho^M=0$, meaning $\lim_{M\to+\infty}S_M=0$.

The second part of the summation in (45), $-\sum_{p=0}^{M}\binom{p+k-1}{k-1}\psi_n^p$, is a partial sum of the negative binomial series. As $M\to\infty$, this yields the Maclaurin series expansion for $-(1-\psi_n)^{-k}$:

Because $\lvert\psi_n\rvert<1$, this series converges absolutely. Substituting these two limits back into the original expression yields

which recovers $E_{\infty,n}^{(k)}$ in (44).

Asymptotic final-point (R-)MGRIT approximate solutions

We now show that for any fixed $k$, both MGRIT with $k$ iterations and R-MGRIT with 1 MGRIT iteration on $k$ consecutive trajectories converge to the same final-point approximate solution, asymptotically for $N_T\to+\infty$.

Proof. We first establish the result for each modal component. For $n=1,\ldots,N$, define $\widetilde{V}_{\infty,n}:=(\phi_n^m-\psi_n)\phi_n^m/(1-\psi_n)$. By Lemma 4.4, we have

Since $\lvert\phi_n\rvert<1$, we also have $\phi_n^{N_T}\to0$ as $N_T\to+\infty$. Thus, using (42) and (43), we obtain

Hence, the MGRIT and R-MGRIT final-point approximations converge to the same limit in every mode.

To recover the full approximate solution vectors, recall that the columns $\{q_n\}_{n=1}^{N}$ of $Q$ form an orthonormal basis. Since $N$ is finite, the modal limits can be combined to give

Finally, by the simultaneous diagonalization of $\Phi$ and $\Psi$,

Therefore, the common limit is $\widetilde{V}_{\infty}^{(k)}u_0$, as claimed. ∎

By the triangle inequality, we can bound the finite-time difference between the MGRIT and R-MGRIT approximate solutions by their differences with respect to their shared asymptotic limit:

In our numerical experiments from § 5, we observe that these right-hand side terms are approximately zero for a sufficiently large $N_t$, which implies the left-hand side difference is also close to zero. Namely, we will infer this in Figure 3 from the fact that the gradient at the C-points stalls for both MGRIT and R-MGRIT, which is consistent with $\widetilde{v}_{N_T}^{(k)}\approx v_{\infty}^{(k)}\approx v_{N_T}^{(k)}$, with similar behavior observed in the nonlinear, nonsmooth MP2 experiments in Figure 5 and Figure 7.

On a truncated R-MGRIT strategy

We now motivate the idea of shortening R-MGRIT trajectories, by showing that the approximate endpoint solution obtained with one MGRIT iteration depends only weakly on the length of the trajectory.

Proof. Under constant initialization, Lemma 4.3 gives

Separating the final term yields the recurrence

Since $\widetilde{v}_{jm}^{(k)}=\Phi^{jm}\widetilde{u}_0^{(k)}-\widetilde{e}_{jm}^{(k)}$, we obtain

Using $\widetilde{v}_m^{(k)}=\Phi^m\widetilde{u}_0^{(k)}$ and iterating gives

Since $\rho(\Psi)<1$, the Neumann series converges, giving

Subtracting this limit from the preceding finite sum and using the geometric series identity gives (58). Finally, since $\Psi$ is unitarily diagonalizable, $\lVert\Psi\rVert_2=\rho(\Psi)<1$, which yields (59). □

Lemma 4.6 provides a motivation for shortening R-MGRIT trajectories. Although the approximate solution along a single trajectory generally converges to $\widetilde{v}_{\infty}^{(k)}\ne0$, its difference from this limit decays geometrically at a rate governed by $\rho(\Psi)$. Thus, sufficiently far along a trajectory, computing additional C-points has little effect on its approximate endpoint. Since R-MGRIT uses only this endpoint to initialize the next trajectory, much of this additional computation may be unnecessary.

More precisely, let $\widetilde{V}_{jm}$ denote the endpoint operator for a trajectory with $j$ coarse-grid intervals. The lemma implies that $\widetilde{V}_{jm}\to\widetilde{V}_{\infty}$ as $j\to\infty$, and hence $(\widetilde{V}_{\infty})^k u_0\to\widetilde{V}_{\infty}^k u_0=v_{\infty}^{(k)}$ for fixed $k$. Therefore, sufficiently long trajectories can be shortened while retaining approximately the same final-point solution after $k$ restarts, potentially reducing the computational cost of R-MGRIT.

Motivated by this observation, we consider truncating each R-MGRIT trajectory to $\overline{N}_t=\gamma N_t$ fine-grid intervals, where $\gamma\in(0,1)$ and $\overline{N}_t$ is rounded to a multiple of $m$ when necessary. For sufficiently large $\overline{N}_t$, the preceding analysis suggests that, for fixed $k$,

The theoretical result does not specify how short the trajectories can be made, nor does it guarantee that truncation preserves the number of restarts required for convergence. We investigate these questions numerically in § 5, using $\gamma=0.2$.

Numerical results for MGRIT vs. R-MGRIT

For our numerical experiments we consider the two model problems, MP1, (7), and MP2, (10), with the optimization methods described in subsections 2.1 and 2.2 for MP1 and MP2, respectively. For each problem, we first determine $N_t$ via a sequential fine-grid solve, then construct our fine-grid mesh of $N_t+1$ points with spacing $\Delta t=1/L$. We choose $10^{-8}$ as the relative tolerance for all iterative methods.

For R-MGRIT, as proposed in subsection 4.5, we consider truncating the fine-grid time domain by a fixed factor $\gamma=0.2$. In this case, we solve all-at-once systems of the form (4) corresponding to a fifth of the number of iterations of the underlying optimization method. We illustrate convergence results for a particular set of choices. for the parameters $N$ (dimension of the Laplacian matrix), $m$ (coarsening factor) and $\hat{\ell}$ (number of levels). Later, we present results for different values of $m$, $\hat{\ell}$ and $N$ to analyze how performance is affected by the choices of these parameters.

The initial MGRIT guess $\mathbf{u}_0 \in\mathbb{R}^N$ is chosen to be random with components in $[0,1]$, after which the condition $\mathbf{u}_0 = \mathbf{w}_0$ is enforced on the all-at-once system (3). For both model problems, we fix a step size $s^{(\ell)} = m^{\ell-1}/L$ for the $\ell$th level operator, $\ell= 1,\ldots,\hat{\ell}$, where $m$ denotes the coarsening factor and $\hat{\ell}$ the maximum number of levels; note that we use the same operator $\Psi$ (with the level-appropriate step size $s^{(\ell)}$) on all levels $\ell\ge2$.

Numerical experiments for gradient descent applied to MP1

We solve MP1, (7), using gradient descent, (8). For this example, the 1D Laplacian matrix $A_1 \in\mathbb{R}^{N \times N}$ in (5) with 2-norm $L$ is taken with $N = n = 80$, the coarsening factor is chosen to be $m = 64$ and the number of levels to be $\hat{\ell} = 2$. We also choose $\mathbf{b} \in\mathbb{R}^N$ to be random with components in $[0,1]$. First, we solve the problem sequentially, with 30,016 iterations required for the convergence of the gradient descent algorithm (8). MGRIT then solves (3) using Algorithm 3.1. We then consider truncating the total number of fine-grid intervals from $N_t$ to $\overline{N}_t = \gamma N_t$. For a truncation factor $\gamma= 0.2$, we have $\overline{N}_t = 6,080$.

Convergence for the sequential and MGRIT solves is illustrated in Figure 2, in which we observe that MGRIT requires 8 iterations for convergence, while R-MGRIT only requires 6 trajectories, so a total of 6 MGRIT iterations. We plot the gradient at the C-point approximate solutions of each iteration of MGRIT and R-MGRIT in Figure 3, in which we observe that the truncated R-MGRIT possesses a similar convergence behavior as the full-time domain R-MGRIT solve. Therefore, by truncating $N_t$ to a sufficiently large (yet relatively small) $\overline{N}_t = \gamma N_t$, we obtain convergence in at most as many iterations as MGRIT.

Figure 2(b) shows that MGRIT and R-MGRIT converge essentially with the same speed and the cost per iteration is the same, so there is no direct benefit in considering multiple inexact trajectories instead of converging MGRIT accurately on a single trajectory. However, the middle panel of Figure 3 suggests that running R-MGRIT with shorter trajectories may provide a substantial cost benefit. This approach is explained in the $\gamma$-truncated numerical results in this section, and gain in computational cost is discussed in § 6.

Although the numerical results illustrated in Figure 3 are only presented for the norm of the gradient, the theory from § 4 allows us to interpret the results in terms of the approximate solutions. Theorem 4.5 predicts that the difference between the MGRIT and R-MGRIT final-point approximate solutions for each fixed iteration index $1 \leq k \leq6$ will be small as the number of C-points becomes sufficiently large for both to approximate their common asymptotic limit of $\mathbf{v}_{\infty}^{(k)}$, and Figure 3 indicates that this happens much earlier than $N_T$. Indeed, Figure 3 suggests that (60) holds for $\overline{N}_T = \gamma N_T$ with $\gamma\leq0.2$ with $\gamma$-truncated R-MGRIT attaining a similar gradient reduction factor to MGRIT and full-time R-MGRIT.

For $k \geq7$, we observe that $N_T$ is no longer large enough for the $\mathbf{v}_{N_T}^{(k)} \approx\mathbf{v}_{\infty}^{(k)}$ approximation to hold and the MGRIT final-point gradient stalls, with a final gradient reduction factor of approximately $10^{-8}$ as predicted by [2 Lemma 5.1]. By approximating new target trajectories rather than the initial one, R-MGRIT iterations of index $k \geq7$ continue to reduce the gradient by approximately the same linear rate, as first illustrated in Figure 2. Thus, R-MGRIT does not require an accurate estimation of $N_t$, the number of iterations required for convergence of gradient descent (8) with the prescribed $10^{-8}$ accuracy threshold.

In Table 1, we present MGRIT and R-MGRIT convergence results for MP1 (7) with different coarsening factors $m$, number of MGRIT levels $\hat{\ell}$ and different sizes for the matrix $A_1$ in (5); we note that all R-MGRIT results hold identically for a truncation factor $\gamma= 0.2$ (not shown). Having established in [2] that the number of iterations does not change much as $\hat{\ell}$ increases, here we only present results for the smallest ($\hat{\ell} = 2$) and largest ($\hat{\ell} = \ell$) possible number of levels given the number of time points $N_t$ and the coarsening factor $m$.

	$\hat{\ell}=2$	$\hat{\ell}=\ell$
$m=4$	9(9)	9(9)
$m=16$	8(7)	8(7)
$m=64$	8(6)	8(6)
$m=256$	7(5)	—
$m=4$	9(9)	9(9)
$m=16$	8(7)	8(7)
$m=64$	8(6)	8(6)
$m=256$	8(6)	—
$m=1024$	7(5)	—

Numerical experiments for proximal gradient applied to MP2

Numerical results for MP2-1D

We now consider solving the EOP, MP2 (10), using proximal gradient method (11) with penalty parameter $\lambda= 900$.

For MGRIT, we choose the alternating proximals operator (13) to be the coarse-grid operator. We present one MP2-1D and one MP2-2D example, similar to those in [2].

For the MP2-1D example, we consider MP2 in (10) with $d=1$ and $N=256$, with MGRIT parameters $m=256$ and $\hat{\ell}=2$. The obstacle $\phi$ is obtained by discretizing $\phi(x)=\max\{0,\sin x\}$ for $x \in\Omega:= [0,3\pi]$ on a uniform mesh of $N+2$ points.

We again start by presenting results for the sequential and (R-)MGRIT convergence in Figure 4. In this example, the proximal gradient method requires 109,888 iterations to reach the prescribed halting tolerance, i.e., that the generalized gradient of $F$, $G_{sF}$ in (21), be reduced by 8 orders of magnitude; we extend the number of time points to $N_t=110,081$ such that $N_T=N_t/m=431$ is an integer. In Figure 4, we observe that MGRIT reaches the prescribed residual tolerance in 10 iterations, while R-MGRIT reaches it in 9 iterations both when using the full $N_t$ number of time points and when truncating the time domain to a fifth of the number of time points. The convergence speed of all 3 methods in Figure 4 is very similar. We present the C-point gradients of each iteration in Figure 5, in which we observe similar behavior in the gradient norm to what we observed in Figure 3 for MP1.

Numerical results for MP2-2D

For the MP2-2D example, we consider MP2 in (10) with $d=2$ and $N^2=64^2$, with MGRIT parameters $m=64$ and $\hat{\ell}=2$. The obstacle $\phi$ is obtained by discretizing $\phi(x,y)=\max{0,\sin x}\max{0,\sin y} for $(x,y) \in\Omega:= [0,3\pi]^2$ on a uniform mesh of $N+2$ points in each spatial dimension. We present the sequential and MGRIT convergence results in Figure 6, noting that the sequential solver required 13,256 iterations to converge; for a coarsening factor $m = 64$, we then work with $N_t = 13,312$, such that $N_T = N_t/m = 208$ denotes the number of coarse time intervals. We present C-point gradients in Figure 7, in which we observe that $\gamma$-truncated R-MGRIT (with $\gamma= 0.2$) converged in 9 iterations, while MGRIT and full-time R-MGRIT converged in 10 iterations. As in the MP2-1D example, we observe that $\gamma$-truncated R-MGRIT required fewer iterations to reduce the gradient by a factor of $10^{-8}$ than MGRIT and full-time R-MGRIT.

While the theory from § 4 is only presented for linear MGRIT problems, we observe generally similar results for the nonlinear MP2, (10), in both one and two dimensions. In particular, in both Figure 4 and Figure 6 we observe that R-MGRIT reduces the gradient by a factor of $10^{-8}$ in the same 9 iterations required for MGRIT to achieve a similar reduction factor, with R-MGRIT continuing to reduce the gradient in subsequent iterations while MGRIT’s final-point generalized gradient stalls. Although the approximation from (60) for $N_T = \gamma N_T$ with $\gamma= 0.2$ is observed to be less accurate in that the convergence to the final gradient is slower in Figure 5 and Figure 7 than in Figure 3, as can be inferred from the rightmost panel of the respective figures, this is advantageous for $\gamma$-truncated R-MGRIT, which is able to reduce the final-point generalized gradient (21) at a slightly faster rate because the gradient undershoots with a slow increase to the final value.

In Table 2 and Table 3, we present convergence results for MP2, (10), with different coarsening factors $m$, number of MGRIT levels $\hat{\ell}$ and different sizes for the matrices $A_1$ in (5) for MP2-1D and $A_2$ in (6) for MP2-2D. As in Table 1, we only present results for the smallest ($\hat{\ell}=2$) and largest ($\hat{\ell}=\ell$) permissible number of levels. Comparing results for MGRIT, R-MGRIT with full time domain, and $\gamma$-truncated R-MGRIT, we observe that the truncated variant achieves comparable gradient reduction in a similar number of iterations.

	$\hat{\ell}=2$	$\hat{\ell}=\ell$		$\hat{\ell}=2$	$\hat{\ell}=\ell$
$m=4$	10(9)[10]	12(13)[11]	$m=4$	10(8)[8]	12(9)[10]
$m=16$	9(8)[8]	10(9)[9]	$m=16$	11(9)[9]	11(10)[9]
$m=64$	10(9)[9]	10(11)[9]	$m=64$	9(7)[9]	10(9)[9]
$m=256$	10(9)[9]	10(9)[9]	$m=256$	9(8)[7]	—
$m=1024$	11(10)[9]	—	$m=1024$	10(9)[9]	—

	$\widehat{\ell}=2$	$\widehat{\ell}=\ell$		$\widehat{\ell}=2$	$\widehat{\ell}=\ell$
$m=4$	$10(10)[10]$	$14(13)[13]$	$m=4$	$9(9)[9]$	$11(11)[10]$
$m=16$	$11(10)[10]$	$12(12)[11]$	$m=16$	$9(8)[8]$	$9(9)[9]$
$m=64$	$10(10)[9]$	$10(10)[9]$	$m=64$	$9(8)[8]$	$9(9)[8]$
$m=256$	$9(9)[7]$	—	$m=256$	$9(9)[8]$	—

R-MGRIT parallel speedup estimates

In § 5, we empirically observed that $\gamma$-truncated R-MGRIT achieves gradient reduction comparable to MGRIT (in the same or fewer number of iterations) with the R-MGRIT time domain truncated to a fifth of the number of time points required for MGRIT convergence to the same accuracy threshold. Because the primary motivation for employing MGRIT is to reduce the wall-clock time required to solve the sequential problem, we now formalize how this temporal truncation translates into practical parallel speedup. In the following, we assume that one R-MGRIT iteration has the same computational cost as one MGRIT iteration as described by the speedup formulas derived in [2]. Thus, in this section we discuss 2-level R-MGRIT with a $\gamma$-truncated temporal domain, with $\gamma=1$ denoting the non-truncated R-MGRIT (or equivalently, MGRIT) speedup estimates.

2-level $\gamma$-truncated R-MGRIT parallel speedup

In the following, recall that $N_t$ denotes the number of time intervals on the fine-grid level $\ell=1$, and let $t_f$ denote the wall-clock time required to apply the fine-grid time-stepping operator $\Phi$ once. Similarly, recall that $N_T$ denotes the total number of time intervals on the coarse-grid level $\ell=2$, and let $t_c$ denote the wall-clock time required to apply the coarse-grid time-stepping operator $\Psi$ once. Let $m$ denote the coarsening factor, given by $m=N_t/N_T$. Let $N_p$ denote the number of processors. By full parallelization we mean that $N_p=N_T=N_t/m$, in which case each coarse-grid step is allocated to its own processor, along with $m$ fine-grid steps. Let $N_{it}$ denote the number of iterations required for convergence of standard MGRIT given the coarsening factor $m$. Finally, we let $\alpha:=t_c/t_f$ denote the ratio of wall-clock time required to compute one coarse-grid step relative to one fine-grid step.

Given the 2-level parallel speedup estimate formula from [2], which is based on similar formulas [17, 24] for the parareal method [15], the speedup formula for this truncated time domain is

where $\overline{N}_{it}$ is the number of iterations required for R-MGRIT convergence on a fine grid with $\overline{N}_T=\overline{N}_T m$ intervals. Note that $\overline{N}_{it}$ still equals both the number of R-MGRIT trajectories and the total number of MGRIT iterations used by R-MGRIT. Given some $\gamma\in(0,1]$, we substitute $\overline{N}_T=\gamma N_t/m=\gamma N_T$ into (61), which yields

Note that if $\overline{N}_{it}=N_{it}$, this formula is identical to the 2-level MGRIT speedup formula from [2], but with the operator cost ratio $\alpha$ scaled down by $\gamma$; letting $\gamma=1$ recovers the 2-level MGRIT formula exactly. By setting the derivative of the denominator of (62) with respect to $m$ equal to zero for a fixed $\gamma$, we maximize $\overline{S}_2(m,\gamma)$ to find the optimal coarsening factor for the $\gamma$-truncated R-MGRIT solve:

Let $m^*:=m_1^*$ denote the MGRIT optimal coarsening factor. Then, the relation between the truncated R-MGRIT and the MGRIT optimal coarsening factors satisfies

By substituting $m^*$ and $m_\gamma^*$ into (62), we find that the ratio between the $\gamma$-truncated R-MGRIT speedup and the standard MGRIT (or non-truncated R-MGRIT) speedup is given by

In practice, we expect that $\overline{N}_{it}\leq N_{it}$ for many small $\gamma$, as seen in Table 1 to Table 3 for $\gamma=0.2$. Therefore, the speedup is increased by a factor of $N_{it}/(\overline{N}_{it}\sqrt{\gamma})$ while the number of processors required for full parallelization is also reduced proportionally to $\sqrt{\gamma}$:

For example, if $\gamma=1/4$ and $\overline{N}_{it}=N_{it}$, we achieve twice the speedup with half of the number of processors required for the standard MGRIT solve.

3-level $\gamma$-truncated R-MGRIT parallel speedup

We now extend the previous 2-level discussion to 3-level methods based on the 3-level speedup estimates from [2]. Following reasoning similar to the 2-level case, the full-parallelization 3-level $\gamma$-truncated speedup formula is given by

By setting the derivative of the denominator in (67) with respect to $m$ equal to zero to maximize $\overline{S}_3(m,\gamma)$, we obtain the optimal coarsening factor $m_\gamma^*$ given by

where $m^*:=m_1^*$ denotes the 3-level MGRIT optimal coarsening factor. As in the 2-level case, truncating R-MGRIT by a factor $\gamma$ increases the speedup and reduces the processor requirement relative to MGRIT, this time proportionally to $\sqrt[3]{\gamma}$ and $\gamma^{2/3}$, respectively:

For example, if $\gamma=1/8$ and $N_{it}=\overline{N}_{it}$, the predicted speedup is doubled while the processor requirement is quartered relative to the full-time 3-level MGRIT solve.

R-MGRIT parallel speedup estimates for MP1 and MP2

We now present the R-MGRIT parallel speedup estimates for MP1 and MP2 in the following tables. First, we compute the optimal $m^{*}$ for $\hat{\ell}=2$ and $\hat{\ell}=3$ by (63) and (68), respectively, given an empirically estimated $\alpha=t_{c}/t_{f}$. We choose $m=m^{*}$ as the coarsening factor, rounding it to the nearest integer if needed, for the MGRIT solves. By running the MGRIT algorithm of our choice, namely Algorithm 3.1 for MP1 and Algorithm 3.2 for MP2, with $m=m^{*}$ as the coarsening factor, we empirically obtain $N_{it}^{*}$ as the number of iterations required for convergence with the theoretically predicted optimal coarsening factors, and we also estimate $\alpha$. We then use $m=m^{*}$ and $N_{it}=N_{it}^{*}$ in the speedup formulas (62) and (67) for 2 and 3 levels, respectively. We then do the same for truncated R-MGRIT, for which we consider a truncation factor $\gamma=0.2$, as in § 5. Note that for R-MGRIT, the $N_{it}^{*}$ reported in the following tables corresponds to $\overline{N}_{it}$ in the discussion from § 6.

In Table 4 to Table 6, we observe that truncating the time domain to a fifth of the number of time points required for standard MGRIT convergence provides a substantial improvement of the standard MGRIT speedup, in accordance with the theory developed in § 6. Notably, our results here suggest that R-MGRIT converges in less parallel wall-clock time by employing only a small fraction of $N_{t}$ on each trajectory, the number of iterations required for the underlying optimization method to converge to a sufficiently accurate approximation of the true optimum.

	$\hat{\ell}$	$m^{*}$	$N_{it}^{*}$	$S_{\hat{\ell}}^{*}$	$N_{p}$	$\hat{\ell}$	$m^{*}$	$N_{it}^{*}$	$S_{\hat{\ell}}^{*}$	$N_{p}$
MGRIT	2	167	8	3.37	109	2	229	8	4.09	132
	3	21	8	6.20	860	3	26	8	7.81	1155
R-MGRIT ($\gamma=0.2$)	2	74	6	10.05	49	2	103	6	12.20	59
	3	13	6	14.12	278	3	15	8	13.36	401

	$\hat{\ell}$	$m^*$	$N_{it}^*$	$S_{\hat{\ell}}^*$	$N_p$	$\hat{\ell}$	$m^*$	$N_{it}^*$	$S_{\hat{\ell}}^*$	$N_p$
MGRIT	2	294	9	1.25	46	2	403	9	3.19	116
	3	20	10	1.06	663	3	30	9	4.91	1548
R-MGRIT ($\gamma= 0.2$)	2	131	8	3.16	21	2	180	8	8.05	52
	3	12	11	1.65	221	3	18	9	8.39	516
	(a) $N^2 = 64^2$, $N_t = 13,256$					(b) $N^2 = 128^2$, $N_t = 46,431$

We note that truncating more aggressively, e.g. choosing $\gamma=0.125$ instead of $\gamma=0.2$, did not significantly change the predicted speedup. We remark that an optimal truncation factor $\gamma$ can not be determined a priori since we cannot estimate $N_{it}^{*}$ as a function of $\gamma$ and the other parameters. However, our numerical results suggest $\gamma=0.2$ is a good baseline parameter choice for our model problems.

Applying Anderson-acceleration to R-MGRIT

In § 4 and § 5, we observed that R-MGRIT is not expected to achieve faster gradient reduction than MGRIT. Thus, in an attempt to obtain faster gradient reduction over MGRIT, we consider applying acceleration techniques to R-MGRIT. Namely, we consider applying Anderson acceleration [1, 8] to accelerate R-MGRIT. To this end, we define $\mathcal{F}$ as the map that takes an initial guess and returns the approximate final point obtained after one MGRIT iteration on the corresponding optimization trajectory. Since R-MGRIT uses this final point to initialize the next trajectory, we can write R-MGRIT as the fixed-point iteration

	$\hat{\ell}$	$m^*$	$N_{it}^*$	$S_{\hat{\ell}}^*$	$N_p$	$\hat{\ell}$	$m^*$	$N_{it}^*$	$S_{\hat{\ell}}^*$	$N_p$
MGRIT	2	363	10	7.56	303	2	849	10	10.88	436
	3	39	10	20.65	2818	3	59	9	33.60	6265
R-MGRIT ($\gamma= 0.2$)	2	162	9	18.80	136	2	380	7	34.79	195
	3	23	10	35.30	956	3	35	8	64.63	2113
	(a) $N = 256$, $N_t = 109,888$					(b) $N = 512$, $N_t = 369,612$

The minimizer $\mathbf{u}^*$ of the underlying optimization problem, $\min_{\mathbf{u}\in\mathbb{R}^N} F(\mathbf{u})$, is a fixed point of this map, i.e., $\mathcal{F}(\mathbf{u}^*) = \mathbf{u}^*$.

The Anderson acceleration iteration for (70) given a window size $w \in\mathbb{N}\cup\{+\infty\}$, denoted by $\mathrm{AA}(w)$, is given by

where for each $k$ the coefficients $\beta_i^{(k)}$ are determined by solving

where $r(\tilde{u}_0^{(k)}) := \tilde{u}_0^{(k)} - \mathcal{F}(\tilde{u}_0^{(k)})$ for all $k$. Put simply, AA($w$) attempts to accelerate (70) by using the $w$ previous initial guesses to generate an improved initial guess for the next all-at-once system. In the tests below we use AA(1), in which case only the previous initial guess is used, and AA($\infty$), corresponding to using the full available history of initial guesses.

In Figure 8, we present numerical experiments similar to those in § 5 in which we apply Anderson acceleration to R-MGRIT with a full time domain and truncated R-MGRIT with $\gamma= 0.2$. We consider the setting of the numerical experiments illustrated in Figure 2 and Figure 6 for MP1-1D, (7), and MP2-2D, (10), respectively; this time, we compare MGRIT with Anderson-accelerated R-MGRIT and Anderson-accelerated truncated R-MGRIT. We observe in Figure 8 that AA($\infty$) reduces the number of R-MGRIT iterations required for convergence in experiments for both model problems. AA(1) achieves a similar improvement in gradient reduction in all instances, except for the $\gamma$-truncated MP1 experiment. Since Anderson acceleration adds trivial overhead cost to the underlying fixed-point iteration, the acceleration provided would improve the parallel speedup estimates presented in § 6 by further reducing $\bar{N}_{it}$, the number of R-MGRIT iterations required for convergence, relative to $N_{it}$, the number of MGRIT iterations required for convergence.

Conclusions

We have extended the parallel-in-iteration MGRIT framework for optimization from [2] by introducing restarted MGRIT (R-MGRIT). Rather than repeatedly improving the approximation of a single optimization trajectory, R-MGRIT performs one MGRIT iteration per trajectory, using its approximate final point to initialize the next. For linear problems under suitable assumptions, we proved that, for any fixed $k$, the final-point approximations obtained from $k$ MGRIT iterations and $k$ R-MGRIT trajectories converge to the same limit as the trajectory length tends to infinity.

We further showed that the approximate endpoint of a single R-MGRIT trajectory approaches its infinite-trajectory limit geometrically, which together with our numerical observations motivates a truncated R-MGRIT strategy in which each trajectory is shortened by a factor $\gamma\in(0,1)$. Our numerical experiments on linear and nonlinear optimization problems showed that truncation generally does not reduce gradient reduction per trajectory, while leading to substantial performance gains: our parallel cost analysis predicts that shortening the trajectories can substantially increase parallel speedup while reducing the number of processors required.

Finally, we investigated the use of Anderson acceleration with R-MGRIT and observed further reductions in the number of trajectories required for convergence in our numerical experiments. Future work will consider extensions of these parallel-in-iteration methods to broader classes of optimization problems and optimization algorithms, along the lines discussed in [2].

Acknowledgements. The authors used Google’s Gemini and OpenAI’s ChatGPT during the final stages of the preparation of this manuscript. All AI suggestions mathematical claims, proofs, computations, citations, and final wording were reviewed, verified, and edited by the authors, who take full responsibility for the manuscript. Los Alamos National Laboratory report number LA-UR-26-28803.