Year: 2022

Abstract:
Anderson mixing (AM) is a powerful acceleration method for fixed-point iterations, but its computation requires storing many historical iterations. The extra memory footprint can be prohibitive when solving high-dimensional problems in a resource-limited machine. To reduce the memory overhead, we propose a novel class of short-term recurrence AM methods (ST-AM). The ST-AM methods only store two previous iterations with cheap corrections. We prove that the basic version of ST-AM is equivalent to the full-memory AM in strongly convex quadratic optimization, and with minor changes it has local linear convergence for solving general nonlinear fixed-point problems. We further analyze the convergence properties of the regularized ST-AM for nonconvex (stochastic) optimization. Finally, we apply ST-AM to several applications including solving root-finding problems and training neural networks. Experimental results show that ST-AM is competitive with the long-memory AM and outperforms many existing optimizers.

heading: 1 INTRODUCTION
text: Anderson mixing (AM) (Anderson, 1965) is a sequence acceleration method (Brezinski et al., 2018) in scientific computing. It is widely used to accelerate the slow convergence of fixed-point iterations (Lin et al., 2019; Fu et al., 2020; An et al., 2017), e.g., the self-consistent field iterations in electronic structure computations (Garza & Scuseria, 2012; Arora et al., 2017). Specifically, we consider a fixed-point iteration xk+1 = g(xk), k = 0, 1, . . . , where g : Rd 7→ Rd is the fixed-point map. By usingm historical iterations, AM(m) aims to extrapolate a new iterate that satisfies certain optimality property. When the function evaluation is costly, the reduction of the number of iterations brought by AM can save a large amount of computation (Fang & Saad, 2009).
AM can be used as a method for solving nonlinear equations (Kelley, 2018) as the fixed-point problem x = g(x) is equivalent to h(x) := x − g(x) = 0. In practice, since computing the Jacobian of h(x) is commonly difficult or even unavailable (Nocedal & Wright, 2006), AM can be seen as a practical alternate for Newton’s method (An et al., 2017). Also, compared with the classical iterative methods such as the nonlinear conjugate gradient (CG) method (Hager & Zhang, 2006), no line-search or trust-region technique is used in AM, which is preferable for large-scale unconstrained optimization. Empirically, it is observable that AM can largely accelerate convergence, though its theoretical analysis is still under-explored. It turns out that in the linear case (Walker & Ni, 2011; Potra & Engler, 2013), the full-memory AM (m = k) is essentially equivalent to GMRES (Saad & Schultz, 1986), a powerful Krylov subspace method that can exhibit superlinear convergence behaviour in solving linear systems (Van der Vorst & Vuik, 1993). For general nonlinear problems, AM is recognized as a multisecant quasi-Newton method (Fang & Saad, 2009; Brezinski et al., 2018). As far as we know, only local linear convergence has been obtained for the limited-memory AM (m < k) in general (Toth & Kelley, 2015; Evans et al., 2020; De Sterck & He, 2021).
For the application of AM, one of the major concerns is the historical length m, a critical factor related to the efficiency of AM (Walker & Ni, 2011). A larger m can incorporate more historical
∗Corresponding author.
information into one extrapolation, but it incurs heavier memory overhead since 2m vectors of dimension d need to be stored in AM(m). The additional memory footprint can be prohibitive for solving high-dimensional problems in a resource-limited machine (Deng, 2019). Using small m can alleviate the memory overhead but may deteriorate the efficacy of AM since much historical information is omitted in the extrapolation (Walker & Ni, 2011; Evans et al., 2020).
To address the memory issue of AM, we deeply investigate the properties of the historical iterations produced by AM and leverage them to develop the short-term recurrence variant, namely ST-AM. The basic version of ST-AM imposes some orthogonalization property on the historical sequence, which is inspired by the CG method (Hestenes & Stiefel, 1952) that enjoys a three-term recurrence. Furthermore, to better suit the more difficult nonconvex optimization, a regularized short-term form is introduced. We highlight the main contributions of our work as follows.
1. We develop a novel class of short-term recurrence AM methods (ST-AM), including the basic ST-AM, the modified ST-AM (MST-AM), and the regularized ST-AM (RST-AM). The basic ST-AM is applicable for linear systems, MST-AM can solve general fixed-point problems, and RST-AM aims for solving stochastic optimization. An important feature of ST-AM is that all methods only need to store two previous iterations with cheap corrections, which significantly reduces the memory requirement compared with the classical AM.
2. A complete theoretical analysis of the ST-AM methods is given. When solving strongly convex quadratic optimization, we prove that the basic ST-AM is equivalent to full-memory AM and the convergence rate is similar to CG method. We also prove that MST-AM has improved local linear convergence for solving fixed-point problems. Besides, we establish the global convergence property and complexity analysis for RST-AM when solving stochastic optimization problems.
3. The numerical results on solving (non)linear equations and cubic-regularized quadratic optimization are consistent with the theoretical results for the basic ST-AM and MST-AM. Furthermore, extensive experiments on training deep neural networks for image classification and language models show that RST-AM is competitive with the long-memory AM and outperforms many existing optimizers such as SGD and ADAM.

heading: 2 RELATED WORK
text: AM is also known as an extrapolation algorithm in scientific computing (Anderson, 2019). A parallel method is Shanks transformation (Shanks, 1955) which transforms an existing sequence to a new sequence for faster convergence. Related classical algorithms include Minimal Polynomial Extrapolation (Cabay & Jackson, 1976) and Reduced Rank Extrapolation (Eddy, 1979), and a framework of these extrapolation algorithms including AM is given in (Brezinski et al., 2018). Note that an elegant recursive algorithm named -algorithm has been discovered for Shanks transformation for scalar sequence (Wynn, 1956), and was later generalized as the vector -algorithm (Wynn, 1962) for handling vector sequences, but this short-term recurrence form is not equivalent to the original Shanks transformation in general (Brezinski & Redivo-Zaglia, 2017). Since AM is closely related to quasi-Newton methods (Fang & Saad, 2009), there are also some works trying to derive equivalent forms of the full-memory quasi-Newton methods using limited memory (Kolda et al., 1998; Berahas et al., 2021), while no short-term recurrence is available. To the best of our knowledge, ST-AM is the first attempt to short-term recurrence quasi-Newton methods.
Recently, there have been growing demands for solving large-scale and high-dimensional fixedpoint problems in scientific computing (Lin et al., 2019) and machine learning (Bottou et al., 2018). For these applications, Newton-like methods (Byrd et al., 2016; Wang et al., 2017; Mokhtari et al., 2018) are less appealing due to the heavy memory and computational cost, especially in nonconvex stochastic optimization, where only sublinear convergence can be expected if only stochastic gradients can be accessed (Nemirovskij & Yudin, 1983). On the other side, first-order methods (Necoara et al., 2019) stand out for their low per-iteration cost, though the convergence can be slow in practice. When training neural networks, SGD with momentum (SGDM) (Qian, 1999), and adaptive learning rate methods, e.g. AdaGrad (Duchi et al., 2011), RMSprop (Tieleman & Hinton, 2012), Adam (Kingma & Ba, 2014), are very popular optimizers. Our methods have the nature of quasiNewton methods while the memory footprint is largely reduced to be close to first-order methods. Thus, ST-AM can be a competitive optimizer from both theoretical and practical perspectives.

heading: 3 METHODOLOGY
text: In this section, we give the details of the proposed ST-AM. We always assume the objective function as f : Rd → R, the fixed-point map g : Rd 7→ Rd. Moreover, we do not distinguish rk = −∇f(xk) and rk = g(xk)− xk in our discussion as∇f(x) = 0 is equivalent to g(x) = −∇f(x) + x = x.

heading: 3.1 ANDERSON MIXING
text: The AM finds the fixed point of g via maintaining two sequences of length m(m ≤ k):
Xk = [∆xk−m,∆xk−m+1, · · · ,∆xk−1], Rk = [∆rk−m,∆rk−m+1, · · · ,∆rk−1] ∈ Rd×m, (1)
where the operator ∆ denotes the forward difference, e.g. ∆xk = xk+1 − xk. Each update of AM can be decoupled into two steps, namely the projection step and the mixing step:
x̄k = xk −XkΓk, (Projection step), xk+1 = x̄k + βkr̄k, (Mixing step), (2)
where r̄k := rk −RkΓk and βk > 0 is the mixing parameter. The Γk is determined by
Γk = arg min Γ∈Rm
‖rk −RkΓ‖2. (3)
Thus, the full form of AM (Fang & Saad, 2009; Walker & Ni, 2011) is
xk+1 = xk + βkrk − (Xk + βkRk) Γk. (4)
Remark 1. To see the rationality of AM, assume g is continuously differentiable, then we have h(xj)−h(xj−1) ≈ h′(xk)(xj−xj−1) around xk, where h′(xk) is the Jacobian of h(x) := x−g(x). So, it is reasonable to assume Rk ≈ −h′(xk)Xk, and we see ‖rk −RkΓ‖2 ≈ ‖rk + h′(xk)XkΓ‖2. Thus, we can recognize (3) as solving h′(xk)dk = h(xk) in a least-squares sense, where dk = XkΓk. The mixing step incorporates rk into the new update xk+1 if βk > 0. Otherwise, xk+1 = x̄k is an interpolation of the previous iterates, leading to a stagnating sequence.

heading: 3.2 THE BASIC SHORT-TERM RECURRENCE ANDERSON MIXING
text: The basic ST-AM is to solve the strongly convex quadratic optimization:
min x∈Rd
f(x) := 1
2 xTAx− bTx, (5)
where A 0. Initialized with p−1 = q−1 = p0 = q0 = 0 ∈ Rd, at the k-th iteration, given the two matrices Pk−1 = (pk−2, pk−1) ∈ Rd×2, Qk−1 = (qk−2, qk−1) ∈ Rd×2 and defining p = xk−xk−1 and q = rk − rk−1, the basic ST-AM constructs
p̃ = p− Pk−1(QTk−1q), q̃ = q −Qk−1(QTk−1q), (6a) pk = p̃/‖q̃‖2, qk = q̃/‖q̃‖2. (6b)
Then, we update Pk = (pk−1, pk), Qk = (qk−1, qk) ∈ Rd×2. Such construction ensures QTkQk = I2 for k ≥ 2 and the storage of Pk and Qk is equal to AM(2). With the corrected Pk and Qk, the ST-AM method modifies the projection and mixing step accordingly, that is,
x̄k = xk − PkΓk, (Projection step), xk+1 = x̄k + βkr̄k, (Mixing step), (7)
where Γk = arg min ‖rk−QkΓ‖2 = QTk rk and r̄k = rk−QkΓk. Thus, the ST-AM replacesXk and Rk in (1) by Pk and Qk respectively and imposes the orthogonalization property on Qk. The details of basic ST-AM are given in Algorithm 2 in Appendix C.1. Define P̄k = (p1, p2, . . . , pk), Q̄k = (q1, q2, . . . , qk), the Krylov subspace Km(A, v) ≡ span{v,Av,A2v, . . . , Am−1v}, the range of X as range(X). We give the properties of the basic ST-AM in Theorem 1. Theorem 1. Let {xk} be the sequence generated by the basic ST-AM. The following relations hold: (i) ‖q̃‖2 > 0, range(P̄k) = range(Xk) = Kk(A, r0), range(Q̄k) = range(Rk) = AKk(A, r0); (ii) Q̄k = −AP̄k, Q̄Tk Q̄k = Ik; (iii) r̄k ⊥ range(Q̄k) and x̄k = x0 + zk, where zk = arg minz∈Kk(A,r0) ‖r0 −Az‖2. If ‖r̄k‖2 = 0, then xk+1 is the exact solution.
The proof is in Appendix C.1. Note that the property (iii) in Theorem 1 exactly describes the relation x̄k = x G k , where x G k is the output of the k-th iteration of GMRES (Saad & Schultz, 1986). Moreover, let x̄AMk be the k-th intermediate iterate in full-memory AM (m = k). We prove that x̄ AM k = x G k (See Proposition 1 in Appendix C.1.), which induces that x̄k = x̄AMk = x G k . This equivalence indicates that ST-AM is more efficient than AM and GMRES since only two historical iterations need to be stored. Moreover, by directly applying the convergence analysis of GMRES (Corollary 6.33 in (Saad, 2003)), we obtain the convergence rate of the basic ST-AM for solving (5): Corollary 1. Suppose the eigenvalues of A lie in [µ,L] with µ > 0, and let {xk} be the sequence generated by the basic ST-AM, then the k-th intermediate residual r̄k satisfies ‖r̄k‖2 ≤
2 (√ L/µ−1√ L/µ+1 )k ‖r0‖2. Moreover, the algorithm finds the exact solution in at most (d+ 1) iterations.
Remark 2. The GMRES can be simplified to an elegant three-term recurrence algorithm called the conjugate residual (CR) method (Algorithm 6.20 in (Saad, 2003)) when solving (5). Thus, a similar simplification for AM is expected to exist. Like CG and Chebyshev acceleration (Algorithm 12.1 in (Saad, 2003)), the convergence rate of ST-AM has the optimal dependence on the condition number, while ST-AM does not form the Hessian-vector products explicitly.

heading: 3.3 THE MODIFIED SHORT-TERM RECURRENCE ANDERSON MIXING
text: For general nonlinear fixed-point problems, global convergence may be unavailable for the basic ST-AM, as a counter-example exists for full-memory AM (Mai & Johansson, 2020). Thus, we propose a modified version of the basic ST-AM (MST-AM) and prove the local linear convergence rate under similar conditions used in (Toth & Kelley, 2015; Evans et al., 2020). Concretely, the MST-AM makes three main changes to the basic ST-AM.
Change 1: Instead of applying the normalization (6b), the MST-AM constructs pk and qk via
ζk = (Q T k−1Qk−1) †QTk−1q, pk = p− Pk−1ζk, qk = q −Qk−1ζk, (8)
where “ † ” is the Moore-Penrose inverse. Accordingly, we choose Γk = arg min ‖rk − QkΓ‖2 = (QTkQk)
†QTk rk. This change relaxes the orthonormality forQk, but keeps the orthogonal condition: QTk−1qk = 0. In fact, Q̄ T k−1qk = 0 in the case of solving (5).
Change 2: MST-AM imposes the boundedness constraints onPk−1ζk andQk−1ζk: If ‖Pk−1ζk‖2 > cp‖p‖2 or ‖Qk−1ζk‖2 > cq‖q‖2, then Pk = Pk−1, Qk = Qk−1, where cp > 0, cq ∈ (0, 1) are predefined constants. It is worth mentioning that the boundedness condition is common in the analysis of AM (Toth & Kelley, 2015; Evans et al., 2020).
Change 3: MST-AM restarts, i.e. setting Pk = Qk = 0 ∈ Rd×2 every m iterations. This restart operation is to limit the number of higher-order terms appeared in the residual expansion in our analysis and we can set m to be a large number in practice.
The detailed description of MST-AM is given in Appendix C.2. In the next theorem, we establish the convergence rate analysis for the MST-AM. Theorem 2. Let {xk} be the sequence generated by MST-AM, x∗ ∈ Rd be a fixed point of g and m be the restarting period for MST-AM. Suppose that in the ball B(ρ) := {x ∈ Rd|‖x−x∗‖2 < ρ} for some ρ > 0, g is Lipschitz continuously differentiable and there are constants κ ∈ (0, 1) and κ̂ > 0 with (i) ‖g(y)− g(x)‖2 ≤ κ‖y−x‖2 for every x, y ∈ B(ρ), and (ii) ‖g′(y)− g′(x)‖2 ≤ κ̂‖y−x‖2 for every x, y ∈ B(ρ), where g′(x) is the Jacobian of g. If ‖x0 − x∗‖2 is sufficiently small, then for the k-th residual rk = g(xk)− xk, the following bound holds:
‖rk+1‖2 ≤ θk(|1− βk|+ κβk)‖rk‖2 + κ̂ mk∑ j=0 O ( ‖rk−j‖22 ) , (9)
where θk = ‖r̄k‖2/‖rk‖2 ≤ 1 and mk = k mod m. Thus, the residuals {rk} converge Q-linearly, and the errors {‖xk − x∗‖2} converge R-linearly. Remark 3. In a local region around x∗, the convergence rate is determined by the first-order term θk(|1−βk|+κβk)‖rk‖2. Since r̄k is the orthogonal projection of rk onto the subspace range(Qk)⊥, θk has the interpretation of the direction-sine between rk and the subspace range(Qk). When θk is
Algorithm 1 RST-AM for stochastic programming
Input: x0 ∈ Rd, βk > 0, αk ∈ [0, 1], δ(1)k > 0, δ (2) k > 0. Output: x ∈ Rd 1: P0, Q0 = 0 ∈ Rd×2, p0, q0 = 0 ∈ Rd 2: for k = 0, 1, . . . , until convergence, do 3: rk = −∇fSk(xk) 4: if k > 0 then 5: p = xk − xk−1, q = rk − rk−1 6: ζk = (QTk−1Qk−1 + δ (1) k P T k−1Pk−1)
†QTk−1q 7: qk = q −Qk−1ζk, pk = p− Pk−1ζk 8: Pk = [pk−1, pk], Qk = [qk−1, qk] 9: end if
10: Check Condition (13) and use smaller αk if (13) is violated 11: Γk = (QTkQk + δ (2) k P T k Pk)
†QTk rk 12: x̄k = xk − αkPkΓk, r̄k = rk − αkQkΓk 13: xk+1 = x̄k + βkr̄k 14: Apply learning rate schedule of αk, βk 15: end for 16: return xk
small, e.g., rk nearly lies in range(Qk), the acceleration by MST-AM is significant. Unlike AM(m) that simply discards the oldest iteration to make space for ∆xk−1 and ∆rk−1, MST-AM incorporates historical information with orthogonalization. In the SPD linear case and without restart, the global orthogonalization property holds, i.e. r̄k ⊥ range(Q̄k), which means there is no loss of historical information, while AM(m) (Evans et al., 2020) does not have such property in this ideal case.

