Year: 2021

Abstract:
Advances in generative modeling and adversarial learning have given rise to renewed interest in smooth games. However, the absence of symmetry in the matrix of second derivatives poses challenges that are not present in the classical minimization framework. While a rich theory of average-case analysis has been developed for minimization problems, little is known in the context of smooth games. In this work we take a first step towards closing this gap by developing average-case optimal first-order methods for a subset of smooth games. We make the following three main contributions. First, we show that for zero-sum bilinear games the average-case optimal method is the optimal method for the minimization of the Hamiltonian. Second, we provide an explicit expression for the optimal method corresponding to normal matrices, potentially non-symmetric. Finally, we specialize it to matrices with eigenvalues located in a disk and show a provable speed-up compared to worst-case optimal algorithms. We illustrate our findings through numerical simulations with a varying degree of mismatch with our assumptions.

heading: 1 INTRODUCTION
text: The traditional analysis of optimization algorithms is a worst-case analysis (Nemirovski, 1995; Nesterov, 2004). This type of analysis provides a complexity bound for any input from a function class, no matter how unlikely. However, since hard-to-solve inputs might rarely occur in practice, the worst-case complexity bounds might not be representative of the observed running time.
A more representative analysis is given by the average-case complexity, averaging the algorithm’s complexity over all possible inputs. This analysis is standard for analyzing, e.g., sorting (Knuth, 1997) and cryptography algorithms (Katz & Lindell, 2014). Recently, a line of work (Berthier et al., 2020; Pedregosa & Scieur, 2020; Lacotte & Pilanci, 2020; Paquette et al., 2020) focused on optimal methods for the optimization of quadratics, specified by a symmetric matrix. While worst-case analysis uses bounds on the matrix eigenvalues to yield upper and lower bounds on convergence, average-case analysis relies on the expected distribution of eigenvalues and provides algorithms with sharp optimal convergence rates. While the algorithms developed in this context have been shown to be efficient for minimization problems, these have not been extended to smooth games.
A different line of work considers algorithms for smooth games but studies worst-case optimal methods (Azizian et al., 2020). In this work, we combine average-case analysis with smooth games, and develop novel average-case optimal algorithms for finding the root of a linear system determined by a (potentially non-symmetric) normal matrix. We make the following main contributions:
1. Inspired by the problem of finding equilibria in smooth games, we develop average-case optimal algorithms for finding the root of a non-symmetric affine operator, both under a normality assumption (Thm. 4.1), and under the extra assumption that eigenvalues of the operator are supported in a disk (Thm. 4.2). The proposed method shows a polynomial speedup compared to the worst-case optimal method, verified by numerical simulations.
2. We make a novel connection between average-case optimal methods for optimization, and average-case optimal methods for bilinear games. In particular, we show that solving the Hamiltonian using an average-case optimal method is optimal (Theorem 3.1) for bilinear games. This result complements (Azizian et al., 2020), who proved that Polyak Heavy Ball algorithm on the Hamiltonian is asymptotically worst-case optimal for bilinear games.

heading: 2 AVERAGE-CASE ANALYSIS FOR NORMAL MATRICES
text: In this paper we consider the following class of problems. Definition 1. Let A ∈ Rd×d be a real matrix and x? ∈ Rd a vector. The non-symmetric (affine) operator (NSO) problem is defined as:
Find x : F (x) def = A(x−x?) = 0 . (NSO)
This problem generalizes that of minimization of a convex quadratic function f , since we can cast the latter in this framework by setting the operator F = ∇f . The set of solutions is an affine subspace that we will denote X ?. We will find convenient to consider the distance to this set, defined as
dist(x, X ?) def= min v∈X? ‖x− v‖2, with X ? = {x ∈ Rd |A(x− x?) = 0} . (1)
In this paper we will develop average-case optimal methods. For this, we consider A and x? to be random vectors, and a random initialization x0. This induces a probability distribution over NSO problems, and we seek to find methods that have an optimal expected suboptimality w.r.t. this distribution. Denoting E(A,x?,x0) the expectation over these random problems, we have that average-case optimal methods they verify the following property at each iteration t
min xt
E(A,x?,x0) dist(xt, X ?) s.t. xi ∈ x0 + span({F (xj)}i−1j=0), ∀i ∈ [1 : t]. (2)
The last condition on xt stems from restricting the class of algorithms to first-order methods. The class of first-order methods encompasses many known schemes such as gradient descent with momentum, or full-matrix AdaGrad. However, methods such as Adam (Kingma & Ba, 2015) or diagonal AdaGrad (Duchi et al., 2011) are not in this class, as the diagonal re-scaling creates iterates xt outside the span of previous gradients. Although we will focus on the distance to the solution, the results can be extended to other convergence criteria such as ‖F (xt)‖2. Finally, note that the expectations in this paper are on the problem instance and not on the randomness of the algorithm.

heading: 2.1 ORTHOGONAL RESIDUAL POLYNOMIALS AND FIRST-ORDER METHODS
text: The analysis of first-order methods simplifies through the use of polynomials. This section provides the tools required to leverage this connection. Definition 2. A residual polynomial is a polynomial P that satisfies P (0) = 1.
Proposition 2.1. (Hestenes et al., 1952) If the sequence (xt)t∈Z+ is generated by a first-order method, then there exist residual polynomials Pt, each one of degree at most t, verifying
xt − x? = Pt(A)(x0 − x?). (3)
As we will see, optimal average-case method are strongly related to orthogonal polynomials. We first define the inner product between polynomials, where we use z∗ for the complex conjugate of z ∈ C.
Definition 3. For P,Q ∈ R[X], we define the inner product 〈·, ·〉µ for a measure µ over C as
〈P,Q〉µ def = ∫ C P (λ)Q(λ)∗ dµ(λ) . (4)
Definition 4. A sequence of polynomials {Pi} is orthogonal (resp. orthonormal) w.r.t. 〈·, ·〉µ if
〈Pi, Pi〉µ > 0 (resp. = 1); 〈Pi, Pj〉µ = 0 if i 6= j.

heading: 2.2 EXPECTED SPECTRAL DISTRIBUTION
text: Following (Pedregosa & Scieur, 2020), we make the following assumption on the problem family.
Assumption 1. x0 − x? is independent ofA, and E(x0,x?)[(x0 − x?)(x0 − x?)>] = R 2 d Id.
We will also require the following definitions to characterize difficulty of a problem class. Let {λ1, . . . , λd} be the eigenvalues of a matrix A ∈ Rd×d. We define the empirical spectral distribution ofA as the probability measure
µ̂A(λ) def = 1d ∑d i=1δλi(λ) , (5)
where δλi is the Dirac delta, a distribution equal to zero everywhere except at λi and whose integral over the entire real line is equal to one. Note that with this definition, ∫ D dµ̂A(λ) corresponds to the proportion of eigenvalues in D. WhenA is a matrix-valued random variable, µA is a measure-valued random variable. As such, we can define its expected spectral distribution
µA def = EA[µ̂A] , (6) which by the Riesz representation theorem is the measure that verifies ∫ f dµ = EA[ ∫ f dµA] for all measureable f . Surprisingly, the expected spectral distribution is the only required characteristic to design optimal algorithms in the average-case.

heading: 2.3 EXPECTED ERROR OF FIRST-ORDER METHODS
text: In this section we provide an expression for the expected convergence in terms of the residual polynomial and the expected spectral distribution introduced in the previous section. To go further in the analysis, we have to assume thatA is a normal matrix.
Assumption 2. The (real) random matrixA is normal, that is, it verifiesAA> = A>A.
Normality is equivalent to A having the spectral decomposition A = UΛU∗, where U is unitary, i.e., U∗U = UU∗ = I. We now have everything to write the expected error of a first-order algorithm applied to (NSO).
Theorem 2.1. Consider the application of a first-order method associated to the sequence of polynomials {Pt} (Proposition 2.1) on the problem (NSO). Let µ be the expected spectral distribution ofA. Under Assumptions 1 and 2, we have
E[dist(xt,X ?)] = R2 ∫ C\{0} |Pt|2 dµ , (7)
Before designing optimal algorithms for certain specific distributions, we compare our setting with the average-case accelerating for minimization problems of Pedregosa & Scieur (2020), who proposed optimal optimization algorithms in the average-case.

heading: 2.4 DIFFICULTIES OF FIRST-ORDER METHODS ON GAMES AND RELATED WORK
text: This section compares our contribution with the existing framework of average-case optimal methods for quadratic minimization problems.
Definition 5. Let H ∈ Rd×d be a random symmetric positive-definite matrix and x? ∈ Rd a random vector. These elements determine the following random quadratic minimization problem
minx∈Rd { f(x) def = 1
2 (x−x?)>H(x−x?)
} . (OPT)
As in our paper, Pedregosa & Scieur (2020) find deterministic optimal first-order algorithms in expectation w.r.t. the matrix H , the solution x?, and the initialization x0. Since they work with problem (OPT), their problem is equivalent to (NSO) with the matrix A = H . However, they have the stronger assumption that the matrix is symmetric, which implies being normal. The normality assumption is restrictive in the case of game theory, as they do not always naturally fit such applications. However, this set is expressive enough to consider interesting cases, such as bilinear games, and our experiments show that our findings are also consistent with non-normal matrices.
Using orthogonal residual polynomials and spectral distributions, they derive the explicit formula of the expected error. Their result is similar to Theorem 2.1, but the major difference is the domain of the integral, a real positive line in convex optimization, but a shape in the complex plane in our case. This shape plays a crucial role in the rate of converge of first-order algorithms, as depicted in the work of Azizian et al. (2020); Bollapragada et al. (2018).
In the case of optimization methods, they show that optimal schemes in the average-case follow a simple three-term recurrence arising from the three-term recurrence for residual orthogonal polynomials for the measure λµ(λ). Indeed, by Theorem 2.1 the optimal method corresponds to the residual polynomials minimizing 〈P, P 〉µ, and the following result holds: Theorem 2.2. (Fischer, 1996, §2.4) When µ is supported in the real line, the residual polynomial of degree t minimizing 〈P, P 〉µ is given by the degree t residual orthogonal polynomial w.r.t. λµ(λ).
However, the analogous result does not hold for general measures in C, and hence our arguments will make use of the following Theorem 2.3 instead, which links the residual polynomial of degree at most t that minimizes 〈P, P 〉µ to the sequence of orthonormal polynomials for µ. Theorem 2.3. [Theorem 1.4 of Assche (1997)] Let µ be a positive Borel measure in the complex plane. The minimum of the integral ∫ C |P (λ)|
2 dµ(λ) over residual polynomials P of degree lower or equal than t is uniquely attained by the polynomial
P ?(λ) =
∑t k=0 φk(λ)φk(0)
∗∑t k=0 |φk(0)|2 , with optimal value ∫ C |P ?(λ)|2 dµ(λ) = 1∑t k=0 |φk(0)|2 , (8)
where (φk)k is the orthonormal sequence of polynomials with respect to the inner product 〈·, ·〉µ.
In the next sections we consider cases where the optimal scheme is identifiable.