heading: 3.4 THE REGULARIZED SHORT-TERM RECURRENCE ANDERSON MIXING
text: Inspired by the recent work on stochastic Anderson mixing (SAM) method (Wei et al., 2021), we develop a regularized ST-AM (RST-AM) for solving nonconvex stochastic optimization problems.
Consider the nonconvex optimization problem minx∈Rd f(x) := 1T ∑T i=1 fξi(x), where fξi : Rd → R is the loss function corresponding to i-th data sample and T is the number of data samples. In mini-batch training, the gradient is evaluated for fSk(xk) := 1 nk ∑ i∈Sk fξi(xk), where Sk ⊆ [T ] := {1, 2, . . . , T} is the sampled mini-batch, and nk := |Sk| is the batch size. In this case, we set rk = −∇fSk(xk) (Line 3 in Algorithm 1), which is an unbiased estimate of the negative gradient.
Recalling from (8), ζk = (QTk−1Qk−1) †QTk−1∆rk−1 = arg min ‖∆rk−1−Qk−1ζ‖2 as q = ∆rk−1 by definition. Since ST-AM is based on a local quadratic approximation (5) in a small region around xk, a large magnitude of ‖Pk−1ζk‖2 tends to make pk = ∆xk−1 − Pk−1ζk overshoot the trust region around xk, thus leading to instability. Consequently, we add a penalty term in the above least squares problem, i.e.
ζk = arg min ‖∆rk−1 −Qk−1ζ‖22 + δ (1) k ‖Pk−1ζ‖ 2 2, (10)
where δ(1)k > 0. The same as SAM (Wei et al., 2021), we also add a regularization term for computing Γk via
Γk = arg min ‖rk −QkΓ‖22 + δ (2) k ‖PkΓ‖ 2 2, (11)
where δ(2)k > 0, and a damping term is used as shown in Line 12 in Algorithm 1. In practice, we choose the two parameters as
δ (1) k = c1‖rk‖22 ‖∆xk−1‖22 + 0 , δ (2) k = max { c2‖rk‖22 ‖pk‖22 + 0 , Cβ−2k } , (12)
where c1, c2, C > 0 are constants, and 0 > 0 is a small constant to bound the denominators away from zero. Assuming ‖pk−1‖2 ≈ ‖pk‖2 = O(‖∆xk−1‖2), the choices of (12) make ‖δ(1)k PTk−1Pk−1‖2 ≈ O(‖rk‖22) and ‖δ (2) k P T k Pk‖2 ≈ O(‖rk‖22) aware of the change of the local curvature: large (small) ‖rk‖2 tends to lead to a large (small) regularization.
Remark 4. One update of xk given by Line 11-13 in Algorithm 1 can be formulated as xk+1 = xk + Hkrk, where Hk = βkI − αkYkZ†kQTk , Yk = Pk + βkQk, Zk = QTkQk + δ (2) k P T k Pk. To guarantee the positive definiteness of Hk, we follow the same procedure in SAM. Let λk be the largest eigenvalue of YkZ † kQ T k +QkZ † kY T k . If αk satisfies αkλk ≤ 2βk(1− µ), (13) then sTkHksk ≥ βkµ‖sk‖22,∀sk ∈ Rd, where µ ∈ (0, 1) is a constant. Note that λk can be cheaply obtained by computing the largest eigenvalue of a matrix of R4×4 (see Appendix C.3.1).
We summarize the RST-AM in Algorithm 1 and establish its convergence properties in the following context. First, we impose the following assumptions on the objective function f . Assumption 1. f : Rd → R is continuously differentiable. f(x) ≥ f low > −∞ for any x ∈ Rd. ∇f is globallyL-Lipschitz continuous; namely ‖∇f(x)−∇f(y)‖2 ≤ L‖x−y‖2 for any x, y ∈ Rd. Assumption 2. For any iteration k, the stochastic gradient ∇fξk(xk) satisfies Eξk [∇fξk(xk)] = ∇f(xk), Eξk [‖∇fξk(xk) − ∇f(xk)‖22] ≤ σ2, where σ > 0, and ξk, k = 0, 1, . . . are independent samples that are independent of {xj}kj=0.
The diminishing condition about βk is +∞∑ k=0 βk = +∞, +∞∑ k=0 β2k < +∞. (14)
We give the convergence properties of RST-AM in nonconvex (stochastic) optimization and proofs are deferred to Appendix C.3.2. Theorem 3. Suppose Assumption 1 hold and {xk} is the sequence generated by full-batch RST-AM, i.e. nk = T . Let βk = β ∈ (0, µ2L(1+C−1) ] be a constant, αk ∈ [0, 1] and satisfies (13), then
1
N N−1∑ k=0 ‖∇f(xk)‖22 ≤ 2(f(x0)− f low) Nµβ , (15)
in the N iterations. To ensure 1N ∑N−1 k=0 ‖∇f(xk)‖22 < , the number of iterations is O(1/ ) . Theorem 4. Suppose Assumptions 1 and 2 hold and {xk} is the sequence generated by RST-AM with batch size nk = n ≤ T . If βk ∈ (0, µ4L(1+C−1) ] and satisfies (14), αk ∈ [0,min{1, β 1 2
k }] and satisfies (13), then
lim inf k→∞
‖∇f(xk)‖2 = 0 with probability 1 and ∃Mf > 0→ E[f(xk)] ≤Mf ,∀k. (16)
If Eξk [‖∇fξk(xk)‖22] ≤Mg,∀k, where Mg > 0 is a constant, we have lim k→∞ ‖∇f(xk)‖2 = 0 with probability 1. (17) Theorem 5. Suppose Assumptions 1 and 2 hold and {xk}N−1k=0 is the first N iterations generated by RST-AM with fixed batch size nk = n. Let βk = min{ µ4L(1+C−1) , D̃ σ √ N }, where D̃ is a problemindependent constant; αk ∈ [0,min{1, β 1 2
k }] and satisfies (13). Let R be a random variable following PR(k) := Prob{R = k} = 1/N , then
E[‖∇f(xR)‖22] ≤ 16DfL(1 + C
−1)
Nµ2 +
σ
µ √ N
( 4Df
D̃ +
4(L+ µ−1)(1 + C−1)D̃
n
) , (18)
where Df := f(x0)− f low and the expectation is taken with respect to R and {Sj}N−1j=0 . To ensure E[‖∇f(xR)‖22] ≤ , the number of iterations is O(1/ 2). Remark 5. The proofs of Theorem 4 and 5 are based on the analysis of SAM (Wei et al., 2021). The theorems show that the convergence of RST-AM is no worse than SGD (Robbins & Monro, 1951). There are two key differences between RST-AM and SAM: RST-AM is based on short-term recurrences while SAM usually maintains longer historical sequences to ensure effectiveness; RSTAM uses additional correction and regularization terms (Line 5-8 in Algorithm 1) to incorporate historical information while SAM simply discards the oldest iteration to make space for ∆xk−1 and ∆rk−1. The reduced memory requirement in RST-AM makes it applicable for solving more challenging problems in machine learning.

heading: 4 EXPERIMENTS
text: We validated the effectiveness of our proposed ST-AM methods in various applications in fixed-point iterations and nonconvex optimization, including linear and nonlinear problems, deterministic and stochastic optimization. Specifically, we first tested ST-AM in linear problems, cubic-regularized quadratic minimization (Carmon & Duchi, 2020) and a multiscale deep equilibrium (MDEQ) model (Bai et al., 2020). Then we applied RST-AM to train neural networks and compared them with several first-order and second-order optimizers. Experimental details are referred to Appendix D.

heading: 4.1 EXPERIMENTS ABOUT THE BASIC ST-AM AND MST-AM
text: We verified the properties of ST-AM declared in Theorem 1 and 2 by solving four problems (details are in Appendix D.1): (I) strongly convex quadratic optimization (corresponding to Theorem 1); (II) solving a nonsymmetric linear system Ax = b (corresponding to κ̂ = 0 in Theorem 2); (III) cubic-regularized quadratic minimization minx∈Rd f(x) := ‖Ax− b‖22 + M3 ‖x‖ 3 2 (corresponding to κ̂ > 0 in Theorem 2); (IV) root-finding problems in MDEQ on CIFAR-10 (Krizhevsky et al., 2009).
The compared methods were gradient descent (GD), fixed-point iteration (FP), conjugate residual method (CR) (Saad, 2003), BFGS (Nocedal & Wright, 2006), Broyden’s method (Broyden, 1965), and the full-memory AM (AM). We used the basic ST-AM to solve Problem I and II, and MST-AM (cp = cq = 1) to solve Problem III and IV.
The numerical results shown in Figure 1 demonstrate the power of ST-AM as a variant of Krylov subspace methods. It significantly accelerates the slow convergence of the GD or FP method, and can outperform AM. Figure 1(a) clearly verifies the correctness of Theorem 1: within the machine precision, the intermediate residual r̄k of ST-AM coincides with the residual rk of CR. Note that AM fails to coincide with CR and ST-AM due to the intrinsic numerical weakness to solve (3), as also pointed out in (Walker & Ni, 2011). Figure 1(b) shows that ST-AM can outperform CR, though both methods enjoy short-term recurrences and are equivalent for solving SPD linear systems. Figure 1(c) also shows MST-AM surpasses BFGS in solving cubic-regularized problems. The tests in MDEQ/CIFAR-10 indicate that MST-AM is comparable to the full-memory methods in the forward root-finding process and converges faster in the backward process. The accuracy is also comparable.

heading: 4.2 EXPERIMENTS ABOUT RST-AM
text: We applied RST-AM to train neural networks, with full-batch training on MNIST (LeCun et al., 1998), and mini-batch training on CIFAR-10/CIFAR-100 and Penn Treebank (Marcus et al., 1993).
Experiments on MNIST. We trained a convolutional neural network (CNN) on MNIST to see the convergence behaviour of RST-AM in nonconvex optimization (cf. Theorem 3), for which we were only concerned about the training loss. Figure 2(a)(b) show that the short-memory SAMs (m = 2, 5) hardly show any improvement over the first-order optimizers SGD and Adam, while RST-AM can close the gap of the long-memory (m = 10) and the short-memory methods. We also considered the effect of preconditioning on RST-AM (see Appendix A.3). The notation “A B” means B method preconditioned by A method. Figure 2(c) indicates that preconditioning also works much better for RST-AM than SAM(2), and RMSprop RST-AM can outperform the non-preconditioned RST-AM.
Experiments on CIFAR. We trained ResNet18/20/32/44/56 (He et al., 2016), WideResNet16-4 (Zagoruyko & Komodakis, 2016) (abbr. WRN16-4) on CIFAR-10, and ResNet18, ResNeXt50 (Xie et al., 2017), DenseNet121 (Huang et al., 2017) on CIFAR-100. The baseline optimizers were SGDM, Adam, AdaBound (Luo et al., 2018), AdaBelief (Zhuang et al., 2020), Lookahead (Zhang et al., 2019), AdaHessian (Yao et al., 2021) and SAM. Table 1(a) reports the final TOP1 test accuracy for training 160 epochs (AdaHessian ran out of memory for training ResNeXt50, DenseNet121). It shows RST-AM improves SAM(2) and has comparable test accuracies to SAM(10). RST-AM also outperforms other baseline optimizers. Table 1(b) reports the memory and computation cost, where
we used SGDM as the baseline and other optimizers were terminated when achieving a comparable or better test accuracy than SGDM. It indicates that RST-AM introduces ≤ 5% extra memory overhead compared with SGDM, and significantly reduces the memory footprint of AM. Since RST-AM needs fewer training epochs, the total running time is less than SGDM.
Table 2: Test perplexity on Penn Treebank for 1,2,3-layer LSTM. Lower is better.
Method 1-Layer 2-Layer 3-Layer
SGDM 83.48±.03 65.89±.18 61.88±.23 Adam 80.33±.15 64.32±.06 59.72±.13 AdaBelief 81.29±.35 64.68±.10 60.46±.07 SAM(2) 80.79±.19 65.52±.29 61.13±.12 SAM(10) 78.78±.14 62.46±.11 58.93±.09 RST-AM 78.41±.18 62.46±.08 58.31±.23
Experiments on Penn Treebank. We trained LSTMs with 1-3 layer(s) on Penn Treebank and report the validation perplexity in Figure 3 and test perplexity in Table 2 (lower is better). The results suggest that RST-AM is comparable to or even better than SAM(10). The improvement of RST-AM over other optimizers is also significant. We report the computation and memory cost in Appendix D.2.4. RST-AM can still surpass Adam while using much fewer epochs, thus reducing the total running time.
Table 4: FID score for SN-GAN.
Method Adam AdaBelief RST-AM
Best FID 13.07±.18 12.80±.09 12.05±.15 Final FID 13.34±.14 13.59±.21 12.50±.29
Adversarial training. We applied RST-AM to adversarial training (Madry et al., 2018) as the outeroptimizer and compared it with SGD by the clean test accuracy and robust test accuracy. The results on CIFAR10/ResNet18 and CIFAR100/DenseNet121 are reported in Table 3. It can be seen that RST-AM can achieve both higher clean test accuracy and robust test accuracy. More results can be found in Appendix D.2.5.
Generative adversarial network (GAN). We tested RST-AM by training a GAN equipped with spectral normalization (SN-GAN) (Miyato et al., 2018), where the generator and discriminator networks were ResNets and the dataset was CIFAR-10. Table 4 shows that RST-AM can achieve lower FID score (better accuracy) than Adam and AdaBelief.

heading: 5 CONCLUSION
text: In this paper, to address the memory issue of Anderson mixing (AM), we develop a novel class of short-term recurrence AM methods (ST-AM) and test it on various applications, including solving linear and nonlinear problems, deterministic and stochastic optimization. We give a complete theoretical analysis of the proposed method. We prove that the basic ST-AM is equivalent to the full-memory AM in strongly convex quadratic optimization. With some minor changes, it has local convergence for nonlinear fixed-point iterations under some common assumptions. We also introduce the regularized form of ST-AM and analyze its asymptotic and non-asymptotic convergence properties. The numerical results show that ST-AM methods are comparable to or even better than the long-memory AM while consuming less memory. The regularized ST-AM also outperforms many existing optimizers in training neural networks in various tasks.

heading: ACKNOWLEDGMENTS
text: This work was supported by the National Key R&D Program of China (No. 2021YFA1001300), National Natural Science Foundation of China (No.61925601), Tsinghua University Initiative Scientific Research Program, National Natural Science Foundation of China (No.11901338), and Huawei Noah’s Ark Lab. We thank all anonymous reviewers for their valuable comments and suggestions on this work.

heading: A ADDITIONAL PRELIMINARIES
text: We provide some additional preliminaries in this section for readers that are not familiar with Anderson mixing, fixed-point iterations and some techniques mentioned in the main paper.

heading: A.1 FIXED-POINT ITERATION
text: The fixed-point problem and the optimization problem are the main application scenarios of our methods. It is worth pointing out that there are some minor differences between these two problems that make algorithm designs different. The key difference is that for a fixed-point problem, the Jacobian (if exists) is generally not symmetric while for optimization, the Hessian is naturally symmetric. In principle, a fixed-point solver can also be applicable for an optimization problem since the first-order necessary condition of minx∈Rd f(x), where f : Rd → R, is ∇f(x) = 0.
Consider a contraction mapping g : Rd 7→ Rd, i.e. for some κ < 1, ‖g(x) − g(y)‖2 ≤ κ‖x − y‖2,∀x, y ∈ Rd. According to the contraction mapping theorem, a unique fixed point x∗ exists for g and the iterates generated by the iteration xk+1 = g(xk) converge to x∗, starting from any x0 ∈ Rd. In practice, the damped fixed-point iteration is also commonly used: given the βk ∈ (0, 2), the update is xk+1 = (1 − βk)xk + βkg(xk) = xk + βkrk, where rk := g(xk) − xk is called the residual.
The fixed-point iteration, also known as Picard iteration in some areas, can converge very slowly in practice. Anderson mixing is a method to improve the convergence.

heading: A.2 ANOTHER FORM OF ANDERSON MIXING
text: The derivation of Anderson mixing (AM) in Section 3.1 explicitly interprets AM as a two-step procedure. In the literature, there is another equivalent form of AM.
Let the projection coefficients Γk = (γ (1) k , . . . , γ (m) k ) T ∈ Rm. Define the auxiliary coefficients {θ(j)k }mj=0 as θ (0) k = γ (1) k , θ (j) k = ∆γ (j) k (j = 1, . . . ,m − 1), θ (m) k = 1 − γ (m) k , then ∑m j=0 θ (j) k =
1 and r̄k = ∑m j=0 θ (j) k rk−m+j . Hence the least squares problem (3) can be reformulated as a
constrained problem min{θ(j)k }mj=0 ‖ ∑m j=0 θ (j) k rk−m+j‖2 s.t. ∑m j=0 θ (j) k = 1, which indicates that as a linear combination of the historical residuals {rj}mj=k−m, r̄k is minimal in terms of theL2-norm. Also, the projection step and the mixing step (2) can be reformulated as x̄k = ∑m j=0 θ (j) k xk−m+j
and xk+1 = (1−βk) ∑m j=0 θ (j) k xk−m+j+βk ∑m j=0 θ (j) k g(xk−m+j), respectively. Such formulation is also adopted in the literature (Toth & Kelley, 2015; Mai & Johansson, 2020; Scieur et al., 2020).
Let Hk be the solution to the constrained optimization problem (Fang & Saad, 2009)
min Hk ‖Hk − βkI‖F subject to HkRk = −Xk, (19)
then iterate (4) is xk+1 = xk + Hkrk. Hence AM is a multisecant quasi-Newton method (Fang & Saad, 2009).

heading: A.3 PRECONDITIONED ANDERSON MIXING
text: Like preconditioning for Krylov subspace methods (Saad, 2003), preconditioning can also be incorporated into Anderson mixing to mitigate the ill-conditioning of the original problem (Wei et al., 2021). The idea is to replace the mixing step in (2) via a preconditioned mixing.
Suppose that there is a basic solver preconditioner(xk, sk), which works as a black-box procedure that updates xk given the residual sk, i.e. xk+1 = preconditioner(xk, sk), then the preconditioned mixing of RST-AM is
xk+1 = preconditioner(x̄k, r̄k).
which substitutes for the Line 13 in Algorithm 1. The simple mixing (Line 13 in Algorithm 1 can be seen as a special case by defining preconditioner(xk, sk) = xk + βksk, i.e. preconditioned by a damped fixed-point iteration. Moreover, if we write the preconditioning operation as
preconditioner(xk, sk) := xk + Gksk, where Gk is the matrix to approximate the inverse Jacobian, then the preconditioned AM is
xk+1 = xk +Gkrk − (Xk +GkRk)(RTkRk)†RTk rk (20)
(cf. the definitions of Xk, Rk in Section 3.1). The matrix Hk = Gk − (Xk + GkRk)(RTkRk)†RTk forms a low-rank updated approximation to the inverse Jacobian: it solves
min Hk ‖Hk −Gk‖F subject to HkRk = −Xk ,
which is a direct extension of (19).
For the stochastic Anderson mixing, using damped projection can be helpful, i.e. x̄k = (1−αk)xk+ αk(xk−XkΓk) = xk−αkXkΓk and r̄k = rk−αkRkΓk correspondingly. Then the preconditioned AM with damped projection is
xk+1 = xk +Gkrk − αk(Xk +GkRk)Γk. (21)

heading: B ADDITIONAL DISCUSSION

heading: B.1 THE MEMORY AND COMPUTATIONAL EFFICIENCY
text: Since the ST-AM methods only need to store two previous iterations, the memory and computational cost can be reduced to be close to first-order methods.
Memory cost. Assume that the iteration number is k and the model parameter size is d. The fullmemory AM stores all previous iterations, thus the additional memory is 2kd. To reduce the memory overhead, the limited-memory AM(m) only maintains the most recent m iterations while discarding the older historical information (cf. (1)). Hence the additional memory of AM(m) is 2md. Choosing a suitable m can be problem-dependent (Walker & Ni, 2011), and it is often suggested that m ≥ 5 (Mai & Johansson, 2020; Fu et al., 2020) to avoid too much historical information being discarded. There is no equivalence between AM(m) and the full-memory AM. On the other side, our ST-AM methods have the same memory footprint as AM(2), i.e. only introducing 4d additional memory. For a large model that contains millions or even billions of parameters, such reduction in memory requirement can be significant.
Computational cost. For ST-AM, besides the gradient evaluations, the main additional computational cost comes from vector corrections, the projection and mixing step. These operations are very cheap: matrix multiplications of R2×d ×Rd×2, R2×d ×Rd×1, and Rd×2 ×R2×1; the pseudoinverse can be exactly solved as the size of the matrix is R2×2. Also, the QTk−1Qk−1 and PTk−1Pk−1 in Line 6 in Algorithm 1 can reuse the results in the previous iteration (cf. Line 11). So the total additional computational cost of RST-AM is O(d).

heading: B.2 APPLICABILITY
text: In the main paper, we develop a class of short-term recurrence Anderson mixing. Among them, the basic ST-AM can serve as a new linear solver for linear systems; MST-AM can be used as a new fixed-point solver for nonlinear equations; RST-AM is a new method for optimization. These methods are based on Anderson mixing and have close relationship between each other. The theoretical analysis and numerical results show the great potential of ST-AM methods for various applications. Here, we give some comparisons between the ST-AM methods and some other classical methods.
ST-AM versus current solvers for linear systems. ST-AM is an iterative solver. It is equivalent to the full-memory AM and GMRES in strongly convex quadratic optimization (or SPD linear systems), and can also have linear convergence rate for general nonsymmetric linear systems. Since it is based on short-term recurrences, like the CG method, the memory and per-iteration cost of ST-AM is economical. On the other side, the LU factorization based methods (Golub & Van Loan, 2013) are direct solvers that can incur overwhelming memory and computation cost for large sparse linear systems. Although sparse direct solvers (Duff et al., 2017) can alleviate overhead, they are often more complicated to implement and difficult for parallel computing. Moreover, iterative solvers can benefit from the preconditioning technique that improves convergence (Saad, 2003). Also, ST-AM
has advantages over other iterative methods since it does not need to directly access the matrix and only the residual is required. For applications that explicitly accessing the matrix is difficult, this property is appealing. Hence, unlike nonlinear CG that relies on line search, it is direct to extend ST-AM to unconstrained optimization where the gradient is commonly available while the Hessian can be too costly to obtain. Moreover, ST-AM is very flexible and any iterative solver can be viewed as a black-box iterative process to be accelerated by ST-AM.
MST-AM versus other quasi-Newton methods for nonlinear equations. MST-AM has the nature of quasi-Newton methods since it is built upon AM that is recognized as a multisecant quasiNewton method (Fang & Saad, 2009). In scientific computing, AM has been successfully applied to solve many difficult nonlinear problems arising from many areas (An et al., 2017; Lin et al., 2019; Fu et al., 2020; Yang, 2021). Since AM only manipulates residuals and does not use line-search or trust-region technique, it is efficient to apply AM to accelerate a slowly convergent black-box iterative process. A comprehensive discussion about the applicability of AM for nonlinear problems and the relation between AM and Broyden’s methods can be found in (Fang & Saad, 2009).
One of the biggest issues of AM and other quasi-Newton methods is the additional memory overhead, because they need to store historical iterations to form the secant equations. To make a compromise, limited-memory quasi-Newton methods such as L-BFGS (Liu & Nocedal, 1989) are proposed in which only limited number of historical iterations are stored. In each update, the oldest iteration is discarded to make space for the up-to-date secant pair. As a result, the limited-memory quasi-Newton methods can lose the local superlinear convergence properties achieved by the fullmemory schemes (Berahas et al., 2021). Also, as stated by Berahas et al. (2021), one does not know a priori what number of historical iterations (i.e. historical length) should be maintained to attain a good performance when solving a particular problem.
Unlike the limited-memory quasi-Newton methods whose performance can be sensitive to the historical length, our MST-AM method only needs to store two corrected historical iterations. MSTAM carefully incorporates historical information through orthogonalization. In the ideal case, i.e. strongly convex quadratic optimization (or SPD linear systems), it is equivalent to the full-memory AM which means there is no loss of historical information. Our experiments verify the theoretical properties of MST-AM and indicate MST-AM can be a competitive method for nonlinear equations.
RST-AM versus first-order methods for stochastic optimization. For stochastic optimization such as training deep neural networks in machine learning, the first-order methods have come to dominate the field due to the low memory and per-iteration cost. The RST-AM is an extension of ST-AM and MST-AM to tackle this challenging problem. RST-AM has theoretical guarantees in deterministic/stochastic optimization and also inherits several advantages of AM and ST-AM:
• Fast convergence in quadratic optimization. It is important to possess such property for an optimizer since a smooth function can be approximated by a quadratic function in the local region around the optima and many techniques such as trust region (Nocedal & Wright, 2006) rely on this local approximation. Adaptive learning rate methods such as AdaGrad (Duchi et al., 2011), Adam (Kingma & Ba, 2014) use a diagonal approximation of the Fisher information matrix that is an approximation of the Hessian. However, in the simple quadratic case, these methods can only roughly match the performance of the Jacobi method for solving linear systems (Saad, 2003), which is better than gradient descent but far inferior to the powerful Krylov subspace methods. The momentum method mimics CG method by incorporating a historical iteration into the search direction. However, the choice of momentum and stepsize can be an art (Sutskever et al., 2013). For RST-AM, there is no need to determine the stepsize and the fast convergence rate of ST-AM can be recovered by simply setting αk = 1 and δ (1) k = δ (2) k = 0. For more difficult functions, the
damping and regularization in RST-AM can be enabled to improve stability.
• Theoretical guarantee in stochastic optimization. If only first-order information can be accessed, the SGD (Ghadimi & Lan, 2013) achieves an optimal convergence rate O(1/ 2) to obtain an -accurate solution (Nemirovskij & Yudin, 1983). In such case, it seems to be a big mismatch for current second-order methods, because there is no theoretical improvement albeit with more memory and computation resource. Such mismatch may also account for the popularity of first-order methods. Since RST-AM has very limited addi-
tional memory overhead, and also achieves the O(1/ 2) complexity, it can be applied to many applications that dominated by first-order methods.
• Flexibility in use. The application of RST-AM can be very flexible. In principle, RST-AM can be applied to improve any slowly convergent black-box iterative process by viewing the latter as a fixed-point iteration. For example, consider accelerating a solver of a commercial software, where we have no access to the underlying codes to provide our custom implementation, or some case rewriting the codes is too cumbersome. Moreover, RSTAM can efficiently incorporates the preconditioning technique. Hence any optimizer, even an optimizer built upon neural networks (Andrychowicz et al., 2016), can be used as a preconditioner for RST-AM. Preconditioning largely enhances the applicability of RSTAM for various applications. For example, RST-AM can be preconditioned by Adam for the language task and SGDM for image classification. So any fine-tuned first-order optimizer can be combined with RST-AM to achieve an overall improvement. As RST-AM is light-weight, this additional cost is marginal and can be largely counteracted by the actual improvement. So in some sense, the purpose of RST-AM is not to totally replace current off-the-shelf optimizers but to achieve collaborative effectiveness: RST-AM is aware of the second-order information while the first-order method can mitigate the ill-conditioning of the problem.
Overall, the ST-AM methods have wide applicability and can be competent methods from both theoretical and practical perspectives.

heading: C PROOFS
text: We give more details about ST-AM and the proofs of the theorems in the main paper.

heading: C.1 THE BASIC ST-AM FOR STRONGLY CONVEX QUADRATIC OPTIMIZATION
text: Recall that the strongly convex quadratic optimization is formulated as
min x∈Rd
f(x) := 1
2 xTAx− bTx, (22)
where A ∈ Rd×d is SPD, b ∈ Rd. Solving (22) is equivalent to solving the SPD linear system
Ax = b. (23)
The detail of the basic ST-AM is given in Algorithm 2.
We first state the relationship of AM with GMRES in the following proposition. Similar results can also be found in (Walker & Ni, 2011).
Let xGk , r G k := b−AxGk denote the k-th GMRES iterate and residual, respectively, and Kk(A, v) := span{v,Av, . . . , Ak−1v} denotes the k-th Krylov subspace generated by A and v. Define ej := (1, 1, . . . , 1)T ∈ Rj for j ≥ 1. Let range(X) denote the linear space spanned by the columns of X . The main results of the full-memory AM are stated in Proposition 1 and Proposition 2.
Proposition 1 (General linear system). For solving a general linear system Ax = b with the fullmemory AM (m = k), suppose that βk > 0 and the fixed-point map is g(x) = (I − A)x + b. If the initial point of AM is x0 = xG0 and rank(Rk) = m, then the intermediate iterate x̄k satisfies x̄k = x G k .
Proof. The definition of the fixed-point map suggests that the residual rk = g(xk)−xk = b−Axk. Since Rk = −AXk and A is nonsingular, we have rank(Xk) = m. We first show
range(Xk) = Kk(A, rG0 ) (24)
by induction. We abbreviate Kk(A, rG0 ) as Kk in this proof.
Algorithm 2 ST-AM for strongly convex quadratic optimization Input: x0 ∈ Rd, βk > 0, 0 < max iter ≤ d. Output: x ∈ Rd
1: P0, Q0 = 0 ∈ Rd×2, p0, q0 = 0 ∈ Rd 2: for k = 0, 1, . . . ,max iter do 3: rk = −∇f(xk) 4: if k > 0 then 5: p = xk − xk−1, q = rk − rk−1 (Compute ∆xk−1,∆rk−1) 6: q̃ = q −Qk−1(QTk−1q), p̃ = p− Pk−1(QTk−1q) (q̃ ⊥ Qk−1) 7: pk = p̃/‖q̃‖2, qk = q̃/‖q̃‖2 (‖qk‖2 = 1) 8: Pk = [pk−1, pk], Qk = [qk−1, qk] (QTkQk = I2) 9: end if
10: Γk = QTk rk 11: x̄k = xk − PkΓk, r̄k = rk −QkΓk (Projection step: r̄k ⊥ Qk) 12: xk+1 = x̄k + βkr̄k (Mixing step) 13: if ‖r̄k‖2 = 0 then 14: break 15: end if 16: end for 17: return xk
First, ∆x0 = β0r0 = β0rG0 since x1 = x0 + β0r0. If k = 1, then the proof is complete. Then, suppose that k > 1 and, as an inductive hypothesis, that range(Xk−1) = Kk−1. With (4) we have
∆xk−1 = xk − xk−1 = βk−1rk−1 − (Xk−1 + βk−1Rk−1)Γk−1 = βk−1(b−Axk−1)− (Xk−1 − βk−1AXk−1)Γk−1 = βk−1b− βk−1A(x0 + ∆x0 + · · ·+ ∆xk−2)− (Xk−1 − βk−1AXk−1)Γk−1 = βk−1r0 − βk−1AXk−1ek−1 − (Xk−1 − βk−1AXk−1)Γk−1. (25)
Since r0 ∈ Kk−1, and by the inductive hypothesis range(Xk−1) ⊆ Kk−1 which also implies range(AXk−1) ⊆ Kk, we know ∆xk−1 ∈ Kk, which implies range(Xk) ⊆ Kk. Since we assume rank(Xk) = m = k which implies dim(range(Xk)) = dim(Kk), we have range(Xk) = Kk, thus completing the induction. As a result, we also have
range(Rk) = range(AXk) = AKk(A, rG0 ). (26)
Recalling that to determine Γk, we solve the least squares problem (3) and Rk = −AXk. We have Γk = arg min
Γ∈Rm ‖rk +AXkΓ‖2. (27)
Since rank(AXk) = rank(Xk) = m, (27) has a unique solution. Also, since rk = b − Axk = b−A(x0 +Xkek) = r0 −AXkek, we have rk +AXkΓ = r0 −AXkek +AXkΓ = r0 −AXkΓ̃, where Γ̃ = ek − Γ. So Γk solves (27) if and only if Γ̃k = ek − Γk solves
min Γ̃∈Rm
‖r0 −AXkΓ̃‖2, (28)
According to (24), (28) is equal to minz∈Km(A,rG0) ‖r0 − Az‖2 which is the GMRES minimization problem. Since the solution of (28) is also unique, we have
x̄k = xk −XkΓk = xk −Xk(ek − Γ̃k) = x0 +XkΓ̃k = xGk .
In Proposition 1, the assumption that Rk has full column rank is critical to ensure no stagnation occurs in AM for solving a general linear system. In fact, for SPD linear systems (23) or strongly convex quadratic optimization (22), when AM breaks down, i.e. Rk is rank deficient, AM obtains the exact solution, as shown in the next proposition.
Proposition 2 (SPD). For applying the full-memory AM to minimize a strongly convex quadratic problem (22), or equivalently, solve a SPD linear system (23), suppose that βk > 0 and the fixedpoint map is g(x) = (I − A)x + b. If rank(Rk) = k holds for 1 ≤ k < s while failing to hold for k = s, where s ≥ 1, then the residual of AM satisfies rs = r̄s−1 = 0.
Proof. The definition of g suggests that the residual rk = g(xk)−xk = b−Axk. The relationRk = −AXk holds during the iterations and the nonsingularity of A implies rank(Xk) = rank(Rk). For s = 1, since the first step of AM is x1 = x0 + β0r0, the assumption rank(R1) = 0 implies that rank(r0) = rank(X1) = 0, i.e. r1 = r̄0 := 0.
For s > 1, because ∆xs−1 = xs − xs−1 = −Xs−1Γs−1 + βs−1r̄s−1, the rank deficiency of Xs implies ∆xs−1 ∈ range(Xs−1), which further implies r̄s−1 ∈ range(Xs−1). So there exists ζ ∈ Rs−1, such that r̄s−1 = Xs−1ζ. Note that according to (3), r̄s−1 ⊥ Rs−1 = −AXs−1, so we have
0 = r̄Ts−1AXs−1 = (Xs−1ζ) TAXs−1 = ζ TXTs−1AXs−1. (29)
Because rank(Xs−1) = s− 1 and A is SPD, we know XTs−1AXs−1 is also SPD. So ζ = 0, which implies r̄s−1 = 0. Hence xs = x̄s−1 and rs = r̄s−1 = 0.
Now we give the proof of Theorem 1.
Proof of Theorem 1. Besides relations (i)-(iii), we add an auxiliary relation here: (iv) rk = r0 + Q̄kΓ̄k ∈ Kk+1(A, r0), where Γ̄k ∈ Rk. We prove the relations (i)-(iv) by induction.
For k = 1, since r̄0 6= 0, according to Proposition 2, rank(∆x0) = rank(X1) = 1, rank(∆r0) = rank(R1) = 1, so q̃ 6= 0, which implies Line 7 in Algorithm 2 is well-defined. The relation (i) holds. Since q̃ = q = ∆x0, p̃ = p = ∆r0, and ∆r0 = −A∆x0, the equality Q̄1 = −AP̄1 also holds. Due to the normalization in Line 7, Q̄T1 Q̄1 = 1. Since r1 = r0 − β0Ar0 and range(Q̄1) = range(Ar0), it is clear that r1 = r0−Q̄1Γ̄1 ∈ K2(A, r0), namely relation (iv). Due to the projection step Line 11, r̄1 ⊥ range(Q1) = range(Q̄1). Also, r̄1 = r1−Q1Γ1 = r0−β0Ar0−Q1Γ1 = r0− Q̄1η1, where the last equality is due to span{Ar0} = range(Q1) = range(Q̄1). For rG1 = r0 − Az1, where z1 = arg minz∈K1(A,r0) ‖r0 − Az‖2, it holds r G 1 ⊥ AK1(A, r0) = range(Q̄1). As a result, both r̄1 and rG1 are the orthogonal projections of r0 onto the subspace range(Q̄1) ⊥, which implies r̄1 = rG1 . So x̄1 = xG1 = x0 + z1 because their residuals are equal and A is nonsingular. Hence relation (iii) holds.
Suppose that k > 1, and as an inductive hypothesis, the relations (i)-(iv) hold for j = 1, . . . , k − 1. Consider the k-th iteration. From Line 6 in Algorithm 2, q̃ ∈ range(∆rk−1, Qk−1), and p̃ ∈ range(∆xk−1, Pk−1). We first prove that q̃ 6= 0 by contradiction. If q̃ = 0, then from Line 6 in Algorithm 2, ∆rk−1 ∈ range(Qk−1) ⊆ range(Q̄k−1), which implies ∆xk−1 ∈ range(Pk−1) ⊆ range(P̄k−1) as ∆rk−1 = −A∆xk−1 and Q̄k−1 = −AP̄k−1 and A is nonsingular. From Line 11 and Line 12, we have
∆xk−1 = xk − xk−1 = −Pk−1Γk−1 + βk−1r̄k−1. (30)
So r̄k−1 ∈ range(Pk−1) ⊆ range(P̄k−1) since ∆xk−1 ∈ range(Pk−1). Hence there exists ζ ∈ Rk−1, such that r̄k−1 = P̄k−1ζ. From the inductive hypothesis, we know r̄k−1 ⊥ Q̄k−1 = −AP̄k−1, so we have
0 = r̄Tk−1AP̄k−1 = (P̄k−1ζ) TAP̄k−1 = ζ TP̄Tk−1AP̄k−1.
Since Q̄Tk−1Q̄k−1 = Ik−1, we know rank(Q̄k−1) = k− 1, which implies rank(P̄k−1) = k− 1 due to Q̄k−1 = −AP̄k−1. Hence P̄Tk−1AP̄k−1 is also SPD. Then ζ = 0 which implies r̄k−1 = 0. It is impossible otherwise Algorithm 2 has terminated in the (k− 1)-th iteration. So q̃ 6= 0 and Line 7 is well-defined.
Since r̄k−1 = rk−1 − Qk−1Γk−1, and rk−1 ∈ Kk(A, r0), range(Qk−1) ⊆ range(Q̄k−1) = AKk−1(A, r0) as the inductive hypothesis, we have r̄k−1 ∈ Kk(A, r0), which together with (30)
and range(Pk−1) ⊆ range(P̄k−1) = Kk−1(A, r0) infers ∆xk−1 ∈ Kk(A, r0). Hence ∆rk−1 = −A∆xk−1 ∈ AKk(A, r0). As a result, qk = q̃/‖q̃‖2 ∈ range(∆rk−1, Qk−1) ⊆ AKk(A, r0). So range(Q̄k) = range(Q̄k−1, qk) ⊆ AKk(A, r0). Moreover, qk /∈ range(Q̄k−1), otherwise q̃ ∈ range(Q̄k−1) that implies ∆rk−1 = q ∈ range(Q̄k−1), which is impossible following the former proof of q̃ 6= 0. So we have range(Q̄k) = AKk(A, r0). Because ∆rk−1 = −A∆xk−1 and Qk−1 = −APk−1 due to Q̄k−1 = −AP̄k−1, Line 6 in Algorithm 2 infers q̃ = −Ap̃, which implies qk = −Apk. So Q̄k = −AP̄k. Since A is nonsingular and range(Q̄k) = AKk(A, r0), we have range(P̄k) = Kk(A, r0). As range(Xk) = Kk(A, r0) and range(Rk) = AKk(A, r0) has been proved in Proposition 1, the relation (i) holds for the k-th iteration.
To prove Q̄Tk Q̄k = Ik, it suffices to show qk ⊥ Q̄k−1, as the equalities Q̄Tk−1Q̄k−1 = Ik−1 and ‖qk‖2 = 1 has already held. It is equivalent to prove q̃ ⊥ Q̄k−1. From the construction of q̃ in Line 6 in Algorithm 2, we know QTk−1q̃ = 0, so q̃ ⊥ span(qk−2, qk−1) (for k = 2, q̃ ⊥ q0 = 0 clearly holds). To further prove q̃ ⊥ range(Q̄k−3)(k ≥ 4), note that
∆rk−1 = −A∆xk−1 = APk−1Γk−1 − βk−1Ar̄k−1 = −Qk−1Γk−1 − βk−1Ar̄k−1, where the second equality is a direct substitution with (30). Therefore,
Q̄Tk−3∆rk−1 = −Q̄Tk−3Qk−1Γk−1 − βk−1Q̄Tk−3Ar̄k−1 = 0− βk−1(AQ̄k−3)Tr̄k−1 = 0, (31)
where the second equality is due to range(Qk−1) = span(qk−2, qk−1) ⊥ range(Q̄k−3) and A is SPD, the third equality is due to r̄k−1 ⊥ range(Q̄k−1) = AKk−1(A, r0) and range(AQ̄k−3) = A2Kk−3(A, r0) ⊆ AKk−1(A, r0). As a result, noting that q = ∆rk−1, we obtain
Q̄Tk−3q̃ = Q̄ T k−3q − Q̄Tk−3Qk−1(QTk−1q) = 0,
which is due to (31) and range(Qk−1) = span{qk−2, qk−1} ⊥ range(Q̄k−3). Therefore, we show that Q̄Tk Q̄ = Ik, which along with Q̄k = −AP̄k proves relation (ii) in the k-th iteration. Next, we prove the relation (iv). We have
rk = r̄k−1 − βk−1Ar̄k−1 = rk−1 −Qk−1Γk−1 − βk−1A(rk−1 −Qk−1Γk−1) = rk−1 −Qk−1Γk−1 − βk−1Ark−1 + βk−1AQk−1Γk−1 = r0 + Q̄k−1Γ̄k−1 −Qk−1Γk−1 − βk−1(Ar0 +AQ̄k−1Γ̄k−1) + βk−1AQk−1Γk−1,
where the last equality is due to rk−1 = r0 + Q̄k−1Γ̄k−1 by the inductive hypothesis. Since range(Q̄k−1) = AKk−1(A, r0) ⊆ AKk(A, r0), range(Qk−1) ⊆ range(Q̄k), span{Ar0} ⊆ AKk(A, r0), range(AQk−1) ⊆ range(AQ̄k−1) ⊆ A2Kk−1(A, r0) ⊆ AKk(A, r0), and range(Q̄k) = AKk(A, r0), it is clear that rk = r0 + Q̄kΓ̄k ∈ Kk+1(A, r0) for some Γ̄k ∈ Rk. The relation (iv) is proved.
Finally, we prove the relation (iii). For proving r̄k ⊥ range(Q̄k), note that r̄k ⊥ span{qk−1, qk} = range(Qk) already holds due to the projection step (Line 10 and Line 11 in Algorithm 2). It suffices to prove r̄k ⊥ range(Q̄k−2). In fact, since we have r̄k = rk − QkΓk, we can prove that rk ⊥ range(Q̄k−2) and QkΓk ⊥ range(Q̄k−2): Since range(Qk) = span{qk−1, qk} ⊥ range(Q̄k−2) as induced from Q̄Tk Q̄k = Ik, it is clear that QkΓk ⊥ range(Q̄k−2). For rk, according to Line 12 in Algorithm 2, we have rk = r̄k−1 − βk−1Ar̄k−1. We have r̄k−1 ⊥ range(Q̄k−1) ⊇ range(Q̄k−2) by the inductive hypothesis. Also, Q̄Tk−2Ar̄k−1 = (AQ̄k−2)
Tr̄k−1 = 0 due to range(AQ̄k−2) = A2Kk−2(A, r0) ⊆ AKk−1(A, r0) = range(Q̄k−1).
Therefore, we obtain r̄k ⊥ range(Q̄k−2), which along with r̄k ⊥ span{qk−1, qk} implies r̄k ⊥ range(Q̄k).
To prove x̄k = xGk := x0 + zk, where zk = arg minz∈Kk(A,r0) ‖r0 − Az‖2, first we have rk = r0 + Q̄kΓ̄k, where Γ̄k ∈ Rk. Hence r̄k = rk −QkΓk = r0 + Q̄kΓ̄k −QkΓk = r0 − Q̄kηk, where ηk ∈ Rk. Since r̄k ⊥ Q̄k, r̄k is the orthogonal projection of r0 onto the subspace range(Q̄k)⊥. On the other side, for GMRES, rGk = r0 − Azk ⊥ AKk(A, r0) = range(Q̄k), so rGk is also the
orthogonal projection of r0 onto the subspace range(Q̄k)⊥. So r̄k = rGk , which further indicates x̄k = x G k . Hence, the relation (iii) holds.
With relations (i)-(iv) being proved in the k-th iteration, we complete the induction.

heading: C.2 MODIFIED ST-AM FOR GENERAL FIXED-POINT ITERATIONS
text: Algorithm 2 is suitable for analysis and implementation in the linear case. For general nonlinear fixed-point iterations, we adopt an alternative form as described in Algorithm 3 which discards the normalization of q̃ in each iteration (Line 7 in Algorithm 2). In Line 7 in Algorithm 3, the orthogonal projection of ∆rk−1 is checked to ensure ∆rk−1 is “less linearly dependent” on range(Qk), which ensures ‖qk‖2 is bounded away from zero; the check of ‖Pk−1ζk‖2 ensures that ‖pk −∆xk−1‖2 ≤ cp‖∆xk−1‖2, which is also important since a large deviation from ∆xk−1 can make ‖pk‖2 > ρ (ρ is the radius introduced in Theorem 2). When this condition cannot be satisfied, the algorithm simply reuses the old Pk−1, Qk−1. The main procedure of MST-AM restarts every m iterations, i.e. Pk, Qk = 0. Such restart mechanism is to restrict the higher-order terms in the residual expansion, as shown in (9). Also, restart can flush out the outdated historical information that may weaken the quality of Pk and Qk that are used to pursue a local first-order approximation of g in MST-AM.
Algorithm 3 MST-AM for nonlinear fixed-point problems Input: x0 ∈ Rd, βk ∈ (0, 1], cp > 0, cq ∈ (0, 1),m > 0. Output: x ∈ Rd
1: P0, Q0 = 0 ∈ Rd×2, p0, q0 = 0 ∈ Rd 2: for k = 0, 1, . . . , until convergence do 3: rk = g(xk)− xk 4: if k mod m 6= 0 then 5: p = xk − xk−1, q = rk − rk−1 6: ζk = (QTk−1Qk−1)
†QTk−1q 7: if ‖Pk−1ζk‖2 ≤ cp‖p‖2 and ‖Qk−1ζk‖2 ≤ cq‖q‖2 then 8: pk = p− Pk−1ζk, qk = q −Qk−1ζk 9: Pk = [pk−1, pk], Qk = [qk−1, qk] (qk ⊥ Qk−1)
10: else 11: Pk = Pk−1, Qk = Qk−1 12: end if 13: else 14: Pk, Qk = 0 ∈ Rd×2, pk, qk = 0 ∈ Rd 15: end if 16: Γk = (QTkQk)
†QTk rk 17: x̄k = xk − PkΓk, r̄k = rk −QkΓk (Projection step: r̄k ⊥ Qk) 18: xk+1 = x̄k + βkr̄k (Mixing step) 19: end for 20: return xk
(The notation “†” is the Moore-Penrose pseudoinverse.)
In the linear case, Algorithm 2 and Algorithm 3 (with m = ∞) are equivalent. Similar to Algorithm 2, we have the following properties held for MST-AM:
Claim 1. In the k-th iteration (k > 0) of Algorithm 1 applied to minimize a strongly convex quadratic problem (5), assuming cp =∞, cq = 1,m =∞, the following relations hold: (i) ‖qk‖2 > 0, range(P̄k) = range(Xk) = Kk(A, r0), range(Q̄k) = range(Rk) = AKk(A, r0); (ii) Q̄k = −AP̄k, qi ⊥ qj(1 ≤ i 6= j ≤ k); (iii) r̄k ⊥ range(Q̄k) and x̄k = x0 + zk, where zk = arg minz∈Kk(A,r0) ‖r0 −Az‖2. If ‖r̄k‖2 = 0, then xk+1 is the exact solution.
The proof of Claim 1 is essentially the same as the proof of Theorem 1, with a special care that Q̄Tk Q̄k = Ik is replaced by the relation that columns of Q̄k are orthogonal to each other, in other words, Q̄Tk Q̄k = diag{‖q1‖22, . . . , ‖qk‖22}.
For minimizing nonlinear functions or accelerating nonlinear fixed-point iterations, the long-term relation that Q̄Tk Q̄k = diag{‖q1‖22, . . . , ‖qk‖22} generally cannot hold, while the orthogonalization procedure in Line 6 still leads to a short-term orthogonalization relation: qi ⊥ qj for |i − j| ≤ 2. Hence, QTkQk = diag{‖qk−1‖22, ‖qk‖22} for k ≥ 1. Thus the pseudoinverse (QTkQk)
† = diag{I(‖qk‖2 6= 0)1/‖qk‖22, I(‖qk‖2 6= 0)1/‖qk‖22}, where I(·) is the indicator function that
I(x) = {
1 x is true, 0 x is false.
Now, we give the proof of Theorem 2.
Proof of Theorem 2. For convenience, we restate the main assumptions of g here:
(i) ‖g(y)− g(x)‖2 ≤ κ‖y − x‖2, κ ∈ (0, 1), for ∀x, y ∈ B(ρ), (32a) (ii) ‖g′(y)− g′(x)‖2 ≤ κ̂‖y − x‖2, κ̂ > 0, for ∀x, y ∈ B(ρ). (32b)
The proof is based on the two lemmas given in Lemma 1 and Lemma 2. With (45b), we know if ‖rk‖2 is small, then ‖xk − x∗‖2 is small as well. Hence we can always ensure xk ∈ B(ρ) if we choose ‖x0 − x∗‖2 to be sufficiently small and ‖rk‖2 is non-increasing during the iterations. At the end of the k-th iteration of Algorithm 3, we have
rk+1 = g(xk+1)− xk+1 = g(xk+1)− g(x̄k) + g(x̄k)− (x̄k + βkr̄k) = (g(xk+1)− g(x̄k)) + (g(x̄k)− x̄k − r̄k) + (1− βk)r̄k. (33)
Let Lk := g(xk+1)− g(x̄k) + (1− βk)r̄k,Hk := g(x̄k)− x̄k − r̄k, then ‖Lk‖2 ≤ κ‖xk+1 − x̄k‖2 + |1− βk|‖r̄k‖2
= κβk‖r̄k‖2 + |1− βk|‖r̄k‖2 = θk(κβk + |1− βk|)‖rk‖2, (34)
which bounds the linear part of the residual rk+1.
For the higher-order termsHk, we have Hk = g(x̄k)− (xk − PkΓk + rk −QkΓk)
= g(x̄k)− g(xk) + (Pk +Qk)Γk.
= g(x̄k)− g(xk) + (pk + qk)γ(2)k + (pk−1 + qk−1)γ (1) k , (35)
where Γk = (γ (1) k , γ (2) k ) T ∈ R2. Since Γk = (QTkQk)†QTk rk, QTkQk = diag{‖qk−1‖22, ‖qk‖22}, it follows that
γ (1) k = I(qk−1 6= 0) rTk qk−1 qTk−1qk−1 , γ (2) k = I(qk 6= 0) rTk qk qTk qk . (36)
Therefore,
|γ(1)k | ≤ I(qk−1 6= 0) ‖rk‖2 ‖qk−1‖2
, |1− γ(1)k | ≤ max { I(qk−1 = 0), I(qk−1 6= 0)
‖rk − qk−1‖2 ‖qk−1‖2
} ,
|γ(2)k | ≤ I(qk 6= 0) ‖rk‖2 ‖qk‖2
, |1− γ(2)k | ≤ max { I(qk−1 = 0), I(qk 6= 0)
‖rk − qk‖2 ‖qk‖2
} . (37)
According to the formula ∫ 1
0 g′(x+ t(y − x))(y − x)dt = g(y)− g(x), we have
g(x̄k)− g(xk) = g(x̄k)− g(xk − pkγ(2)k ) + g(xk − pkγ (2) k )− g(xk)
= ∫ 1 0 −g′(xk − pkγ(2)k − tpk−1γ (1) k )pk−1γ (1) k dt
+ ∫ 1 0 −g′(xk − tpkγ(2)k )pkγ (2) k dt. (38)
Also, by Lemma 2, we have
pk + qk = ∫ 1 0 g′(xk − tpk)pkdt+ κ̂‖∆rπ(k)−1‖2 k−1∑
j=k−mk
O(‖∆rj‖2),
pk−1 + qk−1 = ∫ 1 0 g′(xk−1 − tpk−1)pk−1dt+ κ̂‖∆rυ(k)−1‖2 k−2∑
j=k−mk
O(‖∆rj‖2),
where π(k) denotes that the latest update of qk by Line 8 occurred in the π(k)-th iteration and υ(k) = π(π(k)−1) marks that qk−1 records qυ(k) that is the penultimate update by Line 8 occurring in the υ(k)-th iteration.
By substituting these relations to (35), it follows that
Hk = ∫ 1 0 g′(xk − tpk)pkdt+ κ̂‖∆rπ(k)−1‖2 k−1∑
j=k−mk
O(‖∆rj‖2)  γ(2)k + ∫ 1 0 g′(xk−1 − tpk−1)pk−1dt+ κ̂‖∆rυ(k)−1‖2 k−2∑
j=k−mk
O(‖∆rj‖2)  γ(1)k − ∫ 1
0
g′(xk − pkγ(2)k − tpk−1γ (1) k )pk−1γ (1) k dt− ∫ 1 0 g′(xk − tpkγ(2)k )pkγ (2) k dt
= ∫ 1 0 (g′(xk − tpk)− g′(xk − tpkγ(2)k ))pkγ (2) k dt
+ ∫ 1 0 (g′(xk−1 − tpk−1)− g′(xk − pkγ(2)k − tpk−1γ (1) k ))pk−1γ (1) k dt
+ κ̂‖∆rπ(k)−1‖2 k−1∑
j=k−mk
O(‖∆rj‖2)γ(2)k
+ κ̂‖∆rυ(k)−1‖2 k−2∑
j=k−mk
O(‖∆rj‖2)γ(1)k
= Ak +Bk + Ck +Dk. (39)
Then we can bound each terms of Hk as follows (Here we assume qk 6= 0 and qk−1 6= 0 as qk = 0 leads to γ(2)k = 0 and qk−1 = 0 leads to γ (1) k = 0 where the result is trivial.):
‖Ak‖2 = ∥∥∥∥∫ 1
0
(g′(xk − tpk)− g′(xk − tpkγ(2)k ))pkγ (2) k dt ∥∥∥∥ 2
≤ ∫ 1
0
‖g′(xk − tpk)− g′(xk − tpkγ(2)k )‖2‖pkγ (2) k ‖2dt
≤ κ̂ ∫ 1
0
‖tpk − tpkγ(2)k ‖2‖pk‖2|γ (2) k |dt
= κ̂
2 ‖pk‖22|1− γ (2) k ||γ (2) k |
= κ̂
2 ‖pk‖22 O(‖rk‖2‖rk − qk‖2) ‖qk‖22
= κ̂ k∑ j=k−mk O(‖rj‖22), (40)
where the third equality is due to (37), and the last equality is due to (46) and (47b); ‖Bk‖2 = ∥∥∥∥∫ 1
0
(g′(xk−1 − tpk−1)− g′(xk − pkγ(2)k − tpk−1γ (1) k ))pk−1γ (1) k dt ∥∥∥∥ 2
≤ ∫ 1
0
‖g′(xk−1 − tpk−1)− g′(xk − pkγ(2)k − tpk−1γ (1) k )‖2‖pk−1γ (1) k ‖2dt
≤ κ̂ ∫ 1
0
‖∆xk−1 − pkγ(2)k − tpk−1γ (1) k + tpk−1‖2‖pk−1‖2|γ (1) k |dt
≤ κ̂ ( ‖∆xk−1‖2 + ‖pk‖2|γ(2)k |+ 1
2 ‖pk−1‖2|1− γ(1)k |
) ‖pk−1‖2|γ(1)k |
≤ κ̂ ( ‖∆xk−1‖2 + ‖pk‖2
‖rk‖2 ‖qk‖2 + 1 2 ‖pk−1‖2 ‖rk − qk−1‖2 ‖qk−1‖2
) ‖pk−1‖2
‖rk‖2 ‖qk−1‖2
= κ̂‖rk‖2 (O(‖∆rk−1‖2) +O(‖rk‖2) +O(‖rk − qk−1‖2)
= κ̂ k∑ j=k−mk O(‖rj‖22), (41)
where the last inequality is due to (37), and the second equality is due to (45a), (46), (47b);
‖Ck‖2 = κ̂‖∆rπ(k)−1‖2 k−1∑
j=k−mk
O(‖∆rj‖2)γ(2)k
≤ κ̂‖∆rπ(k)−1‖2 ‖rk‖2 ‖qk‖2 k−1∑ j=k−mk O(‖∆rj‖2)
= κ̂‖∆rπ(k)−1‖2 ‖rk‖2 ‖qπ(k)‖2 k−1∑ j=k−mk O(‖∆rj‖2)
= κ̂ k∑ j=k−mk O(‖rj‖22), (42)
where the first inequality is from (37), and the second equality is due to (47b);
‖Dk‖2 = κ̂‖∆rυ(k)−1‖2 k−2∑
j=k−mk
O(‖∆rj‖2)γ(1)k
≤ κ̂‖∆rυ(k)−1‖2 ‖rk‖2 ‖qk−1‖2 k−2∑ j=k−mk O(‖∆rj‖2)
= κ̂‖∆rυ(k)−1‖2 ‖rk‖2 ‖qυ(k)‖2 k−2∑ j=k−mk O(‖∆rj‖2)
= κ̂ k∑ j=k−mk O(‖rj‖22), (43)
where the first inequality is from (37), and the second equality is due to (47b). Then with the bounds (40), (41), (42) and (43), we obtain
‖Hk‖2 = κ̂ k∑
j=k−mk
O(‖rj‖22). (44)
Combining (34) and (44) to (33), we obtain
‖rk+1‖ ≤ ‖Lk‖2 + ‖Hk‖2 ≤ θk(κβk + |1− βk|)‖rk‖2 + κ̂ k∑
j=k−mk
O(‖rj‖22),
as desired. Since mk ≤ m, the higher-order terms are limited. Then, for ‖x0 − x∗‖2 sufficiently small, residual rk is non-increasing during the iterations.
Remark 6. There may be some concern about if κ̂ can be counteracted by the constant hidden in the Big-O notation. In fact, since κ̂ is the Lipschitz constant of g′ and the constants in O(·) are composed of κ, cp, cq , it follows that κ̂ is unrelated to O(·). Hence a small κ̂ can lead to a small uniform boundedness of the higher-order terms. In the extreme case where κ̂ = 0 , e.g. g is a linear map, the residual only consists of the first-order term Lk. Lemma 1. Under the same assumptions of Theorem 2, for k ≥ 1, we have the following bounds:
(1− κ)‖∆xk−1‖2 ≤ ‖∆rk−1‖2 ≤ (1 + κ)‖∆xk−1‖2, (45a) (1− κ)‖xk − x∗‖2 ≤ ‖rk‖2 ≤ (1 + κ)‖xk − x∗‖2, (45b)
If qk 6= 0, then ‖pk‖2 ‖qk‖2 ≤ 1 + cp (1− κ)(1− cq)
(46)

heading: If the condition in Line 7 is true, then
text: ‖pk‖2 ≤ (1 + cp)‖∆xk−1‖2, (47a) qk 6= 0, (1− cq)‖∆rk−1‖ ≤ ‖qk‖2 ≤ ‖∆rk−1‖2 (47b)
Proof. From the assumption (32a) of g, we have
(1− κ)‖∆xk−1‖2 ≤ ‖xk − xk−1‖2 − ‖g(xk)− g(xk−1)‖2 ≤ ‖g(xk)− g(xk−1)− (xk − xk−1)‖2 = ‖rk − rk−1‖2 = ‖∆rk−1‖2 ≤ ‖g(xk)− g(xk−1)‖2 + ‖xk − xk−1‖ ≤ (1 + κ)‖xk − xk−1‖2,
(1− κ)‖xk − x∗‖2 ≤ ‖xk − x∗‖2 − ‖g(xk)− g(x∗)‖2 ≤ ‖g(xk)− g(x∗)− (xk − x∗)‖2 = ‖rk‖2 ≤ ‖g(xk)− g(x∗)‖2 + ‖xk − x∗‖2 ≤ (1 + κ)‖xk − x∗‖2.
If the condition in Line 7 is true, i.e., ‖Pk−1ζk‖2 ≤ cp‖∆xk−1‖2, ‖Qk−1ζk‖2 ≤ cq‖∆rk−1‖2, we have
‖pk‖2 = ‖∆xk−1 − Pk−1ζk‖2 ≤ ‖∆xk−1‖2 + ‖Pk−1ζk‖2 ≤ (1 + cp)‖∆xk−1‖2, ‖qk‖2 = ‖∆rk−1 −Qk−1ζk‖2 ≥ ‖∆rk−1‖2 − ‖Qk−1ζk‖2 ≥ (1− cq)‖∆rk−1‖2.
The inequality ‖qk‖2 ≤ ‖∆rk−1‖ is due to the fact that qk is the orthogonal projection of ∆rk−1 onto range(Qk−1)⊥. Also, qk 6= 0 must hold otherwise q = Qk−1ζk which violates the condition ‖Qk−1ζk‖2 ≤ cq‖∆rk−1‖2 as cq ∈ (0, 1). If qk 6= 0, then qk must be updated by Line 8 in some previous iteration. We assume the latest update by Line 8 occurred in the j-th iteration, i.e., qk = qj , pk = pj , hence
‖pk‖2 ‖qk‖2 = ‖pj‖2 ‖qj‖2 ≤ (1 + cp)‖∆xj−1‖2 (1− cq)‖∆rj−1‖2 ≤ 1 + cp (1− cq)(1− κ) .
where the first inequality is due to (47a) and (47b) and the second inequality is due to (45a).
Lemma 2. Under the same assumptions of Theorem 2, in the k-th iteration (k ≥ 0) of the restarted MST-AM (Algorithm 3), we have
pk + qk = ∫ 1 0 g′(xk − tpk)pkdt+ κ̂‖∆rπ(k)−1‖2 k−1∑
j=k−mk
O(‖∆rj‖2), (48)
where mk = k mod m, ∆rk−mk−1 := rk−mk , π(k) denotes that the latest update of qk by Line 8 occurred in the π(k)-th iteration and π(k) = k−mk if Line 8 is never executed up to the k-iteration.
Proof. We prove (48) by induction. Denote ζk = (ζ (1) k , ζ (2) k ) T.
For k = 0, the relation trivially holds.
For k = 1, p1 = ∆x0, q1 = ∆r0. It follows that p1 + q1 = g(x1)− g(x0) = ∫ 1
0
g′(x1 − tp1)p1dt.
So (48) holds for k = 1.
Suppose that (48) holds for 0 ≤ j ≤ k − 1, where k ≥ 2. For k ≥ 2, if mk := k mod m = 0 or 1, (48) holds because it is the same as the case of k = 0 or k = 1. Now consider the nontrivial cases. For mk 6= 0 and mk 6= 1, if the condition in Line 7 in Algorithm 3 is true, then π(k) = k. With the convention that ∆rk−mk−1 = rk−mk , we have
pk + qk = ∆xk−1 − Pk−1ζk + ∆rk−1 −Qk−1ζk = g(xk)− g(xk−1)− (pk−1 + qk−1)ζ(2)k − (pk−2 + qk−2)ζ (1) k
= ∫ 1 0 g′(xk − t∆xk−1)∆xk−1dt
− ∫ 1
0
g′(xk−1 − tpk−1)pk−1ζ(2)k dt+ κ̂‖∆rπ(k−1)−1‖2 k−2∑
j=k−1−mk−1
O(‖∆rj‖2)ζ(2)k
− ∫ 1
0
g′(xk−2 − tpk−2)pk−2ζ(1)k dt+ κ̂‖∆rυ(k−1)−1‖2 k−3∑
j=k−2−mk−2
O(‖∆rj‖2)ζ(1)k ,
(49)
where υ(k− 1) = π(π(k− 1)− 1) denotes that qk−2 records qυ(k−1) that is the penultimate update by Line 8 occurring in the υ(k − 1)-th iteration. Considering the terms in pk + qk, we know
g′(xk − t∆xk−1)∆xk−1 − g′(xk−1 − tpk−1)pk−1ζ(2)k − g ′(xk−2 − tpk−2)pk−2ζ(1)k
= g′(xk − t∆xk−1)(∆xk−1 − pk−1ζ(2)k − pk−2ζ (1) k )
+ (g′(xk − t∆xk−1)− g′(xk−1 − tpk−1))pk−1ζ(2)k + (g′(xk − t∆xk−1)− g′(xk−2 − tpk−2))pk−2ζ(1)k . (50)
For ζk = (QTk−1Qk−1) †QTk−1∆rk−1, because Q T k−1Qk−1 = diag{‖qk−2‖22, ‖qk−1‖22}, it follows that
ζ (1) k = I(qk−2 6= 0)
∆rTk−1qk−2
qTk−2qk−2 , ζ
(2) k = I(qk−1 6= 0)
∆rTk−1qk−1
qTk−1qk−1 .
Therefore,
|ζ(1)k | ≤ I(qk−2 6= 0) ‖∆rk−1‖2 ‖qk−2‖2 , |ζ(2)k | ≤ I(qk−1 6= 0) ‖∆rk−1‖2 ‖qk−1‖2 . (51)
Now, to bound the terms in (50), we have
‖(g′(xk − t∆xk−1)− g′(xk−1 − tpk−1))pk−1ζ(2)k ‖2 ≤ κ̂‖xk − xk−1 − t(∆xk−1 − pk−1)‖2‖pk−1ζ(2)k ‖2 ≤ κ̂((1− t)‖∆xk−1‖2 + t‖pk−1‖2)‖pk−1‖2|ζ(2)k | ≤ κ̂((1− t)‖∆xk−1‖2 + t‖pk−1‖2)‖pk−1‖2I(qk−1 6= 0) ‖∆rk−1‖2 ‖qk−1‖2
≤ κ̂‖∆rk−1‖2 k−1∑
j=k−1−mk−1
O(‖∆rj‖2) = κ̂‖∆rk−1‖2 k−1∑
j=k−mk
O(‖∆rj‖2), (52)
where the third inequality is due to (51), and the last inequality is due to (45a), (47a) and (46), where the constants are absorbed into the big-O notation. Similarly, it follows that
‖(g′(xk − t∆xk−1)− g′(xk−2 − tpk−2))pk−2ζ(1)k ‖2 ≤ κ̂‖∆xk−1 + ∆xk−2 − t∆xk−1 + tpk−2‖2‖pk−2ζ(1)k ‖2 ≤ κ̂((1− t)‖∆xk−1‖2 + ‖∆xk−2‖2 + t‖pk−2‖2)‖pk−2‖2|ζ(1)k | ≤ κ̂((1− t)‖∆xk−1‖2 + ‖∆xk−2‖2 + t‖pk−2‖2)‖pk−2‖2I(qk−2 6= 0) ‖∆rk−1‖2 ‖qk−2‖2
≤ κ̂‖∆rk−1‖2 k−1∑
j=k−2−mk−2
O(‖∆rj‖2) = κ̂‖∆rk−1‖2 k−1∑
j=k−mk
O(‖∆rj‖2). (53)
For the remaining terms in (49), we have
κ̂‖∆rπ(k−1)−1‖2 k−2∑
j=k−1−mk−1
O(‖∆rj‖2)ζ(2)k + κ̂‖∆rυ(k−1)−1‖2 k−3∑
k=k−2−mk−2
O(‖∆rj‖2)ζ(1)k
≤ κ̂‖∆rπ(k−1)−1‖2I(qk−1 6= 0) ‖∆rk−1‖2 ‖qk−1‖2 k−2∑ j=k−1−mk−1 O(‖∆rj‖2)
+ κ̂‖∆rυ(k−1)−1‖2I(qk−2 6= 0) ‖∆rk−1‖2 ‖qk−2‖2 k−3∑ k=k−2−mk−2 O(‖∆rj‖2)
≤ κ̂‖∆rk−1‖2 k−2∑
k=k−mk
O(‖∆rj‖2), (54)
where the last inequality is due to (47b) and that qk−1 = qπ(k−1), qk−2 = qυ(k−1).
With relations (49), (50), (52), (53) and (54), and noting that pk = ∆xk−1 − pk−1ζ(2)k − pk−2ζ (1) k , and∥∥∥∥∫ 1 0 (g′(xk − t∆xk−1)∆xk−1 − g′(xk−1 − tpk−1)pk−1ζ(2)k − g ′(xk−2 − tpk−2)pk−2ζ(1)k )dt ∥∥∥∥ 2
≤ ∫ 1
0
‖g′(xk − t∆xk−1)∆xk−1 − g′(xk−1 − tpk−1)pk−1ζ(2)k − g ′(xk−2 − tpk−2)pk−2ζ(1)k ‖2dt,
we can estimate pk + qk as
pk + qk = ∫ 1 0 g′(xk − t∆xk−1)pkdt+ κ̂‖∆rk−1‖2 k−1∑
j=k−mk
O(‖∆rj‖2). (55)
To further obtain (48), notice that the difference∥∥∥∥∫ 1 0 g′(xk − tpk)pkdt− ∫ 1 0 g′(xk − t∆xk−1)pkdt ∥∥∥∥
2
= ∥∥∥∥∫ 1 0 (g′(xk − tpk)− g′(xk − t∆xk−1))pkdt ∥∥∥∥
2 ≤ ∫ 1
0
‖g′(xk − tpk)− g′(xk − t∆xk−1)‖2‖pk‖2dt
≤ κ̂ ∫ 1
0
t‖pk−1ζ(2)k + pk−2ζ (1) k ‖2‖pk‖2dt
≤ κ̂ 2 cp‖∆xk−1‖2‖pk‖2 = κ̂‖∆rk−1‖2O(‖∆rk−1‖2),
where the last inequality is due to the condition ‖Pk−1ζk‖2 ≤ cp‖p‖2 and inequalities (45a) and (47a). Hence the difference can be absorbed in the Big-O notation in (55). Thus (48) holds when the condition in Line 7 is true for the k-th iteration.
We consider the case where the condition in Line 7 is false for the k-th iteration. Then Pk = Pk−1, Qk = Qk−1. In other words, the memory recording pk, qk keeps unchanged from the (π(k)+ 1)-th to k-th iteration. Since π(k) < k and π(π(k)) = π(k), with the inductive hypothesis, we have
pk + qk = pπ(k) + qπ(k)
= ∫ 1 0 g′(xπ(k) − tpπ(k))pπ(k)dt+ κ̂‖∆rπ(k)−1‖2 π(k)−1∑
j=π(k)−mπ(k)
O(‖∆rj‖2)
= ∫ 1 0 g′(xπ(k) − tpk)pkdt+ κ̂‖∆rπ(k)−1‖ k−1∑
j=k−mk
O(‖∆rj‖2). (56)
To further obtain (48), note that the difference∥∥∥∥∫ 1 0 g′(xk − tpk)pkdt− ∫ 1 0 g′(xπ(k) − tpk)pkdt ∥∥∥∥
2 ≤ ∫ 1
0
‖g′(xk − tpk)− g′(xπ(k) − tpk)‖2‖pk‖2dt
≤ κ̂‖xk − xπ(k)‖2‖pk‖2
= κ̂ k−1∑ j=π(k) ‖∆xj‖2‖pk‖2 = κ̂‖∆rπ(k)−1‖2 k−1∑ j=π(k) O(‖∆rj‖2)
can be absorbed into the Big-O notation in (56). Thus (48) also holds when the condition in Line 7 is false for the k-th iteration.
As a result, we complete the induction in the k-th iteration.

heading: C.3 REGULARIZED SHORT-TERM RECURRENCE ANDERSON MIXING

heading: C.3.1 CHECK OF POSITIVE DEFINITENESS
text: We describe the check of positive definiteness in RST-AM, which follows the same procedure of that of SAM (Wei et al., 2021). From Line 11-13 in Algorithm 1, one update of xk RST-AM is xk+1 = xk +Hkrk, where Hk = βkI − αkYkZ†kQTk , Yk = Pk + βkQk, Zk = QTkQk + δkPTk Pk. Hk is generally not symmetric. A critical condition for the convergence analysis of RST-AM is the positive definiteness of Hk, i.e.
sTkHksk ≥ βkµ‖sk‖22, ∀sk ∈ Rd, (57) where µ ∈ (0, 1) is a constant. Next, we give an approach to guarantee it. Let λmin(·) denote the smallest eigenvalue, λmax(·) denote the largest eigenvalue. Since sTkHksk = 1 2s T k (Hk+H T k )sk, Condition (57) is equivalent to λmin ( 1 2 ( Hk +H T k )) ≥ βkµ. With some simple
algebraic operations, we obtain λmin ( 1 2 ( Hk +H T k )) = βk − 12αkλmax(YkZ † kQ T k + QkZ † kY T k ). Let λk def = λmax(YkZ † kQ T k +QkZ † kY
T k ), then Condition (57) is equivalent to αkλk ≤ 2βk(1− µ), (58)
namely, (13) in Remark 4. To check Condition (58), note that
λk = λmax ( (Yk Qk) ( 0 Z†k Z†k 0 )( Y Tk QTk )) = λmax (( Y Tk QTk ) (Yk Qk) ( 0 Z†k Z†k 0 )) . (59)
Since ( Y Tk QTk ) (Yk Qk) , ( 0 Z†k Z†k 0 ) ∈ R4×4, λk can be computed cheaply. This cost is negligible compared with those to form PTk Pk, Q T kQk, which need O(d) flops. After that, to guarantee the positive definiteness, we check if αk satisfies (58) and use a smaller αk if necessary, e.g. αk = 2βk(1− µ)/λk.

heading: C.3.2 PROOFS OF THE THEOREMS
text: We first give the proof of the boundedness of ‖Pk−1ζk‖2 and ‖Qk−1ζk‖2. Lemma 3. For P,Q ∈ Rd×m(d ≥ m), δ > 0, and Z = QTQ+ δPTP , we have
‖PZ†QT‖2 ≤ δ− 1 2 , (60a)
‖QZ†QT‖2 ≤ 1. (60b)
Proof. We first consider the case that Z is nonsingular, then PZ†QT = PZ−1QT, QZ†QT = QZ−1QT. It can be seen that PTP and QTQ are symmetric positive semidefinite, and Z is SPD as it is assumed to be nonsingular. Also, we have δPTP Z and QTQ Z, where the notation “ ” denotes the Loewner partial order, i.e., A B with A,B ∈ Rm×m means that B − A is positive semidefinite. Hence, we have
Z− 1 2 δPTPZ− 1 2 I, Z− 12QTQZ− 12 I,
which implies
‖Z− 12 ( PTP ) Z− 1 2 ‖2 ≤ δ−1,
‖Z− 12 ( QTQ ) Z− 1 2 ‖2 ≤ 1.
Let λmax(·) denote the largest eigenvalue, we have ‖QZ−1QT‖2 = λmax ( QZ−1QT ) = λmax ( QTQZ−1 ) = λmax ( Z− 1 2QTQZ− 1 2 ) ≤ 1,
‖PZ−1QT‖22 = λmax ( PZ−1QTQZ−1PT ) = λmax ( PTPZ−1QTQZ−1
) = λmax ( Z− 1 2 ( PTP ) Z− 1 2 · Z− 12 ( QTQ ) Z− 1 2
) ≤ ‖Z− 12 ( PTP ) Z− 1 2 · Z− 12 ( QTQ ) Z− 1 2 ‖2
≤ ‖Z− 12 ( PTP ) Z− 1 2 ‖2‖Z− 1 2 ( QTQ ) Z− 1 2 ‖2 ≤ δ−1.
Therefore, ‖PZ−1QT‖ ≤ δ− 12 and ‖QZ−1QT‖2 ≤ 1. For the case that Z is singular, it can be proved that
ker(Z) := {x ∈ Rm|Zx = 0} = {x ∈ Rm|Px = 0 and Qx = 0}. (61) In fact, if Px = 0 and Qx = 0, it is obvious that Zx = 0. On the other side, if Zx = 0, then xTZx = 0, i.e. xTQTQx + δxTPTPx = 0. Since 0 PTP, 0 QTQ, δ > 0, it follows that xTQTQx = 0 and δxTPTPx = 0, which further implies Qx = 0 and Px = 0. Hence (61) holds.
Let U1 satisfy UT1 U1 = I and range(U1) = ker(Z), i.e. the orthonormal basis of ker(Z), and U2 satisfy UT2 U2 = I and U T 2 U1 = 0, i.e. the orthonormal basis of ker(Z)
⊥. With the equality (61), we know PU1 = 0, QU1 = 0. Define U = (U1, U2) ∈ Rm×m, then UTU = Im and by direct computation, we have
UTZU = ( 0 0 0 UT2 ZU2 ) ,
where UT2 ZU2 = (QU2) TQU2 + δ(PU2) TPU2 is nonsingular according to the definition of U2. So
Z† = U ( 0 0 0 (UT2 ZU2) −1 ) UT.
As a result, we can further compute PZ†QT and QZ†QT as
PZ†QT = (PU2)(U T 2 ZU2) −1(QU2) T, QZ†QT = (QU2)(U T 2 ZU2) −1(QU2) T.
Then let P̂ = PU2, Q̂ = QU2, and Ẑ = Q̂TQ̂ + δP̂TP̂ , and noting that Ẑ is nonsingular, we can obtain ‖PZ†QT‖2 = ‖P̂ Ẑ−1Q̂T‖2 ≤ δ− 1 2 , ‖QZ†QT‖2 = ‖Q̂Ẑ−1Q̂T‖2 ≤ 1.
As a result of Lemma 3, ‖Pk−1ζk‖2, ‖Qk−1ζk‖2 in Algorithm 1 are bounded by O(‖∆rk−1‖2), as shown in the following corollary. Corollary 2. ‖Pk−1ζk‖2, ‖Qk−1ζ‖2 in Algorithm 1 are bounded, i.e.
‖Pk−1ζk‖2 ≤ (δ(1)k ) − 12 ‖∆rk−1‖2, ‖Qk−1ζk‖2 ≤ ‖∆rk−1‖2, (62)
where ζk = (QTk−1Qk−1 + δ (1) k Pk−1Pk−1) †QTk−1∆rk−1.
Now, we turn to the proofs of the theorems about RST-AM in Section 3.4. For brevity, we use δk to denote δ(2)k , i.e. δk ≡ δ (2) k .
From Assumption 2, for the mini-batch gradient fSk(xk) = 1 nk ∑ i∈Sk fξi(xk), where nk = |Sk|, we have
E[∇fSk(x)|xk] = ∇f(xk), (63a)
E[‖∇fSk(xk)−∇f(xk)‖22|xk] ≤ σ2
nk . (63b)
Note that the update of RST-AM (Line 11-13 in Algorithm 1) can be written as xk+1 = xk +Hkrk, where rk = −∇fSk(xk), and for k ≥ 0,
Hk = βkI − αk (Pk + βkQk) ( QTkQk + δkP T k Pk )† QTk . (64)
To prove the theorems, the critical points are (i) the positive definiteness of the approximate Hessian Hk and (ii) an adequate suppression of the noise in the gradient estimation for the stochastic case.
We first give a lemma related to the projection step. Lemma 4. Suppose that {xk} is generated by RST-AM. If αk ≥ 0, βk > 0, δk > 0, then for any vk ∈ Rd, we have
‖Hkvk‖22 ≤ 2 ( β2k ( 1 + 2α2k − 2αk ) + α2kδ −1 k ) ‖vk‖22. (65)
Proof. The result clearly holds when k = 0 as H0 = β0I . For k ≥ 1, Hkvk = βkvk − (αkPk + αkβkQk)Γk, (66)
where Γk = (QTkQk + δkP T k Pk) †QTk vk solves
min Γ ‖vk −QkΓ‖22 + δk‖PkΓ‖22. (67)
By direct computation, ‖vk−QkΓk‖22 + δk‖PkΓk‖22 = ‖vk‖22−vTkQk(QTkQk+ δkPTk Pk)† ·QTk vk, thus ‖vk −QkΓk‖22 + δk‖PkΓk‖22 ≤ ‖vk‖22. (68) Therefore,
‖Hkvk‖22 = ‖βkvk − (αkPk + αkβkQk)Γk‖22 = ‖βk (vk − αkQkΓk)− αkPkΓk‖22 = ‖βk(1− αk)vk + βkαk(vk −QkΓk)− αkδ − 12 k δ 1 2 k PkΓk‖ 2 2
≤ ( β2k(1− αk)2 + β2kα2k + α2kδ−1k ) · ( ‖vk‖22 + ‖vk −QkΓk‖22 + δk‖PkΓk‖22 ) ≤ ( β2k ( 1 + 2α2k − 2αk ) + α2kδ −1 k ) ( ‖vk‖22 + ‖vk‖22
) = 2 ( β2k ( 1 + 2α2k − 2αk ) + α2kδ −1 k ) ‖vk‖22. (69)
In the above, the first inequality uses the inequality
‖ n∑ i=1 aixi‖22 ≤ ( n∑ i=1 |ai|‖xi‖2 )2 ≤ ( n∑ i=1 a2i )( n∑ i=1 ‖xi‖22 ) , (70)
where ai ∈ R, xi ∈ Rd. The second inequality is based on inequality (68).
With Lemma 4, we can prove the deterministic case of RST-AM.
Proof of Theorem 3. Since 1+2α2k−2αk ≤ 1 for 0 ≤ αk ≤ 1, and δ −1 k ≤ C−1β2, with Lemma 4, we have ‖Hkrk‖22 ≤ 2β2(1 + C−1)‖rk‖22.
Since αk satisfies (58), we have rTkHkrk ≥ βµ‖rk‖22.
Then, under Assumption 1, we have
f(xk+1) ≤ f(xk) +∇f(xk)T(xk+1 − xk) + L
2 ‖xk+1 − xk‖22
= f(xk)− rTkHkrk + L
2 ‖Hkrk‖22
≤ f(xk)− βµ‖rk‖22 + Lβ2(1 + C−1)‖rk‖22 = f(xk)− β ( µ− βL(1 + C−1) ) ‖rk‖22
≤ f(xk)− 1
2 βµ · ‖∇f(xk)‖22, (71)
where the last inequality is due to 0 < β ≤ µ2L(1+C−1) . Thus, f(xk+1) − f(xk) ≤ − 12βµ‖∇f(xk)‖ 2 2.
Summing both sides of this inequality for k ∈ {0, . . . , N − 1} and recalling f(x) > f low in Assumption 1 gives
f low − f(x0) ≤ f(xN )− f(x0) ≤ − 1
2 βµ N−1∑ k=0 ‖∇f(xk)‖22.
Rearranging and dividing further by N yields (15).
The next lemmas and proofs are about the stochastic case. Lemma 5. Suppose that Assumption 2 holds for {xk} generated by RST-AM. In addition, if βk > 0, δk > 0, and αk ≥ 0 and satisfies (13), then
ESk [‖Hkrk‖22] ≤ 2 ( β2k ( 1 + 2α2k − 2αk ) + α2k δk ) · ( ‖∇f(xk)‖22 + σ2 nk ) , (72a)
∇f(xk)TESk [Hkrk] ≤ − 1
2 βkµ‖∇f(xk)‖22 +
1
2
α2k(δ − 12 k + βk) 2
βkµ · σ
2
nk , (72b)
where µ > 0 is the constant introduced in Remark 4.
Proof. (i) From Lemma 4, we have ESk [‖Hkrk‖22] ≤ 2 ( β2k ( 1 + 2α2k − 2αk ) + α2k δk ) ESk [‖rk‖22]. (73)
From Assumption 2, we have
ESk [‖rk‖22] = ESk [‖rk − ESk [rk]‖22] + ‖ESk [rk]‖22 ≤ ‖∇f(xk)‖22 + σ2/nk. (74)
With (73), (74), we obtain (72a).
(ii) Recalling that H0 = β0I , the result holds for k = 0. Define k = ∇fSk(xk) − ∇f(xk) = −rk − ∇f(xk), then Hkrk = Hk (− k −∇f(xk)) . Since αk satisfies (13), we can ensure λmin ( 1 2 ( Hk +H T k )) ≥ βkµ. Thus
∇f(xk)THk∇f(xk) = 1
2 ∇f(xk)T
( Hk +H T k ) ∇f(xk) ≥ βkµ‖∇f(xk)‖22,
which implies ESk [∇f(xk)THk∇f(xk)] ≥ βkµ‖∇f(xk)‖22. (75)
Let Mk = αk (Pk + βkQk) ( QTkQk + δkP T k Pk )† QTk , then Hk = βkI −Mk. With the assumption (63a), i.e. ESk [ k] = 0, we have
ESk [∇f(xk)THk k] = ESk [∇f(xk)T (βk k −Mk k)] = βk∇f(xk)TESk [ k]− ESk [∇f(xk)TMk k] = −ESk [∇f(xk)TMk k].
Using the Cauchy-Schwarz inequality with expectations, we obtain |ESk [∇f(xk)THk k]| = |ESk [∇f(xk)TMk k]| ≤ √ ESk [‖∇f(xk)‖22] √ ESk [‖Mk k‖22]
≤ ‖∇f(xk)‖2 √ ESk [‖Mk k‖22]. (76)
We now bound ‖Mk k‖22. For brevity, let Zk = QTkQk + δkPTk Pk, and N1 = PkZ † kQ T k , N2 = βkQkZ † kQ T k , then
‖Mk‖2 = ‖αk (N1 +N2) ‖2 ≤ αk (‖N1‖2 + ‖N2‖2) ≤ αk(δ − 12 k + βk), (77)
where the last inequality is from Lemma 3.
With (77), we have ‖Mk k‖2 ≤ αk(δ − 12 k + βk)‖ k‖2, which implies
ESk [‖Mk k‖22] ≤ α2k(δ − 12 k + βk) 2ESk [‖ k‖22] ≤ α2k(δ − 12 k + βk)
2 σ 2
nk , (78)
where the last inequality is due to (63b). Now we can obtain the bound of |ESk [∇f(xk)THk k]| as follows (cf. (76)):
|ESk [∇f(xk)THk k]| ≤ ‖∇f(xk)‖2 √ ESk [‖Mk k‖22]
≤ αk(δ − 12 k + βk)‖∇f(xk)‖2 √ ESk [‖ k‖22]
≤ αk(δ − 12 k + βk) σ √ nk ‖∇f(xk)‖2 = √ βkµ‖∇f(xk)‖2 ·
αk(δ − 12 k + βk)√ βkµ σ √ nk
≤ 1 2 βkµ‖∇f(xk)‖22 + 1 2
α2k(δ − 12 k + βk) 2
βkµ · σ
2
nk . (79)
With the inequality (75) and (79), we obtain
∇f(xk)TESk [Hkrk] = −∇f(xk)TESk [Hk ( k +∇f(xk))] = −ESk [∇f(xk)THk∇f(xk)]− ESk [∇f(xk)THk k] ≤ −ESk [∇f(xk)THk∇f(xk)] + |ESk [∇f(xk)THk k]|
≤ −βkµ‖∇f(xk)‖22 + 1
2 βkµ‖∇f(xk)‖22 +
1
2
α2k(δ − 12 k + βk) 2
βkµ · σ
2
nk
= −1 2 βkµ‖∇f(xk)‖22 + 1 2
α2k(δ − 12 k + βk) 2
βkµ · σ
2
nk . (80)
By imposing more restrictions to αk, βk, δk, we can obtain a convenient corollary:
Corollary 3. Suppose that Assumption 2 holds for {xk} generated by RST-AM. C > 0 is the constant in (12). If βk > 0, 0 ≤ αk ≤ min{1, β 1 2 k } and satisfies (58) , then
ESk [‖Hkrk‖22] ≤ 2β2k ( 1 + C−1 ) · ( ‖∇f(xk)‖22 + σ2
nk
) , (81a)
∇f(xk)TESk [Hkrk] ≤ − 1 2 βkµ‖∇f(xk)‖22 + β2k · µ−1
( 1 + C−1 ) σ2 nk . (81b)
Proof. The first result (81a) is easy to obtain by considering (72a) and noticing that 1+2α2k−2αk ≤ 1 when αk ∈ [0, 1] and δ−1k ≤ C−1β2k. Since αk ≤ β 1 2 k , δk ≥ Cβ −2 k and (1+C
− 12 )2 ≤ 2(1+C−1) we have
1
2
α2k(δ − 12 k + βk) 2
βkµ · σ
2
nk ≤ 1 2
βk(C − 12 βk + βk) 2
βkµ · σ
2
nk
= 1
2 µ−1(C− 1 2 + 1)2β2k ·
σ2
nk
≤ β2kµ−1(1 + C−1) σ2
nk .
Substituting it into (72b), we obtain (81b).
Using Corollary 3 we obtain the descent property of RST-AM: Lemma 6. Suppose that Assumptions 1 and 2 hold for {xk} generated by RST-AM. C > 0 is the constant in (12). If 0 < βk ≤ µ4L(1+C−1) , 0 ≤ αk ≤ min{1, β 1 2 k } and satisfies (58), then
ESk [f(xk+1)] ≤ f(xk)− 1
4 βkµ‖∇f(xk)‖22 + β2k
( (L+ µ−1)(1 + C−1) ) σ2 nk . (82)
Proof. According to Assumption 1, we have
f(xk+1) ≤ f(xk) +∇f(xk)T(xk+1 − xk) + L
2 ‖xk+1 − xk‖22
= f(xk) +∇f(xk)THkrk + L
2 ‖Hkrk‖22. (83)
Taking expectation with respect to the mini-batch Sk on both sides of (83) and using Corollary 3 we obtain
ESk [f(xk+1)]
≤ f(xk) +∇f(xk)TESk [Hkrk] + L
2 ESk‖Hkrk‖22
≤ f(xk)− 1
2 βkµ‖∇f(xk)‖22 + β2kµ−1(1 + C−1)
σ2 nk + Lβ2k(1 + C −1)
( ‖∇f(xk)‖22 + σ2
nk ) = f(xk)− βk ( 1
2 µ− βkL(1 + C−1)
) ‖∇f(xk)‖22 + β2k(L+ µ−1)(1 + C−1) σ2
nk . (84)
Then (84) combined with the assumption βk ≤ µ4L(1+C−1) implies (82).
Following the proofs in Wang et al. (2017), we introduce the definition of a supermartingale. Definition 1. Let {Fk} be an increasing sequence of σ-algebras. If {Xk} is a stochastic process satisfying (i) E[|Xk|] <∞, (ii) Xk ∈ Fk for all k, and (iii) E[Xk+1|Fk] ≤ Xk for all k, then {Xk} is called a supermartingale. Proposition 3 (Supermartingale convergence theorem, see, e.g., Theorem 4.2.12 in (Durrett, 2019)). If {Xk} is a nonnegative supermartingale, then limk→∞Xk → X almost surely and E[X] ≤ E[X0].
Now, we prove the convergence theory of RST-AM in the nonconvex stochastic case.
Proof of Theorem 4. (i) Define φk := βkµ4 ‖∇f(xk)‖ 2 2 and L̃ := (L + µ −1)(1 + C−1), γk := f(xk)+ L̃ σ2
n ∑∞ i=k β 2 i . Let Fk be the σ-algebra measuring φk, γk, and xk. From (82) we know that
for any k,
E[γk+1|Fk] = E[f(xk+1)|Fk] + L̃ σ2
n ∞∑ i=k+1 β2i
≤ f(xk)− 1
4 βkµ‖∇f(xk)‖22 + L̃
σ2
n ∞∑ i=k β2i = γk − φk, (85)
which implies that E[γk+1 − f low|Fk] ≤ γk − f low − φk. Since φk ≥ 0, we have 0 ≤ E[γk − f low] ≤ γ0 − f low < +∞. As the diminishing condition (14) holds, we obtain E[f(xk)] ≤Mf for some constant Mf > 0. According to Definition 1, {γk − f low} is a supermartingale. Therefore, Proposition 3 indicates that there exists a γ such that limk→∞ γk = γ with probability 1, and E[γ] ≤ E[γ0]. Note that from (85) we have E[φk] ≤ E[γk]− E[γk+1]. Thus,
E [ ∞∑ k=0 φk ] ≤ ∞∑ k=0 (E[γk]− E[γk+1]) < +∞,
which further yields that ∞∑ k=0 φk = µ 4 ∞∑ k=0 βk‖∇f(xk)‖22 < +∞ with probability 1. (86)
Since ∑∞ k=0 βk = +∞, it follows that (16) holds.
(ii) If the noisy gradient is bounded, i.e.,
Eξk [‖∇fξk(xk)‖22] ≤Mg, (87)
where Mg > 0 is a constant, then a stronger result can be obtained.
For any give > 0, according to (16), there exist infinitely many iterates xk such that ‖∇f(xk)‖2 ≤ . Then if (17) does not hold, there must exist two infinite sequences of indices {si}, {ti} with ti > si, such that for i = 0, 1, . . ., k = si + 1, . . . , ti − 1,
‖∇f(xsi)‖2 ≥ 2 , ‖∇f(xti)‖2 < , ‖∇f(xk)‖2 ≥ . (88)
Then from (86) it follows that
+∞ > ∞∑ k=0 βk‖∇f(xk)‖22 ≥ +∞∑ i=0 ti−1∑ k=si βk‖∇f(xk)‖22 ≥ 2 +∞∑ i=0 ti−1∑ k=si βk with probability 1,
which implies that ti−1∑ k=si βk → 0 with probability 1, as i→ +∞. (89)
According to (74) and (69), we have
E[‖xk+1 − xk‖2|xk] = E[‖Hkrk‖2|xk]
≤ √ 2 ( β2k (1 + 2α 2 k − 2αk) + α2kδ −1 k ) E[‖rk‖2|xk]
≤ βk √ 2(1 + C−1)E[‖rk‖2|xk]
≤ βk √ 2(1 + C−1)(E[‖rk‖22|xk]) 1 2
≤ βk √ 2(1 + C−1)M 1 2 g , (90)
where the last inequalities are due to Cauchy-Schwarz inequality and (87). Then it follows from (90) that
E[‖xti − xsi‖2] ≤ √ 2(1 + C−1)M 1 2 g ti−1∑ k=si βk,
which together with (89) implies that ‖xti − xsi‖2 → 0 with probability 1, as i → +∞. Hence, from the Lipschitz continuity of ∇f , it follows that ‖∇f(xti) − ∇f(xsi)‖2 → 0 with probability 1 as i → +∞. However, this contradicts (88). Therefore, the assumption that (17) does not hold is not true.
Proof of Theorem 5. According to (84) in Lemma 6, we have
N−1∑ k=0 βk ( 1 2 µ− βkL(1 + C−1) ) E‖∇f(xk)‖22
≤ f(x0)− f low + N−1∑ k=0 β2k(L+ µ −1)(1 + C−1) σ2 nk , (91)
where the expectation is taken with respect to {Sj}N−1j=0 . Define
PR(k) def = Prob{R = k} =
βk ( 1 2µ− βkL(1 + C −1) )∑N−1
j=0 βj ( 1 2µ− βjL(1 + C−1) ) , k = 0, . . . , N − 1, (92) then
E [ ‖∇f(xR)‖22 ] =
∑N−1 k=0 βk ( 1 2µ− βkL(1 + C −1) ) E [ ‖∇f(xk)‖22 ]∑N−1 j=0 βj ( 1 2µ− βjL(1 + C−1)
) ≤ Df + σ 2(L+ µ−1)(1 + C−1) ∑N−1 k=0 β 2 k/nk∑N−1
j=0 βj ( 1 2µ− βjL(1 + C−1) ) . (93) Let D̃ be a problem-independent constant. If we choose βk = β := min{ µ4L(1+C−1) , D̃ σ √ N }, and nk = n, then the definition of PR simplifies to PR(k) = 1/N . From (93) we have
E[‖∇f(xR)‖22] ≤ Df + σ
2(L+ µ−1)(1 + C−1)Nβ 2 n∑N−1 j=0 β( 1 2µ− µ 4 )
= Df + σ
2(L+ µ−1)(1 + C−1) · Nβ 2
n
Nβ · 14µ
= 4Df Nβµ
+ σ2(L+ µ−1)(1 + C−1) · β
1 4nµ
≤ 4Df Nµ max
{ 4L(1 + C−1)
µ , σ √ N
D̃
} + 4σ2(L+ µ−1)(1 + C−1)
nµ · D̃ σ √ N
≤ 4Df Nµ
( 4L(1 + C−1)
µ + σ √ N
D̃
) + 4σ(L+ µ−1)(1 + C−1)D̃
nµ √ N
= 16DfL(1 + C
−1)
Nµ2 +
σ
µ √ N
( 4Df
D̃ +
4(L+ µ−1)(1 + C−1)D̃
n
) .
Therefore, to ensure E[‖∇f(xR)‖22] ≤ , the number of iterations is O(1/ 2).

heading: D EXPERIMENTAL DETAILS
text: Our main codes were written based on the PyTorch framework 1 and one GeForce RTX 2080 Ti GPU was used for the tests in training neural networks. Our methods are ST-AM (MST-AM) and the regularized version RST-AM.

heading: D.1 EXPERIMENTAL DETAILS ABOUT ST-AM/MST-AM
text: The experiments about ST-AM were conducted to verify the main theorems, i.e. Theorem 1, Corollary 1 and Theorem 2. Four types of problems were used for the experiments. The conjugate residual (CR) method is a short-term recurrence version of the full-memory GMRES and needs matrix-vector products to fulfill the algorithm. We give the pseudocode of CR method here for readers who are not familiar with this numerical algorithm.
Algorithm 4 CR Algorithm for solving a SPD linear system Ax = b. Input: x0 ∈ Rd. Output: x ∈ Rd
1: r0 = b−Ax0, p0 = r0 2: for k = 0, 1, . . . , until convergence, do 3: αk =
rTkArk (Apk)TApk
4: xk+1 = xk + αkpk 5: rk+1 = rk − αkApk 6: βk = rTk+1Ark+1
rTkArk
7: pk+1 = rk+1 + βkpk 8: Apk+1 = Ark+1 + βkApk 9: end for
10: return xk
It should be pointed out that CR method needs to directly assess the matrixA, and the residual update (Line 5) is based on the linear assumption, which makes it inapplicable for nonlinear problems or the case that A is unavailable. Also, even though finite difference technique can be used to construct the matrix-vector products, the number of gradient evaluations is twice of that of ST-AM.

heading: D.1.1 STRONGLY CONVEX QUADRATIC OPTIMIZATION
text: To construct a case of (5), we considered the least squares problem:
min x∈Rd
f(x) := 1
2 ‖Ax− b‖22, (94)
where A ∈ R`×d, b ∈ R`, which can be reformulated as a form of (5). We generated a dense random matrix A ∈ R500×100 and a dense random vector b ∈ R500 for the test. The gradient descent used a fixed stepsize η ≤ 1‖A‖22 that can guarantee convergence, and the same η was used as the βk for the full-memory AM and ST-AM.
We also conducted additional tests about solving the problem (94) with different condition numbers, where the eigenvalues ofA were set to be uniformly distributed. The results with different condition numbers (cond(ATA) = 102, 104, 106) are shown in Figure 4. It can be observed that the curves of CR, AM and ST-AM nearly coincide, which verify the correctness of Theorem 1. Also, since the eigenvalues are uniformly distributed, the superlinear convergence behaviour may not happen. Nonetheless, CR, AM and ST-AM are still much faster then the GD method.

heading: D.1.2 SOLVING A NONSYMMETRIC LINEAR SYSTEM
text: For the solution of a nonsymmetric linear system Ax = b, (95)
1 Information about this framework is referred to https://pytorch.org.
where A ∈ Rd×d, b ∈ Rd, the fixed-point iteration (FP) is xk+1 = g(xk) := xk + η(b − Axk). Theorem 2 requires g to be a contractive map, i.e. ‖I−ηA‖2 < 1. We used a test matrix “fidap029” from the Matrix Market 2. This matrix is a banded matrix and not symmetric. We used Jacobi preconditioner for all the tested iterative methods. Since solving linear systems is not the main focus of this work, we did not do a thorough test of applying ST-AM to solve various nonsymmetric linear systems.

heading: D.1.3 CUBIC-REGULARIZED QUADRATIC MINIMIZATION
text: The concerned cubic-regularized quadratic minimization is
min x∈Rd
f(x) := 1
2 ‖Ax− b‖22 +
M
3 ‖x‖32, (96)
where A ∈ R`×d, b ∈ R`, and M ≥ 0 is the regularization parameter. It can be computed that
∇f(x) = AT(Ax− b) +M‖x‖2x = (ATA+M‖x‖2I)x−ATb.
Hence, for the gradient descent xk+1 = g(xk) := xk − η∇f(xk), the Jacobian of g is I − ηATA− ηM(‖x‖−12 xxT + ‖x‖2I), which has κ̂ > 0 in the local region B(ρ). We generated a dense random matrix A ∈ R500×100 and a dense vector b ∈ R500 for the test.
We also conducted the tests related to different regularization parameters M , and the cases that M = 0.01, 1, 100 are shown in Figure 5. Figure 5 shows that with the largeM = 100, both AM and MST-AM converge faster. An ablation study was also conducted about the boundedness restriction of ‖Pk−1ζk‖2, ‖Qk−1ζk‖2 in Theorem 2. In Figure 5(b), we show the convergence behaviour of ST-AM without the boundedness check. In the case of the rather large regularization M = 100, STAM does not show a monotone decrease of the residual. To further investigate the cause, we plot the magnitude of ‖Pk−1ζk‖2/‖∆xk−1‖2, ‖Qk−1ζk‖2/‖∆rk−1‖2 in Figure5(c) and Figure 5(d). It can
2 https://math.nist.gov/MatrixMarket/.
be observed that the evolutions of ‖Pk−1ζk‖2/‖∆xk−1‖2, ‖Qk−1ζk‖2/‖∆rk−1‖2 are quite oscillatory in ST-AM, while being roughly bounded below 1 in MST-AM. This phenomenon may accounts for the more stable convergence behaviour of MST-AM and verifies the necessity of the changes of MST-AM compared with ST-AM. In our experiments, we also found the restarting period can be set quite large and has little effect on ST-AM.

heading: D.1.4 ROOT-FINDING PROBLEMS IN THE MULTISCALE DEEP EQUILIBRIUM MODEL
text: The multiscale deep equilibrium (MDEQ) model is a recent extension of the deep equilibrium (DEQ) model (Bai et al., 2019) for computer vision. One of the central engines for these DEQ models is the root-finding problem:
fθ(z;x) = gθ(z;x)− z ⇒ find z∗ s.t. fθ(z∗;x) = 0, (97)
where θ and x are parameters and the input representation, respectively. In (Bai et al., 2020), the Broyden’s method is employed to solve this problem. Since MST-AM is suitable for solving nonlinear systems of equations, we can apply MST-AM to solve (97).
We implemented the MST-AM method and integrated it into the MDEQ framework 3 . The task was image classification on CIFAR-10 and we used the small model for test. We followed the suggested experimental setting of the framework. The optimizer was Adam with learning rate of 0.001 and the weight-decay was 2.5×10−6. The batch size was 128 and the number of epochs was 50. The tested fixed-point solver was used for the forward root-finding process and for the backward root-finding process. The threshold of the steps for the forward process was 18 and the threshold of the steps for the backward process was 20.
We used the built-in AM method and Broyden’s method as the baseline methods. The m for AM was 20, i.e. using the full-memory AM.

heading: D.2 EXPERIMENTAL DETAILS ABOUT RST-AM
text: Our experiments on RST-AM focused on training neural networks. Since ST-AM can be regarded as a special case of RST-AM with δ(1)k = δ (2) k = 0, the basic ST-AM is also covered. In the Line 10 in Algorithm 1, αk should be adjusted to meet the positive definiteness check (13), which can be simplified to (∆xk)Trk ≥ βkµ‖rk‖22 in practice. The adjustment of αk can be (i) αk = min{αk, 2βk(1− µ)/λk}, or (ii) αk = 0. We used the option (ii) since the violation of (13) seldom happened in our tests and option (ii) is more simple to apply.
In the experiments on MNIST, Penn Treebank, we incorporated preconditioning (described in Appendix A.3) into the baseline method SAM and our method. Preconditioning is found to be effective for difficult problems, e.g. mini-batch training with very small batch sizes and the scaling of the model’s parameters being important.

heading: D.2.1 HYPERPARAMETER SETTING OF RST-AM
text: The hyperparameters of RST-AM are easy to tune. The only hyperparameters that need to be carefully tuned are the regularization parameters c1, c2 in δ (1) k and δ (2) k . We found RST-AM is more sensitive to c1, possibly due to the fact that it influences the construction of secant equations. c2 can be quite small in our tests. The hyperparameter C in δ(2)k was set to be very small (C = 1× 10−16) so as to ensure δ(2)k = c2‖rk‖22 ‖pk‖22+ 0
almost all the time. 0 only serves as the safe-guard to prevent the denominator from being zero. We found ‖pk‖2 > 0 and ‖∆xk−1‖2 > 0 always held in the tests, so 0 can be omitted. In the experiments except for deterministic optimization on MNIST, we kept the setting c1 = 1, c2 = 1× 10−7 unchanged and found such setting is quite robust. The other hyperparameters are αk, βk, which can be always initially set as 1. So RST-AM has the same number of hyperparameters to tune as SGDM, since SGDM needs to tune learning rate and the momentum.
3 https://github.com/locuslab/mdeq.

heading: D.2.2 EXPERIMENTS ON MNIST
text: The experiments on MNIST 4 aimed to validate the effectiveness of RST-AM in deterministic optimization, so we were only concerned about the training loss by regarding it as a nonlinear function to be optimized. The training dataset was preprocessed by randomly selecting 10k images from the total 60k images to facilitate the full-batch training.
The baselines were SGDM, Adam, Adagrad, RMSprop and SAM. We tried our best to ensure that the baselines had the best performance in the tests. We tuned the learning rates by log-scale gridsearches from 10−3 to 100. The learning rates of SGDM, Adam, Adagrad, RMSprop were 0.2, 0.001, 0.01, 0.001, respectively. For SAM, we used the hyperparameter setting recommended in (Wei et al., 2021).
For RST-AM, we set c1 = 0.05, c2 = 1 × 10−7, αk = βk = 1. When (∆xk)Trk ≤ 0 occurs, xk+1 = xk + 0.2rk was used as the new update.
For the preconditioned RST-AM, we set c1 = 1 × 10−2, c2 = 1 × 10−7 for the Adagradpreconditioned RST-AM and c1 = 1 × 10−2, c2 = 1 × 10−8 for the RMSprop-preconditioned RST-AM.
For all these tests, the model was trained for 200 epochs.
In Figure 2 in the main paper, we report the RMSprop/Adagrad preconditioned SAM/RST-AM. We give the result about using Adam as the preconditioner in Figure 6. It suggests that Adam does not perform as well as the RMSprop method to serve as a preconditioner for SAM/RST-AM in this task. Nevertheless, Adam RST-AM is still better than Adam and Adam SAM(2), which demonstrates the effect of our proposed short-term recurrence scheme.

heading: D.2.3 EXPERIMENTS ON CIFAR
text: For this group of experiments, we followed the same setting of training ResNet in (He et al., 2016). The batchsize was set to be 128 as commonly suggested. When training forN iterations, the learning rate was decayed at the (bN2 c)-th and the (b 3 4Nc)-th iterations. For SAM and RST-AM, the learning rate decay means decaying αk, βk simultaneously. The experiments were run with 3 random seeds.
The baseline optimizers were SGDM, Adam, AdaBound (Luo et al., 2018), AdaBelief (Zhuang et al., 2020), Lookahead (Zhang et al., 2019) and AdaHessian (Yao et al., 2021). AdaBound and AdaBelief are recently proposed adaptive learning rate methods to improve Adam. Lookahead is a k-step method, which can be seen as a simple sequence interpolation method. In each cycle, Lookahead iterates with an inner-optimizer optim for k steps and then interpolates the first and the last iterates to give the starting point of the next cycle. Compared to Adam, AdaHessian uses Hessian-vector products to construct a diagonal approximation of the Hessian.
Since tuning hyperparameter in large models is very time-consuming, the hyperparameters of each optimizers were tuned through experiments on CIFAR-10/ResNet20. For each optimizer, the hyperparameter setting that has the best final test accuracy on CIFAR-10/ResNet20 was kept unchanged and used for the other tests. We list the hyperparameters of all the tested optimizers here. (Learning rate is abbreviated as lr.)
4Based on the official PyTorch implementation https://github.com/pytorch/examples/blob/master/mnist.
• SGDM: lr = 0.1, momentum = 0.9, weight-decay = 5× 10−4, lr-decay = 0.1. • Adam: lr = 0.001, (β1, β2) = (0.9, 0.999), weight-decay = 5× 10−4, lr-decay = 0.1. • AdaBound: lr = 0.001, (β1, β2) = (0.9, 0.999), final lr = 0.1, gamma = 0.001, weight-
decay = 5× 10−4, lr-decay = 0.1. • AdaBelief: lr = 0.001, (β1, β2) = (0.9, 0.999), eps = 1× 10−8, weight-decay = 5× 10−4,
lr-decay = 0.1.
• Lookahead: optim: SGDM (lr = 0.1, momentum = 0.9, weight-decay = 1 × 10−3), α = 0.8, steps = 10, lr-decay = 0.1.
• AdaHessian: lr = 0.15, (β1, β2) = (0.9, 0.999), eps=1× 10−4, hessian-power: 1, weightdecay = 5× 10−4/0.15, lr-decay = 0.1.
• SAM(2): optim: SGDM (lr = 0.1, momentum = 0, weight-decay = 1.5 × 10−3), αk = 1.0, βk = 1.0, c1 = 0.01, p = 1,m = 2, weight-decay = 1.5× 10−3, lr-decay = 0.06. • SAM(10): optim: SGDM (lr = 0.1, momentum = 0, weight-decay = 1.5 × 10−3), αk = 1.0, βk = 1.0, c1 = 0.01, p = 1,m = 10, weight-decay = 1.5× 10−3, lr-decay = 0.06. • RST-AM: c1 = 1, c2 = 1× 10−7, α0 = β0 = 1, weight-decay = 1× 10−3, lr-decay = 0.1.
For the tests of SGDM, Adam, AdaBelief, Lookahead and SAM(10), we found consistent numerical results as those in (Wei et al., 2021), so we report their results of these methods in the main paper.
Figure 7 shows the test accuracy of different optimizers for training ResNet18/WRN16-4 on CIFAR10 and ResNeXt50/DenseNet121 on CIFAR-100. The full results of final test accuracy are listed in Table 1(a) in the main paper. From Figure 7, we find that the convergence behaviour of RST-AM is less erratic than SAM(10) and SAM(2). In the first 80 epochs, RST-AM converges faster than SAM,
partly due to using a smaller weight decay. The learning rate schedule has a large impact on the convergence of each optimizer. Similar to SAM(10), RST-AM can always climb up and stabilize to a higher test accuracy in the final 40 epochs. It is found that RST-AM is comparable to SAM(10), while improving the short-memory SAM(2).
Since our hyperparameter tuning was conducted on CIFAR10/ResNet20, we give some results about the hyperparameters on this model.
Effect of c2 in RST-AM. Figure 8 shows the effect of c2, where we kept c1 = 1, weight-decay=5× 10−4 fixed, and c2 = 10−1, 10−3, 10−5, 10−7, 10−9. It indicates that RST-AM is not sensitive to c2.
Effect of c1 in RST-AM. Figure 9 shows the effect of c1, where we kept c2 = 10−7, weightdecay=5×10−4 fixed, and c1 = 0.01, 0.1, 1, 10. It suggests that with smaller c1, RST-AM converges faster in terms of train accuracy, but the test accuracy is worse.
Effect of weight decay in RST-AM. Weight-decay is a common hyperparameter that can affect the generalization of each optimizer. In Figure 10, we show the convergence behaviour of RST-AM with different weight-decays. It suggests that with a too small weight-decay, RST-AM tends to be overfitting in the test dataset.
Table 1(a) in the main paper reports the test accuracy of each optimizer when training for the same epochs. In fact, as shown in Figure 11, within 100 epochs, RST-AM can achieve a better test accuracy than SGD. So if the running time of the training process matters, it is expected that RSTAM can use less total running time due to the large number of reduction in training epochs. In
Table 1(b), we set SGD as the baseline, and compare the computation and memory cost with SGD. In Table 5, we give more details about the saving in training epochs and the final test accuracy. It can be seen that RST-AM can achieve a comparable or better test accuracy than SGDM with less computation time, while reducing the memory footprint of SAM.
We also tried using Adam as a preconditioner for RST-AM but found the final test accuracy was often worse. For example, the test accuracy of Adam RST-AM for CIFAR-10/ResNet20 is only 90.79%. We suppose it is the worse generalization ability of Adam (Luo et al., 2018) that makes Adam not suitable as a preconditioner for RST-AM in the image classification task.

heading: D.2.4 EXPERIMENTS ON PENN TREEBANK
text: Our experimental setting on training LSTM models on Penn Treebank dataset was based on the official implementation of AdaBelief (Zhuang et al., 2020). Results in Table 2 were measured with 3 random seeds. The parameter setting of the LSTM models were the same as that of AdaBelief. The baseline optimizers were SGDM, Adam, AdaBelief, and SAM. We tuned hyperparameters on the validation dataset for each optimizer.
For SGDM, we tuned the learning rate (abbr. lr) via grid-search in {1, 10, 30, 100} and found that lr=10 performs best on 2,3-layer LSTM. For 1-layer LSTM, we set lr=30 and momentum=0 as that in AdaBelief because we found such setting is better.
For Adam, we tuned the learning rate via grid-search in {1×10−3, 2×10−3, 5×10−3, 8×10−3, 1× 10−2, 2× 10−2} and found 5× 10−3 performs best. For AdaBelief, we tuned the learning rate and found 5×10−3 is better than 1×10−2 used in (Zhuang et al., 2020).
For SAM, we used the recommended setting in (Wei et al., 2021) and used the baseline Adam as the preconditioner for SAM. Historical lengths of m = 2 and m = 10 were tested for SAM, i.e. denoted as SAM(2), SAM(10), respectively.
For our method RST-AM, we kept c1 = 1, c2 = 1 × 10−7 unchanged and used the same preconditioner (Adam) as SAM.
The batch size was 20. We trained for 500 epochs and decayed the learning rate by 0.1 at the 250th and 375th epoch.
In Table 6, we report the memory footprint and per-epoch running time of SAM(10) and RST-AM compared with Adam. It indicates that RST-AM also largely reduces the memory overhead of the long-memory SAM(10) in the language model task. Since the batch size is very small, the cost of
stochastic gradient evaluation is quite cheap, which makes the additional cost of matrix computation in RST-AM and SAM(10) become considerable. However, if we consider achieving a comparable validation/test perplexity to that of Adam, RST-AM does not need to train for the same number of epochs as Adam. Table 7 shows the test perplexity of Adam, SAM(2), SAM(10), and RST-AM for training 200 epochs. By comparing the results of Table 2 and Table 7, we see RST-AM is better than Adam with much fewer training epochs, thus RST-AM can save a large amount of training time, as shown in Table 6.

heading: D.2.5 EXPERIMENTS ON ADVERSARIAL TRAINING
text: The adversarial training considers the empirical adversarial risk minimization (EARM) problem:
min x∈Rd
1
T T∑ i=1 max ‖ξ̄i−ξi‖2≤ fξ̄i(x), (98)
where ξ̄i is the i-th adversarial data within the -ball centered at ξi. We followed the standard PGD adversarial training in (Madry et al., 2018), using projection gradient descent (PGD) to solve the inner maximization problem and the tested optimizers (SGD and RST-AM) to solve the outer minimization problem. The experiments were conducted on CIFAR10/ResNet18, CIFAR10/WRN16-4, CIFAR100/ResNet18, and CIFAR100/DenseNet121. We trained the neural networks for 200 epochs and decayed the learning rate at the 100th and 150th epoch.
Since adversarial training is much more time consuming than the ordinary training in Section D.2.3, we tuned the hyperparameters in CIFAR10/ResNet20, and applied the same hyperparameters to other models. We found the setting that c1 = 1, c2 = 1 × 10−7 and weight-decay = 0.001 is still suitable for this task.
The CIFAR-10 (CIFAR-100) dataset contains 50K images for training and 10K images for testing. Since it is found that the phenomenon of overfitting is severer (Rice et al., 2020) in adversarial training, we randomly selected 5K images from the total 50K training dataset as the validation dataset (the other 45K images remained as the training dataset), and chose the best checkpoint model in the validation set to evaluate on the test dataset. We consider two types of test accuracy: (i) the clean test accuracy, where clean data was used for model evaluation; (ii) the robust test accuracy, where corrupted data was used for model evaluation. The attacking methods are FGSM (Goodfellow et al., 2014), PGD-20 (Madry et al., 2018), and C&W∞ attack (Carlini & Wagner, 2017).
In Tabel 8, we report the average results of tests with three different random seeds. It can be found that RST-AM can achieve both higher clean test accuracies and robust test accuracies across various models on CIFAR10/CIFAR100. To see the convergence behaviour of SGD and RST-AM, we plot the curves of the clean validation accuracy and the PGD-10 attacked validation accuracy in Figure 12. We can observe the phenomenon of robust overfitting (Rice et al., 2020) from these curves, which justifies our experimental setting with validation set for the checkpoint selection. Nonetheless, the numerical results suggest that RST-AM can still be better than SGDM in adversarial training.
It is also found that due to the heavy cost of gradient evaluations in PGD adversarial training, the additional computational cost incurred by RST-AM is negligible and the per-epoch running time of SGD and RST-AM is roughly the same. So we do not report it here.

heading: D.2.6 EXPERIMENTS ON TRAINING A GENERATIVE ADVERSARIAL NETWORK
text: We describe our setting of training a generative adversarial network (GAN) here. Like the adversarial training, the GAN training process is also a min-max problem. The stability of an optimizer is critical for the training process. To demonstrate the applicability of RST-AM, we conducted experiments on a GAN which was equipped with spectral normalization (Miyato et al., 2018) (SNGAN). The experimental setting was the same as that of AdaBelief (Zhuang et al., 2020): the dataset was CIFAR-10; ResNets were used as the generator and the discriminator networks; the steps for optimization in the discriminator and the generator per iteration were 5 and 1, respectively; the minbatch size was 64 and the total iteration number was 100000. The Frechet Inception Distance (FID) (Heusel et al., 2017) was used as the evaluation metric: lower FID score means better accuracy.
The baseline optimizer were Adam and AdaBelief as they perform well in this task (Zhuang et al., 2020). We also used the recommended hyperparameter settings for the two optimizers. For our RST-AM method, due to the ill-conditioning of the min-max problem, we used the AdaBelief as the preconditioner and set the damping parameter αk = 0.6. The regularization parameters c1 = 1 and c2 = 1× 10−7 were still kept unchanged just as the previous experiments.
In Figure 13, we show the curve of FID score for each optimizer, which is the average of three independent runs. It indicates that the RST-AM method is stable for this min-max optimization problem.
Since we only tuned the damping parameter αk for RST-AM, we report the FID scores of other choices of αk in Table 9 during our experiment. It shows that even with the suboptimal choices of αk, e.g. αk = 0.1, 0.5, 0.8, RST-AM can still outperform the baselines.

heading: D.3 ADDITIONAL EXPERIMENTS
text: To further test the performance of our method in training neural networks on larger datasets or different models, we conducted additional experiments of the image classification task in ImageNet (Deng et al., 2009) and the Transformer (Vaswani et al., 2017) based neural machine translation task in the IWSTL14 DE-EN (Cettolo et al., 2014) dataset.

heading: D.3.1 EXPERIMENTS ON IMAGENET
text: We trained ResNet50 on ImageNet with SGDM and RST-AM. We used the built-in ResNet50 model in PyTorch. We ran the tests of each optimizer with three random seeds and four GeForce RTX 2080
Ti GPUs were used for each test. The hyperparameters of SGDM were set as the recommended setting in PyTorch. For RST-AM, the hyperparameters were kept the same as those in the CIFAR experiments. The weight-decay was 1 × 10−4. The number of the total training epochs is 90. The learning rate decay of SGD was at the 30th and 60th epochs. For RST-AM, since the experiments on CIFAR show it can often converge faster to an acceptable solution than SGDM, we adopted the early learning rate decay strategy recommended in (Zhang et al., 2019): decay the αk and βk for RST-AM at the 30th, 50th and 70th epochs.
In Table 10, we report the TOP1 accuracy in the validation dataset during training. Note that we also report the epoch number for each optimizer to achieve the accuracy equal or exceeding 75.90%. It shows RST-AM needs 72 epochs while SGDM needs 88 epochs. The curves of the training accuracy and test accuracy are shown in Figure 14. The results suggest that RST-AM is still a competitive optimizer in training a larger model in a larger dataset.

heading: D.3.2 TRAINING TRANSFORMER
text: We conducted the neural machine translation task with Transformer. We implemented our RST-AM method and integrated it into the fairseq framework 5. The basic experimental setting was set as the recommended setting in (Yao et al., 2021). We trained the model for 50 epochs and the BLEU (Papineni et al., 2002) score was calculated using the average model of the last five checkpoints. We added a baseline optimizer RAdam (Liu et al., 2019) which was inspired by the warmup procedure in training Transformer.
The numerical results reported in Table 11 show that Adam is well-suited for this neural machine translation task, though it does not perform well in the image classification task. Also, the results demonstrate that RST-AM can still outperform Adam in this task.
In Table 12, we report the BLEU scores evaluated in the test dataset at the 40th, 45th and 50th epochs. The results show the better performance of RST-AM over Adam and the per-epoch computational cost is nearly the same. To achieve a comparable solution to Adam, RST-AM can save 10% training time.

Review summaries:
The paper introduces 3 new variants of Anderson Mixing methods relying on very limited short term memory. This makes these methods more attractive for usual machine learning workloads. The paper also provides detailed analysis of the performance of each of these AM algorithms, as well as thorough experiments, showcasing improved performance of select neural network training.

The paper proposes a new class of memory-efficient Anderson mixing (AM) methods. Compared with classical Anderson mixing which requires saving m historical iterates, the new variants require only storing two historical iterates while keeping good performance. Convergence analyses are given to the proposed variants showing a variant for strongly convex quadratic problem enjoys the same guarantee as full-memory AM, and another variant converges with a similar convergence rate as SGD (O(1/\sqrt{number of iterations})) on nonconvex problems. Experiments on MNIST, CIFAR-10, and PENN TREEBANK with some popular deep models validates the superior performance of one of the proposed method.