heading: 3 AVERAGE-CASE OPTIMAL METHODS FOR BILINEAR GAMES
text: We consider the problem of finding a Nash equilibrium of the zero-sum minimax game given by
min θ1 max θ2
`(θ1, θ2) def = (θ1 − θ?1)>M(θ2 − θ?2) . (9)
Let θ1,θ?1 ∈ Rd1 ,θ2,θ?2 ∈ Rd2 ,M ∈ Rd1×d2 and d def = d1 + d2. The vector field of the game (Balduzzi et al., 2018) is defined as F (x) = A(x− x?), where
F (θ1, θ2) = [ ∇θ1`(θ1, θ2) −∇θ2`(θ1, θ2) ] = [ 0 M −M> 0 ] ︸ ︷︷ ︸
=A
([ θ1 θ2 ] ︸︷︷︸
=x
− [ θ?1 θ?2 ] ︸ ︷︷ ︸ =x? ) = A(x− x?) . (10)
As before, X ? denotes the set of points x such that F (x) = 0, which is equivalent to the set of Nash equilibrium. If M is sampled independently from x0,x? and x0 − x? has covariance R 2
d Id, Assumption 1 is fulfilled. Since A is skew-symmetric, it is in particular normal and Assumption 2 is also satisfied.
We now show that the optimal average-case algorithm to solve bilinear problems is Hamiltonian gradient descent with momentum, described below in its general form. Contrary to the methods in Azizian et al. (2020), the method we propose is anytime (and not only asymptotically) average-case optimal.
Optimal average-case algorithm for bilinear games.
Initialization. x−1 = x0 = ( θ1,0, θ2,0 ) , sequence {ht,mt} given by Theorem 3.1. Main loop. For t ≥ 0,
gt = F (xt − F (xt))− F (xt) ( = 12∇‖F (xt)‖ 2 by (12) )
xt+1 = xt − ht+1gt +mt+1(xt−1 − xt) (11)
The quantity 12‖F (x)‖ 2 is commonly known as the Hamiltonian of the game (Balduzzi et al., 2018), hence the name Hamiltonian gradient descent. Indeed, gt = ∇ ( 1 2‖F (x)‖ 2 ) when F is affine:
F (x− F (x))− F (x) = A(x−A(x− x?)− x?)−A(x− x?) = −A(A(x− x?)) = A>(A(x− x?)) = ∇ ( 1
2 ‖A(x− x?)‖2
) = ∇ ( 1
2 ‖F (x)‖2
) . (12)
The following theorem shows that (11) is indeeed the optimal average-case method associated to the minimization problem minx ( 1 2‖F (x)‖ 2 ) , as the following theorem shows.
Theorem 3.1. Suppose that Assumption 1 holds and that the expected spectral distribution of MM> is absolutely continuous with respect to the Lebesgue measure. Then, the method (11) is average-case optimal for bilinear games when ht, mt are chosen to be the coefficients of the average-case optimal minimization of 12‖F (x)‖ 2.
How to find optimal coefficients? Since 12‖F (x)‖ 2 is a quadratic problem, the coefficients {ht, mt} can be found using the average-case framework for quadratic minimization problems of (Pedregosa & Scieur, 2020, Theorem 3.1).
Proof sketch. When computing the optimal polynomial xt = Pt(A)(x0 − x?), we have that the residual orthogonal polynomial Pt behaves differently if t is even or odd.
• Case 1: t is even. In this case, we observe that the polynomial Pt(A) can be expressed as Qt/2(−A2), where (Qt)t≥0 is the sequence of orthogonal polynomials w.r.t. the expected spectral density of −A2, whose eigenvalues are real and positive. This gives the recursion in (11).
• Case 2: t is odd. There is no residual orthogonal polynomial of degree t for t odd. Instead, odd iterations do correspond to the intermediate computation of gt in (11), but not to an actual iterate.

heading: 3.1 PARTICULAR CASE: M WITH I.I.D. COMPONENTS
text: We now show the optimal method when the entries of M are i.i.d. sampled. For simplicity, we order the players such that d1 ≤ d2. Assumption 3. Assume that each component ofM is sampled iid from a distribution of mean 0 and variance σ2, and we take d1, d2 →∞ with d1d2 → r < 1.
In such case, the spectral distribution of 1d2MM > tends to the Marchenko-Pastur law, supported in [`, L] and with density:
ρMP (λ) def =
√ (L− λ)(λ− `)
2πσ2rλ , where L def= σ2(1 +
√ r)2, ` def = σ2(1− √ r)2. (13)
Proposition 3.1. WhenM satisfies Assumption 3, the optimal parameter of scheme (11) are
ht = − δtσ2√r , mt = 1 + ρδt, where ρ = 1+r√ r , δt = (−ρ− δt−1)−1, δ0 = 0. (14)
Proof. By Theorem 3.1, the problem reduces to finding the optimal average-case algorithm for the problem minx 12‖F (x)‖ 2. Since the expected spectral distribution of 1d2MM > is the MarchenkoPastur law, we can use the optimal algorithm from (Pedregosa & Scieur, 2020, Section 5).

heading: 4 GENERAL AVERAGE-CASE OPTIMAL METHOD FOR NORMAL OPERATORS
text: In this section we derive general average-case optimal first-order methods for normal operators. First, we need to assume the existence of a three-term recurrence for residual orthogonal polynomials (Assumption 4). As mentioned in subsection 2.4, for general measures in the complex plane, the existence of a three-term recurrence of orthogonal polynomials is not ensured. In Proposition B.3 in Appendix B we give a sufficient condition for its existence, and in the next subsection we will show specific examples where the residual orthogonal polynomials satisfy the three-term recurrence. Assumption 4 (Simplifying assumption). The sequence of residual polynomials {ψt}t≥0 orthogonal w.r.t. the measure µ, defined on the complex plane, admits the three-term recurrence
ψ−1 = 0, ψ0 = 1, ψt(λ) = (at + btλ)ψt−1(λ) + (1− at)ψt−2(λ). (15)
Under Assumption 4, Theorem 4.1 shows that the optimal algorithm can also be written as an average of iterates following a simple three-terms recurrence.
Theorem 4.1. Under Assumption 4 and the assumptions of Theorem 2.1, the following algorithm is optimal in the average case, with y−1 = y0 = x0:
yt = atyt−1 + (1− at)yt−2 + btF (yt−1)
xt = Bt
Bt + βt xt−1 + βt Bt + βt yt , βt = φ 2 t (0), Bt = Bt−1 + βt−1, B0 = 0 . (16)
where (φk(0))k≥0 can be computed using the three-term recurrence (upon normalization). Moreover, E (A,x?,x0) dist(xt,X ?) converges to zero at rate 1/Bt.
Remark. Notice that it is not immediate that (16) fulfills the definition of first-order algorithms stated in (2), as yt is clearly a first-order method but xt is an average of the iterates yt. Using that F is an affine function we see that xt indeed fulfills (2).
Remark. Assumption 4 is needed for the sequence (yt)t≥0 to be computable using a three-term recurrence. However, for some distribution, the associated sequence of orthogonal polynomials may admit another recurrence that may not satisfy Assumption 4.

heading: 4.1 CIRCULAR SPECTRAL DISTRIBUTIONS
text: In random matrix theory, the circular law states that ifA is an n×nmatrix with i.i.d. entries of mean C and variance R2/n, as n → ∞ the spectral distribution of A tends to the uniform distribution on DC,R. In this subsection we apply Theorem 4.1 to a class of spectral distributions specified by Assumption 5, which includes the uniform distribution on DC,R. Even though the random matrices with i.i.d entries are not normal, we will see in section 6 that the empirical results for such matrices are consistent with our theoretical results under the normality assumption. Assumption 5. Assume that the expected spectral distribution µA is supported in the complex plane on the disk DC,R of center C ∈ R, C > 0 and radius R < C. Moreover, assume that the spectral density is circularly symmetric, i.e. there exists a probability measure µR supported on [0, R] such for all f measurable and r ∈ [0, R], dµA(C + reiθ) = 12π dθ dµR(r).
Proposition 4.1. If µ satisfies Assumption 5, the sequence of orthonormal polynomials is (φt)t≥0,
φt(λ) = (λ− C)t
Kt,R , where Kt,R =
√∫ R 0 r2t dµR(r) . (17)
Example. The uniform distribution in DC,R is to dµR = 2rR2 dr, and Kt,R = R t/ √ t+ 1.
From Proposition 4.1, the sequence of residual polynomials is given by φt(λ)/φt(0) = ( 1− λC )t , which implies that Assumption 4 is fulfilled with at = 1, bt = − 1C . Thus, by Theorem 4.1 we have
Theorem 4.2. Given an initialization x0(y0 = x0), if Assumption 5 is fulfilled with R < C and the assumptions of Theorem 2.1 hold, then the average-case optimal first-order method is
yt = yt−1 − 1CF (yt−1), βt = C 2t/K2t,R, Bt = Bt−1 + βt−1,
xt = Bt
Bt + βt xt−1 + βt Bt + βt yt. (18)
Moreover, E (A,x?,x0) dist(xt,X ?) converges to zero at rate 1/Bt.
We now compare Theorem 4.2 with worst-case methods studied in Azizian et al. (2020). They give a worst-case convergence lower bound of (R/C)2t on the quantity dist(zt,X ?) for first-order methods (zt)t≥0 on matrices with eigenvalues in the disk DC,R. By the classical analysis of first-order methods, this rate is achievable by gradient descent with stepsize 1/C, i.e. the iterates yt defined in (18). However, by equation (79) in Proposition D.3 we have that under slight additional assumptions (those of Proposition 5.2), limt→∞ E [dist(xt,X ?)]/E [dist(yt,X ?)] = 1− R 2
C2 holds. That is, the average-case optimal algorithm outperforms gradient descent by a constant factor depending on the conditioning R/C.

heading: 5 ASYMPTOTIC BEHAVIOR
text: The recurrence coefficients of the average-case optimal method typically converges to limiting values when t → ∞, which gives an ”average-case asymptotically optimal first-order method” with constant coefficients. For the case of symmetric operators with spectrum in [`, L], Scieur & Pedregosa (2020) show that under mild conditions, the asymptotically optimal algorithm is the Polyak momentum method with coefficients depending only on ` and L. For bilinear games, since the average-case optimal algorithm is the average-case optimal algorithm of an optimization algorithm, we can make use of their framework to obtain the asymptotic algorithm (see Theorem 3 of Scieur & Pedregosa (2020)). Proposition 5.1. Assume that the expected spectral density µMM> ofMM> is supported in [`, L] for 0 < ` < L, and strictly positive in this interval. Then, the asymptotically optimal algorithm for bilinear games is the following version of Polyak momentum:
gt = F (xt − F (xt))− F (xt) xt+1 = xt + (√ L− √ `√
L+ √ `
)2 (xt−1 − xt)− ( 2√
L+ √ `
)2 gt
(19)
Notice that the algorithm in (19) is the worst-case optimal algorithm from Proposition 4 of Azizian et al. (2020). For the case of circularly symmetric spectral densities with support on disks, we can also compute the asymptotically optimal algorithm. Proposition 5.2. Suppose that the assumptions of Theorem 4.2 hold with µR ∈ P([0, R]) fulfilling µR([r,R]) = Ω((R − r)κ) for r in [r0, R] for some r0 ∈ [0, R) and for some κ ∈ Z. Then, the average-case asymptotically optimal algorithm is, with y0 = x0:
yt = yt−1 − 1CF (yt−1), xt = ( R C )2 xt−1 + ( 1− ( R C )2) yt.
(20)
Moreover, the convergence rate for this algorithm is asymptotically the same one as for the optimal algorithm in Theorem 4.2. Namely, limt→∞ E [dist(xt,X ?)]Bt = 1.
The condition on µR simply rules out cases in which the spectral density has exponentially small mass around 1. It is remarkable that in algorithm (20) the averaging coefficients can be expressed so simply in terms of the quantity R/C. Notice also that while the convergence rate of the algorithm
is slower than the convergence rate for the optimal algorithm by definition, both rates match in the limit, meaning that the asymptotically optimal algorithm also outperforms gradient descent by a constant factor 1− R 2
C2 in the limit t→∞.

heading: 6 EXPERIMENTS
text: We compare some of the proposed methods on settings with varying degrees of mismatch with our assumptions.
Bilinear Games. We consider min-max bilinear problems of the form (10), where the entries of M are generated i.i.d. from a standard Gaussian distribution. We vary the ratio r = d/n parameter for d = 1000 and compare the average-case optimal method of Theorems 3.1 and 5.1, the asymptotic worst-case optimal method of (Azizian et al., 2020) and extragradient (Korpelevich, 1976). In all cases, we use the convergence-rate optimal step-size assuming knowledge of the edges of the spectral distribution.
The spectral density for these problems is displayed in the first row of Figure 1 and the benchmark results on the second row. Average-case optimal methods always outperform other methods, and the largest gain is in the ill-conditioned regime (r ≈ 1).
Circular Distribution. For our second experiment we choose A as a matrix with iid Gaussian random entries, therefore the support of the distribution of its eigenvalue is a disk. Note that A does not satisfy the normality assumption of Assumption 2. Figure 1 (third row) compares the average-case optimal methods from Theorems 4.2 and 5.2 on two datasets with different levels of conditioning. Note that the methods converge despite the violation of Assumption 2, suggesting a broader applicability than the one proven in this paper. We leave this investigation for future work.

heading: 7 DISCUSSION AND FUTURE RESEARCH DIRECTIONS
text: In this paper, we presented a general framework for the design of optimal algorithms in the averagecase for affine operators F , whose underlying matrix is possibly non-symmetric. However, our approach presents some limitations, the major one being the restriction to normal matrices. Fortunately, the numerical experiments above suggests that this assumption can be relaxed. Developing a theory without that assumption is left for future work. Another avenue for future work is to analyze the nonlinear-case in which the non-symmetric operatorA is non-linear, as well as the case in which it is accessed through a stochastic estimator (as done by (Loizou et al., 2020) for the worst-case analysis).

heading: ACKNOWLEDGEMENTS
text: C. Domingo-Enrich has been partially funded by “la Caixa” Foundation (ID 100010434), under agreement LCF/BQ/AA18/11680094, and partially funded by the NYU Computer Science Department.
Bilinear Problems

heading: A PROOF OF THEOREM 2.1
text: A.1 PRELIMINARIES
Before proving Theorem 2.1, we quickly analyze the distance function (1), recalled below,
dist(x, X ?) def= min v∈X? ‖x− v‖2.
The definition of the distance function is not practical for the theoretical analysis. Fortunately, it is possible to find a simple expression that uses the orthogonal projection matrix Π to the kernel Ker(A). Since Π is an orthogonal projection matrix to the kernel of a linear transformation, it satisfies Π = ΠT , Π2 = Π, and AΠ = 0. (21) The normality assumption onA implies also that
ΠA = 0. (22)
Indeed, the spectral decomposition ofA is A = [U1|U2] [ Λ 0 0 0 ] [U1|U2]∗,
and then Π = U2U∗2 . The next proposition uses Π to derive the explicit solution of the (1). Proposition A.1. We have that
dist(y, X ?) = ‖(I −Π)(y − x?)‖2 ∀x? ∈ X ?.
Proof. We first parametrize the set of solution X ?. By definition we have
X ? = {x : A(x− x?) = 0}.
Which can be written in terms of the kernel ofA as
X ? = {x? + Πw : w ∈ Rd}.
From this, we can rewrite the distance function (1) as
dist(y, X ?) = min w∈Rd ‖y − (x? + Πw)‖2.
The minimum can be attained at different points, but in particular atw = −(y−x?), which proves the statement.
We now simplifies further the result of the previous proposition in the case where xt is generated by a first order method. Proposition A.2. For every iterate xt of a first-order methods, i.e., xt satisfies
xt − x? = Pt(A)(x0 − x?), deg(Pt) ≤ t, P (0) = I,
we have that
dist(xt, X ?) = ‖xt − x?‖2 − ‖Π(x0 − x?)‖2.
Proof. We start with the result of Proposition A.1,
dist(xt, X ?) = ‖(I −Π)(xt − x?)‖2.
The norm can be split into
‖(I −Π)(xt − x?)‖2 = ‖xt − x?‖2 + ‖ Π2︸︷︷︸ =Π by (21) (xt − x?)‖2 − 2‖Π(xt − x?)‖2
= ‖xt − x?‖2 − ‖Π(xt − x?)‖2.
Since xt is generated by a first order method, we have
xt − x? = Pt(A)(x0 − x?), Pt(0) = 1.
Since P (0) = 1, the polynomial can be factorized as P (A) = I + AQt−1(A), Qt−1 being a polynomial of degree t− 1. Therefore, ‖Π(xt − x?)‖2 reads
‖Π(xt − x?)‖2 = ‖Π (I +AQt−1(A)) (x0 − x?)‖2
= ‖Π(x0 − x?) + ΠA︸︷︷︸ =0 by (22) Qt−1(A)(x0 − x?)‖2
= ‖Π(x0 − x?)‖2 ,
which prove the statement.
A.2 PROOF OF THE THEOREM
We are now ready to prove the main result.
Theorem 2.1. Consider the application of a first-order method associated to the sequence of polynomials {Pt} (Proposition 2.1) on the problem (NSO). Let µ be the expected spectral distribution of A. Under Assumptions 1 and 2, we have
E[dist(xt,X ?)] = R2 ∫ C\{0} |Pt|2 dµ , (7)
Proof. We start with the result of Proposition A.2,
dist(xt, X ?) = ‖xt − x?‖2 − ‖Π(x0 − x?)‖2.
We now write the expectation of the distance function, E[dist(xt, X ?)] = E [ ‖xt − x?‖2 − ‖Π(x0 − x?)‖2 ] = E [ ‖Pt(A)(x0 − x?)‖2 − ‖Π(x0 − x?)‖2
] = E [ trPt(A)Pt(A) T (x0 − x?) (x0 − x?)T − tr Π2(x0 − x?)(x0 − x?)T ]
= EA [ trPt(A)Pt(A) TE [ (x0 − x?) (x0 − x?)T |A ] − tr ΠE [ (x0 − x?)(x0 − x?)T |A ]] = REA [ trPt(A)Pt(A) T − tr Π ]
= RE [ d∑ i=1 |P (λi)|2 − tr Π ]
= RE [∫ C\{0} |P (λ)|2δλi(λ) + |P (0)|2 · [# zero eigenvalues]− tr Π ]
However, |P (0)|2 = 1 and tr Π corresponds to the number of zero eigenvalues ofA, therefore,
E[dist(xt, X ?)] = RE [∫ C\{0} |P (λ)|2δλi(λ) ] = R ∫ C\{0} P (λ)µ(λ).

heading: B PROOFS OF THEOREM 3.1 AND PROPOSITION 3.1
text: Proposition B.1. [Block determinant formula] If A,B,C,D are (not necessarily square) matrices, det [ A B C D ] = det(D)det(A−BD−1C), (23)
if D is invertible.
Definition 6 (Pushforward of a measure). Recall that the pushforward f∗µ of a measure µ by a function f is defined as the measure such that for all measurable g,∫
g(λ) d(f∗µ)(λ) =
∫ g(f(λ)) dµ(λ). (24)
Equivalently, if X is a random variable with distribution µ, then f(X) has distribution f∗µ.
Proposition B.2. Assume that the dimensions of M ∈ Rdx×dy fulfill dx ≤ dy and let r = dx/dy . Let µMM> be the expected spectral distribution of the random matrix MM> ∈ Rdx×dx , and assume that it is absolutely continuous with respect to the Lebesgue measure. The expected spectral distribution ofA is contained in the imaginary line and is given by
µA(iλ) =
( 1− 2
1 + 1r
) δ0(λ) +
2|λ| 1 + 1r µMM>(λ 2) . (25)
for λ ∈ R. If dx ≥ dy , then (25) holds with µM>M in place of µMM> and 1/r in place of r.
Proof. By the block determinant formula, we have that for s 6= 0, det (sId1+d2 −A) = ∣∣∣∣sId1 −MM> sId2 ∣∣∣∣ = det(sId2)det(sId1 +Ms−1Id2M>) = sd2−d1det(s2Id1 +MM >)
(26)
Thus, for every eigenvalue −λ ≤ 0 of −MM>, both i √ λ and −i √ λ are eigenvalues of A. Since rank(MM>) = rank(M), we have rank(A) = 2rank(M). Thus, the rest of the eigenvalues of A are 0 and there is a total of d− 2d1 = d2 − d1 of them. Notice that
d1 d1 + d2 = 1
d1+d2 d1
= 1
1 + 1r (27)
Let f+(λ) = i √ λ, f−(λ) = −i √ λ, and let (f+)∗µMM> (resp., (f−)∗µMM> ) be the pushforward measure of µMM> by the function f+ (resp., f−). Thus, by the definition of the pushforward measure (Definition 6),
µA(iλ) =
( 1− 2
1 + 1r
) δ0(λ) + 1
1 + 1r (f+)∗µMM>(λ) +
1
1 + 1r (f−)∗µMM>(λ) (28)
We compute the pushforwards (f+)∗µMM> , (f−)∗µMM> performing the change of variables y = ±i √ λ under the assumption that µMM>(λ) = ρMM>(λ)dλ:∫
R≥0 g ( ±i √ λ ) dµMM>(λ) = ∫ R≥0 g ( ±i √ λ ) ρMM>(λ)dλ = ∫ ±iR≥0 g (y) ρMM>(|y|2)2|y|d|y|,
(29)
which means that the density of (f+)∗µMM> at y ∈ iR≥0 is 2|y|ρMM>(|y|2) and the density of (f−)∗µMM> at y ∈ −iR≥0 is also 2|y|ρMM>(|y|2).
Proposition B.3. The condition
∀P,Q polynomials 〈P (λ), λQ(λ)〉 = 0 =⇒ 〈λP (λ), Q(λ)〉 = 0 (30)
is sufficient for any sequence (Pk)k≥0 of orthogonal polynomials of increasing degrees to satisfy a three-term recurrence of the form
γkPk(λ) = (λ− αk)Pk−1(λ)− βkPk−2(λ), (31)
where
γk = 〈λPk−1(λ), Pk(λ)〉 〈Pk(λ), Pk(λ)〉 , αk = 〈λPk−1(λ), Pk−1(λ)〉 〈Pk−1(λ), Pk−1(λ)〉 , βk = 〈λPk−1(λ), Pk−2(λ)〉 〈Pk−2(λ), Pk−2(λ)〉 (32)
Proof. Since λPk−1(λ) is a polynomial of degree k, and (Pj)0≤j≤k is a basis of the polynomials of degree up to k, we can write
λPk−1(λ) = k∑ j=0 〈λPk−1, Pj〉 〈Pj , Pj〉 Pj(λ) (33)
Now, remark that for all j < k − 2, 〈Pk−1, λPj〉 = 0 because the inner product of Pk−1 with a polynomial of degree at most k − 2. If we make use of the condition (30), this implies that 〈λPk−1, Pj〉 = 0 for all j < k − 2. Plugging this into (33), we obtain (31).
Proposition B.4. Let ΠRt be the set of polynomials with real coefficients and degree at most t. For t ≥ 0 even, the minimum of the problem
min Pt∈ΠRt ,Pt(0)=1 ∫ iR\{0} |Pt(λ)|2|λ|ρMM>(|λ|2) d|λ| (34)
is attained by an even polynomial with real coefficients.
Proof. Since dµ(iλ) def= |λ|ρMM>(|λ|2) d|λ| is supported in the imaginary axis and is symmetric with respect to 0, for all polynomials P,Q,
〈λP (λ), Q(λ)〉 = ∫ iR λP (λ)Q(λ)∗dµ(λ) = − ∫ iR P (λ)λ∗Q(λ)∗dµ(λ) = −〈P (λ), λQ(λ)〉.
(35)
Hence, 〈P (λ), λQ(λ)〉 = 0 implies 〈λP (λ), Q(λ)〉 = 0. By Proposition B.3, a three-term recurrence (31) and (32) for the orthonormal sequence (φt)t≥0 of polynomials holds.
By Proposition B.5, the orthonormal polynomials (φt)t≥0 of even (resp. odd) degree are even (resp. odd) and have real coefficients. Hence, for all t ≥ 0 even∑t
k=0 φk(λ)φk(0) ∗∑t
k=0 |φk(0)|2 =
∑t/2 k=0 φ2k(λ)φ2k(0)
∗∑t/2 k=0 |φ2k(0)|2
(36)
is an even polynomial with real coefficients. By Theorem 2.3, this polynomial attains the minimum of the problem
min Pt∈ΠCt ,Pt(0)=1 ∫ iR\{0} |Pt(λ)|2|λ|ρMM>(|λ|2) d|λ| (37)
and, a fortiori, the minimum of the problem in (34), in which the minimization is restricted polynomials with real coefficients instead of complex coefficients.
Proposition B.5. The polynomials (φt)t≥0 of the orthonormal sequence corresponding to the measure µ(iλ) = |λ|ρMM>(|λ|2)d|λ| have real coefficients and are even (resp. odd) for even (resp. odd) k.
Proof. The proof is by induction. The base case follows from the choice φ0 = 1. Assuming that φk−1 ∈ R[X] by the induction hypothesis, we show that αk = 0 (where αk is the coefficient from (31) and (32)):
〈λφk−1(λ), φk−1(λ)〉 = ∫ iR λ|φk−1(λ)|2|λ|ρMM>(|λ|2)d|λ|
= ∫ R≥0 iλ(|φk−1(iλ)|2 − |φk−1(−iλ)|2)λρMM>(λ2)dλ = 0 (38)
The last equality follows from |φk−1(iλ)|2 = |φk−1(−iλ)|2, which holds because φk−1(iλ)∗ = φk−1(−iλ), and in turn this is true because φk−1 ∈ R[X] by the induction hypothesis. Once we have seen that αk = 0, it is straightforward to apply the induction hypothesis once again to show that φk also satisfies the even/odd property. Namely, for k even (resp. odd), γkPk = λPk−1 − βkPk−2, and the two polynomials in the right-hand side have even (resp. odd) degrees.
Finally, φk must have real coefficients because φk−1 and φk−2 have real coefficients by the induction hypothesis, and the recurrence coefficient βk is real, as
〈λPk−1(λ), Pk−2(λ)〉 = ∫ iR λφk−1(λ)φk−2(λ) ∗|λ|ρMM>(|λ|2)d|λ|
= ∫ R≥0 iλ(φk−1(iλ)φk−2(iλ) ∗ − φk−1(iλ)∗φk−2(iλ))λρMM>(λ2)dλ
= − ∫ R≥0 2λIm(φk−1(iλ)φk−2(iλ)∗)λρMM>(λ 2)dλ ∈ R.
(39)
Proposition B.6. Let t ≥ 0 even. Assume that on R>0, the expected spectral density µMM> has Radon-Nikodym derivative ρMM> with respect to the Lebesgue measure. If
Q?t/2 def = arg min
Pt/2∈ΠRt/2, Pt/2(0)=1
∫ R>0 Pt/2(λ) 2 dµ−A2(λ), (40)
and
P ?t def = arg min
Pt∈ΠRt , Pt(0)=1
∫ iR\{0} |Pt(λ)|2|λ|ρMM>(|λ|2) d|λ|, (41)
then P ?t (λ) = Q ? t/2(−λ 2).
Proof. First, remark that the equalities in (40) and (41) are well defined because the arg min are unique by Theorem 2.3. Without loss of generality, assume that dx ≤ dy (otherwise switch the players), and let r def= dx/dy < 1. Since,
−A2 = [ MM> 0
0 M>M
] , (42)
each eigenvalue of MM> ∈ Rdx×dx is an eigenvalue of −A2 with doubled duplicity, and the rest of eigenvalues are zero. Hence, we have µ−A2 = ( 1− 2/(1 + 1r ) ) δ0 + 2µMM>/(1 + 1 r ). Thus, for all t ≥ 0,
Q?t = arg min Pt∈ΠRt , Pt(0)=1
∫ R>0 Pt(λ) 2 dµ−A2(λ) = arg min
Pt∈ΠRt , Pt(0)=1
∫ R>0 Pt(λ) 2ρMM>(λ) dλ (43)
By Proposition B.4, for an even t ≥ 0 the minimum in (41) is attained by an even polynomial with real coefficients. Hence,
min Pt∈ΠRt , Pt(0)=1
∫ iR\{0}
|Pt(λ)|2|λ|ρMM>(|λ|2) d|λ| = min Pt/2∈ΠRt/2, Pt/2(0)=1
∫ iR\{0} |Pt/2(λ2)|2|λ|ρMM>(|λ|2) d|λ|
= 2 min Pt/2∈ΠRt/2, Pt/2(0)=1
∫ R>0 |Pt/2((iλ)2)|2λρMM>(λ2) dλ = 2 min
Pt/2∈ΠRt/2, Pt/2(0)=1
∫ R>0 Pt/2(λ 2)2λρMM>(λ 2) dλ
= min Pt/2∈ΠRt/2, Pt/2(0)=1
∫ R>0 Pt/2(λ) 2ρMM>(λ) dλ
(44)
Moreover, for any polynomialQt/2 that attains the minimum on the right-most term, the polynomial Pt(λ) = Qt/2(−λ2) attains the minimum on the left-most term. In particular, using (43), P ?t (λ) def = Q?t/2(−λ 2) attains the minimum on the left-most term.
Theorem 3.1. Suppose that Assumption 1 holds and that the expected spectral distribution of MM> is absolutely continuous with respect to the Lebesgue measure. Then, the method (11) is average-case optimal for bilinear games when ht, mt are chosen to be the coefficients of the average-case optimal minimization of 12‖F (x)‖ 2.
Proof. Making use of Theorem 2.1 and Proposition B.2, we obtain that for any first-order method using the vector field F ,
E[dist(xt,X ?)] = R2 ∫ C\{0} |Pt(λ)|2 dµA(λ) = 2R2 1 + 1r ∫ iR\{0} |Pt(λ)|2|λ|ρMM>(|λ|2) d|λ|
(45)
Let Q?t/2, P ? t be as defined in (41) and (40). For t ≥ 0 even the iteration t of the average-case optimal method for the bilinear game must satisfy
xt − PX?(x0) = P ?t (A)(x0 − PX?(x0)) = Q?t/2(−A 2)(x0 − PX?(x0)) (46)
On the other hand, the first-order methods for the minimization of the function 12‖F (x)‖ 2 make use of the vector field ∇ (
1 2‖F (x)‖
2 )
= A>(Ax + b) = −A2(x − x?). Let µ−A2 be the spectral density of −A2. By Theorem 2.1, the average-case optimal first-order method for the minimization problem is the one for which the residual polynomial Pt (Proposition 2.1) minimizes the functional∫ R P 2 t dµ−A2 . That is, the residual polynomial is Q ? t . From (46), we see that the t-th iterate of the average-case optimal method for F is equal to the t/2-th iterator of the average-case optimal method for ∇ ( 1 2‖F (x)‖ 2 ) .

heading: C PROOFS OF THEOREM 4.1 AND THEOREM 4.2
text: Theorem 4.1. Under Assumption 4 and the assumptions of Theorem 2.1, the following algorithm is optimal in the average case, with y−1 = y0 = x0:
yt = atyt−1 + (1− at)yt−2 + btF (yt−1)
xt = Bt
Bt + βt xt−1 + βt Bt + βt yt , βt = φ 2 t (0), Bt = Bt−1 + βt−1, B0 = 0 . (16)
where (φk(0))k≥0 can be computed using the three-term recurrence (upon normalization). Moreover, E (A,x?,x0) dist(xt,X ?) converges to zero at rate 1/Bt.
Proof. We prove by induction that xt − x? = ∑t k=0 φk(A)φk(0)
∗∑t k=0 φk(0) 2 (x0 − x?) (47)
The base step t = 0 holds trivially because φ0 = 1. Assume that (47) holds for t − 1. Subtracting x? from (16), we have
xt − x? = ∑t−1 k=0 φk(0)
2∑t k=0 φk(0) 2 (xt−1 − x?) + φt(0) 2∑t k=0 φk(0) 2 (yt − x?) (48)
If
φt(0) 2(yt − x?) = φt(0)φt(A)(x0 − x?), (49)
by the induction hypothesis for t− 1 and (48), we have xt − x? = ∑t−1 k=0 φt(0)φt(A)∑t k=0 φk(0) 2 (x0 − x?) + φt(0)φt(A)∑t k=0 φk(0) 2 (x0 − x∗)
= ∑t k=0 φt(0)φt(A)∑t k=0 φk(0) 2 (x0 − x∗),
(50)
which concludes the proof of (47). The only thing left is to show (49), again by induction. The base case follows readily from y0 = x0 in (16). Dividing by φt(0)2, we rewrite (49) as
yt − x? = φt(A)
φt(0) (x0 − x?) = ψt(A)(x0 − x?), (51)
where ψt is the t-th orthogonal residual polynomial of sequence. By Assumption 4, ψt must satisfy the recurrence in (15). If we subtract x∗ from the second line of (16), we apply the induction hypothesis and then the recurrence in (15), we obtain
yt − x? = at(yt−1 − x?) + (1− at)(yt−2 − x?) + btF (yt−1) = at(yt−1 − x?) + (1− at)(yt−2 − x?) + btA(yt−1 − x∗) = atψt−1(A)(x0 − x?) + (1− at)ψt−2(A)(x0 − x?) + btAψt−1(A)(x0 − x?) = ψt(A)(x0 − x?), (52)
thus concluding the proof of (49).
Proposition C.1. Suppose that Assumption 5 holds with C = 0, that is, the circular support of µ is centered at 0. Then, the basis of orthonormal polynomials for the scalar product
〈P,Q〉 = ∫ DR,0 P (λ)Q(λ)∗ dµ(λ) is φk(λ) = λk Dk,R , ∀k ≥ 0, (53)
where Kk,R = √ 2π ∫ R
0 r2kdµR(r).
Proof. First, we will show that if µ satisfies Assumption 5 with C = 0, then 〈λi, λj〉 = 0 if j, k ≥ 0 with j 6= k (without loss of generality, suppose that j > k).
〈λj , λk〉 = ∫ DR,0 λj(λ∗)k dµ(λ) = ∫ DR,0 λj−k|λ|2k dµ(λ)
= ∫ R 0 1 2π ∫ 2π 0 (reiθ)j−kr2k dθ dµR(r) = 1 2π ∫ 2π 0 eiθ(j−k) dθ ∫ R 0 rj+k dµR(r)
= ei2π − 1
2πi(j − k) ∫ R 0 rj+k dµR(r) = 0
(54)
And for all k ≥ 0, 〈λk, λk〉 = ∫ DR,0 |λk|2 dµ(λ) = ∫ R 0 1 2π ∫ 2π 0 r2k dθ dµR(r) = ∫ 2π 0 r2k dµR(r). (55)
Proposition 4.1. If µ satisfies Assumption 5, the sequence of orthonormal polynomials is (φt)t≥0,
φt(λ) = (λ− C)t
Kt,R , where Kt,R =
√∫ R 0 r2t dµR(r) . (17)
Proof. The result follows from Proposition C.1 using the change of variables z → z + C. To compute the measure µR for the uniform measure on DC,R, we perform a change of variables to circular coordinates:∫
DC,R
f(λ) dµ(λ) = 1
πR2 ∫ R 0 ∫ 2π 0 f(C + reiθ)r dθ dr = ∫ R 0 ∫ 2π 0 f(C + reiθ) dθ dµR(r).
=⇒ dµR(r) = r
πR2 dr
(56)
And ∫ R 0 r2t dµR(r) = 1 πR2 ∫ R 0 r2t+1 dr = 1 π R2t 2t+ 2 =⇒ Kt,R = Rt/ √ t+ 1. (57)
Theorem 4.2. Given an initialization x0(y0 = x0), if Assumption 5 is fulfilled with R < C and the assumptions of Theorem 2.1 hold, then the average-case optimal first-order method is
yt = yt−1 − 1CF (yt−1), βt = C 2t/K2t,R, Bt = Bt−1 + βt−1,
xt = Bt
Bt + βt xt−1 + βt Bt + βt yt. (18)
Moreover, E (A,x?,x0) dist(xt,X ?) converges to zero at rate 1/Bt.
Proof. By Proposition 4.1, the sequence of residual orthogonal polynomials is given by ψt(λ) = φt(λ)/φt(0) = ( 1− λC )t . Hence, Assumption 4 is fulfilled with at = 1, bt = − 1C , as ψt(λ) = ψt−1(λ) − λCψt−1(λ). We apply Theorem 4.1 and make use of the fact that φk(0) 2 = C 2k
K2t,R . See
Proposition D.3 for the rate on dist(xt,X ?).

heading: D PROOF OF PROPOSITION 5.2
text: Proposition D.1. Suppose that the assumptions of Theorem 4.2 hold with the probability measure µR fulfilling µR([r,R]) = Ω((R − r)κ) for r in [r0, R] for some r0 ∈ [0, R) and for some κ ∈ Z. Then,
lim t→∞
C2t
K2t,R∑t k=0 C2k
K2k,R
= 1− R 2
C2 . (58)
Proof. Given > 0, let c ∈ Z≥0 be the minimum such that
1∑c i=0 ( R2 C2 )i ≤ (1 + ) 1∑∞ i=0 ( R2 C2
)i = (1 + )(1− R2C2 )
(59)
Define Qt,R def = R 2t
K2t,R . Then,
C2t
K2t,R∑t k=0 C2k
K2k,R
= C2t R2tQt,R∑t k=0 C2k R2k Qk,R = Qt,R∑t
k=0
( R2
C2
)t−k Qk,R
(60)
Now, on one hand, using that Qt,R is an increasing sequence on t,
Qt,R∑t k=0 ( R2 C2 )t−k Qk,R ≥ 1∑t k=0 ( R2 C2 )t−k ≥ 1∑∞ k=0 ( R2 C2 )k = 1− R2C2 (61) On the other hand, for t ≥ c ,
Qt,R∑t k=0 ( R2 C2 )t−k Qk,R ≤ Qt,R∑t k=t−c ( R2 C2 )t−k Qk,R = Qt,R∑t k=t−c ( R2 C2 )t−k ( Qt,R − ∫ t k d dsQs,R ds ) (62)
Thus, we want to upper-bound ∫ t k d dsQs,R ds. First, notice that
d
ds Qs,R =
d
ds (∫ R 0 ( r R )2s dµR(r) )−1 = ∫ R 0 ( r R )2s (− log( rR )) dµR(r)(∫ R 0 ( r R )2s dµR(r)
)2 (63) By concavity of the logarithm function we obtain log(Rr ) ≤ R r0 − 1 for r ∈ [r0, R]. Choose r0 close enough to R so that Rr0 − 1 ≤ /c . We obtain that∫ R 0 ( r R )2s log ( R r ) dµR(r) ≤ ∫ r0 0 ( r R )2s log ( R r ) dµR(r) + ∫ R r0 ( r R )2s(R r0 − 1 ) dµR(r).
(64)
Thus,
∫ t k d ds Qs,R ds ≤ ∫ t k ∫ r0 0 ( r R )2s log ( R r ) dµR(r)(∫ R
0 ( r R )2s dµR(r)
)2 ds+ ∫ t k ∫ R r0 ( r R )2s ( R r0 − 1 )
dµR(r)(∫ R 0 ( r R )2s dµR(r) )2 ds. (65) Using that log x ≤ x, for k ∈ [t− c , t] we can bound the first term of (65) as∫ t
k
∫ r0 0 ( r R )2s log ( R r ) dµR(r)(∫ R
0 ( r R )2s dµR(r)
)2 ds ≤ ∫ t k ∫ r0 0 ( r R )2s−1 dµR(r)(∫ R
0 ( r R )2s dµR(r) )2 ds ≤ (t− k) ( r0 R
)2k−1(∫ R 0 ( r R )2t dµR(r)
)2 ≤ c (r0 R )2(t−c )−1 Q2t,R
≤ c (r0 R )2(t−c )−1 1 (c1)2 (2t+ 1)2κ t→∞−−−→ 0.
(66)
In the last inequality we use that by Proposition D.2, for t large enough, Qt,R = R 2t
K2t,R ≤ (2t +
1)k/c1. For k ∈ [t− c , t], the second term of (65) can be bounded as∫ t k ∫ R r0 ( r R )2s R r0
dµR(r)(∫ R 0 ( r R )2s dµR(r)
)2 ds ≤ (t− k)(Rr0 − 1 ) 1∫ R 0 ( r R )2t dµR(r)
≤ c ( R r0 − 1 ) 1∫ R 0 ( r R )2t dµR(r) ≤ Qt,R.
(67)
From (65), (66) and (67), we obtain that for t large enough, for k ∈ [t− c , t],∫ t k d ds Qs,R ds ≤ 2 Qt,R. (68)
Hence, we can bound the right-hand side of (62):
Qt,R∑t k=t−c ( R2 C2 )t−k ( Qt,R − ∫ t k d dsQs,R ds ) ≤ Qt,R∑t k=t−c ( R2 C2 )t−k (Qt,R − 2 Qt,R)
= 1 (1− 2 ) ∑t k=t−c ( R2 C2 )t−k = 1 (1− 2 ) ∑c k=0 ( R2 C2
)k ≤ 1 + 1− 2 ( 1− R 2 C2 ) .
(69)
The last inequality follows from the definition of c in (59). Since is arbitrary, by the sandwich theorem applied on (60), (61) and (69),
lim t→∞
C2t
K2t,R∑t k=0 C2k
K2k,R
= 1− R 2
C2 . (70)
Proposition D.2. Under the assumptions of Theorem 4.2, we have that there exists c1 > 0 such that for t large enough,
K2t,R ≥ c1R2t(2t+ 1)−κ. (71)
Proof. By the assumption on µR, there exist r0, c1, κ > 0 such that
K2t,R def = 2π ∫ R 0 r2t dµR(r) = 2π ∫ r0 0 r2t dµR(r) + 2π ∫ R r0 r2t dµR(r)
≥ 2πc1 ∫ R r0 r2t(R− r)κ−1 dr = −2πc1 ∫ r0 0 r2t(R− r)κ−1 dr + 2πc1 ∫ R 0 r2t(R− r)κ−1 dr ≥ −2πc1Rr2t0 + 2πc1R2t+κB(2t+ 1, κ). (72)
where the beta function B(x, y) is defined as
B(x, y) def = ∫ 1 0 rx+1(1− r)y+1 dr. (73)
Using the link between the beta function and the gamma function B(x, y) = Γ(x)Γ(y)/Γ(x + y), and Stirling’s approximation, we obtain that for fixed y and large x,
B(x, y) ∼ Γ(y)x−y. (74)
Hence, for t large enough, B(2t+ 1, κ) ∼ Γ(κ)(2t+ 1)−κ = (κ− 1)!(2t+ 1)−κ. Hence, from (72) we obtain that there exist c′1 depending only on κ and r0 such that for t large enough
K2t,R ≥ −2πc1Rr2t0 + 2πc1R2t+κ(k − 1)!(2t+ 1)−κ ≥ c′1R2t(2t+ 1)−κ. (75)
Proposition 5.2. Suppose that the assumptions of Theorem 4.2 hold with µR ∈ P([0, R]) fulfilling µR([r,R]) = Ω((R − r)κ) for r in [r0, R] for some r0 ∈ [0, R) and for some κ ∈ Z. Then, the average-case asymptotically optimal algorithm is, with y0 = x0:
yt = yt−1 − 1CF (yt−1), xt = ( R C )2 xt−1 + ( 1− ( R C )2) yt.
(20)
Moreover, the convergence rate for this algorithm is asymptotically the same one as for the optimal algorithm in Theorem 4.2. Namely, limt→∞ E [dist(xt,X ?)]Bt = 1.
Proof. The proof follows directly from Theorem 4.2 and Proposition D.1. See (77) and (79) in Proposition D.3 for the statement regarding the convergence rate.
Proposition D.3. For the average-case optimal algorithm (18),
Edist(xt,X ?) = ξopt(t) def = 1∑t k=0 C2k
K2k,R
(76)
For the average-case asymptotically optimal algorithm (20),
Edist(xt,X ?) = ξasymp(t) def = ( 1− ( R
C )2)2 t∑ k=1 K2k,R C2k ( R C )4(t−k) + ( R C )4t (77)
For the iterates yt in (18), i.e. gradient descent with stepsize 1/C, we have
Edist(yt,X ?) = ξGD(t) def = K2t,R C2t
(78)
Moreover, for all t ≥ 0, we have ξopt(t) ≤ ξasymp(t), and under the assumptions of (5.1),
lim t→∞
ξopt(t)
ξasymp(t) = 1, lim t→∞
ξopt(t) ξGD(t) = ξasymp(t) ξGD(t) = 1−
( R
C
)2 (79)
Proof. To show (76), (77), (78), we use the expression xt−x? = Pt(A)(x0−x?) (Proposition 2.1) and then evaluate ‖Pt‖2µ = ∫ C\{0} |Pt| 2 dµ (Theorem 2.1).
For (76), the value of ‖Pt‖2µ follows directly from Theorem 2.3, which states that the value for the optimal residual polynomial Pt is
1∑t k=0 |φk(0)|2 = 1∑t k=0 C2k
K2k,R
. (80)
A simple proof by induction shows that for the asymptotically optimal algorithm (20), the following expression holds for all t ≥ 0:
xt − x? =
(( R
C
)2t + ( 1− ( R
C )2) t∑ k=1 ( 1− A C )k ( R C )2(t−k)) (x0 − x?) (81)
Thus,
Pt(λ) =
( R
C
)2t + ( 1− ( R
C )2) t∑ k=1 ( 1− λ C )k ( R C )2(t−k)
=
( R
C
)2t φ0(λ) + ( 1− ( R
C )2) t∑ k=1 Kk,R Ck φk(λ) ( R C )2(t−k) ,
(82)
which concludes the proof of (77), as
‖Pt‖2µ =
( 1− ( R
C )2)2 t∑ k=1 K2k,R C2k ( R C )4(t−k) + ( R C )4t . (83)
By equation (52), yt − x? = (
1− A C
)t (y0 − x?) =
Kt,R Ct φk(A)(y0 − x?) (84)
Thus, for the yt iterates, ‖Pt‖2µ = K2t,R C2t , and (78) follows.
Now, ξopt(t) ≤ ξasymp(t),∀t ≥ 0 is a consequence of ξopt(t) being the rate of the optimal algorithm. And
lim t→∞
ξopt(t) ξGD(t) = lim t→∞
C2t
K2t,R∑t k=0 C2k
K2k,R
= 1− R 2
C2 (85)
follows from Proposition D.1. To show limt→∞ ξopt(t) ξGD(t) = 1 − R 2
C2 , which concludes the proof, we rewrite
ξasymp(t) =
( R
C )2t(1− (R C )2)2 t∑ k=1 1 Qk,R ( R C )2(t−k) + ( R C )2t , (86) using that by definition, Qk,R = R2k/K2k,R. Now, let c ∈ Z≥0 such that
∞∑ k=c ( R C )2k ≤ . (87)
Using the same argument as in Proposition D.1 (see (68)), for t large enough and k ∈ [t− c , t],∫ t k d ds Qs,R ds ≤ 2 Qt,R. (88)
Hence, for t large enough,( 1− ( R
C )2)2 t∑ k=1 1 Qk,R ( R C )2(t−k) + ( R C )2t
= ( 1− ( R
C )2)2( t∑ k=t−c
1 Qt,R − ∫ t k d dsQs,R
( R
C
)2(t−k) + t−c ∑ k=1 1 Qk,R ( R C )2(t−k)) + ( R C )2t
≤ ( 1− ( R
C
)2)2( 1
(1− 2 )Qt,R t∑ k=t−c ( R C )2(t−k) + t−c ∑ k=1 ( R C )2(t−k)) +
≤ ( 1− ( R
C
)2)( 1
(1− 2 )Qt,R +
( 1− ( R
C
)2) ) + ,
(89) which can be made arbitrarily close to ( 1− ( R C )2) 1 Qt,R by taking > 0 small enough. Plugging
this into (86), we obtain that we can make ξasymp(t) arbitrarily close to ( 1− ( R C )2) (R C )2t 1 Qt,R
=( 1− ( R C )2) ξGD(t) by taking t large enough.

Review summaries:
In this submission, first-order methods for solving smooth games are studied in the average case. In particular, first-order methods are derived and studied that are average-case optimal for certain optimization problems. In particular average-optimal first-order methods for solving zero-sum minimax games are presented. Also for finding the root of non-symmetric affine operators average-case optimal operators are derived if either the relevant matrix is normal or the eigenvalues are supported in a disk. Some experiments with the derived methods are conducted but the focus lies clearly on the theoretical results.

The paper combines two recent areas of interest in optimization for machine learning: average-case analysis for acceleration and learning in bilinear games.  Average-case optimal methods are proposed for some distributions on eigenvalues for the Jacobian of the games vector-field. Also, a connection is made between average-case optimal methods in bilinear games and optimally solving the Hamiltonian.