Science, Faculty of
Computer Science, Department of
DSpace
UBCV
Fang, Huang
Doctor of Philosophy - PhD
University of British Columbia
First-order methods are gaining substantial interest in the past two decades because of their superior performance in solving today's large-scale problems. In this thesis, we study some widely used first-order methods for problems that satisfy certain structures. Specifically, in the first part, we contribute to coordinate optimization and show that greedy coordinate descent (GCD) has an implicit screening ability that usually selects coordinates that are nonzero at the solution, which explains why GCD works exceptionally well for problems that admit sparse solutions. We also extend the elegant safe-screening rule that depends on duality gap to atomic-norm regularized problems. In the second part, we study online mirror descent (OMD) with unknown time horizon and unbounded domain, which is known to suffer from linear regret. We provide a stabilization technique and show that the stabilized-OMD can achieve sublinear regret. We also build the connection between stabilized-OMD and dual averaging. In the third part, we derive improved iteration complexity of the stochastic subgradient method for over-parameterized models that satisfy an interpolation condition. The obtained iteration complexity matches the rate of the stochastic gradient method applied to smooth problems that also satisfy an interpolation condition. Our analysis partially explains the empirical observation that nonsmoothness in modern machine learning models sometimes does not slow down the training process.
First-order methods for structured optimizationbyHuang FangB.Sc, The Central University of Finance and Economics, China, 2011M.Sc, The University of California, Davis, US, 2015A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFDoctor of PhilosophyinTHE FACULTY OF GRADUATE AND POSTDOCTORALSTUDIES(Computer Science)The University of British Columbia(Vancouver)October 2021© Huang Fang, 2021The following individuals certify that they have read, and recommend to the Facultyof Graduate and Postdoctoral Studies for acceptance, the thesis entitled:First-order methods for structured optimizationsubmitted by Huang Fang in partial fulfillment of the requirements for the degreeof Doctor of Philosophy in Computer Science.Examining Committee:Michael P. Friedlander, Computer ScienceSupervisorNicolas J.A. Harvey, Computer ScienceSupervisory Committee MemberMark Schmidt, Computer ScienceSupervisory Committee MemberBrian Wetton, MathematicsUniversity ExaminerBruce Shepherd, Computer ScienceUniversity ExaminerChih-Jen Lin, Computer Science, National Taiwan UniversityExternal ExamineriiAbstractFirst-order methods are gaining substantial interest in the past two decades becauseof their superior performance in solving today’s large-scale problems. In this thesis,we study some widely used first-order methods for problems that satisfy certainstructures. Specifically, in the first part, we contribute to coordinate optimizationand show that greedy coordinate descent (GCD) has an implicit screening abilitythat usually selects coordinates that are nonzero at the solution, which explainswhy GCD works exceptionally well for problems that admit sparse solutions. Wealso extend the elegant safe-screening rule that depends on duality gap to atomic-norm regularized problems. In the second part, we study online mirror descent(OMD) with unknown time horizon and unbounded domain, which is known tosuffer from linear regret. We provide a stabilization technique and show that thestabilized-OMD can achieve sublinear regret. We also build the connection betweenstabilized-OMD and dual averaging. In the third part, we derive improved iterationcomplexity of the stochastic subgradient method for over-parameterized models thatsatisfy an interpolation condition. The obtained iteration complexity matches therate of the stochastic gradient method applied to smooth problems that also satisfyan interpolation condition. Our analysis partially explains the empirical observationthat nonsmoothness in modern machine learning models sometimes does not slowdown the training process.iiiLay SummaryOptimization algorithms are prevalent in modern data-driven applications. Therecent immense growth of data volumes has spurred the rapid development of newvarieties of algorithms efficient for large problems. Despite their success in practice,the theory behind these algorithms is incomplete and often fails to explain empiricalperformance. This thesis aims to fill some of these gaps, contribute to a betterunderstanding of some widely used optimization algorithms, and shed light on thedesign of new algorithms.ivPrefaceThe main body of this thesis is based several collaborative papers that are eitherpublished or currently under review.• The work presented in the first part of Chapter 2 (Section 2.3) was publishedin AISTATS, 2020 (Fang et al., 2020a). I am the primary contributor tothis paper. The convergence analysis and experiments are developed by me.Zhenan Fan participated in some discussions and was involved in polishingthe proofs. Yifan Sun improved the writing and presentation of the paperand provided some useful suggestions. Michael Friedlander provided overallsupervision, technical guidance, writing assistance, and financial support.• The material in the second part of Chapter 2 (Section 2.4) was submitted andcurrently under review. Zhenan Fan and me contributed equally to this paper,the main theorem of this paper is based on our discussion and the case study onnuclear-norm regularized problem is developed by me. Michael Friedlanderprovided overall supervision, technical guidance, writing assistance, andfinancial support.• The material in Chapter 3 was published in ICML, 2020 (Fang et al., 2020b).This paper originates from a conjecture in Nick Harvey’s graduate courseCPSC531H. I disproved the conjecture when taking Nick Harvey’s course andextended this result to a research paper. The stabilization technique owns toNick Harvey. I proposed to apply the stabilization technique to online mirrordescent (OMD) and developed the sublinear regret bound. The first-orderregret bound is based on my discussion with Nick Harvey. Victor Portellaprovided the analysis to compare stabilized-OMD and dual averaging. NickvHarvey, Victor Portella and me were heavily involved in writing up the proofs.Michael Friedlander provided overall supervision, technical guidance, writingassistance, and financial support.• The material in Chapter 4 was published in ICLR, 2021 (Fang et al., 2021). Iam the primary contributor to this paper. The improved iteration complexities,lower bound analysis, experiments and initial draft are from me. Zhenan Fanparticipated in some discussions and checked the proofs. Michael Friedlanderprovided overall supervision, technical guidance, writing assistance, andfinancial support.viTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Notation and preliminaries . . . . . . . . . . . . . . . . . . . . . 21.2 Gradient and subgradient descent . . . . . . . . . . . . . . . . . . 31.3 Coordinate descent . . . . . . . . . . . . . . . . . . . . . . . . . 51.4 Stochastic gradient descent . . . . . . . . . . . . . . . . . . . . . 71.5 Mirror descent . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91.6 Summary of contributions . . . . . . . . . . . . . . . . . . . . . 102 Coordinate descent and sparse optimization . . . . . . . . . . . . . . 122.1 Different selection rules . . . . . . . . . . . . . . . . . . . . . . . 122.2 Greedy coordinate descent (GCD) . . . . . . . . . . . . . . . . . 152.3 GCD for sparse optimization . . . . . . . . . . . . . . . . . . . . 17vii2.3.1 Problem setup . . . . . . . . . . . . . . . . . . . . . . . 182.3.2 Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . 192.3.3 Improved selection rule . . . . . . . . . . . . . . . . . . . 242.3.4 Numerical experiments . . . . . . . . . . . . . . . . . . . 272.3.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . 302.4 Gap-based safe-screening rules for atomic-norm regularized problem 312.4.1 From coordinates to atoms . . . . . . . . . . . . . . . . . 312.4.2 Some technical tools . . . . . . . . . . . . . . . . . . . . 322.4.3 Problem setup . . . . . . . . . . . . . . . . . . . . . . . 332.4.4 The gap-based safe-screening rule . . . . . . . . . . . . . 342.4.5 Gap-based safe-screening rule for nuclear norm . . . . . . 372.4.6 Approximation with partial SVD . . . . . . . . . . . . . . 382.4.7 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . 433 Online mirror descent with unknown time horizon . . . . . . . . . . 443.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 453.1.1 Definitions and notations . . . . . . . . . . . . . . . . . . 453.1.2 OMD and DA with constant stepsize . . . . . . . . . . . . 483.1.3 OMD and DA with dynamic stepsize . . . . . . . . . . . 483.2 Stabilized OMD . . . . . . . . . . . . . . . . . . . . . . . . . . . 493.2.1 Dual-stabilized OMD . . . . . . . . . . . . . . . . . . . . 503.2.2 Primal-stabilized OMD . . . . . . . . . . . . . . . . . . . 533.2.3 Dual averaging . . . . . . . . . . . . . . . . . . . . . . . 563.2.4 Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . 603.3 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 603.3.1 Strongly-convex mirror maps . . . . . . . . . . . . . . . 603.3.2 Prediction with expert advice . . . . . . . . . . . . . . . . 613.4 Comparing DS-OMD and DA . . . . . . . . . . . . . . . . . . . 663.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 694 Fast convergence of stochastic subgradient descent under interpolation 704.1 Background and motivation . . . . . . . . . . . . . . . . . . . . . 704.1.1 Practical algorithms based on SGD . . . . . . . . . . . . 71viii4.1.2 Parallel and distributed SGD . . . . . . . . . . . . . . . . 724.1.3 Variance reduction . . . . . . . . . . . . . . . . . . . . . 734.1.4 SGD with the interpolation condition . . . . . . . . . . . 734.2 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . 744.3 Main results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 764.3.1 Bounds and Lipschitz properties of the generalized gradient 764.3.2 Convergence rate of stochastic subgradient descent . . . . 794.3.3 Lower bounds . . . . . . . . . . . . . . . . . . . . . . . . 834.4 Numerical experiments . . . . . . . . . . . . . . . . . . . . . . . 854.4.1 Teacher-student setup . . . . . . . . . . . . . . . . . . . . 864.4.2 Classify 4’s and 9’s on MNIST dataset . . . . . . . . . . . 864.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 875 Conclusion and future work . . . . . . . . . . . . . . . . . . . . . . 885.1 Coordinate optimization . . . . . . . . . . . . . . . . . . . . . . 885.2 Mirror descent . . . . . . . . . . . . . . . . . . . . . . . . . . . . 895.3 Stochastic subgradient descent . . . . . . . . . . . . . . . . . . . 89Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91A Appendix for Chapter 2 . . . . . . . . . . . . . . . . . . . . . . . . . 104A.1 Proofs for Section 2.3 . . . . . . . . . . . . . . . . . . . . . . . . 104A.2 Proofs for Section 2.4 . . . . . . . . . . . . . . . . . . . . . . . . 113B Appendix for Chapter 3 . . . . . . . . . . . . . . . . . . . . . . . . . 115B.1 Standard facts . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115B.1.1 Scalar inequalities . . . . . . . . . . . . . . . . . . . . . 115B.1.2 Bregman divergence properties . . . . . . . . . . . . . . . 117B.2 Proofs for Section 3.3.1 . . . . . . . . . . . . . . . . . . . . . . . 119B.3 Proofs for Section 3.3.2 . . . . . . . . . . . . . . . . . . . . . . . 119B.4 Proofs for Section 3.4 . . . . . . . . . . . . . . . . . . . . . . . . 120C Appendix for Chapter 4 . . . . . . . . . . . . . . . . . . . . . . . . . 125C.1 Proofs for Section 4.3 . . . . . . . . . . . . . . . . . . . . . . . . 125ixList of TablesTable 2.1 Properties of the experimental data. Here, d denotes the numberof features and n denotes the number of samples. . . . . . . . . 27Table 2.2 Commonly used sets of atoms and their gauge and support func-tion representations . . . . . . . . . . . . . . . . . . . . . . . 33Table 4.1 Iteration complexity of deterministic gradient and stochasticgradient methods. . . . . . . . . . . . . . . . . . . . . . . . . 71Table 4.2 Iteration complexity of SGD with and without the interpolationcondition. SC stands for strongly convex. . . . . . . . . . . . . 74xList of FiguresFigure 1.1 Illustration of gradient descent . . . . . . . . . . . . . . . . . 3Figure 1.2 Illustration of cooridnate descent . . . . . . . . . . . . . . . . 5Figure 1.3 Illustration of stochastic gradient descent . . . . . . . . . . . 7Figure 2.1 Exploratory investigations. . . . . . . . . . . . . . . . . . . . 20Figure 2.2 Illustrations for Theorem 2.3.1 and 2.3.2 . . . . . . . . . . . . 23Figure 2.3 Comparison between different kinds of initialization . . . . . 28Figure 2.4 Compare ∆-GCD with different choices of ∆. . . . . . . . . . 29Figure 3.1 Illustration of the t-th iteration of DS-OMD. . . . . . . . . . 50Figure 4.1 The performance of SSGD with smooth and nonsmooth lossfunctions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87xiAcknowledgmentsThis thesis would not have been possible without the help and kindness of manypeople throughout the last several years. First, I would like to thank my supervisorMichael Friedlander for providing constant support and giving me the freedom toexplore a diverse set of projects. I have learned very much from Michael’s broadknowledge in optimization and computation. His rigorous research attitude andcheerful spirit have given me an example to pursue.I would also like to thank Cho-Jui Hsieh, who led me on the path of researchand helped me building my research skills. He has been a role model for me topursue throughout my Ph.D. study. I am also grateful to Nick Harvey, with whomI have had the privilege to collaborative. I learned how to think as a theoreticianfrom him. Next, I would like to thank the university examiners, Bruce Shepherdand Brian Wetton, and the external examiner, Chih-Jen Lin for reading my thesisand providing valuable feedback.I am thankful to my other collaborators Zhenan Fan, Victor Portella, Yifan Sunfor helping me to be more productive than I would have been able to on my own.I would like to thank Naomi Graham and Emma Hansen for proofreading somematerials of the thesis and providing helpful feedback to improve the presentation. Iwould like to thank my friends Zhenan Fan, Wen Xiao, Liran Li, Xin Ding, QiuyanLiu and Anyi Bao for the company during my Ph.D. study, they make my life inVancouver more joyful.Last but not least, my academic journey is not possible without the understandingand encouragement from my parents. I would like to thank them as they alwayssupport me to pursue what I want. I would like to dedicate this thesis to them.xiiChapter 1IntroductionOptimization is a fundamental aspect of many fields, including machine learning,data mining, signal processing, and bioinformatics. In modern applications suchas recommender systems, computer vision, and natural language processing, boththe amount of data and the model’s size can be very large. For example, someof today’s social networks could have hundreds of millions of users, and somerecent neural network architectures could have billions of parameters. As a result,efficient optimization algorithms that can handle both big data and complex modelsare of great interest in recent years. Many optimization algorithms were proposedor re-discovered in the last two decades to solve today’s large-scale optimizationproblems. Due to the low numerical accuracy required by most machine learningand data mining tasks, first-order optimization methods with cheap per-iterationcomputational cost dominate certain fields in machine learning. In particular,stochastic gradient descent (SGD) and coordinate descent (CD) are the two mostimportant representatives. SGD and its variants are dominant for the big-N problems— i.e., problems with a large number of samples, while CD and its variants arehighly efficient in handling the structured big-p problems — i.e., models with alarge number of parameters.The empirical success of first-order methods has driven extensive researchin recent years on understanding the effectiveness of these methods as well asdeveloping new variants of first-order methods for emerging applications. However,despite some notable progress made by some pioneering researchers in recent years,1there are still many open problems left unsolved in this area. Given the importanceand popularity of first-order optimization algorithms, there is a great need to solvethese open problems. This thesis aims to solve some of these open problems andcontribute to a better understanding of first-order optimization methods in differentscenarios. In this chapter, we summarize some classical results of some importantfirst-order algorithms and describe the contributions of this thesis.1.1 Notation and preliminariesThroughout this thesis, unless otherwise specified, we use capital letters A,B, . . .to denote matrices, lowercase letters x,y,w, . . . to denote vectors, except L whichusually denote the smoothness parameter, and m,n,d which are commonly usedas the dimension of matrices or vectors. We use Greek letters α,β ,γ, . . . to denotescalars. The ith entry of a given vector x is denoted as xi, and we use x(t) to representthe t-th iterate in an algorithm. The norm ‖ ·‖ stands for Euclidean-norm. We define[n] to be the set {1,2, . . . ,n}. We denote a solution of an optimization problem tobe x∗ and the optimal function value f (x∗) to be f ∗, we assume that x∗ exist andf ∗ >−∞ throughout this thesis.Before proceeding to the classical results of first-order methods, we review somestandard definitions that are commonly used in the literature of convex optimization.Definition 1.1.1. A function f : Rd → R is convex iff (λx+(1−λ )y)≤ λ f (x)+(1−λ ) f (y)for any x,y ∈ Rd , and any λ ∈ [0,1].Definition 1.1.2. For a closed, convex function f : Rd → R, the subdifferential off at x is defined as∂ f (x) :={v ∈ Rd∣∣∣ f (y)≥ f (x)+ 〈v,y− x〉, ∀y ∈ Rd} .Definition 1.1.3. The dual norm of a given norm ‖ · ‖ is defined as‖z‖∗ := sup{〈z,x〉 | ‖x‖ ≤ 1,x ∈ Rd}.2x(0)Figure 1.1: Illustration of gradient descentDefinition 1.1.4. A function f : Rd → R is γ-Lipschitz continuous with respect toa norm ‖ · ‖ for some γ > 0 if ∀x,y ∈ Rd ,| f (x)− f (y)| ≤ γ‖x− y‖.Definition 1.1.5. A function f : Rd → R is L-smooth with respect to a norm ‖ · ‖for some L> 0 if it is differentiable and ∀x,y ∈ Rd ,‖∇ f (x)−∇ f (y)‖ ≤ L‖x− y‖∗,where ‖ · ‖∗ is the dual norm paired with ‖ · ‖.Definition 1.1.6. A function f : Rd → R is µ-strongly convex with respect to anorm ‖ · ‖ for some µ ≥ 0 if ∀x,y ∈ Rd ,f (x)≥ f (y)+ 〈g,x− y〉+ µ2‖x− y‖2, ∀g ∈ ∂ f (y).1.2 Gradient and subgradient descentGradient and subgradient descent (GD, subGD) are perhaps the oldest first-orderoptimization methods that can be traced back to 1847 (Cauchy, 1847). Consider theunconstrained optimization problemminx∈Rdf (x), (1.1)3where f is a continuous function. In the t-th iteration, (sub)GD updates the iterateby moving it to the negative (sub)gradient directionx(t+1) = x(t)−ηtg, g ∈ ∂ f (x(t)),where ηt is the stepsize used in the t-th iteration. An illustration of GD is shown inFigure 1.1. It is known that GD has the following convergence properties:• When f is convex and L-smooth, GD with constant stepsize ηt = 1/L isguaranteed to achieve an ε-accurate solution1 inO(Lε−1) iterations (Nesterov,2004). This rate can be improved to O(Lε−1/2) by adopting Nesterov’sacceleration technique (Nesterov, 1983) and was proven to be optimal forfirst-order methods in the convex and smooth setting;• When f is µ-strongly convex and L-smooth, GD with constant stepsizeηt = 1/L achieves a linear convergence rate O((L/µ) log(ε−1))(Nesterov,2004);• When f is convex but not necessarily smooth, GD with decaying stepsizeηt ∝ t−1/2 converges at the rateO(Lε−2) (Nemirovski and Yudin, 1983; Shor,1984). If we treat L as a constant, then the rate O(ε−2) has been proven to beoptimal for first-order methods (Nemirovski and Yudin, 1983) in the convexand nonsmooth setting;• When f is L-smooth but not necessarily convex, GD may not converge to aglobal solution. Instead, GD with constant stepsize ηt = 1/L can convergeto an ε-stationary point2 in O(Lε−1) iterations. This is a well-known resultthat can be found in optimization textbooks, but its origin is not known to theauthor’s knowledge.Note that many results in the literature treat the smoothness or strongly convexparameters as constants and omit them in the big-O notation. However, as wewill show in Section 1.3, these parameters play an important role in understandingthe convergence of coordinate descent. In order to compare the convergence rates1A point x is an ε-accurate solution if f (x)− f ∗ ≤ ε .2A point x is an ε-stationary point if ‖∇ f (x)‖2 ≤ ε .4x(0)Figure 1.2: Illustration of cooridnate descentbetween different first-order algorithms, we will not treat them as constants in thischapter.1.3 Coordinate descentCoordinate descent (CD) is a simple extension of GD. Instead of updating thewhole vector x in each iteration, CD selects a single coordinate (or a small setof coordinates) according to some coordinate selection rules and apply standardgradient step to them. Consider the single coordinate selection case with the updaterulex(t+1) = x(t)−ηt∇it f (x(t))eit ,where it is the coordinate selected by some selection rules in the t-th iteration andeit is a zero vector with a one in the it-th entry. An illustration of CD is shown inFigure 1.2.CD is an old algorithm. Its history can be traced back to 1874 as an iterativemethod due to Seidel to solve linear systems, also known as the Gauss-Seidelmethod. Empirically, CD is simple and effective. Its superior performance againstother classical algorithms on some machine learning problems makes it a popularoptimizer for a wide range of applications, including clustering (Lloyd, 1982),support vector machines (Chang and Lin, 2011; Joachims, 1999; Platt, 1998) andLASSO (Hastie et al., 2008; Sylvain Sardy and Tseng, 2000). On the theory side,pioneering researchers made important contributions to understand the convergence5behavour of CD in late 20th century (Bertsekas and Tsitsiklis, 1989; Luo and Tseng,1992, 1993; Tseng, 2001). However, the global convergence rate of CD was notclear until Nesterov’s seminal work (Nesterov, 2012), in which Nesterov explicitlyanswered why CD could be significantly faster than GD.Consider problem 1.1, where f is L-smooth. Let Li denote the coordinate-wisesmoothness parameter of f , namely for any i ∈ [d], i.e.,|∇i f (x+αei)−∇i f (x)| ≤ Li|α|, ∀x ∈ Rd ,α ∈ R.Let Lmax := maxi∈[d]Li. It is easy to verify that Lmax ≤ L ≤ dLmax. With a simplerandom coordinate selection rule, Nesterov demonstrated that one could adopt alonger stepsize 1/Lmax compared with the stepsize 1/L used in GD when updatinga single coordinate. Based on this observation, Nesterov established the followingconvergence rates for randomized CD (RCD):• When f is convex, RCD is guaranteed to obtain an ε-accurate solution withinO(dLmaxε−1) iterations in expectation;• When f is µ-strongly convex, RCD has the rate O ((dLmax/µ) log(ε−1)) inexpectation.Recall that Lmax ≤ L≤ dLmax, the above rates indicate that RCD is slower than GD(see the rates in Section 1.2). However, Nesterov pointed out that for a wide rangeof problems, RCD can be implemented cleverly such that running d RCD stepscost roughly the same time as one GD iteration. Note that the convergence rate ofRCD in terms of epoch (d RCD iterations) is O(Lmaxε−1) for convex objective andO ((Lmax/µ) log(ε−1)) for strongly convex objective. Therefore RCD is faster thanGD in terms of runtime. We refer interested readers to Nesterov (2012) for moredetails.Extensive research has been carried on since Nesterov’s 2012 work. Here wejust list some notable results: convergence analysis for proximal RCD (Richta´rikand Taka´c, 2014), different coordinate selection rules (Beck and Tetruashvili, 2013;Lee and Wright, 2018; Nutini et al., 2015; Recht and Re´, 2012; Saha and Tewari,2013; Sun and Ye, 2021), accelerated RCD (Allen Zhu et al., 2016; Lee and Sidford,2013; Lin et al., 2015), and parallel CD with multiple computing cores (Bradley6x(0)Figure 1.3: Illustration of stochastic gradient descentet al., 2011; Hannah et al., 2019; Jaggi et al., 2014; Liu et al., 2014a; Richta´rik andTaka´c, 2011; You et al., 2016).1.4 Stochastic gradient descentFor many machine learning and data mining applications, we are interested insolving an empirical-risk minimization problem that can be written as a finite summinimization problemminx∈Rdf (x) :=1nn∑i=1fi(x),where n is the number of training samples that could be very large in modernapplications, for example 107 for the imagenet dataset3.GD or subGD are expensive in this setting since they need to call n gradientoracles for a single iteration. Instead of calculating the exact gradient in eachiteration, stochastic gradient descent (SGD) (Robbins and Monro, 1951) methodcalculates an unbiased estimator of the true gradient from a small subset of thewhole training samples. Consider sampling one training sample in each iteration,the update rule of SGD isx(t+1) = x(t)−ηtg, g ∈ ∂ fit (x(t)),for some it uniform randomly chosen from {1,2, . . . ,n} in each iteration. An3http://www.image-net.org/7illustration of SGD is shown in Figure 1.3. Obviously, the per iteration cost ofSGD is n time cheaper than GD. At the cost of using cheap and noisy gradient forupdate, SGD converges slower than GD in terms of the number of iterations. SGDhas the following convergence guarantee under the bounded variance assumption(supx∈Rd E[‖∇ fi(x)−∇ f (x)‖]≤ σ for some σ > 0):• When f is convex and γ-Lipschitz continuous, SGD with decaying stepsizeηt ∝ t−1/2 could obtain an ε-solution within O((γ2+σ2)ε−2) iterations (Ne-mirovski et al., 2009) in expectation;• When f is µ-strongly convex and γ-Lipschitz continuous, SGD convergesto global minimum with a rate O((γ2 +σ2)µ−1ε−1) (Lacoste-Julien et al.,2012; Nemirovski et al., 2009) in expectation;• When f is L-smooth but not necessarily convex, SGD with decaying stepsizeηt ∝ t−1/2 could approach an ε-stationary point within O(Lε−1 +Lσ2ε−2)iterations (Ghadimi and Lan, 2013) in expectation.Compared with the rates of GD stated in Section 1.2, the convergence rates ofSGD described above are worse. However, an extraordinary fact of SGD is that itsconvergence rates are independent of the total number of samples n. Therefore, SGDfits perfectly to applications with a large number of training samples and requirelow numerical accuracy. The superior empirical performance of SGD makes it andits variants such as ADAM (Kingma and Ba, 2015) the default choice for trainingmodern machine learning models. Note that modern machine learning modelsoften yield nonconvex and nonsmooth objective, but SGD still usually performsunreasonably well in these tasks.An important advance of SGD in the last decade is the variance reductiontechnique proposed by Schmidt et al. (2017), who proposed the first linear conver-gent SGD (under strong convexity) based algorithm called SAG (Schmidt et al.,2017). Following this work, there is a line of research that tries to improve SGDwith variance reduction under different scenarios. Examples include the well-known SVRG (Johnson and Zhang, 2013), SAGA (Defazio et al., 2014), andSARAH (Nguyen et al., 2017) algorithms.81.5 Mirror descentThe mirror descent (MD) algorithm (Beck and Teboulle, 2003; Nemirovski andYudin, 1983) is a generalization of the classical projected subgradient descent(PGD) method. The major difference between MD and PGD is that MD explicitlydistinguishes the “primal” and “dual” spaces and applies the gradient update in thedual space. Consider problem 1.1 with the constraint x ∈ X , the update rule of MDcan be written asx(t+1) = argminx∈X{〈x,g〉+ 1ηtBΦ(x,x(t))}, g ∈ ∂ f (x(t)),where BΦ is the Bregman divergence induced by the mirror map Φ : Rd → R,detailed explanation of these technical terms will be described in Chapter 3.We can recover the update rule of PGD by choosing the mirror map Φ to be2-norm squared. Therefore PGD is a special case of MD. MD’s main advantage overPGD is that MD, for certain applications, can enjoy a better Lipschitz constant andsimpler projection than PGD by choosing an appropriate mirror map. A typical ex-ample that favours MD than PGD is when the constraint set is a probability simplex,i.e., X = {∑di=1 xi = 1,xi ≥ 0}. In this case, MD has the following advantages:• When f is γ-Lipschitz continuous with respect to ‖ · ‖∞. It follows that f is√dγ-Lipschitz continuous with respect to ‖ · ‖2 and the number of iterationsrequired to obtain an ε-solution would be O(√dγε−2). By using negativeentropy as the mirror map, MD can improve the rate toO(Lε−2), which couldbe significant when d is large.• The projection in Euclidean space requires us to do a sorting while theprojection used in MD is a simple scaling operation by using negative entropyas the mirror map.Given the above observations, MD with the negative entropy mirror map (alsoknown as the exponentiated gradient descent) perfectly fits problems with Lipschitzconstant related to infinity norm and probability simplex constraints. A classicexample of this type is the prediction with expert advice problem in the onlinelearning community. Therefore MD is an useful algorithmic template for online9learning problems and has been heavily studied in recent years.1.6 Summary of contributionsAs we described in previous sections, the theory of first-order optimization algo-rithms is heavily studied. However, there are some deficiencies with the currentanalysis:• For some algorithms, the current analyses work for quite general problems.At the cost of generality, the analyses are sometimes too pessimistic and failto explain some common empirical observations for problems that satisfycertain structures.• For some other algorithms, the analyses work only for a restricted class ofproblems and cannot be applied to problems that satisfy some other importantstructures, such as problems with atomic-norm regularization instead of theclassic 1-norm regularization.This thesis aims to improve and extend the classical analysis of some first-orderoptimization algorithms for problems that satisfy certain structures.In Chapter 2, we focus on coordinate descent and sparse optimization. First,we consider applying coordinate descent with greedy selection rule (GCD) forcomposite objectives in form f +g. When the composite problem admits sparsesolutions, empirical evidence in the literature suggests that GCD, when initializedas zero vector, has an implicit variable selection ability. Therefore, GCD canusually select coordinates that are nonzero at the solution and converge significantlyfaster than randomized CD. By leveraging the composite problem structure andsparse solution, we present an improved analysis of GCD for sparse optimizationand theoretically explained why GCD has the implicit variable selection property.Second, we consider problems with general atomic-sparsity and extend the analysisof the gap safe screening rule (Ndiaye et al., 2017) to general atomic sets. Wealso study the effectiveness of gap safe screening rule for problems with low-ranksolutions.In Chapter 3, we try to improve mirror descent in the context of online learning.In the literature of online convex optimization, online mirror descent (OMD) and10dual averaging (DA) are two fundamental algorithmic templates. They are knownto have a very similar (or even identical) performance guarantee in most scenarioswhen a fixed stepsize is used. However, for dynamic stepsize, OMD is provablyinferior to dual averaging. It is known that OMD with a dynamic stepsize schedulingcan suffer from linear regret. We modify the OMD algorithm through a simpletechnique called stabilization and give essentially the same abstract regret bound forstabilized-OMD and DA by modifying the classical OMD convergence analysis in acareful and modular way. Simple corollaries of these bounds show that OMD withstabilization and DA enjoy the same performance guarantees in many applicationseven under dynamic stepsize scheduling.In Chapter 4, we study the convergence behaviour of the stochastic subgradientdescent (SSGD) method applied to over-parameterized nonsmooth optimizationproblems that satisfy an interpolation condition. By leveraging the composite struc-ture of the empirical risk minimization problems, we prove that SSGD converges,respectively, with rates O(ε−1) and O(log(ε−1)) for convex and strongly convexobjectives when interpolation holds. These rates coincide with established ratesfor the stochastic gradient descent (SGD) method applied to smooth problems thatalso satisfy an interpolation condition. Our analysis provides a partial explanationfor the empirical observation that sometimes SGD and SSGD behave similarly fortraining smooth and nonsmooth machine learning models. We also prove that therate O(ε−1) is optimal for the subgradient method in the convex and interpolationsetting.11Chapter 2Coordinate descent and sparseoptimizationCoordinate descent (CD) is gaining increasing interest due to its simplicity andeffectiveness in the last two decades. It is the state-of-the-art optimization algorithmfor problems that satisfy the coordinate-friendly structure (Nutini et al., 2015; Penget al., 2016), that is the computation time required to update all coordinates isroughly the same as the time of one gradient descent step. Problems of this typeincluding the `1 regularized least squares (LASSO) (Sylvain Sardy and Tseng, 2000)and the support vector machines (SVMs) (Platt, 1998).2.1 Different selection rulesThe coordinate selection rule used in the CD algorithm plays an important rolein the convergence of CD. Five selection rules are commonly used in practice —uniform random, non-uniform random, cyclic, random permuted cyclic and greedyselection rules.• Randomized CD (RCD): RCD uniform randomly selects a coordinate from{1,2, . . . ,d} to update in each iteration. It is the version that Nesterov an-alyzed in his seminal work (Nesterov, 2012). The convergence analysis ofrandomized CD is simple since the coordinates selected in different iterationsare independent from each other. RCD has good theoretical properties, it is12provably faster (in terms of run time) than GD for problems that satisfy thecoordinate-friendly structure. However, the implementation of RCD is cacheunfriendly due to random access in memory.• Non-uniform random CD (Non-uniform RCD): Non-uniform RCD is sim-ilar to RCD except that it adopts non-uniform sampling to select coordinateto update. A typical example of non-uniform sampling is the Lipschitz sam-pling, this sampling strategy sample a coordinate to update with probabilitythat proportional to the coordinate’s Lipschitz constant. This non-uniformsampling strategy was first studied by Nesterov (2012). Nesterov showed thatthe Lipschitz sampling could improve the convergence rate of RCD when theLipschitz constants among coordinates are highly imbalanced.• Cyclic CD (CCD): CCD is the oldest version of CD and can be traced backto the Gauss-Seidel algorithm for solving linear systems. CCD predefine apermutation pi of {1,2, . . . ,d} and selects coordinate pi((i−1) mod d+1) in theith iteration. CCD is easy to implement and more cache friendly than RCDbecause of the nature of sequential access. Though CCD usually converges ata similar rate as (or sometimes even faster than) RCD in practice, its conver-gence analysis is hard due to the inherent dependency among the coordinatesselected. In fact, Nesterov stated that “it is almost impossible to estimatethe rate of convergence” for CCD (Nesterov, 2012). The non-asymptoticconvergence rate of CCD was later investigated by Saha and Tewari (2013),who obtained a non-asymptotic rate by imposing the isotonicity assumption,which is a strong assumption that may not hold in practice. Later Beck andTetruashvili (2013) obtained the first global convergence rate of CCD withoutany further assumptions (and actually with simple and clean proofs). However,the rate that Beck and Tetruashvili obtained for CCD is much worse than RCDand this result seems to be inconsistent with empirical observations, Beckand Tetruashvili stated that the rate they derived may be improvable. This gapbetween theory and practice was filled by Sun and Ye (2021), who formallyproved that there exist instances such that CCD is provably Ω(d2) slowerthan RCD (Sun and Ye, 2021). However, their counterexample does not applyto random permuted cyclic CD (Lee and Wright, 2018), which we describe13next.• Random permuted cyclic CD (RPCD): RPCD can be viewed as a midwaypoint between RCD and CCD. It produces a random permutation at thebeginning of each epoch and runs CCD with the permutation generated. Thecoordinates being selected are independent between epochs but dependentwithin each epoch. In terms of cache-friendliness, RPCD is better than RCDbut worse than CCD. RPCD is easy to implement and performs similarlytp (or even better than) RCD empirically in many cases. Many efficientsolvers for machine learning problems such as LIBLINEAR (Fan et al.,2008) are actually built on RPCD. The convergence rate of RPCD can bedirectly obtained from Beck and Tetruashvili’s analyses for CCD (Beckand Tetruashvili, 2013). But again, this rate cannot explain the empiricalperformance of RPCD. Provably showing that RPCD outperforms RCD forleast square problems turns out to be a hard mathematical problem: Recht andRe´ (2012) abstracted the RPCD versus RCD problem into a formal conjecturecalled the “matrix AMGM” inequality (Recht and Re´, 2012)1, which is thematrix version of the arithmetic-geometric mean inequality. The matrixAMGM inequality was recently proven to be false (Lai and Lim, 2020; Sa,2020) and there exist instances such that RPCD is provably inferior to RCD.• Greedy CD (GCD): GCD is also known as the Gauss-Southwell (GS) rulefor solving linear system. Different from RCD, CCD and RPCD, GCDselects the coordinate from argmaxi∈[d] |∇i f (x(t))|, which has the largestmarginal change at each iteration. It is worth noting that one may not beable to implement the greedy rule efficiently even if the problem satisfy thecoordinate-friendly structure. That is, the run time of one epoch of GCDis usually much more expensive than one gradient descent step. However,exceptions do exist, for example the dual problem of kernel SVM (a quadraticprogramming problem with box constraints) whose kernel matrix has similarnumber of nonzeros among its columns, in these cases GCD can be imple-1The original matrix AMGM inequality from Recht and Re´ (2012) was originally motivated by theproblem of random permuted SGD versus SGD, but its argument can be easily incorporated to CD.14mented efficiently2 and one epoch of GCD has roughly the same computationcost as one gradient descent iteration. The state-of-the-art solvers for kernelSVM — LIBSVM (Chang and Lin, 2011) and SVMLight (Joachims, 1999)both are built on GCD.The computational costs of the above selection rules are different. Random,cyclic and random permuted cyclic selection rule shares similar property. They canall benefit from the coordinate-friendly structure, that is, the runtime of d iterationswith these three selection rules cost roughly the same time as one gradient descentstep. For GCD, as mentioned above, the coordinate-friendly structure alone is notsufficient for it to have efficient implementation. To implement the greedy selectionrule efficiently, stronger problem structures (such as the Hessian matrix need tobe both row- and column-wise sparse, see a more comprehensive discussion inNutini et al. (2015)) are required. The computational cost of non-uniform RCD issimilar to GCD, and it depends on the sparsity structure of problems. For problemswith a balanced sparsity structure (for example, the least-square problem with abalanced number of non-zeros for different columns), the time required to updatedifferent coordinate costs similar time. Therefore, d steps of non-uniform RCD costsimilar time as one gradient descent step in this case. For problems with a highlyunbalanced sparsity structure (for example, least-square problem with one columncompletely dense and other columns very sparse), the computational cost of RCDwith the Lipschitz sampling could be very high; one coordinate update could be asexpensive as one gradient descent step in the worst case.2.2 Greedy coordinate descent (GCD)GCD and GD are closely related. Both of these two methods can be interpretedas steepest descent method on the first-order Taylor expansion of the objectivefunction, where GD is the steepest descent with respect to 2-norm and GCD isthe steepest descent with respect to 1-norm (Boyd and Vandenberghe, 2004, §9.4).Different from the convergence rate of RCD, the convergence of GCD is more“similar” to GD and both of them enjoy dimension free convergence rate. For convex2But require to store the kernel matrix in memory.15and smooth (but not necessarily strongly convex) objective function f , GCD has theconvergence rate (see, for example, in work from Dhillon et al. (2011)3)f (x(t))− f ∗ ≤ LmaxR1(x(0))t,where R1(x(0)) = supx∈X ∗ ‖x− x(0)‖21, X ∗ is the solution set, Lmax is the maximumcoordinate-wise smoothness parameter, see Chapter 1.3 for its definition. Note thatthe convergence rate of GD (Nesterov, 2004) isf (x(t))− f ∗ ≤ LR(x(0))2t.where L is the smoothness parameter with respect to 2-norm and R(x(0)) is the initialcondition under 2-norm. Therefore convergence rates of both GCD and GD areindependent from the problem size d, their only difference is how the smoothnessparameter and distance to solution is measured, where the parameters of GCDdepend on 1-norm and the parameters of GD depend on 2-norm.For strongly convex objective, GCD has the convergencef (x(t))− f ∗ ≤(1− µ1Lmax)t (f (x(0))− f ∗).This rate is from the refined analysis of GCD from Nutini et al. (2015), where µ1 isthe strongly convex parameter with respect to 1-norm. Compared to the convergenceof GD for strongly convex objectivef (x(t))− f ∗ ≤(1− µ2L)t (f (x(0))− f ∗),which is built on the smoothness and strongly convex parameters (L and µ2) withrespect to 2-norm. By norm inequality, we know that µ2/d ≤ µ1 ≤ µ2 and Lmax ≤L≤ dLmax. The above rates implies that GCD converges faster than RCD (in termsof iterations) since µ1 ≤ µ2. Moreover, these rates demonstrate that GCD could befaster than GD when µ1/Lmax > µ2/L and slower than GD otherwise.3The original proof from Dhillon et al. (2011) used the smoothness parameter with respect to1-norm. However, their proof can be modified to use the coordinate-wise smoothness parameter Lmax.162.3 GCD for sparse optimizationWe consider the composite problemminx∈RdF(x) := f (x)+g(x), (2.1)where f is a strongly convex and L-smooth function, and the regularizer g is afunction that is separable and convex, but not necessarily smooth. Note that thedual formulation of the SVM problem with a bias term has an additional linearconstraint and does not satisfy (2.1) because in this case g is not separable, theformulation in (2.1) only applies to dual SVM without the bias term. For sparsityinducing regularizer including 1-norm regularization and non-negative constraints,it is observed in practice that GCD could converge significantly faster than RCD.However, the refined analysis from Nutini et al. (2015) does not apply to compositeproblems and cannot explain why GCD is faster than RCD in this scenario. Thisgap was filled by Karimireddy et al. (2019), who proved that the boundF(x(t))−F∗ ≤(1− µ1Lmax)dt/2e(F(0)−F∗) (2.2)holds for 1-norm regularization and box-constraints by adding a post-processing stepin each iteration. This result extends Nutini et al.’s refined analysis for compositeproblems and illustrates the theoretical advantage of GCD over RCD.A weakness of the above result is that the rate (2.2) suggests GCD appliedto strongly convex composite problems, with either 1-norm regularization or non-negative constraints, has the same rate as it does for problems without regularizers— i.e., minimizing only f (x) instead of the sum f (x)+g(x). However, GCD appliedto composite problems with sparsity inducing regularization is usually significantlyfaster than its non-regularized counterpart and empirically exhibit screening ability,that is the greedy selection rule can mostly focus on coordinates that are nonzeroat solution. Therefore an improved convergence analysis is needed to explain thisphenomenon. We present our work on GCD for sparse optimization in this sectionto fill this gap between theory and practice.17Algorithm 1 A generic template for GCD1: Input: functions f and gi ∀i ∈ [d].2: W0 = /03: x(0) = 04: for t = 0,1,2, . . . do5: Coordinate selection: select i according to the GS-s rule6: Gradient step: x(t+12 ) = x(t)− (1/Li)∇i f (x(t))ei7: Prox step: x(t+1) = prox(1/Li)gi(x(t+12 ))8: Post-processing; see (2.4)9: Update working set: Wt+1 =Wt ∪{i}; see definition 2.3.110: end for2.3.1 Problem setupFormally, we make the following assumptions on f and g in problem (2.1):• f is Li coordinate-wise smooth ∀i ∈ [d], and we let Lmax = maxi∈[d]Li.• f is L∞-smooth with respect to the ∞-norm.• f is µp strongly convex with respect to the p-norm, p ∈ {1,2}.• g = λ‖ · ‖1 or g = δ≥0, where δ≥0 is the indicator function on the set {xi ≥0 | x ∈ Rd} that vanishes on the nonnegative orthant, and is +∞ otherwise.(It is an open problem if our analysis could be extended to general separableregularizers that are nonsmooth at 0 for all i ∈ [d]).We consider proximal CD algorithm with the GS-s selection rule (Nutini et al.,2015):Selection rule 2.1 (GS-s rule). Select coordinate i ∈ argmaxj∈[d]Q j(x(t)), whereQ j(x) = mins∈∂g j(x j)|∇ j f (x)+ s|. (2.3)The GS-s rule that appeared in Nutini et al. (2015) is a natural extension of thevanilla greedy selection rule to composite problems, it has been widely used in theliterature (Bertsekas, 1999; Li and Osher, 2009; Shevade and Keerthi, 2003; Wuand Lange, 2008). GS-r, GS-q rules are other variants that are commonly being18used, see the discussion from Nutini et al. (2015) and the references therein. Herewe focus on the GS-s rule. The detailed algorithm of GCD for composite problemis shown in Algorithm 1, where the post-processing step is defined asx(t+1)i := 0 if x(t+1)i x(t)i < 0, (2.4)and the proximal operator is defined asproxηgi(x) := argminu∈Rd{12‖u− x‖22+ηgi(ui)}.We also make the following definitions to help our analysis.Definition 2.3.1. The working set Wt is the set of indices selected up to and includ-ing iteration t. Define also W :=⋃∞t=0Wt as the overall working set.Definition 2.3.2. The support of a vector x is the set supp(x) = { i | xi 6= 0}.2.3.2 AnalysisWe now study the screening ability of GCD and develop a bound on the size of theworking set W , i.e., focus on coordinates that are nonzero at solution. We requirethe following quantity, often used in sparsity pattern identification:δi := min{−∇i f (x∗)− `i, ui+∇i f (x∗)} , (2.5)where ∂gi(x∗i ) = [`i,ui]; see Hare and Lewis (2007); Lewis and Wright (2011);Nutini et al. (2017); Sun et al. (2019). The constant δi is closely related to thedistance to the relative interior of the sparse manifold, and is an important quantityin sparse manifold identification (Lewis and Wright, 2011). Optimality conditionsfor (2.1) imply that δi = 0 if x∗i 6= 0, and δi ≥ 0 if x∗i = 0. Because x∗ is unique(by strong convexity), these quantities are problem-specific and algorithmicallyinvariant. The definition in (2.5) leads to the following identification result.Lemma 2.3.1 (Nutini et al., 2017). If for some t > 0 and some coordinate i ∈ [d],x∗i = 0 and x(t) satisfies|∇i f (x(t))−∇i f (x∗)| ≤ δi and x(t)i = 0,19(a) The evolution of the number of inac-tive variables.(b) Objective progress, where ρt is de-fined in (2.8).Figure 2.1: Exploratory investigations.then after one coordinate proximal gradient step x(t+1)i = 0.This lemma suggests that if ∇i f (x(t)) is close to ∇i f (x∗) and x∗i = 0 , then theith entry of x(t) will be correctly identified as 0.Numerical motivationWe present some numerical observations to motivate our analysis. Consider aLASSO problem from random synthetic data. Definef (x) =12‖Ax−b‖22 and g(x) = λ‖x‖1, (2.6)where A ∈ R50×104 and b := Ax]+ ε . The elements Ai j,εi, and non-zeros in thesolution x] are distributed as standard Gaussians. We randomly select 10 elementsfrom x] to be non-zeros and set λ to be 2.Figures 2.1a and 2.1b show the evolution of the number of “inactive” variablesand objective progress for Algorithm 1.In Figure 2.1a, we define#inactive :=d∑i=11{|∇i f (x(t))−∇i f (x∗)| ≤ δi}. (2.7)According to Lemma 2.3.1, this quantity measures how many variables are staying“inactive”, i.e., will not move away from 0 in the next iteration. From Figure 2.1a,20we observe that most variables are initially incorrectly labeled as “active”, i.e.,|∇i f (x(t))−∇i f (x∗)|> δi, but a large number of them quickly switch to “inactive”within first few iterations.In Figure 2.1b, we illustrate the objective progress at each step by plotting ρt ,defined to satisfyF(x(t+1))−F∗ = (1−ρt)(F(x(t))−F∗). (2.8)These experiments illustrate the fact that the initial convergence of GCD, which, forsparse solutions, may be sufficient to quickly identify the few non-zeros. From thisexperiment we observe that• GCD converges fast initially and ∇ f (x(t)) quickly approaches ∇ f (x∗) whenx(t) is still sparse; and• before supp(x(t)) has grown significantly, the coordinates i where x∗i = 0have mostly become inactive, and thus future coordinates that enter W areconstrained to supp(x(t)).We then rigorously characterize these observations. Before proceeding to ourresults, we first introduce some needed concepts.Definition 2.3.3. The function f is µ(τ)p -strongly convex with respect to ‖ · ‖p andsparse vectors if ∀x,y ∈ Rd such that whenever |supp(x)∪ supp(y)| ≤ τ ,f (x)≥ f (y)+ 〈∇ f (y),x− y〉+ µ(τ)p2‖x− y‖2p,where p ∈ {1,2}.The concept of strongly convexity with respect to sparse vectors is not new.For example, see the idea of restricted strongly convex (Negahban and Wainwright,2012) in the literature.21It can be easily verified that µ(τ)1 and µ(τ)2 satisfy the following conditions:µ1 = µ(d)1 ≤ µ(d−1)1 ≤ ·· · ≤ µ(1)1 ,µ2 = µ(d)2 ≤ µ(d−1)2 ≤ ·· · ≤ µ(1)2 ,µ(τ)2 /τ ≤ µ(τ)1 ≤ µ(τ)2 ∀τ ∈ [d]. (2.9)Next, we present a formal analysis to answer why GCD may converge fastinitially, and give a bound on the size of the working set W .Theorem 2.3.1 (Fast initial convergence). Let τ = |supp(x∗)| and let {x(i)}∞i=1 bethe iterates generated by Algorithm 1 with the GS-s rule (selection rule 2.1). Thenfor any t < d− τ ,F(x(t))−F∗ ≤dt/2e∏i=1(1− µ(τ+i−1)1Lmax)(F(0)−F∗) (2.10)≤dt/2e∏i=1(1− µ2(τ+ i−1)Lmax)(F(0)−F∗) . (2.11)Proof is placed in Appendix. The bound in (2.11) follows from (2.9). InTheorem 2.3.1 we show two different bounds, which allow us to draw comparisonsto existing results, below.Bound (2.11). Nesterov (2012) and Richta´rik and Taka´c (2014) established thatRCD exhibits the rateE[F(x(t))−F∗]≤(1− µ2dLmax)(F(0)−F∗) .Compared to (2.11), we see that the dimension d is replaced by the quantity (τ+i−1), which, if τ = |supp(x∗)| is small and we are in the first few iterations, maybe much smaller than d. This reflects the fast initial convergence often observed inpractice; cf. Figure 2.1b.Bound (2.10). Nutini et al. (2015) and Karimireddy et al. (2019) established forthe GCD method the linear convergence rate described by (2.2). Compared to (2.10),22(a) The trend of approximate µ(τ+i)1 ,where τ = 10 in this example.(b) The curve illustrates the upper boundof |W | in Theorem 2.3.2.Figure 2.2: Illustrations for Theorem 2.3.1 and 2.3.2we see that µ1 is replaced with the quantity µ(τ+i−1)1 , which is potentially muchlarger in the early stages (i is small), particularly if τ is small (x∗ is sparse). Thisis confirmed by Figure 2.2a, which shows how in practice this quantity can bemuch larger than µ1. In particular, when τ and i are small, µ2/d µ(τ+i−1)1 andthe convergence rate for GCD is initially significantly faster than RCD, even in theworst case.The rate we derived here is based on two important ingredients: zero initializa-tion and sparse solution. The screening ability of GCD does not hold without eitherof these properties.To better understand the effect of the quantity µ(τ+i−1)1 on the convergencerate, we conduct a simulation using the LASSO problem described in eq. (2.6).The term µ(τ+i−1)1 is hard to compute in general, and thus here we set τ = 10, andfor each i ∈ {1,2, . . . ,30} we generate 103 random (τ+ i−1)-sparse vectors andapproximate µ(τ+i−1)1 as the minimum of ‖Ax‖22/‖x‖21 over the sample vectors. Theplot of approximate µ(τ+i−1)1 against i is shown in Figure 2.2a, and its pattern clearlysupports our previous argument.To develop an upper bound for the size of working set, we define the errormeasurepδ (α) =d∑i=11{α ≤ δi},which we use to quantify the number of inactive elements in the iterates, as ineq. (2.7).23Theorem 2.3.2 (Working set bound). Let τ = |supp(x∗)| and {x(i)}∞i=1 be the iteratesgenerated by Algorithm 1 with the GS-s rule (selection rule 2.1). Then|W | ≤mint∈[d]{Bt + t} , (2.12)whereBt := d+ τ− pδ(L∞ supi≥t{‖x(i)− x∗‖1}).The term L∞ is the smoothness parameter with respect to infinity norm. Thedetailed proof is placed in Appendix. We give a short interpretation to betterunderstand this bound. Note that Bt is a decreasing function of t. Thus, |W | isbounded by the infimum of the sum of a decreasing and increasing function. (SeeFigure 2.2b.) Theorem 2.3.2 implies that if x(t) converges quickly to x∗ (i.e., witht d), then the bound (2.12) will be far less than d.Again, consider the synthetic LASSO problem (A ∈ R50×104 ,λ = 2) as a con-crete example to illustrate this bound. In this example, the curve of Bt + t is shownin Figure 2.2b and the infimum of Bt + t is about 1000 in this case. This experimentdemonstrates that the bound we derived in Theorem 2.3.2 is non-trivial, especiallyfor problems whose d is large.Furthermore, we can use Theorem 2.3.1 and 2.3.2 to derive an alternative boundthat depends only on the constant µ(τ+i)1 , i ∈ [d− τ], instead of the iterates x(i)’s.Corollary 2.3.1. Let τ = |supp(x∗)| and let {x(i)}∞i=1 be the sequence of iteratesgenerated by Algorithm 1 with the GS-s rule (selection rule 2.1). Then Bt inbound (2.12) can be replaced byBt := d+ τ− pδ[2L2∞µ1t−1∏i=0(1− µ(τ+i)1Lmax)R]1/2 ,where R = F(0)−F∗ is the initial objective gap.2.3.3 Improved selection ruleOur analysis in previous section provides an upper bound of |W |. In this section,we propose a variant of the GS-s rule that could favour an even smaller working set.24The resulting algorithm, which we call ∆-GCD, is Algorithm 1 with the followingmodified selection rule.Selection rule 2.2 (∆-GS-s rule). Given the fixed parameter ∆ ∈ (0,1], select coor-dinatei ∈argmaxi∈[d]Qi(x(t)), if ∆maxi∈[d]Qi(x(t))2 ≥maxi∈WtQi(x(t))2argmaxi∈WtQi(x(t)), if ∆maxi∈[d]Qi(x(t))2 <maxi∈WtQi(x(t))2where Wt denotes the set of indices accrued thus far and Qi is defined by (2.3).Note that when ∆= 1, the ∆-GS-s and GS-s rules are equivalent. Intuitively, the∆-GS-s rule, with small ∆, is more likely to focus on the current working set; on theother hand, a large ∆ encourages the algorithm to include unexplored coordinatesand expand the current working set. Thus ∆ controls the trade-off between the sizeof working set and the progress we can make when staying in the current workingset. This is similar to the exploration/exploitation trade-off in the context of onlinelearning (Auer et al., 1995).Theorem 2.3.3. Let {x(i)}∞i=1 be the iterates generated by Algorithm 1 with the∆-GS-s rule (selection rule 2.2) and let W∆ be the final working set. Then for allt > 0,F(x(t))−F∗ ≤(1− ∆µ(|W∆|)1Lmax)dt/2e(F(0)−F∗) (2.13)≤(1− ∆µ2|W∆|Lmax)dt/2e(F(0)−F∗) . (2.14)Theorem 2.3.3 explicitly described the trade-off between the convergence rateand the size of working set. Similar to Theorem 2.3.1, we provide two boundsfor easier interpretation: eq. (2.13) can be viewed as a refinement of the strongconvexity parameter in Karimireddy et al. (2019) and Nutini et al. (2015); andeq. (2.14) where the variable dimension dependency that appears in Nesterov (2012)and Richta´rik and Taka´c (2014) is replaced by the size of the final working set.The ∆-GCD variant is expected to outperform standard GCD when the latter has25a comparatively large working set, and ∆-GCD can reduce the size of working setwith an appropriate value of ∆.To better understand the relationship between W∆ and ∆, we present an descrip-tion of W∆ as ∆→ 0. First, consider the standard GCD algorithm where at eachiteration we additionally minimize the objective over the current working set, i.e.,the next iterate is obtained asx(t+1) := argminsupp(x)⊆Wt+1f (x)+g(x), (2.15)where all variables not in the working set are held fixed at 0. The algorithmterminates when the iterate x(t+1) is optimal for (2.15). The resulting methodis known as the totally corrective greedy algorithm, which is closely related toorthogonal matching pursuit for sparse least squares; see works from Pati et al.(1993); Davis et al. (1997); and Foucart and Rauhut (2013). Note that when∆→ 0, the result selection rule is close to the selection rule used in Algorithm 1from Shevade and Keerthi (2003). We denote the final working set from this schemeas W ].Next, note that as ∆→ 0, the ∆-GS-s selection rule tends to select indicesfrom the current working set, and thus the ∆-GS algorithm converges to a solutionof (2.15). However, when the ∆-GCD iterate is close to the exact minimizer, the∆-GS-s rule must eventually expand the working set. As the following result shows,W∆ converges to W ].Theorem 2.3.4. Let k = |W ]| and let {x(t)}kt=0 be the iterates generated by thetotally corrective greedy algorithm. Assume that argmaxi∈[d]Qi(x(t)) are singletonsfor t = 0,1, . . . ,k−1 and δi > 0 ∀x∗i = 0. Then ∃ ε > 0, such thatW∆ =W ]for any ∆< ε .If the totally corrective greedy algorithm can yield a small working set, then weexpect that a sufficiently small value of ∆ would also yield a small working set (butcould probably slow down the convergence according to Theorem 2.3.3). Hence,26Table 2.1: Properties of the experimental data. Here, d denotes the number offeatures and n denotes the number of samples.Datasets colon leukemia make circle ijcnn1d 2,000 7,129 2 22n 62 72 1,000 35,000our new algorithm ∆-GCD can be viewed as a flexible greedy algorithm between thetwo extreme cases — standard GS-GCD and the totally corrective greedy algorithm.2.3.4 Numerical experimentsWe conduct experiments on both real world data and synthetic data to illustratethe importance of zero initialization and exam the effectiveness of our proposed∆-GCD.The statistics of our experimental data are shown in Table 2.1, where the datasetscolon, leukemia, and ijcnn1 were obtained from the LIBSVM website4 (Chang andLin, 2011). The make circle dataset were generated from the scikit-learn package(Pedregosa et al., 2011). We solve the LASSO problem over the colon and leukemiadatasets, and the dual RBF kernel SVM over the make circle and ijcnn1 datasets.For ijcnn1, we follow the parameter settings described by Hsieh et al. (2014), andthus set γ = 2 and C = 32, where γ is the free parameter in the RBF kernel and 1/Cis the hinge-loss weight parameter. All experiments are conducted on a machinewith 4 CPUs and 16GB memory.CodeThe code to reproduce our experimental results is publicly available at https://github.com/fanghgit/Greed Meets Sparsity.Zero v.s. other initializationsWe compare the convergence of standard GCD for solving LASSO problems overdifferent initialization strategies, we include the following initialization strategies4https://www.csie.ntu.edu.tw/ cjlin/libsvmtools/datasets/27(a) LASSO, data: leukemia, λ = 0.1. (b) LASSO, data: colon, λ = 0.1.Figure 2.3: Comparison between different kinds of initializationinto our comparisons:• Zero initialization: x(0) = 0.• Random initialization: x(0) is generated from Gaussian distributionsN (0,σ Id),for σ ∈ {1,0.1,0.01}.• Least-squares initialization: x(0) = (AT A+ λ Id)−1AT b. This initializationstarts the method at a low objective value, but is not sparse.In Figure 2.3, we can see that zero initialization clearly outperforms otherinitialization strategies, and random initialization tends to be the worst. In particular,GCD with zero initialization is able to get close to a solution even before one pass ofall coordinates. On the other hand, although GCD with least-squares initializationhas a better initial objective value than zero initialization, it suffers from slowconvergence and requires at least a full pass of all coordinates before reaching thesame low error, which is consistent with our intuition. Random initialization withdifferent standard deviations also vary in their performance and random initializationwith smaller variance tends to converge faster; however they are still outperformedby the zero initialization for the same reasons as the least-squares initialization.Evaluation of ∆-GCDWe evaluate the proposed ∆-GCD algorithm on LASSO, 1-norm regularized logisticregression, and kernel SVM problems. As shown in Figure 2.4, the value of ∆28(a) LASSO, data: leukemia, λ = 0.1.Solid line is objective value, dashedline is # non-zeros in x(t).(b) LASSO, data: leukemia, λ = 0.01.Solid line is objective value, dashedline is # non-zeros in x(t).(c) L1-regularized logistic regression,data: colon, λ = 0.1.(d) L1-regularized logistic regression,data: leukemia, λ = 0.1.(e) kernel SVM, data: make circle, C =10,γ = 0.5.(f) kernel SVM, data: ijcnn1, C = 32,γ =2.Figure 2.4: Compare ∆-GCD with different choices of ∆.has a clear impact on the size of the working set, where smaller values of ∆ tendto promote sparser iterates for all the test problems. This trend is more obviouswhen the underlying solution is less sparse, as shown in Figures 2.4b and 2.4f. Thisis because vanilla GCD produces a working set that is much larger than needed.Sometimes this stronger screening ability of a smaller ∆ can lead to slightly faster29convergence compared to standard GCD (i.e., ∆ = 1) as shown in Figures 2.4aand 2.4b. A by-product of ∆-GCD is early identification of the final sparsity pattern,which can be leveraged in two-stage methods (Bertsekas, 1976; Daniilidis et al.,2009; Ko et al., 1994; Wright, 2012). However, the acceleration functionality of∆-GCD is not present for all test problems, since vanilla GCD already has a strongscreening ability for constraining the size of the working set.2.3.5 DiscussionBy bringing techniques from sparsity pattern identification and convergence analysisof GCD, we formally analyze the screening ability of GCD and explicitly answeredwhy GCD is usually fast for sparse optimization. We also propose an improvedselection rule with a stronger ability to encourage sparse iterates and connect toexisting algorithms.For future work, it would be interesting to generalize our analysis and relaxthe strong-convex assumption on the function f . In particular, one could considerproblems where x∗ may not be unique (but supp(x∗) may be). The core of ouranalysis relies on understanding the convergence of the iterates themselves, and notjust the function values. Thus, the challenge in generalizing our analysis to moregeneral smooth objectives requires a different proof technique. Another directionis to tighten the working set bound in Theorem 2.3.2. The bound illustrated inFigure 2.2b is still about 10 times worse than the actual size of the working set,and for small value of λ , the actual working set become larger and our bound canbe trivially larger than d. Our analysis also require the regularizer to be 1-normor nonnegative constraint, it is an open problem to extend our analysis to generalregularizer that are non-smooth at origin.302.4 Gap-based safe-screening rules for atomic-normregularized problemThe safe-screening rules, originally proposed by Ghaoui et al. (2012), generally referto approaches that correctly identify the coordinates that can be safely discardedwithout changing the solution or hindering the optimization process. It has beensuccessfully applied to coordinate optimization algorithms to reduce the overallcomputational effort and achieved promising result for 1-norm regularized problem.Due to the empirical success of safe-screening rule on sparse optimization, variousscreening strategies (Atamtu¯rk and Go´mez, 2020; Bao et al., 2020; Bonnefoy et al.,2015; Kuang et al., 2017; Liu et al., 2014b; Ndiaye et al., 2017; Raj et al., 2015;Wang et al., 2014, 2013; Xiang et al., 2017; Zhang et al., 2017) have been proposedfor difference tasks in recent years. In particular, Ndiaye et al. (2017) proposedan effective and elegant screening framework based on duality gap called the gapsafe-screening rule. The gap-based safe-screening rule has been shown to be usefulfor 1- and group-norm regularized problems.2.4.1 From coordinates to atomsSparse optimization problems are characterized by solutions that are sparse vectors,i.e., few nonzero entries. The idea of vector sparsity can be generalized to “atomic”sparsity. That is, the problem attains a solution that is sparse with respect tosome atomic set A, where A could potentially contain infinite number of elements.Formally, for a given atomic set A⊆ Rd , an optimal solution x∗ can be representedasx∗ = ∑a∈Acaa, ca ≥ 0, (2.16)where only a small number of coefficients ca are nonzero. The archetypal exampleis a sparse vector, which is sparse with respect to the set of signed canonical unit vec-torsA= {±e1, . . . ,±ed}. 1-norm regularization is the standard approach to producesparse solution. The atoms that participate nontrivially in the decomposition (2.16)represent latent structure in the solution. The notion of atomic sparsity is prevalentin machine learning (Argyriou et al., 2007; Meinshausen and Bu¨hlmann, 2006;Tibshirani, 1996; Yuan and Lin, 2006) and signal processing (Candes et al., 2015),31and has been formalized in the context of inverse problems by Chandrasekaranet al. (2012).The gap-based safe-screening rules have been shown to be effective for 1-normregularized problems. In the following sections, we will extend the gap-based safe-screening rule to general atomic sets and investigate whether it is possible to usethe generalized gap-based safe-screening rule to save computation for nuclear-normregularized problems.2.4.2 Some technical toolsWe introduce in this section the basic tools of convex analysis and atomic sparsitythat serve as the cornerstone of our analysis, most of the material for sections 2.4.1to 2.4.3 is a summary of the framework described by Fan et al. (2020). We makethe blanket assumption that the atomic set A⊆ Rd is compact, and that the originis contained in its convex hull (we do not assume that A is convex). The gaugefunction to the set A measure the magnitude of a function relative to that set.Definition 2.4.1 (Gauge function). The gauge function with respect to A is definedasγA(x) = inf{∑a∈Aca∣∣∣∣∣ x = ∑a∈Acaa, ca ≥ 0,∀a ∈ A}. (2.17)The gauge function is always convex, nonnegative, and positively homogeneous.However, it is not necessarily a norm because it may not be symmetric (unless Ais centrosymmetric) and may vanish or take infinite value at points off the origin(unless the origin is in the relative interior of A). Definition 2.4.1 makes explicitthe role of a gauge function as a convex penalty for atomic sparsity. The supportof a vector x describes the atoms that contribute positively in the decompositiondescribed by (2.17).Definition 2.4.2 (Atomic support). The atomic support for a point x ∈ Rd withrespect to the set A is defined to be the set SA(x) that satisfiesγA(x) =∑a∈SA(x)ca, x =∑a∈SA(x)caa, and ca > 0 ∀a ∈ SA(x).The atomic set of signed 1-hot unit vectors A = {±ei | i = 1,2, . . . ,d}, for32Atomic sparsity A γA(x) SA(x) σA(z)non-negative cone({e1, . . . ,ed }) δ≥0 cone({ei | xi > 0}) δ≤0element-wise {±e1, . . . ,±ed } ‖ · ‖1 {sign(xi)ei | xi 6= 0} ‖ · ‖∞low rank {uvT | ‖u‖2 = ‖v‖2 = 1} nuclear-norm singular vectors of x spectral normPSD & low rank {uuT | ‖u‖2 = 1} tr+δ0 eigenvectors of x max{λmax,0}Table 2.2: Commonly used sets atom sets and the corresponding gauge and supportfunctions. The indicator function δC(x) is zero if x is in the set C and+∞ otherwise. The commonly used group-norm is also an atomic norm;see Fan et al. (2020, Example 4.7).example, the support SA(x) coincides with the nonzero elements of x with thecorresponding sign. The support function, defined below, is dual to the gaugefunction, and provides a key tool for identifying atoms associated with the supportof a vector.Definition 2.4.3 (Exposed faces and ε-exposed faces). The exposed face and ε-exposed face for a point z ∈ Rd with respect to the set A is defined byFA(z) = {a ∈ A | 〈a,z〉= σA(z)} , FA(z,ε) = {a ∈ A | 〈a,z〉 ≥ σA(z)− ε }(2.18)where σA(z) = supa∈A〈a,z〉 is the support function with respect to A.Note that when ε = 0, the ε-exposed face coincides with the exposed face. Welist in Table 2.2 some commonly used atomic sets, their corresponding gauge andsupport functions, and atomic supports.2.4.3 Problem setupWe consider the gauge regularized problemminimizex∈Rdf (Mx) subject to γA(x)≤ τ, (P)where f is an L-smooth and convex function, M : Rd → Rm is a linear operatorand τ > 0 is a given constant can control the sparsity level. The gauge functionconstraint in (P) can generate a solution that is sparse with respect to the atomic setA, and a wide range of application including LASSO and matrix completion can beformulated as (P) (with an appropriate choice of the atomic set).33The dual problem of (P) isminimizey∈Rmd(y) := f ∗(y)+ τσA(M∗y), (D)where f ∗(y) = supw∈Rm 〈y,w〉− f (w) is the convex conjugate function of f , andM∗ : Rm→ Rd is the adjoint operator of M, which satisfies 〈Mx,y〉= 〈x,M∗y〉 forall x ∈ Rd and y ∈ Rm. The dual problem plays an important role in safe-screeningrules. Most safe-screening rules including the very first one proposed by Ghaouiet al. (2012) rely on dual variables to discard redundant features.2.4.4 The gap-based safe-screening ruleNext, we extend the gap-based safe-screening rule from Ndiaye et al. (2017) togeneral atomic sets with the language of atomic sparsity. Note that a similar analysiscan be found in a recent work from Sun and Bach (2020).Given a feasible solution x of (P), safe-screening rule aims to infer the atomswith nonzero coefficients at solution based on x. We know that SA(x) could bedifferent from SA(x∗) even for feasible solution x that is arbitrarily close to x∗.Fortunately, the following support-face relationship from Fan et al. (2020) allow usto constrain SA(x∗) by constructing a superset of FA(M∗y∗) from the dual problem.Lemma 2.4.1 (Fan et al., 2020, Proposition 4.5 and Theorem 5.1). Let x∗ and y∗ beoptimal primal-dual solutions for problem (P) and (D). ThenSA(x∗)⊆FA(M∗y∗).Define the atomic operator norm by ‖M‖A := maxa∈A ‖Ma‖2, we can get thefollowing result based on Lemma 2.4.1.Theorem 2.4.1 (Generalized gap-based safe-screening rule). Given a primal feasi-ble point x and a dual feasible point y, we haveSA(x∗)⊆FA(M∗y,ε), (2.19)where ε = 2‖M‖A√2L( f (Mx)+d(y)), and f (Mx)+d(y) is the duality gap.34Proof of Theorem 2.4.1. We require the following standard technical tool for theproof of Theorem 2.4.1Lemma 2.4.2 (Kakade et al., 2009, Theorem 6). If f is L-smooth, then f ∗ is 1/L-strongly convex.Let y∗ denote the optimal dual variable for (D). We show that ‖y− y∗‖2 canbe bounded by the duality gap. By Lemma 2.4.2, we know that f ∗ is 1/L-stronglyconvex, and it follows that d(y) is also 1/L-strongly convex. Then by the definitionof strongly convexity, we have∀s ∈ ∂d(y∗), d(y)≥ d(y∗)+ 〈s,y− y∗〉+ 12L‖y− y∗‖22.By the optimality condition, we know that 0 ∈ ∂d(y∗). Therefore, by reordering theinequality, we can get‖y− y∗‖2 ≤√2L(d(y)−d(y∗))≤√2L( f (Mx)+d(y)) ∀x s.t. γA(x)≤ τ. (2.20)Next, we show that FA(M∗y∗)⊆FA(M∗y,ε). For any a ∈ FA(M∗y∗),〈a,M∗y〉= 〈a,M∗y∗〉+ 〈a,M∗y−M∗y∗〉≥ σA(M∗y∗)−(maxa∈A‖Ma‖)‖y− y∗‖2= σA(M∗y)− (σA(M∗y)−σA(M∗y∗))−‖M‖A‖y− y∗‖2≥ σA(M∗y)−2‖M‖A‖y− y∗‖2≥ σA(M∗y)−2‖M‖A√2L( f (Mx)+d(y))≥ σA(M∗y)− ε.The last line is true by the definition of ε . Combining above result with Lemma 2.4.1,the proof is finished.The proof of Theorem 2.4.1 mirrors the original proof of the gap-based safe-screening rule for 1-norm regularization (Ndiaye et al., 2017) — redundant atoms(or coordinates) can be identified by bounding ‖y− y∗‖ with the duality gap. By35using the support-face relationship described in Lemma 2.4.1 instead of the first-order optimality condition with 1-norm regularization, we can generalize their proofto atomic sets.Next we characterize the atom identification property of the generalized gap-based safe-screening rule.Proposition 2.4.1 (Atomic identification). Let {x(t)}∞t=1 and {y(t)}∞t=1 be sequencesthat converge to optimal primal and dual solutions x∗ and y∗ respectively. Let{ε(t)}∞t=1 be the gaps defined in Theorem 2.4.1 evaluated at x(t) and y(t). Then the setA(t) := ∩tj=1FA(M∗y( j),ε( j)) has the Painlevee´-Kuratowski set limit (Rockafellarand Wets, 2009, p. 111)limt→∞A(t) = FA(M∗y∗). (2.21)Proposition 2.4.1 ensures that the safe-screening rule (2.19) is guaranteed toeventually discard superfluous atoms as long as we have available an iterativesolver that can generate primal iterates that converge to a solution. For polyhedralatomic set, e.g., atomic set with finite elements, it is straightforward to verify thatProposition 2.4.1 implies the following finite-time atom identification property:∃ T > 0 such that A(t) = FA(M∗y∗) ∀t > T.The implementation of the generalized gap-based safe-screening rule for poly-hedral atomic sets is also straightforward. One can store all atoms in memory duringcomputation, and the gap-based safe-screening rule offers a computable way todiscard redundant atoms periodically during the optimization. When FA(M∗y,ε)is small enough, let Â := FA(M∗y,ε) := {aˆi}ri=1, we can solve the reduced low-dimensional problemminimizex∈Rr,x≥0f(Mr∑i=1aˆixi)subject tor∑i=1xi ≤ τinstead of the original high-dimensional problem efficiently using algorithm such asaccelerated projected gradient descent.A remarkable aspect of the generalized gap-based safe-screening rule (and also36the original gap safe-screening rule) is that it depends solely on the duality gap, andhence is algorithm agnostic. As long as we have an algorithm that guarantees dualitygap converges to 0, the gap-based safe-screening rule will recover FA(M∗y∗) withina finite number of iterations (for finite atomic sets). The gap-based safe-screeningrule has been successfully applied to algorithms such as conditional gradient descentand projected coordinate descent to achieve promising performance.2.4.5 Gap-based safe-screening rule for nuclear normWe explore a key question whether the generalized safe-screening rule can provideany computational advantage for atomic sets A with infinite number of atoms. Inparticular we consider the nuclear-norm regularized problems whose atomic set isthe set of rank-one matrices, i.e.,A= {uvT | u ∈ Rn,v ∈ Rm,‖u‖2 = ‖v‖2 = 1} .In the following proposition, we show that the ε-exposed face of M∗y contains allthe singular vectors of M∗y when ε is strictly positive.Proposition 2.4.2 (Limitation of ε-Face). Let M∗y=UΣV T be the full SVD decom-position of M∗y, where Σ= diag(σ1,σ2, . . . ,σmin{n,m}),σi ≥ 0 ∀ 0< i≤min{n,m}.For any ε ≥ 0, the ε-exposed face can be explicitly expressed asFA(M∗y,ε) ={U p(V q)T∣∣∣∣∣ min{n,m}∑i=1 σi piqi ≥ σ1− ε,‖p‖2 = ‖q‖2 = 1}.Then for any ε > 0, there exist p,q with all entries being nonzero, such thatU p(V q)T ∈ FA(M∗y,ε).Proof of Proposition 2.4.2. By the definition of FA(M∗y,ε),FA(M∗y,ε) = {uvT | 〈uvT ,M∗y〉 ≥ σ1− ε, ‖u‖2 = ‖v‖2 = 1}= {uvT | 〈uvT ,UΣV T 〉 ≥ σ1− ε, ‖u‖2 = ‖v‖2 = 1} .We know that U,V are orthonormal matrices, by setting u =U p and v = V p for37some p,q ∈ Rmin{n,m} such that ‖p‖2 = ‖q‖2 = 1, we obtainFA(M∗y,ε) = {U p(V q)T | 〈U p(V q)T ,UΣV T 〉 ≥ σ1− ε, ‖p‖2 = ‖q‖2 = 1}= {U p(V q)T | pTΣq≥ σ1− ε, ‖p‖2 = ‖q‖2 = 1}={U p(V q)T∣∣∣∣∣ min{n,m}∑i=1 σi piqi ≥ σ1− ε,‖p‖2 = ‖q‖2 = 1}.The above finishes the proof.Proposition 2.4.2 indicates that FA(M∗y,ε) contains not only the top singularvectors of M∗y but also the bottom singular vectors—even if ε is arbitrarily closeto 0. This result is unfortunate since the gap-based safe-screening rule stated inTheorem 2.4.1 does not allow us to discard any singular vectors of M∗y and thusrequire a full SVD decomposition of M∗y even if the duality gap is arbitrarily closeto zero.2.4.6 Approximation with partial SVDThe face of A exposed by the vector M∗y∗ is given byFA(M∗y∗) ={uvT | uT (M∗y∗)v = σ1(M∗y∗)},where σ1(M∗y∗) is the largest singular value of M∗y∗. Therefore, when there arefew singular vectors associated with the largest singular value, only the top fewsingular vectors of M∗y∗ are actual useful atoms. This property motivates us touse the partial SVD decomposition of M∗y to extract the reduced atomic set. Thishard-thresholding technique has been widely used as a heuristic. Formally, given adual feasible solution y with partial SVD decompositionM∗y =UrΣrV Tr Ur ∈ Rn×r,Vr ∈ Rm×r,r min{n,m},we construct the corresponding reduced atomic setÂ= {Ur pqTV Tr | ‖p‖2 = ‖q‖2 = 1} ,38and solve the reduced problem over Â.First, we give a concrete example showing that the partial SVD of M∗y is notable to give us a safe cover of FA(M∗y∗) even when FA(M∗y∗) is a singleton and yarbitrarily close to y∗.Example 2.4.1 (Limitation of partial SVD). Consider the problemminimizeX∈Rn×n12‖X−Z‖2F subject to ‖X‖∗ ≤ 1, (2.22)whereZ =U diag(2,0.1, . . . ,0.1)V T and U =V =√1− ε 0 . . . −√ε0 1 . . ..... . .√ε 0√1− εn×nfor some ε ∈ (0,1). The dual problem isminimizeY∈Rn×n12‖Y −Z‖2F −12‖Z‖2F +‖Y‖2. (2.23)The solution for problem (2.22) and problem (2.23) areX∗ =U diag(1,0, . . . ,0)V T and Y ∗ = Z−X∗ =U diag(1,0.1, . . . ,0.1)V T .Let U := [u1,u2, . . . ,un],V := [v1,v2, . . . ,vn], then obviously SA(X∗) = FA(Y ∗) =u1vT1 is a singleton. We construct the following dual feasible solutionŶ = diag(1,0.1, . . . ,0.1).Let Û ,V̂ be the singular vectors of Ŷ , then Û = V̂ = [e1,e2, . . . ,en], where ei is then-dimensional vector with 1 at the ith entry and 0 at other entries. In order tocover u1 = [√1− ε,0 . . . ,√ε]T , we need both the top singular and bottom singularvectors e1 and en, and therefore any top-r SVD decomposition of Ŷ with r < n willend up to be “unsafe”. It is also easy to verify that ‖Ŷ −Y ∗‖F = O(√ε). Notethat our argument holds for any ε ∈ (0,1). Therefore, ∀ε ∈ (0,1), ∃y ∈ B(y∗,ε)39such that only full SVD decomposition of M∗y can guarantee a safe coverage ofFA(M∗y∗).This result shows that the screening rule with partial SVD is not safe. Therefore,we can only resort to an approximate screening rule. We use the one-sided Hausdorffdistanceρ(A1,A2) := supa1∈A1infa2∈A2‖a1−a2‖Fto measure the similarity between any two subsets A1 and A2 of the atomic set A.The next result shows that there is a set Â that is close to SA(x∗), then there mustexist a point in Â that is close to x∗.Proposition 2.4.3 (Hausdorff error bound). Given Â ⊆ A, there exist x ∈ cone(Â)such that‖x− x∗‖F ≤ ρ(SA(x∗),Â)√|SA(x∗)|‖x∗‖F .Proof. Proof of Proposition 2.4.3 Let x∗ = ∑a∈SA(x∗) caa,ca > 0. By the definitionof the one-sided Hausdorff distance ρ , for any a ∈ SA(x∗), there exist a correspond-ing aˆ ∈ Â such that‖aˆ−a‖F ≤ ρ(SA(x∗),Â).Let xˆ = ∑a∈SA(x∗) caaˆ, then it is easy to verity that xˆ ∈ cone(Â) and‖x− x∗‖F ≤ ρ(SA(x∗),Â) ∑a∈SA(x∗)ca(i)≤ ρ(SA(x∗),Â)√|SA(x∗)|‖x∗‖F ,(i) is true since the decomposition x∗ = ∑a∈SA(x∗) caa,ca > 0 is an orthonormaldecomposition and ‖x∗‖2F =∑c2a when our atomic-set is the set of rank-one matrices.Next, we study the approximation ability of the partial SVD decomposition of agiven feasible dual solution M∗y to FA(M∗y∗).Theorem 2.4.2. Let y be a dual feasible vector. Let M∗y=UrΣrV Tr ,Ur ∈Rn×r,Vr ∈Rm×r be the top-r SVD decomposition where r <min{n,m}. Denote {σi}min{n,m}i=1as the singular values and Â= {Ur pqTV Tr | ‖p‖2 = ‖q‖2 = 1} be the our reduced40atomic set, assume σ1 > σr+1, thenρ(FA(M∗y∗),Â)≤ ρ(FA(M∗y,ε),Â) =√2min{εσ1−σr+1 ,1},where ε is the same as defined in Theorem 2.4.1.Proof. Proof for Theorem 2.4.2 By the definition of ρ(·, ·), it is straightforward thatρ(A,C)≤ ρ(B,C) ∀A,B,C ⊆ Rn×m such that A⊆ B.We know that FA(M∗y∗) ⊆ FA(M∗y,ε) from Theorem 2.4.1, then obviously wehaveρ(FA(M∗y∗),Â)≤ ρ(FA(M∗y,ε),Â).For any A1,A2 ⊆A,ρ(A1,A2) =√supa1∈A1infa2∈A2‖a1−a2‖2F =√√√√2−2( infa1∈A1supa2∈A2〈a1,a2〉), (2.24)where the second equality holds since ‖a1‖F = ‖a2‖F = 1 by the definition of A.Define A1 = FA(M∗y,ε) and A2 = Â ={Ur pqTV Tr∣∣ ‖p‖2 = ‖q‖2 = 1}, whereUr,Vr are the top-r singular vectors of M∗y. Let k := min{n,m}, C1 = {(p,q) |∑ki=1σi piqi ≥ σ1 − ε, ‖p‖2 = ‖q‖2 = 1, p,q ∈ Rk} and C2 = {(pˆ, qˆ) | ‖ pˆ‖2 =‖qˆ‖2 = 1, pˆ, qˆ ∈ Rr}, thenρ(A1,A2) =√2−2(minp,q∈C1maxpˆ,qˆ∈C2〈U pqTV T ,Ur pˆqˆTV Tr 〉)=√√√√2−2( minp,q∈C1maxpˆ,qˆ∈C2(r∑i=1pi pˆi)(r∑i=1qiqˆi))=√2−2(minp,q∈C1‖p1:r‖2‖q1:r‖2)(2.25)41Now we consider the subproblem in (2.25):minp,q‖p1:r‖2‖q1:r‖2 (P1)subject tok∑i=1σi piqi ≥ σ1− ε, ‖p‖2 = ‖q‖2 = 1, p,q ∈ Rk.If p∗ and q∗ is a solution of the problem (P1), then it is easy to verify thatp˜ =[‖p∗1:r‖2,0, . . . ,‖p∗r+1:k‖2,0, . . . ,0]and q˜ =[‖q∗1:r‖2,0, . . . ,‖q∗r+1:k‖2,0, . . . ,0]is also a valid solution. Therefore there must exist solution p∗,q∗ such that pi =qi = 0 ∀i /∈ {1,r+1}, that is only p∗1,q∗1 and p∗r+1 and q∗r+1 are greater or equal than0. This allow us to further reduce the problem tominp1,q1,pr+1,qr+1p1q1subject to σ1 p1q1+σr+1 pr+1qr+1 ≥ σ1− ε,p21+ p2r+1 = q21+q2r+1 = 1, p1,q1, pr+1,qr+1 ≥ 0.It is easy to verify that when σ1−σr+1 ≥ ε , the above problem attains solution atp1 = q1 =√σ1−σr+1− εσ1−σr+1 and pr+1 = qr+1 =√1− p21.When σ1−σr+1 < ε , the solution is simply p1 = q1 = 0, pr+1 = qr+1 = 1. Thereforethe optimal value of (P1) is max{1−ε/(σ1−σr+1),0}, plug this into eq. (2.25) andthe proof is finished.Oustry developed a related result based on the two-sided Hausdorff distance (Ous-try, 2000, Theorem 2.11). Directly applying his Theorem to our context wouldend up with a bound O(√ε/(σr−σr+1)), which is looser than the bound shown inTheorem 2.4.2 because σ1 ≥ σr ≥ σr+1.422.4.7 DiscussionOur extension of gap-based safe-screening rules to the various forms of atomic-norm regularization is based on the convex calculus of sublinear functions. Ourproposed screening rules can provide practical computational advantages when theatomic sets are polyhedral. As demonstrated by Example 2.4.1, however, thereare limitations of the rule when used for non-polyhedral atomic sets. In that case,Theorem 2.4.2 provides an error bound based on the truncated SVD.Further research opportunities remain, particularly for designing meaningfulsafe-screening rules for non-polyhedral sets. For example, it seems possible todesign safe-screening rules for nuclear-norm regularized problems that are particularto the search directions generated by the conditional-gradient method.43Chapter 3Online mirror descent withunknown time horizonIn modern big data applications such as recommendation and advertisement, datausually comes in a stream and one is required to make decisions in an online manner.From the perspective of optimization, such problems are usually formulated underthe online convex optimization (OCO) framework. In OCO, a player is required tomake a sequence of online decisions over discrete time steps. Each decision incursa cost given by a convex function that is only revealed to the player before theymake that decision. The goal of the player is to minimize the total cost.Formally, let T denote the number of decisions required to make. For each timestep t ∈ {1,2, . . . ,T}, our algorithm proposes a point x(t) from a closed convex setX ⊆ Rd , and an adversary simultaneously picks a convex cost function ft . Thisfunction penalizes the proposal x(t) by the amount ft(x(t)), and the cost of theiteration at time t is defined as ft(x(t)). The goal of the algorithm is to minimize theregret against an unknown comparison point z ∈ X at time TRegret(T,z) :=T∑t=1ft(x(t))−T∑t=1ft(z).The Regret(T,z) measures the difference between the total cost of the algorithmand the cost of the competitor z up to time T . Our goal is to develop algorithm thatensures its regret is sublinear in T against any competitor. In this way, the average44cost incurred by the algorithm will be guaranteed to be smaller or equal to the costfrom the best competitor as T → ∞.Online mirror descent (OMD) and dual averaging (DA) are two important algo-rithm templates for OCO from which many classical online learning algorithms canbe derived as special cases; see works from Shalev-Shwartz (2012) and McMahan(2017) for examples. For Lipschitz continuous functions { ft}t≥1, when the numberof decisions to be made T is known in advance, the algorithm could use T as aninput and the performance of OMD and DA (with properly chosen stepsize) areshown to be very similar (Hazan, 2016). That is, they achieve essentially the sameregret bound when using the same stepsize scheduling. However, in the morechallenging setting when T is not known a priori (unknown time horizon), thereis a fundamental difference in the regret rates of OMD and DA with a similartime-varying stepsize scheduling. While DA can guarantee sublinear regret boundO(√T )1 for any T > 0 (Nesterov, 2009),there are instances such that OMD couldsuffer asymptotically linear regret, i.e., Ω(T ) (Orabona and Pa´l, 2018).In this chapter, we introduce a stabilization technique to fix OMD in the un-known time horizon setting and give essentially the same regret bound for stabilized-OMD and DA. We present our convergence analysis in a careful and modular waythat allows for straightforward and flexible proofs. We also adapt our stabilized-OMD for composite objective setting.3.1 BackgroundWe review some standard technical tools used by OMD and briefly describe someknown properties of OMD and DA in the literature of OCO.3.1.1 Definitions and notationsBoth OMD and DA are parameterized by a special convex functionΦ, often referredas a regularizer or a mirror map (for X ), which among other properties needs2 tobe of Legendre type (Rockafellar, 1970, Chapter 26). Formally, throughout the1This is rate equivalent to the iteration complexity O(ε−2). To keep our notation consistent withthe literature, we describe the rate in terms of T instead of ε in this chapter.2One may relax this condition in some cases. For a detailed discussion on the conditions neededon the mirror map, see (Bubeck, 2011, § 5.2).45chapter we assume that the function Φ : D¯ → R is a closed convex function suchthat intD¯ ∩ riX 6= /0 (where riX denotes the relative interior of X ), and whoseconjugate is differentiable on Rd . Moreover, we also assume that Φ is of Legendretype, which means that Φ is strictly convex on its domain3 and essentially smooth,that is, for D := intD¯ we have• D is nonempty,• Φ is differentiable on D, and• limx→∂D ‖∇Φ(x)‖=+∞, where ∂D is the boundary of D, i.e., ∂D := clD\D.The gradient of the mirror map ∇Φ : D → Rd and the gradient of its conjugate∇Φ∗ : Rd →D are mutually inverse bijections between the primal space D and thedual space Rd . We will adopt the following notational convention. Any vector inthe primal space will be written without a hat, such as x ∈D. The same letter with ahat, namely xˆ, will denote the corresponding dual vector:xˆ := ∇Φ(x) and x := ∇Φ∗(xˆ) for all letters x .Essential smoothness ensures not only that Φ is differentiable on the interior ofits domain, but also that the slope of Φ increases arbitrarily fast near the boundaryof its domain. The latter guarantees, at least intuitively, that the function to beincreasing near and in the direction of the boundary of its domain. This propertyis fundamental for mirror descent to be well-defined (although not essential fordual averaging) since it ensures that the Bregman Projection onto X is attained bya point on D where Φ is differentiable, and uniqueness is a consequence of thestrict convexity of Φ. Some mirror maps we shall look are classical cases of theOCO literature such as the negative entropy x ∈ Rd+ 7→ ∑di=1 xi lnxi and the squared2-norm 12‖ · ‖22, and details on the reasons they are mirror maps can be found inthe literature (Bubeck, 2011, 2015; Shalev-Shwartz, 2012). In particular, Bubeck(2011, Section 5.2) discussed the properties of functions of Legendre type and why3In fact we only need Φ to be strictly convex on some convex subsets of the domain (Rockafellar,1970, Chapter 26), but for the sake of simplicity we assume that Φ is strictly convex on its entiredomain.46Algorithm 2 Pseudocode for both online mirror descent and dual averaging withdynamic stepsize given by ηt on iteration t. These methods differ only in how theiterate yˆ(t+1) is updated.Input: x(1) ∈ X ∩D,η : N→ R>0.for t = 1,2, . . . doIncur cost ft(x(t)) and receive gt ∈ ∂ ft(x(t))xˆ(t) = ∇Φ(x(t))[OMD update] yˆ(t+1) = xˆ(t)−ηtgt[DA update] yˆ(t+1) = xˆ(1)−ηt ∑i≤t giy(t+1) = ∇Φ∗(yˆ(t+1))x(t+1) =ΠΦX (y(t+1))end forrequiring the conjugate of the mirror map to be differentiable on the whole space isnot necessary for mirror descent to be well-defined if one restricts the gradient stepsin the dual space in some ways.Given a mirror map Φ, the Bregman divergence with respect to Φ is defined byDΦ(x,y) := Φ(x)−Φ(y)−〈∇Φ(y), x− y〉, ∀x ∈ D¯,∀y ∈ D. (3.1)Throughout this chapter it will be convenient to use the notationDΦ(ab ;c) := DΦ(a,c)−DΦ(b,c) = Φ(a)−Φ(b)−〈∇Φ(c), a−b〉. (3.2)In the important special case where Φ(x) = 12‖x‖22, the Bregman divergence relatesto the Euclidean distance, i.e., DΦ(x,y) = 12‖x− y‖22. When D = Rd+ and Φ(x) =∑di=1 xi logxi, the Bregman divergence becomes the generalized Kullback-Leibler(KL) divergence. The projection operator induced by the Bregman divergence iswritten as ΠΦX (y) := argmin{DΦ(x,y) | x ∈ X}.A general template for optimization in the mirror descent framework is shown inAlgorithm 2. The two classical algorithms, online mirror descent and dual averaging,are incarnations of this, differing only in how the dual variable yˆ(t) is updated.473.1.2 OMD and DA with constant stepsizeWhen the time horizon T is known in advance, constant stepsize that depends onT can be adopted to achieve sublinear regret. In particular, as characterized by thefollowing Theorem, OMD and DA with the same stepsize can obtain exactly thesame regret bound.Theorem 3.1.1 (Nesterov, 2009, Theorem 1, Hazan, 2016, Theorem 5.6). Supposethat Φ is ρ-strongly convex with respect to a norm ‖ · ‖ and pick a constant stepsizeηt := η > 0 for all t ≥ 1. Let {x(t)}t≥1 be the sequence of iterates generatedby Algorithm 2. Then for any sequence of convex functions { ft}t≥1 with eachft : X → R, the following bound holds for both OMD and DA updates,Regret(T,z) ≤T∑t=1η‖gt‖2∗2ρ+DΦ(z,x(1))η. (3.3)When T is known as a priori, O(√T ) regret can be obtained by setting η =1/√T in eq. (3.7). Interestingly, though OMD and DA with constant stepsize havesame regret bound, the proofs used to derive this bound tend to be quite different.3.1.3 OMD and DA with dynamic stepsizeIn the unknown time horizon scenario, a dynamic stepsize with ηt ∝ 1/√t is usuallyadopted in the literature of online learning (Beck and Teboulle, 2003; Zinkevich,2003). Moreover, when the domain X is bounded, both OMD and DA with stepsizeηt ∝ 1/√t have O(√T ) regret bounds (with differing constants). However, whenwe allow the domain X to be unbounded, OMD is shown to be to provably worsethan DA:Theorem 3.1.2 (Linear regret for OMD, Orabona and Pa´l, 2018, Theorem 3). Setηt = 1/√t. Let {x(t)}t≥1 denote the sequence of iterates generated by Algorithm 2with OMD update. For any T ≥ 3 there exists a sequence of convex 1-Lipschitzfunctions { ft}Tt=1 and an initial point x(1) ∈ X such thatsupz∈XDΦ(z,x(1)) is bounded and Regret(T,z) = Ω(T ),48while Algorithm 2 with DA update can always guarantee sublinear regret boundO(√T ) using a similar stepsizes (which differ only by constants).Moreover, there are examples showing that for offline 1-dimensional gradientdescent (i.e., mirror descent with 2 norm square regularization), a stepsize that iseither asymptotically o(1/√t) or ω(1/√t) cannot achieve regret bound O(√t) forall t > 0. So OMD with stepsizes of the form t−α with α > 0 cannot obtain optimalregret when X is unbounded. A natural question is if we can improve OMD to makeit provably work with dynamic stepsizes. In the next section we provide a fix forOMD with dynamic stepsizes through a stabilization technique and later we showits connection with dynamic DA.3.2 Stabilized OMDThe intuition for the idea is as follows. Suppose Z ⊆X is a set of comparison pointswith respect to which we wish our algorithm to have low regret. Usually, we assumesupz∈Z DΦ(z,x(1)) is bounded, that is, the initial point is not too far (with respect tothe Bregman divergence) from any comparison point. Since supz∈Z DΦ(z,x(1)) isbounded (but not necessarily supz∈Z,x∈X DΦ(z,x)), the point x(1) is the only pointin X that is known to be somewhat close (with respect to the Bregman divergence)to all the other points in X . Thus, iterates computed by the algorithm should remainreasonably close to x(1) so that no other point z ∈ Z is too far from the iterates. Ifthere were such a point z, an adversary could later chose functions so that pickingz every round would incur low loss. At the same time, OMD would take manyiterations to converge to z since consecutive OMD iterates tend to be close withrespect to the Bregman divergence. That is, the algorithm would have high regretagainst z. To prevent this, the stabilization technique modifies each iterate x(t) tomix in a small fraction of x(1). This idea is not entirely new: it appears, for example,in the original Exp3 algorithm (Auer et al., 2002a), although for different reasons.There are two ways to realize the stabilization idea.Primal Stabilization. Replace x(t) with a convex combination of x(t) and x(1).Dual Stabilization. Replace yˆ(t) with a convex combination of yˆ(t) and xˆ(1)(Recall from Algorithm 2 that yˆ(t) is the dual iterate computed by taking a gradientstep). An illustration for dual stabilization is shown in Figure 3.1.49PrimalDual}Figure 3.1: Illustration of the t-th iteration of DS-OMD.After a draft of this chapter was made publicly available (the preliminary versionreleased in ICML 2020), we were informed that an idea similar to primal stabiliza-tion had appeared in the Robust Optimistic Mirror Descent algorithm (Kangarshahiet al., 2018). Their setting is somewhat different since they perform optimistic steps.Furthermore, their results are somewhat weaker in terms of constant factors andsince they cannot handle Bregman projections.Additionally, after the ICML version of this paper was published we were toldabout ideas similar to primal stabilization were involved in the Twisted MirrorDescent(TMD) algorithm (Gyo¨rgy and Szepesva´ri, 2016). More specifically, TMDis a meta-algorithm adds a step at each iteration controlled by some sequence offunctions that may depend on previous iterates and has, as a special case, primalstabilization. For the case of prediction with expert advice, they showed (Gyo¨rgyand Szepesva´ri, 2016, Example 6) how primal stabilization yields good bounds onthe shifting regret. In our work, we extended the idea of primal stabilization forcases beyond the one of prediction with expert’s advice.3.2.1 Dual-stabilized OMDAlgorithm 3 gives pseudocode showing our modification of OMD to incorporatedual stabilization.Theorem 3.2.1 (Regret bound for dual-stabilized OMD). Assume that ηt ≥ηt+1 > 050Algorithm 3 Dual-stabilized online mirror descent, with dynamic stepsize ηt . Theparameters γt control the amount of stabilization.Input: x(1) ∈ X , η : N→ R+, γ : N→ (0,1]for t = 1,2, . . . doIncur cost ft(x(t)) and receive gt ∈ ∂ ft(x(t))xˆ(t) = ∇Φ(x(t)) B map primal iterate to dual spacewˆ(t+1) = xˆ(t)−ηtgt B gradient step in dual space (3.4)yˆ(t+1) = γtwˆ(t+1)+(1− γt)xˆ(1) B stabilization in dual space (3.5)y(t+1) = ∇Φ∗(yˆ(t+1)) B map dual iterate to primal spacex(t+1) =ΠΦX (y(t+1)) B project onto feasible region (3.6)end forfor all t ≥ 1. Define γt = ηt+1/ηt ∈ (0,1] for all t ≥ 1. Let {x(t)}t≥1 be the sequenceof iterates generated by Algorithm 3. Then for any sequence of convex functions{ ft}t≥1 with ft : X → R for each t ≥ 1,Regret(T,z) ≤T∑t=1DΦ( x(t)x(t+1) ;w(t+1))ηt+DΦ(z,x(1))ηT+1∀T > 0. (3.7)Note that strong convexity of Φ is not assumed. As we will see in Section 3.3.1,the term DΦ( x(t)x(t+1) ;w(t+1)) can be easily bounded when the mirror map is stronglyconvex. This yields sublinear regret for ηt ∝ 1/√t, which is not the case for theclassical OMD when supx,y∈X DΦ(x,y) = +∞.Proof (of Theorem 3.2.1).The first step is the same as in the standard OMD proof. For all z ∈ X ,ft(x(t))− ft(z)(i)≤ 〈gt , x(t)− z〉(ii)=1ηt〈 xˆ(t)− wˆ(t+1), x(t)− z〉(iii)=1ηt(DΦ(x(t),w(t+1))−DΦ(z,w(t+1))+DΦ(z,x(t))),(3.8)51where (i) is from the subgradient inequality, (ii) follows from eq. (3.4) and (iii) isby Proposition B.1.6.The next step exhibits the main point of stabilization. Without stabilization wewould havex(t+1) =ΠΦX (w(t+1)) and DΦ(z,w(t+1))≥ DΦ(z,x(t+1))+DΦ(x(t+1),w(t+1))by Proposition B.1.8, so eq. (3.8) would lead to a telescoping sum involving DΦ(z, ·)if the stepsize were fixed. With a dynamic stepsize the analysis is trickier: we need aclaim that leads to telescoping terms by relating DΦ(z,w(t+1)) to DΦ(z,x(t+1)).Claim 3.2.1. Assume that γt = ηt+1/ηt ∈ (0,1]. Then(3.8)≤ DΦ(x(t)x(t+1) ;w(t+1))ηt+( 1ηt+1− 1ηt)︸ ︷︷ ︸telescopesDΦ(z,x(1))+DΦ(z,x(t))ηt− DΦ(z,x(t+1))ηt+1︸ ︷︷ ︸telescopes.Proof. First we derive the inequalityγt(DΦ(z,w(t+1))−DΦ(x(t+1),w(t+1)))+ (1− γt)DΦ(z,x(1))(i)≥ γtDΦ( zx(t+1) ;w(t+1)) + (1− γt)DΦ( zx(t+1) ;x(1))(ii)= DΦ(zx(t+1) ;y(t+1))(iii)≥ DΦ(z,x(t+1))where (i) is from the fact that DΦ(x(t+1),x(1)) ≥ 0 and γt ≤ 1, (ii) follows fromProposition B.1.7 and eq. (3.5) and (iii) is by Proposition B.1.8 and eq. (3.6).Rearranging and using γt > 0 yieldsDΦ(z,w(t+1)) ≥ DΦ(x(t+1),w(t+1))−( 1γt−1)DΦ(z,x(1))+1γtDΦ(z,x(t+1)).(3.9)52Plugging this into eq. (3.8) yields(3.8) =1ηt(DΦ(x(t),w(t+1))−DΦ(z,w(t+1))+DΦ(z,x(t)))≤ 1ηt(DΦ(x(t),w(t+1))−DΦ(x(t+1),w(t+1))+( 1γt−1)DΦ(z,x(1))− 1γtDΦ(z,x(t+1))+DΦ(z,x(t))),by eq. (3.9). The claim follows by the definition of γt .The final step is very similar to the standard OMD proof. Summing eq. (3.8)over t and using Claim 3.2.1 leads to the desired telescoping sum.T∑t=1(ft(x(t))− ft(z))≤T∑t=1(DΦ( x(t)x(t+1) ;w(t+1))ηt+( 1ηt+1− 1ηt)DΦ(z,x(1))+DΦ(z,x(t))ηt− DΦ(z,x(t+1))ηt+1)≤T∑t=1DΦ( x(t)x(t+1) ;w(t+1))ηt+(1η1+T∑t=1(1ηt+1− 1ηt))DΦ(z,x(1))=T∑t=1DΦ( x(t)x(t+1) ;w(t+1))ηt+DΦ(z,x(1))ηT+1.3.2.2 Primal-stabilized OMDAlgorithm 4 gives pseudocode showing our modification of OMD to incorporateprimal stabilization.The algorithm is analyzed in the following Theorem.53Algorithm 4 Online mirror descent with primal stabilization.Input: x(1) ∈ Rd , η : N→ R, γ : N→ R.for t = 1,2, . . . doIncur cost ft(x(t)) and receive gt ∈∂ ft(x(t))xˆ(t) = ∇Φ(x(t)) B map primal iterate to dual spacewˆ(t+1) = xˆ(t)−ηtgt B gradient step in dual space (3.12)w(t+1) = ∇Φ∗(wˆ(t+1)) B map dual iterate to primal space (3.13)y(t+1) =ΠΦX (w(t+1)) B project onto feasible region (3.14)x(t+1) = γty(t+1)+(1− γt)x(1) B stabilization in primal space (3.15)end forTheorem 3.2.2 (Regret bound for primal-stabilized OMD). Assume that ηt ≥ηt+1 > 0 for all t ≥ 1. Define γt = ηt+1/ηt ∈ (0,1] for all t ≥ 1. Let {x(t)}t≥1 bethe sequence of iterates generated by Algorithm 4. Furthermore, assume thatfor all z ∈ X , the map x 7→ DΦ(z,x) is convex on X . (3.10)Then for any sequence of convex functions { ft}t≥1 with each ft : X → R,Regret(T,z) ≤T∑t=1DΦ( x(t)y(t+1) ;w(t+1))ηt+DΦ(z,x(1))ηT+1∀T > 0. (3.11)Proof (of Theorem 3.2.2).Let z ∈ X . The first step is identical to the proof of Theorem 3.2.1 since theupdate rule in (3.12) is exactly the same as (3.4). Therefore, we have that (3.8)holds, that is,ft(x(t))− ft(z) ≤ 1ηt(DΦ(x(t),w(t+1))−DΦ(z,w(t+1))+DΦ(z,x(t))).54Claim 3.2.2. Assume that γt = ηt+1/ηt ∈ (0,1]. Then(3.8)≤DΦ( x(t)y(t+1) ;w(t+1))ηt+( 1ηt+1− 1ηt)︸ ︷︷ ︸telescopesDΦ(z,x(1))+DΦ(z,x(t))ηt− DΦ(z,x(t+1))ηt+1︸ ︷︷ ︸telescopes.Proof. First, we derive the inequalityγt(DΦ(z,w(t+1))−DΦ(y(t+1),w(t+1)))+ (1− γt)DΦ(z,x(1))= γtDΦ( zy(t+1) ;w(t+1)) + (1− γt)DΦ(z,x(1))≥ γtDΦ(z,y(t+1))+(1− γt)DΦ(z,x(1)) (by Proposition B.1.8 and eq. (3.14))≥ DΦ(z,x(t+1)) (by eq. (3.15), (3.10) and γt ∈ (0,1]) .Rearranging and using γt > 0 yieldsDΦ(z,w(t+1)) ≥ DΦ(y(t+1),w(t+1))−( 1γt−1)DΦ(z,x(1))+1γtDΦ(z,x(t+1)).(3.16)Plugging this into eq. (3.8) yields(3.8) =1ηt(DΦ(x(t),w(t+1))−DΦ(z,w(t+1))+DΦ(z,x(t)))≤ 1ηt(DΦ(x(t),w(t+1))−DΦ(y(t+1),w(t+1))+( 1γt−1)DΦ(z,x(1))− 1γtDΦ(z,x(t+1))+DΦ(z,x(t))),by eq. (3.16). The claim follows by the definition of γt .The final step is very similar to the proof of Theorem 3.2.1. The only differ-ence is that we are using Claim 3.2.2 instead of Claim 3.2.1 and we replaceDΦ( x(t)x(t+1) ;w(t+1)) with DΦ( x(t)y(t+1) ;w(t+1)). Formally, summing (3.8) over t and55using Claim 3.2.2 leads to the desired telescoping sum, that is,T∑t=1(ft(x(t))− ft(z))≤T∑t=1(DΦ( x(t)y(t+1) ;w(t+1))ηt+( 1ηt+1− 1ηt)DΦ(z,x(1))+DΦ(z,x(t))ηt− DΦ(z,xt+1)ηt+1)≤T∑t=1DΦ( x(t)y(t+1) ;w(t+1))ηt+(1η1+T∑t=1(1ηt+1− 1ηt))DΦ(z,x(1))=T∑t=1DΦ( x(t)y(t+1) ;w(t+1))ηt+DΦ(z,x(1))ηT+1.3.2.3 Dual averagingIn this section, we show that Nesterov’s dual averaging algorithm can be obtainedfrom a small modification to dual-stabilized online mirror descent. Furthermore ourproof of Theorem 3.2.1 can be adapted to analyze this algorithm.The main difference between DS-OMD and dual averaging is in the gradientstep, as we now explain. In iteration t+1 of DS-OMD, the gradient step is takenfrom xˆ(t), the dual counterpart of the iterate x(t):DS-OMD gradient step: wˆ(t+1) = xˆ(t)−ηtgt .Suppose that the algorithm is modified so that the gradient step is taken from yˆ(t),the dual point from iteration t before projection onto the feasible region. (Here yˆ(1)is defined to be xˆ(1).) The resulting gradient step is:Lazy gradient step: wˆ(t+1) = yˆ(t)−ηtgt . (3.17)As before, we setyˆ(t+1) = γtwˆ(t+1)+(1− γt)xˆ(1) (3.18)56Algorithm 5 Dual averaging with stepsize re-indexed as η2,η3, . . .Input: x(1) ∈ X , η : N→ R+, γ : N→ (0,1]yˆ(1) = ∇Φ(x(1))for t = 1,2, . . . doIncur cost ft(x(t)) and receive gt ∈ ∂ ft(x(t))yˆ(t+1) = xˆ(1)−ηt+1∑i≤t gi B dual averaging updatey(t+1) = ∇Φ∗(yˆ(t+1)) B map dual iterate to primal spacex(t+1) =ΠΦX (y(t+1)) B project onto feasible regionend forwhere γt = ηt+1/ηt . Then a simple inductive proof yields the following claim.Claim 3.2.3. wˆ(t) = xˆ(1)−ηt−1∑i<t gi and yˆ(t) = xˆ(1)−ηt ∑i<t gi for all t > 1.Thus, the algorithm with the lazy gradient step can be written as in Algorithm 5.This is equivalent to Algorithm 2 with the DA update, except that ηt in Algorithm 2corresponds to ηt+1 in Algorithm 5.Theorem 3.2.3 (Regret bound for dual averaging). Assume that ηt ≥ ηt+1 > 0 forall t > 1. Let {x(t)}t≥1 be the sequence of iterates generated by Algorithm 5. Thenfor any sequence of convex functions { ft}t≥1 with each ft : X → R,Regret(T,z) ≤T∑t=1DΦ( x(t)x(t+1) ;∇Φ∗(xˆ(t)−ηtgt))ηt+DΦ(z,x(1))ηT+1∀T > 0.(3.19)The proof parallels the proof of Theorem 3.2.1.Proof (of Theorem 3.2.3).The first step is very similar to the proof of Theorem 3.2.1. For all z ∈ X ,ft(x(t))− ft(z) (3.20)≤ 〈gt , x(t)− z〉 (subgradient inequality)=1ηt〈 yˆ(t)− wˆ(t+1), x(t)− z〉 (by eq. (3.17))=1ηt(DΦ(x(t),w(t+1))−DΦ(z,w(t+1))+DΦ( zx(t) ;y(t))), (3.21)57where we have used Proposition B.1.5 instead of Proposition B.1.6.As in the proof of Theorem 3.2.1, the next step is to relate DΦ(z,w(t+1)) toDΦ(z,y(t+1)) so that eq. (3.21) can be bounded using a telescoping sum. Thefollowing claim is similar to Claim 3.2.1.Claim 3.2.4. Assume that γt = ηt+1/ηt ∈ (0,1]. Then(3.21)≤DΦ(x(t)x(t+1) ;w(t+1))ηt+( 1ηt+1− 1ηt)︸ ︷︷ ︸telescopesDΦ(z,x(1))+DΦ(zx(t) ;y(t))ηt− DΦ(zx(t+1) ;y(t+1))ηt+1︸ ︷︷ ︸telescopes.Proof. The first two steps are identical to the proof of Claim 3.2.1.γt(DΦ(z,w(t+1))−DΦ(x(t+1),w(t+1)))+ (1− γt)DΦ(z,x(1))(i)≥ γtDΦ( zx(t+1) ;w(t+1)) + (1− γt)DΦ( zx(t+1) ;x(1))(ii)= DΦ(zx(t+1) ;y(t+1)).where (i) is from the fact that DΦ(x(t+1),x(1)) ≥ 0 and γt ≤ 1, (ii) is by Proposi-tion B.1.7 and eq. (3.18). Rearranging and using γt > 0 yieldsDΦ(z,w(t+1)) ≥ DΦ(x(t+1),w(t+1))−( 1γt−1)DΦ(z,x(1))+DΦ(zx(t+1) ;y(t+1))γt.(3.22)Plugging this into eq. (3.21) yields(3.21) =1ηt(DΦ(x(t),w(t+1))−DΦ(z,w(t+1))+DΦ( zx(t) ;y(t)))≤ 1ηt(DΦ(x(t),w(t+1))−DΦ(x(t+1),w(t+1))+( 1γt−1)DΦ(z,x(1))− DΦ(zx(t+1) ;y(t+1))γt+DΦ(zx(t) ;y(t))),by eq. (3.22). The claim follows by the definition of γt .58The final step is very similar to the proof of Theorem 3.2.1. Summing eq. (3.21)over t and using Claim 3.2.4 leads to the desired telescoping sum.T∑t=1(ft(x(t))− ft(z))≤T∑t=1(DΦ( x(t)x(t+1) ;w(t+1))ηt+( 1ηt+1− 1ηt)DΦ(z,x(1))+DΦ(zx(t) ;y(t))ηt− DΦ(zx(t+1) ;y(t+1))ηt+1)≤T∑t=1DΦ( x(t)x(t+1) ;w(t+1))ηt+(1η1+T∑t=1(1ηt+1− 1ηt))DΦ(z,x(1))=T∑t=1DΦ( x(t)x(t+1) ;w(t+1))ηt+DΦ(z,x(1))ηT+1. (3.23)For the second inequality we have also used that DΦ( zx(1) ;y(1)) = DΦ(z,x(1))since x(1) = y(1).Notice that eq. (3.23) is syntactically identical to eq. (3.7); the only differenceis the definition of w(t+1) in these two settings. It turns out to be more usefulfor applications to provide a convenient upper bound on eq. (3.23), which isthe conclusion of this theorem. Referring to eq. (3.17), we see that w(t+1) isclosely related to y(t), but less closely related to x(t). To control DΦ( x(t)x(t+1) ;w(t+1)), itturns out to be convenient to apply Proposition B.1.9 as follows. Taking p = y(t),pi = x(t) =ΠΦX (y(t)), v = x(t+1) and qˆ = ηtgt , we obtainDΦ( x(t)x(t+1) ;w(t+1))(i)= −DΦ( vpi ;∇Φ∗(pˆ− qˆ))(ii)≤ −DΦ( vpi ;∇Φ∗(pˆi− qˆ))= DΦ( x(t)x(t+1) ;∇Φ∗(xˆ(t)−ηtgt)),where (i) is from the fact that wˆ(t+1) = yˆ(t)−ηtgt = pˆ− qˆ and (ii) is by Proposi-tion B.1.9. Plugging this into (3.23) completes the proof.593.2.4 RemarksInterestingly, the doubling trick (Shalev-Shwartz, 2012) on OMD can be viewedas an incarnation of stabilization. To see this, set ηt := 1/√2blg tc and γt :=1{t is a power of 2}. Then, for each dyadic interval of length 2`, the first iterate is x(1)and a fixed learning rate 1/√2` is used. Thus, with these parameters, Algorithm 3reduces to the doubling trick.One should note that in Theorem 3.2.1 the stabilization parameter γt used inround t ≥ 1 depends on the stepsizes for rounds t and t + 1, need to “peek intothe future”. Thus, to use stabilization as in Theorem 3.2.1 the stepsize for round tcan depend on information available only up to round t−1. This will come intoplay, for example, when we derive first-order regret bounds on Section 3.3.2 wherethe stepsize is based on the subgradients of the past functions (instead of simplydepending on the count of rounds). Reindexing the stepsizes could fix the problem,but then the proof of Theorem 3.2.1 would look syntactically odd. Although thisdependence on the future may seem unnatural, in Section 3.4 we shall see that undersome mild conditions, stabilized OMD coincides exactly with DA with dynamicstepsizes after reindexing. This extends the same behavior observed between OMDand DA when the stepsizes are fixed. In this sense, stabilization may seem as anatural way to fix OMD for dynamic stepsizes.3.3 Applications3.3.1 Strongly-convex mirror mapsWe now analyze the algorithms of the previous section in the scenario that the mirrormap is strongly convex. Let ηt ,γt , ft be as above. The following result is a corollaryof Theorems 3.2.1, 3.2.2 and 3.2.3.Corollary 3.3.1 (Regret bound for dual-stabilized OMD). Suppose that Φ is ρ-strongly convex on X with respect to a norm ‖ · ‖. Let {x(t)}∞t=1 be the iteratesproduced by Algorithms 3, 4 or 5. (For Algorithm 4, the additional assumption60(3.10) is required.) ThenRegret(T,z) ≤T∑t=1ηt‖gt‖2∗2ρ+DΦ(z,x(1))ηT+1∀T > 0.This is identical to Nesterov’s bound for dual averaging (Nesterov, 2009,eq. 2.15) (taking his λi = 1 and his βi = 1/ηi). The proof is based on the fol-lowing simple proposition, which bounds the Bregman divergence when Φ isstrongly convex. See, e.g., Bubeck (Bubeck, 2015, pp. 300). A proof is given inAppendix B.2.Proposition 3.3.1. Suppose that Φ is ρ-strongly convex on X with respect to ‖ · ‖.Consider any x,x′ ∈ X and qˆ ∈ Rd . ThenDΦ( xx′ ;∇Φ∗(xˆ− qˆ)) ≤ ‖qˆ‖2∗/2ρ.Proof (of Corollary 3.3.1). The regret bounds proven by Theorems 3.2.1, 3.2.2and 3.2.3 all involve a summation with terms of the form3.2.1 : DΦ( x(t)x(t+1) ;w(t+1))3.2.2 : DΦ( x(t)y(t+1) ;w(t+1))3.2.3 : DΦ( x(t)x(t+1) ;∇Φ∗(xˆ(t)−ηtgt)).For Theorems 3.2.1 and 3.2.3, we have x(t+1) ∈ X , whereas for Theorem 3.2.2we also have y(t+1) ∈ X by eq. (3.14). For Theorems 3.2.1 and 3.2.2 we havew(t+1) = ∇Φ∗(xˆ(t)−ηtgt) by eq. (3.4) and eq. (3.10). Therefore all of these termsmay be bounded using Proposition 3.3.1 with x = x(t) and qˆ = ηtgt . This yields theclaimed bound.3.3.2 Prediction with expert adviceNext, we consider the setting of “prediction with expert advice”. In this setting, Dis Rd>0, X is the simplex ∆d ⊂ Rd , and the mirror map is Φ(x) = ∑di=1 xi logxi. (OnX , Φ is the negative of the entropy function.) The gradient of the mirror map and61its conjugate are∇Φ(x)i = ln(xi)+1 and ∇Φ∗(xˆ)i = exp(xˆi−1). (3.24)For any two points a ∈ D¯, b ∈ D, a small calculation shows that DΦ(a,b) is thegeneralized KL-divergenceDKL(a,b) =d∑i=1ai ln(ai/bi)−‖a‖1+‖b‖1.Note that the KL-divergence is convex on its second argument for any a ∈ D¯ =Rd≥0since the functions − ln(·) and absolute value are both convex. This means that thegeneral regret bounds for all algorithms discussed in Section 3.2, including primalstabilized OMD, hold in this setting.In this section, we will use the theorems in Section 3.2 to derive regret boundsfor this setting without much extra-work. As an intermediate step, we will derivebounds that use the following function:Λ(a,b) := DKL(a,b)+‖a‖1−‖b‖1+ ln‖b‖1 =d∑i=1ai ln(ai/bi)+ ln‖b‖1,which is a useful tool in the analysis of algorithms for the experts’ problem. Forexamples, see works from de Rooij et al. (2014, §2.1) and Cesa-Bianchi et al. (2007,Lemma 4). An initial observation shows that Λ is non-negative in the experts’setting.Proposition 3.3.2. Λ(a,b)≥ 0 for all a ∈ X , b ∈ D.Proof. Let us write Λ(a,b) =−∑di=1 ai ln biai + ln(∑di=1 bi). Since a is a probabilitydistribution, we may apply Jensen’s inequality to show that this expression isnon-negative.The following result is a corollary of Theorems 3.2.1, 3.2.2 and 3.2.3.Corollary 3.3.2. Assume that ηt ≥ ηt+1 > 0 for all t ≥ 1. Define γt = ηt+1/ηt ∈(0,1]. Let x(1) be the uniform distribution 1/d and let x(2),x(3), . . . be the iterates62produced by Algorithms 3, 4 or 5. ThenRegret(T,z) ≤T∑t=1Λ(x(t),∇Φ∗(xˆ(t)−ηtgt))ηt+lndηT+1∀T > 0. (3.25)The proof is a direct consequence of the following proposition, which is provenin Appendix B.3.Proposition 3.3.3. Let a,b ∈ X and c ∈ D. Then DΦ(ab ;c)≤ Λ(a,c).Proof (of Corollary 3.3.2). First of all, recall that the DKL is convex on its secondargument, which allows us to use the bound from eq. (3.11) for primal stabilizedOMD. As in the proof of Corollary 3.3.1, we first observe that the regret bounds (3.7),(3.11) and (3.19) all have sums with terms of the form DΦ( x(t)u(t) ;∇Φ∗(xˆ(t)−ηtgt))for some u(t) ∈ X . These terms may be bounded using Proposition 3.3.3. Finally,the standard inequality supz∈X DKL(z,x(1))≤ lnd completes the proof.Anytime regretAs another corollary of Corollary 3.3.2 we now derive an anytime regret bound inthe case of bounded costs. This matches the best known bound appearing in theliterature (Bubeck, 2011, Theorem 2.4) (Gerchinovitz, 2011, Proposition 2.1).Corollary 3.3.3. Suppose that gt ∈ [0,1]d for all t. Define ηt = 2√ln(d)/t andγt = ηt+1/ηt . Let x(1) be the uniform distribution 1/d and let x(2),x(3), . . . be theiterates produced by Algorithms 3, 4 or 5. ThenRegret(T,z) ≤√T lnd ∀T ≥ 1, z ∈ X . (3.26)The proof follows from Corollary 3.3.2 and Hoeffding’s Lemma, as shownbelow.Lemma 3.3.1 (Hoeffding’s Lemma (Cesa-Bianchi and Lugosi, 2006, Lemma 2.2)).Let X be a random variable with a≤ X ≤ b. Then for any s ∈ R,lnE[esX ]− sEX ≤ s2(b−a)28.63Proof (of Corollary 3.3.3). Denote gti to be the ith entry of the vector gt . Byeq. (3.24) we have ∇Φ∗(xˆ(t)−ηtgt)i = x(t)i exp(−ηtgti) for each i ∈ [d]. This to-gether with Lemma 3.3.1 for s =−ηt yieldsΛ(x(t),∇Φ∗(xˆ(t)−ηtgt)) = ηt〈x(t), gt 〉+ ln(d∑j=1x(t)j e−ηt gt j)≤ η2t8. (3.27)Plugging this and ηt = 2√lnd/t into eq. (3.25), we obtain the boundRegret(T,z) ≤√lnd(14T∑t=11√t+√T +12)(3.28)≤√lnd(2√T −14+√T +0.52)≤√T lnd (3.29)by Fact B.1.3 and sub-additivity of square root.First-order regret boundThe regret bound presented in previous section depends on√T ; this is known asthe “zeroth-order” regret bound. In some scenarios the cost of the best expert up totime T can be far less than T . This makes the problem somewhat easier, and it ispossible to improve the regret bound. Formally, let L∗T denote the cost of the bestexpert at time T . Then L∗T ≤ T due to our assumption that all costs are at most 1. A“first-order” regret bound depends on√L∗T instead of√T . See, for example, thebook from Cesa-Bianchi and Lugosi (2006, §2.4).The only modification to the algorithm is to change the stepsize. If the costsare “smaller than expected”, then intuitively time is progressing “slower thanexpected”. We will adopt an elegant idea from Auer et al. (2002b), which is to usethe algorithm’s cost itself as a measure of the progression of time, and to incorporatethis into the stepsize. They call this a “self-confident” stepsize.Corollary 3.3.4. Assume that gt ∈ [0,1]d . Set ηt =√ln(d)/(1+∑i<t〈gi, xi 〉) andγt = ηt+1/ηt . Denote the minimum total cost of any expert up to time T as L∗T :=64min j∈[d]∑Tt=1 gt j. ThenRegret(T,z) ≤ 2√ln(d)L∗T +8lnd ∀T ≥ 1.The main ingredients are the following alternative bound on Λ, which is provenin Appendix B.3, and some standard scalar inequalities.Proposition 3.3.4. Let a ∈ X , qˆ ∈ [0,1]d and η > 0. Then Λ(a,∇Φ∗(aˆ−η qˆ)) ≤η2〈a, qˆ〉/2.Proof (of Corollary 3.3.4). From Corollary 3.3.2 and Proposition 3.3.4, we haveT∑t=1(〈gt , x(t) 〉−〈gt , z〉) ≤ T∑t=1ηt〈gt , x(t) 〉2+lndηT+1. (3.30)The algorithm’s total cost at time t is denoted At = ∑i≤t〈gi, x(i) 〉. Recall that thetotal cost of the best expert at time T is L∗T = minz∈∆d ∑Tt=1〈gt , z〉 and the stepsizeis ηt =√ln(d)/(1+At−1). Substituting these into eq. (3.30),AT −L∗T ≤√lnd(12T∑t=1〈gt , x(t) 〉√1+At−1+√1+AT)≤√lnd(√AT +√AT +1)by Proposition B.1.1 with ai = 〈gi, x(i) 〉 and u = 1. Rewriting the previous inequal-ity, we have shown thatAT −L∗T ≤ 2√ln(d)AT +√lnd.By Proposition B.1.2 we obtainAT −L∗T ≤ 2√ln(d)L∗T +√lnd+2(lnd)3/4+4lnd.Since the left-hand side equals Regret(T,z), the result follows.Now we compare our bound with some existing results in the literature: ourconstant term of 2 obtained in Corollary 3.3.4 is better than the constant (√2/(√2−651)) obtained by the doubling trick (Cesa-Bianchi and Lugosi, 2006, Exercise 2.8),and the constant (2√2) in Auer et al. (2002b)’s work but worse than the constant(√2) of the best known first-order regret bound (Yaroshinsky et al., 2004), which isobtained by a sophisticated algorithm. We also match the constant 2 of the Hedgealgorithm from de Rooij et al. (2014, Theorem 8). Their result is actually moregeneral; we could similarly generalize our analysis, but that would deviate too farfrom the main purpose of this chapter.3.4 Comparing DS-OMD and DAIn this section we shall write the iterates of dual-stabilized OMD in two equivalentforms. First we shall write it in a proximal-like formulation similar to the mirrordescent formulation in Beck and Teboulle (2003), shedding some light into theintuition behind dual-stabilization. We then write the iterates from DS-OMD in aform very similar to the original definition of DA in Nesterov (2009). The laterwill allow us to intuitively understand why OMD does not play well with dynamicstep-size and to derive simple sufficient conditions under which DS-OMD and DAgenerate the same iterates, mimicking the relation between OMD and DA for a fixedstepsize.Beck and Teboulle (2003) showed that the iterate x(t+1) for round t +1 fromOMD is the unique minimizer over X of ηt〈gt , · 〉+DΦ(·,x(t)), where gt ∈ ∂ ft(x(t)).The next proposition extends this formulation to DS-OMD, recovering the resultfrom Beck and Teboulle when γt = 1. The proof, which can be found in Ap-pendix B.4, is a simple application of the optimality conditions of eq. (3.31).Proposition 3.4.1. Let { ft}t≥1 be a sequence of convex functions with ft : X → Rfor each t ≥ 1. Let η : N→R>0 and γ : N→ [0,1]. Let {x(t)}t≥1 and {gt}t≥1 be asin Algorithm 3. Then, for any t ≥ 1,{x(t+1)}= argminx∈X(γt(ηt〈gt , x〉+DΦ(x,x(t)))+(1− γt)DΦ(x,x(1))). (3.31)In spite of their similar descriptions, Orabona and Pa´l (2018) showed that OMDand DA may behave in extremely different ways even on the well-studied experts’problem with similar choices of stepsizes. This extreme difference in behavior is not66clear from the classical algorithmic description of these methods as in Algorithm 2.In the case of DA, it is well-known that DA can be seen as an instance of the FTRLalgorithm; see Bubeck (2015, §4.4) or Hazan (2016, §5.3.1). More specifically,if {x(t)}t≥1 and {gt}t≥1 are as in Algorithm 2 with the DA update, then for everyt ≥ 0, 4{x(t+1)}= argminx∈X(ηt+1t∑i=1〈gi, x〉−〈 xˆ(1), x〉+Φ(x)). (3.32)In the next theorem, proven in Appendix B.4, we write DS-OMD in a similar form,but with vectors from the normal cone of X creeping into the formula due to theback and forth between the primal and dual spaces. Recall that the normal cone ofX at a point x ∈ X is the setNX (x) := { p ∈ Rd | 〈 p, z− x〉 ≤ 0 ∀ z ∈ X } .The result in McMahan (2017, Theorem 11) is similar but slightly more intricatedue to the use of time-varying mirror maps. Moreover, this result does not directlyapply when we have stabilization.Theorem 3.4.1. Let { ft}t≥1 with ft :X →R be a sequence of convex functions andlet η : N→ R>0 be non-increasing. Let {x(t)}t≥1 and {gt}t≥1 be as in Algorithm 3.Then, there are {pt}t≥1 with pt ∈ NX (x(t)) for all t ≥ 1 such that, if γi = 1 for alli≥ 1, then for all t ≥ 0{x(t+1)}= argminx∈X( t∑i=1〈ηigi+ pi, x〉−〈 xˆ(1), x〉+Φ(x))(3.33)and if γi = ηi+1ηi for all i≥ 1, then for all t ≥ 0{x(t+1)}= argminx∈X(ηt+1t∑i=1〈gi+ p′i, x〉−〈 xˆ(1), x〉+Φ(x)). (3.34)where p′t :=1ηt pt ∈ NX (x(t)) for every t ≥ 1.With the above theorem, we may compare the iterates generated from DA, OMD,and DS-OMD by comparing the formulas eq. (3.32), eq. (3.33), and eq. (3.34). For4The 〈∇Φ(x(1)), x〉 term disappears if x(1) minimizes Φ on X .67the simple unconstrained case where X = Rd we have NX (x(t)) = {0} for eacht ≥ 1 and DA and DS-OMD are identical. However, if the stepsize is not constant,OMD is not equivalent to the latter methods. In particular, if ηt ∝ 1/√t, eq. (3.33)shows that the subgradients of the earlier-seen functions have a bigger weight onthe iterates if compared to the subgradients of functions from later rounds. In otherwords, OMD may be sensitive to the ordering of the functions, and adversarialorderings may affect its performance.When X is an arbitrary convex set, DA and DS-OMD are not necessarilyequivalent anymore due to the vectors from the normal cone of X . If we know thatthe iterates live in the relative interior of X , the next lemma shows that these vectorsdo not affect the set of minimizers from eq. (3.34).Lemma 3.4.1. For any x˚∈ riX , we have NX (x˚) = (−NX (x˚))∩NX (x˚). In particular,for any p ∈ NX (x˚) we have 〈 p, x〉= 〈 p, x˚〉 for every x ∈ X .With this lemma, we can easily derive simple and intuitive conditions underwhich DS-OMD and DA are equivalent.Corollary 3.4.1. Let D ⊆ Rd be the interior of the domain of Φ, let {x(t)}t≥1 bethe DS-OMD iterates as in Algorithm 3 and let {x(t)′}t≥1 be the DA iterates as inAlgorithm 2 with DA updates. If D∩X ⊆ riX , then x(t) = x(t)′ for each t ≥ 1.Proof. Let t ≥ 1. Since x(t)=ΠΦX (y(t)), where y(t) is as in Algorithm 3, Lemma B.4.2implies x(t) ∈D∩X ⊆ riX . By Lemma 3.4.1 we have that the vectors on the normalcone in eq. (3.34) do not affect the set of minimizers, which implies that eq. (3.32)and eq. (3.34) are equivalent.An important special case of the above corollary is the prediction with expertadvice setting as in Section 3.3.2, where D = Rd>0 and X is the simplex ∆d . Inthis setting, X ∩D = {x ∈ (0,1)d | ∑di=1 xi = 1}= riX . By the previous corollaryDS-OMD and DA produce the same iterates in this case even for dynamic stepsizes.Classical OMD and DA were already known to be equivalent in the experts’ settingfor a fixed learning rate (Hazan, 2016, §5.4.2). In contrast, with a dynamic stepsize,the DA and OMD iterates are certainly different, since OMD with a dynamiclearning rate may have linear regret (Orabona and Pa´l, 2018), whereas DA hassublinear regret.683.5 DiscussionIn this chapter we modified OMD via stabilization in order to guarantee sublinearregret even when using the method with a dynamic stepsize. We showed that (primaland dual) stabilized-OMD recover the regret bounds enjoyed by DA in the anytimesetting, presented some applications of our results, and analyzed the similarities anddifferences between DS-OMD, OMD, and DA.Our bounds for the problem of prediction with expert advice nearly matchthe current state-of-the-art. A distinctive feature of our proofs are their relativesimplicity if compared to other results from the literature. It is our hope that thesimplicity of our analysis framework allows it to be extended to other problems.Moreover, the modularity of our proofs allowed us to extend this analysis for DA, afact interesting on its own since drastically different analysis techniques are usuallyused to analyze DA in the literature (such as the Follow the Leader-Be the LeaderLemma and optimality conditions of eq. (3.32), see Section 2.3 from Shalev-Shwartz(2012) for an example). This together with our analysis from Section 3.4 helpsdemystify the connections between DA and OMD, since in spite of having similardescriptions they had extremely different analyses and behaved wildly differentlyin some scenarios. We believe that a better understanding between the differencesbetween DA and OMD will be helpful in future applications and in the design ofnew algorithms.69Chapter 4Fast convergence of stochasticsubgradient descent underinterpolation4.1 Background and motivationWe consider the empirical-risk minimization problem:minx∈Rdf (x) :=1nn∑i=1fi(x), (P)where n is the number of data points, and fi’s are continuous functions that measureshow well our model x fit the data. The empirical-risk minimization formulation isprevalent in data-fitting problems—from problems that as simple as linear regressionto complicated modern applications such as using deep neural networks for imageclassification and language modeling.Stochastic (sub)gradient descent (SGD) is an simple first-order algorithm withthe update rule1:x(t+1) = x(t)−ηtgt gt ∈ ∂ fi(x(t)),1One could also update the iterate based on a small batch of data points, which is known as themini-batch SGD. In this thesis, we consider the single sample case for simplicity.70objective f (x) smooth nonsmooth strongly cvxGD O(ε−1) O(ε−2) O(log(ε−1))SGD O(ε−2) O(ε−2) O(ε−1)Table 4.1: Iteration complexity of deterministic gradient and stochastic gradi-ent methods.where i is uniform randomly sampled from {1,2, . . . ,n} in each iteration. The itera-tion complexities of SGD and GD for convex objective functions are summarized inTable 4.1. Although in general SGD converges slower than GD, a remarkable prop-erty of SGD is that both the per iteration cost and the convergence rate of SGD areindependent from n. This property makes SGD particularly effective for (P) whenn is large. In fact, SGD and its variants are indeed the most popular algorithms formodern big data applications including image classification and language modeling.Given the empirical success and popularity of SGD, a huge line of works havebeen devoted into• designing new variants of SGD to improve its convergence either theoreticallyor empirically;• refining the classical analysis of SGD in various settings e.g., convex andnonconvex , and thus better understand SGD from a theoretical point of view.This chapter belongs to the latter category. Next, we will briefly review somemilestones in the recent development of SGD and describe the motivation for thischapter.4.1.1 Practical algorithms based on SGDSGD is the state-of-the-art algorithm for big data problems. Many practical algo-rithms developed for a wide range of applications are related to SGD. For example,the randomized Kaczmarz algorithm (Kaczmarz, 1937; Needell et al., 2014) forsolving linear systems and the Pegasos algorithm (Shalev-Shwartz et al., 2007)for solving linear support vector machines. In the field of deep learning and re-inforcement learning, specialized algorithms such as ADAM (Kingma and Ba,2015), RMSprop (Tieleman and Hinton, 2012), TRPO (Schulman et al., 2015) and71PPO (Schulman et al., 2017), etc. are developed to adapt SGD to different problemstructures. Beside efficiency, SGD and its variants are believed to exhibit implicitregularization (Gunasekar et al., 2017), which is another factor that makes SGD thedominate algorithm for training machine learning models.4.1.2 Parallel and distributed SGDThe computation performance of single processor is reaching its limit due to en-ergy dissipation. Designing scalable algorithms that obtain near-linear speed upwith multiple computing cores to address real-world huge scale problems becomesan important research topic. Under such background, many parallel SGD algo-rithms (Goyal et al., 2017; Li et al., 2014; Lian et al., 2015; Liu et al., 2015; Rechtet al., 2011; You et al., 2018) are proposed to adapt SGD to multi-CPUs, multi-GPUs or multi-machines setting. Both synchronous and asynchronous parallelSGD algorithms have been studied heavily in both theory and practice. On thetheory side, Bertsekas and Tsitsiklis first formally brought optimization and parallelcomputation together, their seminal work (Bertsekas and Tsitsiklis, 1989) estab-lished the theoretical foundations of parallel optimization and the theoretical toolstherein are still being used nowadays (the book from Bertsekas and Tsitsiklis is veryforward-thinking in 1989 since parallel computing is less of interest at that time).More recently, HOGWILD! (Recht et al., 2011), an asynchronous parallel SGD withtheoretical guarantees and empirical success, brought resurged interest in parallelSGD for machine learning tasks; the convergence of asynchronous parallel SGD invarious settings are also studied by subsequent works (Lian et al., 2015; Liu et al.,2015). On the practical side, researchers are keeping decreasing the training time ofneural networks on the imagenet dataset by scaling SGD to hundreds or even thou-sands of GPUs. The training time is reduced from the original 6 days (Krizhevskyet al., 2012) to 1 hour in 2017 (Goyal et al., 2017), around 10 minutes in 2018 (Youet al., 2018) and less than 5 minutes by the time of the completion of this thesis2.2Benchmark on imagenet training: https://dawn.cs.stanford.edu/benchmark/ImageNet/train.html,last access time 07/02/2021.724.1.3 Variance reductionThe variance reduction technique that originates from Schmidt et al. (2017) iscertainly an exciting advancement of SGD in the last decade. The resulting variance-reduced SGD algorithms including SAG (Schmidt et al., 2017), SVRG (Johnsonand Zhang, 2013), SAGA (Defazio et al., 2014) and SARAH (Nguyen et al., 2017)improved the iteration complexities of SGD from O(ε−2) to O(ε−1) for convexand smooth objectives and O(ε−1) to O(log(ε−1)) for strongly convex and smoothobjectives while maintaining the cheap per iteration cost. These rates matches therates of deterministic GD.4.1.4 SGD with the interpolation conditionAlthough variance-reduced SGD theoretically attains a faster convergence rate thanvanilla SGD, however practitioners find that it cannot outperform vanilla SGDin the practice of training modern machine learning models. This discrepancybetween theory and practice was recently filled by Ma et al. (2018), who showedthat for over-parameterized models that satisfy the interpolation condition, vanillaSGD has an inherent variance reduction functionality as iterates converging to thesolution. The interpolation condition means that our model has the ability to fitall training data perfectly; we will give a more formal description of this conceptin Section 4.2. Therefore, in the practice of training over-parameterized models,vanilla SGD already has the fast convergence rate that similar to the deterministicGD, and further explicitly adding the variance reduction step to SGD will not giveus a better convergence.The detailed convergence rate of SGD for convex objective functions under theinterpolation condition is summarized in Table 4.2. As we can see from the table,there is a gap between the convergence rate of SGD for smooth and nonsmoothobjectives under interpolation. However, in the practice of training machine learningmodels, the nonsmootheness from the model does not cause much trouble (Glorotet al., 2011; Goodfellow et al., 2016). Neural networks with nonsmooth activationfunction such as ReLU activation can usually be trained as fast as (or even fasterthan) the one with smooth activation function such as softplus. Therefore there is adiscrepancy between theory and practice in the nonsmooth world.73objective f (x) smooth nonsmooth SC + smooth SC + nonsmoothSGD O(ε−2) O(ε−2) O(ε−1) O(ε−1)SGD + interpolation O(ε−1) O(ε−2) O(log(ε−1)) O(ε−1)Table 4.2: Iteration complexity of SGD with and without the interpolationcondition. SC stands for strongly convex.In this chapter, we aim to fill the gap between the convergence rates of SGDfor smooth and nonsmooth objectives under the interpolation condition. We willdescribe some semi-smoothness properties of the empirical-risk minimization prob-lem. These properties, together with the interpolation condition, allow us to provethat stochastic subgradient method has iteration complexity O(ε−1) for convexobjectives, and O(log(ε−1)) for strongly-convex objectives. These rates improvedthe classic bounds O(ε−2) and O(ε−1) for convex and strongly-convex objectivesand match the convergence rates of SGD for convex and smooth objectives underinterpolation. We also prove that the iteration bound O(ε−1) is optimal in theconvex and interpolation setting. In contrast to the case with a smooth objectivefunction, subgradient-based methods cannot be further accelerated for nonsmoothmodel—even with the interpolation assumption.4.2 PreliminariesFirst, we introduce a terminology used in this chapter. We use SSGD to denote thestochastic subgradient descent, which is widely used when the objective function isnonsmooth.We consider our objective to be the unconstrained empirical risk-minimizationproblemminimizex∈Rdf (x) :=1nn∑i=1fi(x) where each fi(x) := `(hi(x)) (4.1)and n is the number of data points. Throughout the paper, we use x∗ to denote anysolution of eq. (4.1), and thus f ∗ := f (x∗) is the optimal objective value. We assumethat the loss function ` : R→ R≥0 is 1-dimensional function that is nonnegative74with inf`= 0, convex and 1-smooth, i.e.,|`′(α)− `′(β )| ≤ |α−β | ∀α,β ∈ R.Without loss of generality, we also assume `(0) = 0. Common examples for theloss function include• 2-norm loss: `(x) = x2;• logistic loss: `(x) = log(1+ exp(x));• 2-norm hinge loss: `(x) = (max{0,x})2.The n functions hi’s are Lipschitz continuous with respect to a fixed parameter L:|hi(x)−hi(y)| ≤ L‖x− y‖ ∀x,y ∈ Rd , ∀i ∈ [n].We make no assumption on their smoothness properties. Many key machine learningtasks can be formulated as (4.1), including training deep neural networks withnonsmooth activations, such as the ReLU function. Here and throughout, thefunction ‖ · ‖ is the 2-norm of a vector, unless otherwise specified.Our analysis relies on the Clark generalized gradient (Clarke, 1990) of thenonconvex and nonsmooth function h, defined as the convex hull of all valid limitinggradients:∂h(x) := conv{u | u = limk→∞∇h(xk), xk→ x} ,This definition additionally requires h to be almost everywhere differentiable byRademacher’s theorem. We refer readers to Clarke (1990), and more recently, toZhang et al. (2020), who use this generalized gradient in a related analysis. Theseproperties of the generalized gradient are needed for our analysis:(chain rule) ∂ fi(x) = `′(h(x)) ·∂h(x),(gradient bound) ‖∂hi(x)‖ ≤ L, ∀x ∈ Rd , ∀i ∈ [n].The second property follows from the L-Lipschitz continuity of each function hi.We define the norm of the generalized gradient at a vector x as‖∂hi(x)‖= sup{‖z‖ | z ∈ ∂hi(x)}.75Algorithm 6 Stochastic subgradient descent. The learning rate function ηt :N→R+returns the learning rate at iteration t.1: Initialize: x(1) ∈ Rd2: for t = 1,2, . . . do3: select i ∈ {1,2, . . . ,n} uniformly at random4: compute g(t) ∈ ∂ fi(x(t))5: x(t+1) = x(t)−ηtg(t)6: end forAlgorithm 6 describes the SSGD method.The interpolation conditionBy our assumption on fi’s, we can immediately conclude that f ≥ 0. The inter-polation condition means that our model has the ability to fit all training samplesperfectly and therefore achieve zero training loss at the solution. Formally, theinterpolation condition is defined as f ∗ = 0. It has been shown that the interpolationcondition holds for overparameterized neural networks (Jacot et al., 2018), and theinterpolation condition is gaining increasing interest in recent years. Note that wedo not make interpolation as an assumption here and most of our analysis in thefollowing sections holds in general without assuming interpolation. We will high-light the interpolation condition when the presented results require the interpolationcondition.4.3 Main resultsWe present some semi-smoothness properties and the convergence analysis in thissection.4.3.1 Bounds and Lipschitz properties of the generalized gradientOur analysis hinges on establishing that the objective function f is “almost” dif-ferentiable at points with small objective value. This is implied by the followingproposition, which holds even without the interpolation condition.Proposition 4.3.1 (Generalized growth condition). Assume that the assumptions76stated in Section 4.2 hold. Then for all d-vectors x,‖∂ fi(x)‖2 ≤ 2L2 fi(x) ∀i ∈ [n]. (4.2)Consequently,‖∂ f (x)‖2 ≤ 1nn∑i=1‖∂ fi(x)‖2 ≤ 2L2 f (x). (4.3)Eq. (4.3) implies that if the objective value f (x) = 0, then the generalizedgradient contains only the origin: ∂ f (x) = {0}. It thus follows from (Clarke, 1981,Property 10) that f is differentiable at any point with zero objective value. Thismeans that if there is a solution with a zero value, then it must be a fixed point of thesubgradient method. This property does not hold for functions that are nonsmooth atsolution (for example the absolute value function f (x) = |x|) and makes it possiblefor the subgradient method to converge to solution with constant learning rate.Proof of Proposition 4.3.1. Fix an arbitrary d-vector x. The subdifferential ∂ fi(x)=`′(hi(x)) ·∂hi(x). Because each function hi is L-Lipschitz continuous, we can thendeduce that‖∂ fi(x)‖2 ≤ L2 [`′(hi(x)) ]2(i)= L2(`′(hi(x))− `′(0))2(ii)≤ 2L2(`(hi(x))− `(0)) (4.4)= 2L2 fi(x).Step (i) follows from the assumption that `(0) = minλ∈R `(λ ) = 0, which implies`′(0) = 0. Step (ii) follows from the fact that any convex L-smooth function u :Rn→ R satisfies the boundu(a)−u(b)−〈∇u(b),a−b〉 ≥ 12L‖∇u(a)−∇u(b)‖2;77see Nesterov (2014, Theorem 2.1.5). Make the identificationsu = `, a = hi(x), b = 0to immediately obtain eq. (4.4) and thus the proof for eq. 4.2.Eq. (4.3) can be obtained directly from Jensen’s inequality:‖∂ f (x)‖2 =∥∥∥1nn∑i=1∂ fi(x)∥∥∥2 ≤ 1nn∑i=1‖∂ fi(x)‖2 ≤ 1nn∑i=12L2 fi(x) = 2L2 f (x).The composite structure of the functions fi allows us to develop a semi-Lipschitzbound on their generalized gradients. Moreover, we show that it is possible to obtaina global convex majorant for each function fi, which holds without assumingconvexity. Interestingly, eq. (4.6) overlaps with the “semi-smoothness” property ofover-parameterized neural networks derived by Allen-Zhu et al. (2019, Theorem 4).In contrast to the result from Allen-Zhu et al. (2019), however, our result does notuse any special properties of over-parameterized neural networks. Proposition 4.3.2may thus be of more general interest.Proposition 4.3.2 (Semi-smoothness). Assume that the assumptions stated in Sec-tion 4.2 hold. Then for all vectors x1 and x2, and each i ∈ [n],‖∂ fi(x2)−∂ fi(x1)‖ ≤ L2‖x2− x1‖+2L√2min{ fi(x1), fi(x2)}, (4.5)andfi(x2)≤ fi(x1)+ 〈∂ fi(x1),x2− x1〉+ L22‖x2− x1‖2+2L‖x2− x1‖√2 fi(x1).(4.6)The proof of Proposition 4.3.2 is rather tedious, interested readers can find it inAppendix C.784.3.2 Convergence rate of stochastic subgradient descentWe now present a global convergence analysis for the SSGD algorithm under theadditional assumption that the objective f of eq. (4.1) is convex. We develop abound on the expected progress of the objective value that depends on the optimalvalue f ∗, rather than on the Lipschitz bound on the function itself, which is theusual bound in the literature. Our proof is based on a simple modification of theclassical proof of subgradient descent method.Theorem 4.3.1 (Global convergence rate of SSGD). Assume f is convex. Then forany positive integer T and any learning rate function ηt that satisfies ∑Tt=1(ηt −L2η2t )> 0,mint∈[T ]E[ f (x(t))− f ∗]≤ ‖x(1)− x∗‖2+2L2 f ∗∑Tt=1η2t2∑Tt=1(ηt −L2η2t ). (4.7)Proof. Let g(t)i ∈ ∂ fi(x(t)). Then each iterate x(t) satisfies the boundE[‖x(t+1)− x∗‖2∣∣∣ x(t)]=1nn∑i=1‖x(t)−ηtg(t)i − x∗‖2= ‖x(t)− x∗‖2−2ηt〈1nn∑i=1g(t)i ,x(t)− x∗〉+1nn∑i=1η2t ‖g(t)i ‖2(i)≤ ‖x(t)− x∗‖2−2ηt( f (x(t))− f (x∗))+ 1nn∑i=1η2t ‖g(t)i ‖2 (4.8)(ii)≤ ‖x(t)− x∗‖2−2ηt( f (x(t))− f (x∗))+2η2t L2 f (x(t)), (4.9)where (i) follows from the convexity of f , and (ii) follows from Proposition 4.3.1.Take expectations on both sides of the inequality (4.9) and rearrange to obtain(2ηt −2L2η2t ) ·E[( f (x(t))− f ∗)]≤ E[‖x(t)− x∗‖2]−E[‖x(t+1)− x∗‖2]+2L2η2t f ∗.(4.10)79Summing inequality (4.10) over t ∈ {1,2, . . . ,T} yieldsT∑t=1(2ηt −2L2η2t )E[ f (x(t))− f ∗]≤ ‖x(1)− x∗‖2+2L2 f ∗T∑t=1η2t .Divide both sides by 2∑Tt=1(ηt −L2η2t )> 0, we obtain the desired result.This proof mirrors closely the classical proof of SSGD, which assumes that thesubdifferential of each fi is bounded, i.e., ‖∂ fi(x)‖≤G for some constant G. We useProposition 4.3.1 in eq. (4.8) to avoid the bounded subgradient assumption and toexpress the convergence rate using the minimal value f ∗. This modification allowsus to leverage the interpolation assumption that f ∗ = 0. We use Theorem 4.3.1 toimmediately deduce the following convergence rate result.Corollary 4.3.1 (Global convergence rate of SSGD with constant learning rate).Assume that f is convex and that the learning rate ηt = 1/(2L2) is constant for allt > 0. Then for any positive integer T , the SSGD iterates x(t) satisfymint∈[T ]E[ f (x(t))− f ∗]≤ (2L2/T )‖x(1)− x∗‖2+ f ∗.Furthermore, when f ∗ = 0 (interpolation holds),mint∈[T ]E[ f (x(t))]− f ∗]≤ (2L2/T )‖x(1)− x∗‖2.Theorem 4.3.1 indicates that the SSGD method converges at rate O(ε−1) wheninterpolation holds. This rate matches the convergence rate of SGD for smoothobjective functions under interpolation (Schmidt and Le Roux, 2013). When inter-polation does not hold, Theorem 4.3.1 implies that SSGD, on expectation, couldobtain objective value lower than 2 f ∗+ ε in O(ε−1) time, which could be close tof ∗ when f ∗ is close to 0 (interpolation nearly holds).Next, we derive convergence rate under the stronger assumption that f is stronglyconvex, which means thatf (x1)≥ f (x2)+ 〈g,x1− x2〉+ µ2 ‖x1− x2‖2 ∀x1,x2 ∈ Rd , ∀g ∈ ∂ f (x2) (4.11)80for some constant µ > 0. Note that recent works indicated that in order to provelinear convergence rate, the strong convexity assumption can be relaxed to otherweaker assumptions that could hold even for some nonconvex functions (Karimiet al., 2016; Qian et al., 2019). For simplicity, we assume strong convexity in ouranalysis.Theorem 4.3.2 (Global convergence rate of SSGD under strong convexity). Assumethat f is µ-strongly convex and that the learning rate ηt = 1/L2 is constant for allt > 0. Then for any positive integer T , the SSGD iterates x(t) satisfyE[‖x(T )− x∗‖2]≤(1− µL2)T−1‖x(1)− x∗‖2+ 2µf ∗. (4.12)Proof. Use the definition of the SSGD iteration to obtainE[‖x(t+1)− x∗‖2∣∣∣ x(t)](i)= ‖x(t)− x∗‖2−2ηt〈1nn∑i=1g(t)i ,x(t)− x∗〉+1nn∑i=1η2t ‖g(t)i ‖2(ii)≤ ‖x(t)− x∗‖2−2ηt( f (x(t))− f ∗)−µηt‖x∗− x(t)‖2+ 1nn∑i=1η2t ‖g(t)i ‖2(iii)≤ (1−µηt)‖x(t)− x∗‖2−2ηt( f (x(t))− f ∗)+2η2t L2 f (x(t)) (4.13)(iv)=(1− µL2)‖x(t)− x∗‖2+ 2L2f ∗, (4.14)where (i) follows from the same argument of the proof of Theorem 4.3.1; (ii) followsfrom the µ-strong convexity of f (see eq. (4.12)); (iii) follows from eq. (4.2); and(iv) follows from the definition of the learning rate ηt = 1/L2.Taking expectation to both sides of eq. (4.14) and recursively apply it to t ∈{1,2, . . . ,T} to deduceE[‖x(T )− x∗‖2]≤(1− µL2)T−1‖x(1)− x∗‖2+T−2∑t=0(1− µL2)t 2L2f ∗(i)≤(1− µL2)T−1‖x(1)− x∗‖2+ 2µf ∗,81where (i) follows from the fact that ∑∞t=0(1− µL2)t= L2/µ .Theorem 4.3.2 indicates that the SSGD converges to the ball centered at x∗with radius√2 f ∗/µ at a linear rate. If interpolation also holds, SSGD convergesto the solution linearly. Again, this rate matches the convergence rate of SGD forsmooth and strongly convex objectives in the interpolation setting (Ma et al., 2018;Schmidt and Le Roux, 2013). Similar to Corollary 4.3.1, when interpolation isnearly satisfied, SSGD converges linearly to an√2 f ∗/µ-approximate solution.Corollary 4.3.1 and Theorem 4.3.2 also provide insight into the effect of learningrate schedules on the performance of SSGD. A learning rate schedule ηt that decaysas t−1/2 is optimal because it causes the algorithm to exhibit a complexity boundof O(ε−2), which is the theoretical lower bound (Nemirovski and Yudin, 1983).However, the learning rate scheduleO(t−1/2) is slow in practice. Corollary 4.3.1 andTheorem 4.3.2 partially explain the discrepancy between the theory and practice:many machine learning models exhibit interpolation or near interpolation, andan aggressive constant learning-rate schedule works better than the conservativeworst-case optimal learning rate ηt =O(t−1/2).Other than the learning rate scheduling ηt = 1/L2, we can adopt a more carefullytuned learning rate scheduling (Stich, 2019) to obtain a refined convergence rate.We start with eq. (4.13), by setting ηt ≤ 1/(2L2), we obtainE[‖x(T )− x∗‖2]≤ (1−µηt)‖x(t)− x∗‖2−2ηt( f (x(t))− f ∗)+2η2t L2 f (x(t))≤ (1−µηt)‖x(t)− x∗‖2−2(ηt −η2t L2)( f (x(t))− f ∗)+2η2t L2 f ∗≤ (1−µηt)‖x(t)− x∗‖2−ηt( f (x(t))− f ∗)+2η2t L2 f ∗. (4.15)The last line is true by the condition ηt ≤ 1/(2L2). Eq. (4.15) coincides with Stich(2019, Lemma 1) by setting his γt = ηt , σ2 = 2L2 f ∗. Then we can directly applyStich (2019, Theorem 5) to get the following corollary.Corollary 4.3.2. Assume that f is µ-strongly convex for some µ ≥ 0. For allpositive integer T , there exist a constant learning rate scheduling ηt :=α ≤ 1/(2L2),82such that the SSGD iterates x(t) satisfymint∈[T+1]E[ f (x(t))− f ∗]≤ 64L2‖x(1)− x∗‖2 exp(− µT4L2)+72L2 f ∗µTfor µ > 0, andmint∈[T ]E[ f (x(t))− f ∗]≤ 4L2‖x(1)− x∗‖2T+4L2 f ∗‖x(1)− x∗‖√Tfor µ = 0.4.3.3 Lower boundsWe have proven that, under interpolation, SGD for smooth problems and SSGD fornonsmooth problems exhibit the same convergence rates. This causes us to considerthe following questions:• Is it possible to induce momentum-type acceleration for SSGD under interpo-lation? Vaswani et al. (2019) and Liu and Belkin (2020) showed recently thatacceleration is possible for SGD under interpolation. It seems plausible thatsimilar techniques could be used for nonsmooth problems.• Can the composite structure fi = `◦hi, which is central to our analysis, be usedto establish an improved convergence rate for SSGD without interpolation?In other words, we know the lower bound on the iteration complexity for anysubgradient method for nonsmooth functions is Ω(ε−2). What, then, is thelower bound for minimizing the structured nonsmooth function `◦h?Unfortunately, the answer to these questions is “no”, as we show below. First, wedisprove the first conjecture by deriving the lower iteration bound for any algorithmwith access only to a subgradient oracle in the interpolation setting.Theorem 4.3.3 (Lower bound with interpolation). Given t < d and positive con-stants L and R, and an initial vector x(1). Let ` = 12(·)2. Then there exists an L-Lipschitz function h such that f := l ◦h is convex, f ∗ =minx f (x) = 0, ‖x(1)−x∗‖ ≤83R such thatmin1≤s≤tf (x(s))− f ∗ ≥ L2R22(t+1). (4.16)Proof. With out loss of generality, we assume the initial point x(1) = 0. Let h(x) =L‖x− x∗‖∞, andx∗ =[ R√t+1+ ε,R√t+1+ε2, . . . ,R√t+1+ε2t︸ ︷︷ ︸t+1 entries, 0, . . . ,0︸ ︷︷ ︸(d− t−1) entries],where ε is some arbitrary constant greater than 0. It is easy to check that ` is1-smooth, h is L-Lipschitz, f = `◦h is convex, and f ∗ = 0.We follow the proof template from Nemirovski and Yudin (1983) and assumethat any algorithm that uses only subgradient information generates iteratesx(s) ∈ span{∂ f (x(1)), . . . ,∂ f (x(t))}for all s≤ t. By our construction, we can further obtain x(s) ∈ span{e1,e2, . . . ,et}for all s≤ t. Therefore, for all s≤ t,f (x(s))− f ∗ ≥ 12(LR√t+1+Lε2s−1)2.Note that the above holds ∀ε > 0. Taking ε → 0+, we completes the proof.Theorem 4.3.3 established Ω(ε−1) iteration complexity for subgradient methodunder interpolation. Combining this result with the O(ε−1) iteration complexity inCorollary 4.3.1, we can conclude that the rateO(ε−1) is optimal in the interpolationsetting and no acceleration is possible for subgradient-based methods. Then weproceed to the lower bound without assuming interpolation.Theorem 4.3.4 (Lower bound without interpolation). Given t < d,L,R > 0 andan initial point x(1). Let `(·) = 12(·)2. Then there exist an L-Lipschitz function h84satisfying f := l ◦h is convex, ‖x(1)− x∗‖ ≤ R such thatmin1≤s≤tf (x(s))− f ∗ ≥ L2R2√t+1. (4.17)Proof. Similar to the proof of Theorem 4.3.3, we set x(1) = 0. But we change theconstruction of h(x) to h(x) = L(‖x− x∗‖∞+R), and x∗ is defined in the same wayas in the proof of Theorem 4.3.3:x∗ =[ R√t+1+ ε,R√t+1+ε2, . . . ,R√t+1+ε2t︸ ︷︷ ︸t+1 entries, 0, . . . ,0︸ ︷︷ ︸(d− t−1) entries],where ε is some arbitrary constant greater than 0. Following the same analysis ofthe proof of Theorem 4.3.3. x(s) ∈ span{e1,e2, . . . ,et} ∀s≤ t andf (x(s))− f ∗ ≥ L22(R√t+1+ε2s−1+R)2− L2R22≥ L2R2√t+1∀s≤ t,By taking ε → 0+, we have ‖x(1)− x∗‖ ≤ R and completes the proof.Theorem 4.3.4 provides an Ω(ε−2) iteration complexity for subgradient-basedmethods for solving the structured function ` ◦ h. This matches the lower boundfor solving general nonsmooth objective function with subgradient-based method,and implies that the structure ` ◦ h itself without interpolation cannot give us animproved iteration complexity of subgradient-based methods.4.4 Numerical experimentsThe smoothness of the loss function ` is crucial to our analysis. We now present somenumerical experiments to compare the convergence of SSGD for training ReLUneural networks with smooth and nonsmooth loss functions. We denote {xi}ni=1 astraining samples and {yi}ni=1 as training labels. yˆi stands for the prediction of thei-th sample xi from the trained model. Note that we use the letter x to represent themodel in previous sections, to be consistent to convention, we use it to denote datapoints in this section.854.4.1 Teacher-student setupWe randomly generate a small neural network with one hidden layer (16 neurons andReLU activation) as the teacher network. The network takes 16 dimensional vectorsas inputs and outputs a scalar. Then we generate 128 random vectors from theGaussian distribution as our training data {xi}128i=1 ⊂R16 and get their correspondinglabels {yi}128i=1 as the output of the teacher network. In order to ease the trainingand satisfy the interpolation assumption, we overparameterize the student neuralnetwork and set it to be a one hidden layer network with 512 neurons and ReLUactivation. We train the student network with different loss functions: squared losse.g., (yi− yˆi)2 and absolute loss e.g., |yi− yˆi| with different learning rates. Thetraining curves are shown in Figure 4.1a. We can observe that the training curve withsquared loss is smoother than the curve with absolute loss. For absolute loss, theperformance of SSGD is more sensitive to the change of learning rate and we needto decrease the learning rate to obtain a lower objective value. These observationsvalidate the importance of the smoothness of loss function under the interpolationsetting.4.4.2 Classify 4’s and 9’s on MNIST datasetWe train the LeNet (Lecun et al., 1998) on the MNIST dataset to classify 4’sand 9’s. To convert this task to a binary classification problem, we transformthe labels {yi}ni=1 to {−1,+1}n. Then we run SSGD to train the model withdifference loss functions: logistic loss e.g., log(1+ exp(−yiyˆi)) and L1-hinge losse.g., max{0,1− yiyˆi} and with different learning rates. The training curves arepresented in Figure 4.1b. Different from the teacher-student experiment, SSGDin this task perform similarly with smooth and nonsmooth loss functions. Weconjecture that this is because the objective function with L1-hinge loss almostsatisfies Proposition 4.3.1 locally at the solution, namely eq. 4.2 holds for mosttraining samples in a neighbourhood of the solution. Our observation supports thisconjecture. We observe that more that 95% of our final predictions yˆi’s satisfy thecondition |yˆi|> 2. Since the L1-hinge loss is locally smooth when |yˆi|> 2, we canthus say that most training samples satisfy eq. 4.2 locally at the solution. While forthe teacher-student training problem, its objective is nonsmooth at solution (when86(a) Teacher-student training (b) Classify 4’s and 9’s on MINST.Figure 4.1: The performance of SSGD with smooth and nonsmooth loss func-tions.zero residual is attained) since the absolute value function `(x) := |x| is nonsmoothat 0. Therefore the objective of the teacher-student training problem does not satisfyProposition 4.3.1 locally as solution and running SSGD to solve it could suffer fromslow convergence.4.5 DiscussionAn empirical-risk minimization problem based on composite functions has sufficientstructure to allow for a tight convergence analysis that explains the effectiveness ofstochastic subgradient descent methods on nonsmooth problems with interpolation.Surprisingly, the complexity bounds O(ε−1) and O(log(ε−1)) that we prove underinterpolation match those of stochastic gradient descent for smooth functions.87Chapter 5Conclusion and future workThis thesis contributes to a better understanding of the theory behind some widelyused first-order optimization algorithms for problems that satisfy specific structures.We must admit that the progress described in this thesis is just a small step towardsunderstanding first-order methods thoroughly. There are still a lot of interesting andimportant open problems on first-order optimization algorithms unsolved. Next, wesummarize some of the results obtained in this thesis and discuss possible futuredirections.5.1 Coordinate optimizationIn Chapter 2, we provided theoretical justification for the implicit screening func-tionality of GCD. However, our analysis only works for composite problems with1-norm regularization or non-negative constraints. Whether it is possible to extendour analysis to general regularizers that are nonsmooth at origin is left unanswered.In Chapter 2, we also established the fast convergence rate of GCD, but the ex-pensive greedy selection rule is still the bottleneck for GCD’s implementation inpractice. Developing a reliable approximate greedy selection rule that could reducethe overhead while preserving the fast convergence rate is a possible direction toexplore.In the literature of coordinate descent, it is often assumed that the objective isunconstrained or the constraint set is separable. Specialized analysis of CD with88simple nonseparable constraint set such as the linear constraint {x | ∑di=1 aixi = b}has also appeared in the literature (Tseng and Yun, 2010). Developing a generalframework to analyze CD’s convergence (if possible) for a more general nonsepara-ble constraint set is an important topic from the perspective of theory and algorithmdesign.5.2 Mirror descentIn Chapter 3, we provided a careful study of OMD when using dynamic stepsize.By modifying the OMD via a stabilization technique, we obtain the O(√T ) regretand therefore fixed the divergence issue of OMD under unknown time horizon andunbounded domain. Through the stabilization technique, we are also able to analyzethe similarities and difference between stabilized-OMD, OMD and DA.Here we post one possible future direction. We know that the theory of OMDand DA works for any mirror map that satisfies the basic assumptions listed inSection 3.1.1. However, to the author’s knowledge, the only examples that candemonstrate the advantage of OMD and DA over vanilla projected subgradientdescent (PGD) are the experts’ problem and the bandit problem (when negativeentropy mirror map is used). It is of interest to identify new applications and thecorresponding new mirror maps that could enjoy faster convergence rate with OMDand DA than PGD.5.3 Stochastic subgradient descentIn Chapter 4, we identified a convex-composite structure and develop some semi-smoothness properties of the empirical-risk minimization problem. The semi-smoothness properties allow us to derive improved convergence rates of SSGDwhen the interpolation holds.Our convergence analysis is based on the convexity properties of f . As wemention in connection with Theorem 4.3.2, the strong convexity assumption can berelaxed to weaker conditions, but even these exclude important nonconvex modelsthat appear in neural networks. It is still an open problem as to whether a linearconvergence rate for subgradient methods under a weaker assumptions, such as therestricted secant inequality (RSI). In Section 4.3.3 we proved that the rate O(ε−1)89is optimal for subgradient-based methods under interpolation. However, this lowerbound holds only for subgradient-based algorithms. Smoothing techniques basedon Moreau envelopes (Nesterov, 2005) can sometimes lead to acceleration fornonsmooth optimization, which may be further avenue to explore for obtaining anaccelerated SSGD method.90BibliographyAllen-Zhu, Z., Li, Y., and Song, Z. (2019). A convergence theory for deep learningvia over-parameterization. In Proceedings of ICML, pages 242–252. → page 78Allen Zhu, Z., Qu, Z., Richta´rik, P., and Yuan, Y. (2016). Even faster acceleratedcoordinate descent using non-uniform sampling. In Proceedings of ICML, pages1110–1119. → page 6Argyriou, A., Evgeniou, T., and Pontil, M. (2007). Multi-task feature learning. InProceesings of NeurIPS. → page 31Atamtu¯rk, A. and Go´mez, A. (2020). Safe screening rules for `0-regression. InProceesings of ICML. → page 31Auer, P., Cesa-Bianchi, N., Freund, Y., and Schapire, R. E. (1995). Gambling in arigged casino: The adversarial multi-arm bandit problem. In 36th AnnualSymposium on Foundations of Computer Science, pages 322–331. → page 25Auer, P., Cesa-Bianchi, N., Freund, Y., and Schapire, R. E. (2002a). Thenonstochastic multiarmed bandit problem. SIAM Journal on Computing, 32(1).→ page 49Auer, P., Cesa-Bianchi, N., and Gentile, C. (2002b). Adaptive and self-confidenton-line learning algorithms. Journal of Computer and System Sciences,64(1):48–75. → pages 64, 66, 115Bao, R., Gu, B., and Huang, H. (2020). Fast oscar and owl with safe screeningrules. In Proceesings of ICML. → page 31Beck, A. and Teboulle, M. (2003). Mirror descent and nonlinear projectedsubgradient methods for convex optimization. Operations Research Letters,31(3):167–175. → pages 9, 48, 66, 11891Beck, A. and Tetruashvili, L. (2013). On the convergence of block coordinatedescent type methods. SIAM Journal on Optimization, 23(4):2037–2060. →pages 6, 13, 14Bertsekas, D. P. (1976). On the goldstein-levitin-polyak gradient projection method.IEEE Transactions on automatic control, 21(2):174–184. → page 30Bertsekas, D. P. (1999). Nonlinear Programming. Athena Scientific. → page 18Bertsekas, D. P. and Tsitsiklis, J. N. (1989). Parallel and Distributed Computation:Numerical Methods. Prentice-Hall, Inc., USA. → pages 6, 72Bonnefoy, A., Emiya, V., Ralaivola, L., and Gribonval, R. (2015). Dynamicscreening: Accelerating first-order algorithms for the lasso and group-lasso.IEEE Trans. Signal Process., 63(19):5121–5132. → page 31Boyd, S. and Vandenberghe, L. (2004). Convex Optimization. CambridgeUniversity Press, New York, NY, USA. → pages 15, 107Bradley, J. K., Kyrola, A., Bickson, D., and Guestrin, C. (2011). Parallel coordinatedescent for l1-regularized loss minimization. In Proceedings of ICML, pages321–328. → page 6Bubeck, S. (2011). Introduction to online optimization. unpublished. → pages45, 46, 63Bubeck, S. (2015). Convex optimization: Algorithms and complexity. Foundationsand Trends in Machine Learning, 8(3-4):231–357. → pages 46, 61, 67, 118Candes, E. J., Eldar, Y. C., Strohmer, T., and Voroninski, V. (2015). Phase retrievalvia matrix completion. SIAM review, 57(2):225–251. → page 31Cauchy, A.-L. (1847). Methode ge´ne´rale pour la re´solution des systemes d’e´qationssimultane´es. Comptes Rendus, 25(2):536–538. → page 3Cesa-Bianchi, N. and Lugosi, G. (2006). Prediction, learning, and games.Cambridge University Press. → pages 63, 64, 66Cesa-Bianchi, N., Mansour, Y., and Stoltz, G. (2007). Improved second-orderbounds for prediction with expert advice. Machine Learning, 66(2-3):321–352.→ page 62Chandrasekaran, V., Recht, B., Parrilo, P. A., and Willsky, A. S. (2012). The convexgeometry of linear inverse problems. Foundations of Computationalmathematics, 12(6):805–849. → page 3292Chang, C.-C. and Lin, C.-J. (2011). LIBSVM: A library for support vectormachines. ACM Transactions on Intelligent Systems and Technology,2:27:1–27:27. Software available at http://www.csie.ntu.edu.tw/∼cjlin/libsvm. →pages 5, 15, 27Clarke, F. H. (1981). Generalized gradients of lipschitz functionals. Advances inMathematics, 40(1):52–67. → page 77Clarke, F. H. (1990). Optimization and nonsmooth analysis, volume 5. SIAM. →page 75Daniilidis, A., Sagastiza´bal, C., and Solodov, M. (2009). Identifying structure ofnonsmooth convex functions by the bundle technique. SIAM Journal onOptimization, 20(2):820–840. → page 30Davis, G., Mallat, S., and Avellaneda, M. (1997). Adaptive greedy approximations.Constructive approximation, 13(1):57–98. → page 26de Rooij, S., van Erven, T., Gru¨nwald, P. D., and Koolen, W. M. (2014). Follow theleader if you can, hedge if you must. Journal of Machine Learning Research(JMLR), 15:1281–1316. → pages 62, 66Defazio, A., Bach, F., and Lacoste-Julien, S. (2014). Saga: A fast incrementalgradient method with support for non-strongly convex composite objectives. InProceedings of NeurIPS, volume 27, pages 1646–1654. → pages 8, 73Dhillon, I. S., Ravikumar, P., and Tewari, A. (2011). Nearest neighbor based greedycoordinate descent. In Proceedings of NeuraIPS, pages 2160–2168. → page 16Fan, R.-E., Chang, K.-W., Hsieh, C.-J., Wang, X.-R., and Lin, C.-J. (2008).LIBLINEAR: A library for large linear classification. Journal of MachineLearning Research, 9:1871–1874. → page 14Fan, Z., Jeong, H., Sun, Y., and Friedlander, M. P. (2020). Atomic decompositionvia polar alignment: The geometry of structured optimization. Foundations andTrends in Optimization, 3(4):280–366. → pages 32, 33, 34Fang, H., Fan, Z., and Friedlander, M. P. (2021). Fast convergence of stochasticsubgradient method under interpolation. In Proceedings of ICLR. → page viFang, H., Fan, Z., Sun, Y., and Friedlander, M. (2020a). Greed meets sparsity:Understanding and improving greedy coordinate descent for sparse optimization.In Proceeding of AISTATS. → page v93Fang, H., Harvey, N. J., Portella, V. S., and Friedlander, M. P. (2020b). Onlinemirror descent and dual averaging: keeping pace in the dynamic case. InProceedings of ICML, volume 119. → page vFoucart, S. and Rauhut, H. (2013). A Mathematical Introduction to CompressiveSensing. Birkha¨user Basel. → page 26Gerchinovitz, S. (2011). Prediction of individual sequences and prediction in thestatistical framework: some links around sparse regression and aggregationtechniques. PhD thesis, Universite´ Paris-Sud. → page 63Ghadimi, S. and Lan, G. (2013). Stochastic first- and zeroth-order methods fornonconvex stochastic programming. SIAM Journal on Optimization,23(4):2341–2368. → page 8Ghaoui, L. E., Viallon, V., and Rabbani, T. (2012). Safe feature elimination insparse supervised learning. Pacific Journal of Optimization, 8(4):667–698. →pages 31, 34Glorot, X., Bordes, A., and Bengio, Y. (2011). Deep sparse rectifier neuralnetworks. In Proceedings of AISTATS, pages 315–323. → page 73Goodfellow, I., Bengio, Y., and Courville, A. (2016). Deep Learning. The MITfPress. → page 73Goyal, P., Dolla´r, P., Girshick, R., Noordhuis, P., Wesolowski, L., Kyrola, A.,Tulloch, A., Jia, Y., and He, K. (2017). Accurate, large minibatch sgd: Trainingimagenet in 1 hour. arXiv e-prints. → page 72Gunasekar, S., Woodworth, B. E., Bhojanapalli, S., Neyshabur, B., and Srebro, N.(2017). Implicit regularization in matrix factorization. In Advances in NeuralInformation Processing Systems, volume 30. → page 72Gyo¨rgy, A. and Szepesva´ri, C. (2016). Shifting regret, mirror descent, and matrices.In Proceedings of ICML’16, pages 2943–2951. → page 50Hannah, R., Feng, F., and Yin, W. (2019). A2BCD: asynchronous acceleration withoptimal complexity. In International Conference on Learning Representations,ICLR. → page 7Hare, W. L. and Lewis, A. S. (2007). Identifying active manifolds. AlgorithmicOperations Research, 2(2):75. → page 1994Hastie, T., Friedman, J. H., and Tibshirani, R. (2008). Regularization paths andcoordinate descent. In Proceedings of the 14th ACM SIGKDD InternationalConference on Knowledge Discovery and Data Mining, Las Vegas, Nevada, USA,August 24-27, 2008, page 3. → page 5Hazan, E. (2016). Introduction to online convex optimization. Foundations andTrends in Optimization, 2(3-4):157–325. → pages 45, 48, 67, 68Hsieh, C.-J., Si, S., and Dhillon, I. (2014). A divide-and-conquer solver for kernelsupport vector machines. In Proceedings of ICML, pages 566–574. → page 27Jacot, A., Hongler, C., and Gabriel, F. (2018). Neural tangent kernel: Convergenceand generalization in neural networks. In Proceedings of NeuraIPS, pages8580–8589. → page 76Jaggi, M., Smith, V., Taka´c, M., Terhorst, J., Krishnan, S., Hofmann, T., and Jordan,M. I. (2014). Communication-efficient distributed dual coordinate ascent. InProceedings of NeuraIPS, pages 3068–3076. → page 7Joachims, T. (1999). Making large-scale SVM learning practical. In Advances inKernel Methods - Support Vector Learning, chapter 11, pages 169–184. MITPress, Cambridge, MA. → pages 5, 15Johnson, R. and Zhang, T. (2013). Accelerating stochastic gradient descent usingpredictive variance reduction. In Proceedings of NeurIPS, volume 26, pages315–323. → pages 8, 73Kaczmarz, S. (1937). Angenaherte aufl ¨ osung von systemen linearer gleichungen.Bulletin International de l’Acade´mie Polonaise des Sciences et des Lettres,35:335–357. → page 71Kakade, S., Shalev-Shwartz, S., and Tewari, A. (2009). On the duality of strongconvexity and strong smoothness: Learning applications and matrixregularization. Unpublished Manuscript, http://ttic. uchicago.edu/shai/papers/KakadeShalevTewari09. pdf, 2(1). → page 35Kangarshahi, E. A., Hsieh, Y., Sahin, M. F., and Cevher, V. (2018). Let’s be honest:An optimal no-regret framework for zero-sum games. In Proceedings ofICML’18, pages 2493–2501. → page 50Karimi, H., Nutini, J., and Schmidt, M. (2016). Linear convergence of gradient andproximal-gradient methods under the polyak-łojasiewicz condition. InProceedings of ECML PKDD, pages 795–811. → page 8195Karimireddy, S. P., Koloskova, A., Stich, S. U., and Jaggi, M. (2019). Efficientgreedy coordinate descent for composite problems. In Proceedings of AISTATS,pages 2887–2896. → pages 17, 22, 25, 104, 105Kingma, D. P. and Ba, J. (2015). Adam: A method for stochastic optimization. InProceedings of ICLR. → pages 8, 71Ko, M., Zowe, J., et al. (1994). An iterative two-step algorithm for linearcomplementarity problems. Numerische Mathematik, 68(1):95–106. → page 30Krizhevsky, A., Sutskever, I., and Hinton, G. E. (2012). Imagenet classificationwith deep convolutional neural networks. In Proceeding of NeurIPS, volume 25.→ page 72Kuang, Z., Geng, S., and Page, D. (2017). A screening rule for l1-regularized isingmodel estimation. In Proceesings of NeurIPS. → page 31Lacoste-Julien, S., Schmidt, M., and Bach, F. R. (2012). A simpler approach toobtaining an o(1/t) convergence rate for the projected stochastic subgradientmethod. ArXiv, abs/1212.2002. → page 8Lai, Z. and Lim, L. (2020). Recht-re noncommutative arithmetic-geometric meanconjecture is false. In Proceedings of ICML, volume 119, pages 5608–5617. →page 14Lecun, Y., Bottou, L., Bengio, Y., and Haffner, P. (1998). Gradient-based learningapplied to document recognition. In Proceedings of the IEEE, pages 2278–2324.→ page 86Lee, C.-P. and Wright, S. J. (2018). Random permutations fix a worst case forcyclic coordinate descent. IMA Journal of Numerical Analysis, 39(3):1246–1275.→ pages 6, 13Lee, Y. T. and Sidford, A. (2013). Efficient accelerated coordinate descent methodsand faster algorithms for solving linear systems. In Proceedings of FOCS, pages147–156. → page 6Lewis, A. S. and Wright, S. J. (2011). Identifying activity. SIAM Journal onOptimization, 21(2):597–614. → page 19Li, M., Andersen, D. G., Smola, A. J., and Yu, K. (2014). Communication efficientdistributed machine learning with the parameter server. In Proceeding ofNeurIPS, volume 27. → page 7296Li, Y. and Osher, S. (2009). Coordinate descent optimization for l1 minimizationwith application to compressed sensing; a greedy algorithm. Inverse Problemsand Imaging, 3(3):487—-503. → page 18Lian, X., Huang, Y., Li, Y., and Liu, J. (2015). Asynchronous parallel stochasticgradient for nonconvex optimization. In Proceeding of NeurIPS, volume 28. →page 72Lin, Q., Lu, Z., and Xiao, L. (2015). An accelerated randomized proximalcoordinate gradient method and its application to regularized empirical riskminimization. SIAM Journal on Optimization, 25(4):2244–2273. → page 6Liu, C. and Belkin, M. (2020). Accelerating SGD with momentum forover-parameterized learning. In Proceedings of ICLR. → page 83Liu, J., Wright, S. J., Re´, C., Bittorf, V., and Sridhar, S. (2014a). An asynchronousparallel stochastic coordinate descent algorithm. In Proceedings of ICML, pages469–477. → page 7Liu, J., Wright, S. J., Re´, C., Bittorf, V., and Sridhar, S. (2015). An asynchronousparallel stochastic coordinate descent algorithm. Journal of Machine LearningResearch, 16(10):285–322. → page 72Liu, J., Zhao, Z., Wang, J., and Ye, J. (2014b). Safe screening with variationalinequalities and its application to lasso. In Proceesings of ICML. → page 31Lloyd, S. P. (1982). Least squares quantization in PCM. IEEE Trans. InformationTheory, 28(2):129–136. → page 5Luo, Z. and Tseng, P. (1992). On the convergence of the coordinate descent methodfor convex differentiable minimization. Journal of Optimization Theory andApplications, 72(1):361–379. → page 6Luo, Z.-Q. and Tseng, P. (1993). Error bounds and convergence analysis of feasibledescent methods: a general approach. Annals of Operations Research,46(1):157–178. → page 6Ma, S., Bassily, R., and Belkin, M. (2018). The power of interpolation:Understanding the effectiveness of SGD in modern over-parametrized learning.In Proceedings of ICML, pages 3331–3340. → pages 73, 82McMahan, H. B. (2017). A survey of algorithms and analysis for adaptive onlinelearning. Journal of Machine Learning Research, 18:90:1–90:50. → pages45, 67, 11597Meinshausen, N. and Bu¨hlmann, P. (2006). High dimensional graphs and variableselection with the lasso. The Annals of Statistics, 34. → page 31Ndiaye, E., Fercoq, O., Gramfort, A., and Salmon, J. (2017). Gap safe screeningrules for sparsity enforcing penalties. Journal of Machine Learning Research,18:128:1–128:33. → pages 10, 31, 34, 35Needell, D., Ward, R., and Srebro, N. (2014). Stochastic gradient descent, weightedsampling, and the randomized kaczmarz algorithm. In Proceeding of NeurIPS,volume 27. → page 71Negahban, S. N. and Wainwright, M. J. (2012). Restricted strong convexity andweighted matrix completion: Optimal bounds with noise. Journal of MachineLearning Research, 13:1665–1697. → page 21Nemirovski, A., Juditsky, A. B., Lan, G., and Shapiro, A. (2009). Robust stochasticapproximation approach to stochastic programming. SIAM Journal onOptimization, 19(4):1574–1609. → page 8Nemirovski, A. and Yudin, D. (1983). Problem Complexity and Method Efficiencyin Optimization. Wiley Interscience. → pages 4, 9, 82, 84Nesterov, Y. (2009). Primal-dual subgradient methods for convex problems.Mathematical Programming, 120(1):221–259. → pages 45, 48, 61, 66Nesterov, Y. (2012). Efficiency of coordinate descent methods on huge-scaleoptimization problems. SIAM Journal on Optimization, 22(2):341–362. → pages6, 12, 13, 22, 25Nesterov, Y. (2014). Introductory Lectures on Convex Optimization: A BasicCourse. Springer Publishing Company, Incorporated, 1 edition. → page 78Nesterov, Y. E. (1983). A method for solving the convex programming problemwith convergence rate o(1/k2). Dokl. Akad. Nauk SSSR, 269:543–547. → page 4Nesterov, Y. E. (2004). Introductory Lectures on Convex Optimization - A BasicCourse, volume 87 of Applied Optimization. Springer. → pages 4, 16Nesterov, Y. E. (2005). Smooth minimization of non-smooth functions.Mathematical Programming, 103(1):127–152. → page 90Nguyen, L. M., Liu, J., Scheinberg, K., and Taka´cˇ, M. (2017). SARAH: A novelmethod for machine learning problems using stochastic recursive gradient. InProceeding of ICML, volume 70, pages 2613–2621. → pages 8, 7398Nutini, J., Schmidt, M., and Hare, W. (2017). ”active-set complexity” of proximalgradient: How long does it take to find the sparsity pattern? Optimization Letter.→ pages 19, 109Nutini, J., Schmidt, M. W., Laradji, I. H., Friedlander, M. P., and Koepke, H. A.(2015). Coordinate descent converges faster with the gauss-southwell rule thanrandom selection. In Proceedings of the International Conference on MachineLearning, pages 1632–1641. → pages 6, 12, 15, 16, 17, 18, 19, 22, 25, 104Orabona, F. and Pa´l, D. (2018). Scale-free online learning. Theor. Comput. Sci.,716:50–69. → pages 45, 48, 66, 68Oustry, F. (2000). A second-order bundle method to minimize the maximumeigenvalue function. Mathematical Programming, 89(1):1–33. → page 42Pati, Y. C., Rezaiifar, R., and Krishnaprasad, P. S. (1993). Orthogonal matchingpursuit: Recursive function approximation with applications to waveletdecomposition. In Proceedings of 27th Asilomar conference on signals, systemsand computers, pages 40–44. IEEE. → page 26Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O.,Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A.,Cournapeau, D., Brucher, M., Perrot, M., and Duchesnay, E. (2011). Scikit-learn:Machine learning in Python. J. Mach. Learn. Res., 12:2825–2830. → page 27Peng, Z., Wu, T., Xu, Y., Yan, M., and Yin, W. (2016). Coordinate friendlystructures, algorithms and applications. Annals of Mathematical Sciences andApplications, 1:57–119. → page 12Platt, J. (1998). Sequential minimal optimization: A fast algorithm for trainingsupport vector machines. Technical Report MSR-TR-98-14. → pages 5, 12Qian, X., Richta´rik, P., Gower, R. M., Sailanbayev, A., Loizou, N., and Shulgin, E.(2019). SGD with arbitrary sampling: General analysis and improved rates. InProceedings of ICML, pages 5200–5209. → page 81Raj, A., Olbrich, J., Ga¨rtner, B., Scho¨lkopf, B., and Jaggi, M. (2015). Screeningrules for convex problems. ArXiv, abs/1609.07478. → page 31Recht, B. and Re´, C. (2012). Toward a noncommutative arithmetic-geometric meaninequality: Conjectures, case-studies, and consequences. In The 25th AnnualConference on Learning Theory, pages 11.1–11.24. → pages 6, 1499Recht, B., Re, C., Wright, S., and Niu, F. (2011). Hogwild!: A lock-free approachto parallelizing stochastic gradient descent. In Proceeding of NeurIPS,volume 24. → page 72Richta´rik, P. and Taka´c, M. (2011). Efficient serial and parallel coordinate descentmethods for huge-scale truss topology design. In Proceedings of theInternational Conference on Operations Research, pages 27–32. → page 7Richta´rik, P. and Taka´c, M. (2014). Iteration complexity of randomizedblock-coordinate descent methods for minimizing a composite function.Mathematical Programming, 144(1-2):1–38. → pages 6, 22, 25Robbins, H. and Monro, S. (1951). A stochastic approximation method. Annals ofMathematical Statistics, 22(3):400–407. → page 7Rockafellar, R. T. (1970). Convex Analysis. Princeton University Press, Princeton.→ pages 45, 46, 121Rockafellar, R. T. and Wets, R. J.-B. (2009). Variational analysis, volume 317.Springer Science & Business Media. → pages 36, 114Sa, C. D. (2020). Random reshuffling is not always better. In Proceedings ofNeurIPS. → page 14Saha, A. and Tewari, A. (2013). On the nonasymptotic convergence of cycliccoordinate descent methods. SIAM Journal on Optimization, 23(1):576–601. →pages 6, 13Schmidt, M. and Le Roux, N. (2013). Fast convergence of stochastic gradientdescent under a strong growth condition. arXiv e-prints. → pages 80, 82Schmidt, M., Roux, N. L., and Bach, F. R. (2017). Minimizing finite sums with thestochastic average gradient. Mathematical Programming, 162(1-2):83–112. →pages 8, 73Schulman, J., Levine, S., Abbeel, P., Jordan, M., and Moritz, P. (2015). Trust regionpolicy optimization. In Proceedings of ICML, volume 37, pages 1889–1897. →page 71Schulman, J., Wolski, F., Dhariwal, P., Radford, A., and Klimov, O. (2017).Proximal policy optimization algorithms. arXiv e-prints, page arXiv:1707.06347.→ page 72100Shalev-Shwartz, S. (2012). Online learning and online convex optimization.Foundations and Trends in Machine Learning, 4(2):107–194. → pages45, 46, 60, 69Shalev-Shwartz, S., Singer, Y., and Srebro, N. (2007). Pegasos: Primal estimatedsub-gradient solver for svm. In Proceedings of ICML, page 807–814. → page 71Shevade, S. K. and Keerthi, S. S. (2003). A simple and efficient algorithm for geneselection using sparse logistic regression. Bioinformatics, 19(17):2246–2253. →pages 18, 26Shor, N. Z. (1984). Minimization methods for non-differentiable functions.Springer Series in Computational Mathematics, Springer. → page 4Stich, S. U. (2019). Unified optimal analysis of the (stochastic) gradient method.arXiv e-prints, page arXiv:1907.04232. → page 82Sun, R. and Ye, Y. (2021). Worst-case complexity of cyclic coordinate descent:o(n2) gap with randomized version. Mathematical Programming,185(1):487–520. → pages 6, 13Sun, Y. and Bach, F. R. (2020). Safe screening for the generalized conditionalgradient method. ArXiv, abs/2002.09718. → page 34Sun, Y., Jeong, H., Nutini, J., and Schmidt, M. W. (2019). Are we there yet?manifold identification of gradient-related proximal methods. In Proceedings ofAISTATS, pages 1110–1119. → page 19Sylvain Sardy, A. G. B. and Tseng, P. (2000). Block coordinate relaxation methodsfor nonparametric wavelet denoising. Journal of Computational and GraphicalStatistics, 9(2):361–379. → pages 5, 12Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. Journal ofthe Royal Statistical Society (Series B), 58:267–288. → page 31Tieleman, T. and Hinton, G. (2012). Lecture 6.5—RmsProp: Divide the gradient bya running average of its recent magnitude. COURSERA: Neural Networks forMachine Learning. → page 71Tseng, P. (2001). Convergence of a block coordinate descent method fornondifferentiable minimization. Journal of Optimization Theory andApplications, 109(3):475–494. → page 6101Tseng, P. and Yun, S. (2010). A coordinate gradient descent method for linearlyconstrained smooth optimization and support vector machines training.Computational Optimization and Applications, 47. → page 89Vaswani, S., Bach, F., and Schmidt, M. (2019). Fast and faster convergence of SGDfor over-parameterized models and an accelerated perceptron. In Proceedings ofAISTATS, pages 1195–1204. → page 83Wang, J., Zhou, J., Liu, J., Wonka, P., and Ye, J. (2014). A safe screening rule forsparse logistic regression. In Proceesing of NeurIPS, pages 1053–1061. → page31Wang, J., Zhou, J., Wonka, P., and Ye, J. (2013). Lasso screening rules via dualpolytope projection. In Advances in Neural Information Processing Systems 26,pages 1070–1078. → page 31Wright, S. J. (2012). Accelerated block-coordinate relaxation for regularizedoptimization. SIAM Journal on Optimization, 22(1):159–186. → page 30Wu, T. T. and Lange, K. (2008). Coordinate descent algorithms for lasso penalizedregression. The Annals of Applied Statistics, 2(1):224—-244. → page 18Xiang, Z. J., Wang, Y., and Ramadge, P. J. (2017). Screening tests for lassoproblems. IEEE Trans. Pattern Anal. Mach. Intell., 39(5):1008–1027. → page 31Yaroshinsky, R., El-Yaniv, R., and Seiden, S. S. (2004). How to better use expertadvice. Machine Learning, 55(3):271–309. → page 66You, Y., Lian, X., Liu, J., Yu, H., Dhillon, I. S., Demmel, J., and Hsieh, C. (2016).Asynchronous parallel greedy coordinate descent. In Proceedings of NeuraIPS,pages 4682–4690. → page 7You, Y., Zhang, Z., Hsieh, C.-J., Demmel, J., and Keutzer, K. (2018). Imagenettraining in minutes. In Proceedings of the International Conference on ParallelProcessing, page 1–10. → page 72Yuan, M. and Lin, Y. (2006). Model selection and estimation in regression withgrouped variables. Journal of the Royal Statistical Society (Series B), 68:49–67.→ page 31Zhang, J., Lin, H., Sra, S., and Jadbabaie, A. (2020). On complexity of findingstationary points of nonsmooth nonconvex functions. In Proceedings of ICML.→ page 75102Zhang, W., Hong, B., Liu, W., Ye, J., Cai, D., He, X., and Wang, J. (2017). Scalingup sparse support vector machines by simultaneous feature and sample reduction.In Proceesings of ICML. → page 31Zhang, X. (2013). Bregman divergence and mirror descent lecture notes. Availableat http://users.cecs.anu.edu.au/∼xzhang/teaching/bregman.pdf. Last accessedJuly. 10, 2021. → pages 117, 118Zinkevich, M. (2003). Online convex programming and generalized infinitesimalgradient ascent. In Proceedings of ICML, pages 928–936. → page 48103Appendix AAppendix for Chapter 2A.1 Proofs for Section 2.3PreliminariesWe introduce some notations and Lemmas that appear in works from Nutini et al.(2015) and Karimireddy et al. (2019).We say that a coordinate gradient step is bad if the iterate crosses origin or endat origin, i.e., x(t+1)i x(t)i < 0 or x(t+1)i = 0, otherwise we call this step a good step,see more details in Karimireddy et al. (2019). We denote the set of good steps untilthe t-th iteration as Gt . Because a bad step always follows by a good step, it is easyto verify that|Gt | ≤⌈ t2⌉. (A.1)Recall the selection rule in Section 2.3:Selection rule A.1 (GS-s rule). Select i ∈ argmax j Q j(x(t)) whereQi(x) = mins∈∂gi|∇i f (x)+ s|. (A.2)Lemma A.1.1 (Karimireddy et al., 2019, Theorem 1). Assume f is µ1-stronglyconvex with respect to 1-norm, then the iterates generated from Algorithm 1 with104GS-s rule (selection rule A.2) satisfyF(x(t))−F(x∗)≤(1− µ1Lmax)dt/2e(F(0)−F(x∗)) .Lemma A.1.2 (Karimireddy et al., 2019, Lemma 2). Consider g to be 1-normregularization or non-negative constraint. If the t-th iteration is a good step, thenwe haveF(x(t+1))≤ F(x(t))− 12Lmaxmaxi∈[d]Qi(x(t))2, (A.3)where Qi is defined in the GS-s rule (selection rule A.2).The above two Lemmas due to the work from Karimireddy et al. (2019). The sec-ond Lemma comes from the ‘proof sketch’ of (Karimireddy et al., 2019, Lemma 2,page 5).Proof of Lemma 2.3.1Proof. If i is not select by Algorithm 1 at the t-th iteration, then x(t+1)i = 0 triviallyremains at 0.If i is selected at the t-th iteration, by assuming |∇i f (x(t))−∇i f (x∗)| ≤ δi, weknow that−δi+∇i f (x∗)≤ ∇i f (x(t))≤ δi+∇i f (x∗)(i)⇒−ui ≤ ∇i f (x(t))≤−li, (A.4)where (i) follows directly from the definition of δi :=min{−∇i f (x∗)− li,ui+∇i f (x∗)}.Next we show that proxgi/Li(0− 1Li∇i f (x(t)))= 0. By the definition of theproximal operatorproxgi/Li(0− 1Li∇i f (x(t)))= argminy∈R{12(y−(− 1Li∇ fi(x(t))))2+1Ligi(y)}.This minimization problem is strongly convex and thus attains a unique solution105that satisfies0 ∈ y+ 1Li∇i f (x(t))+1Li∂gi(y). (A.5)Knowing that −ui ≤ ∇i f (x(t))≤−li from (A.4) and int∂gi(0) = (li,ui) by thedefinition of li and ui. We can easily conclude that y= 0 satisfies (A.5) and thereforex(t+1)i = proxgi/Li(0− 1Li∇i f (x(t)))= 0.Proof of Theorem 2.3.1Proof. Let t ≤ d− τ and recall the definition of good steps until the t-th iterationfrom section A.1: let Gt = {i1, i2, . . . , ik}, where k = |Gt | ≥ dt/2e.At iteration im for m ∈ [k], we know x(im) is guaranteed to be m− 1–sparsebecause the number of non-zeros of the iterate at most can increase by one forone good step and will not increase for bad steps. By assuming f is µ(τ+m−1)1strongly convex w.r.t. 1-norm and τ+m−1–sparse vectors, we know that F is alsoµ(τ+m−1)1 strongly convex w.r.t. 1-norm and τ+m−1–sparse vectors. Moreover|supp(y)∪ supp(x(im))| ≤ τ +m− 1 is true ∀y ∈ Rd that is τ-sparse. Then by thedefinition of µ(τ+m−1)1 , we haveF(y)≥ F(x(im))+ 〈∂F(x(im)),y− x(im)〉+ µ(τ+m−1)12‖y− x(im)‖21 (A.6)for any τ-sparse vectors y (with a little bit abuse of notation, we use ∂F(x(t)) todenote any vectors in the subdifferential of F(x(t))). Minimize both sides of (A.6)106w.r.t. y that is τ-sparse, we getF(x∗)≥ F(x(im))− sup‖y‖0≤τ(〈−∂F(x(im)),y− x(im)〉− µ(τ+m−1)12‖y− x(im)‖21)≥ F(x(im))− supy∈Rd(〈−∂F(x(im)),y− x(im)〉− µ(τ+m−1)12‖y− x(im)‖21)(i)= F(x(im))−(µ(τ+m−1)12‖ · ‖21)∗(−∂F(x(im)))(ii)= F(x(im))− 12µ(τ+m−1)1‖∂F(x(im))‖2∞,where (i) is from the definition of conjugate function, and (ii) is from the fact that(12‖ · ‖21)∗= 12‖ · ‖2∞ (Boyd and Vandenberghe, 2004).More specifically,F(x∗)≥ F(x(im))− 12µ(τ+m−1)1‖∇ f (x(im))+u‖2∞ ∀u ∈ ∂g(x(im)).By the definition of Qi(·) in the GS-s rule (selection rule 2.3), we further haveF(x∗)≥ F(x(im))− 12µ(τ+m−1)1maxi∈[d]Qi(x(im))2. (A.7)Recall Lemma A.1.2, we haveF(x(im+1))≤ F(x(im))− 12Lmaxmaxi∈[d]Qi(x(im))2.Plug the above equation into (A.7)F(x∗)≥ F(x(im))− Lmaxµ(τ+m−1)1(F(x(im+1))−F(x(im)))⇒ F(x(im+1))−F∗ ≤(1− µ(τ+m)1Lmax)(F(x(im))−F∗).107By applying the above inequality recursively, we getF(x(t))−F∗ ≤k∏m=1(1− µ(τ+m−1)1Lmax)(F(0)−F∗)≤d t2e∏i=1(1− µ(τ+i−1)1Lmax)(F(0)−F∗),which completes the proof.Proof of Theorem 2.3.3Proof. This proof is essentially the same as Theorem 2.3.1. The difference is that,by the definition of the ∆-GS-s rule (selection rule 2.2), the Lemma A.1.2 becomesF(x(t+1))−F(x(t))≤− ∆2Lmaxmaxi∈[d]Qi(x(t))2at each good step t.Knowing that supp(x(t))⊂W∆, we have |supp(x∗)∪ supp(x(t))| ≤ |W∆| ∀t > 0.Then we can incorporate the new Lemma into the analysis of Theorem 2.3.1 and getF(x(t))−F∗ ≤(1− ∆µ(|W∆|)1Lmax)d t2 e(F(0)−F∗)≤(1− ∆µ2|W∆|Lmax)d t2 e(F(0)−F∗) .Proposition A.1.1. Letx˜( j) := argminsupp(x)⊆Wjf (x)+g(x),and Qi defined as (A.2) (the GS-s rule). Then Qi’s are continuous at x˜( j) on thesupport Wj for all i, j ∈ [d].108Proof. Given i, j ∈ [d]. By the optimality condition, we know thatQi(x˜( j)) = 0 ∀i ∈Wj.We consider 2 cases — i∈Wj and i /∈Wj. For i∈Wj, following Nutini et al. (2017)’sanalysis, we know that as x→ x˜( j), Qi(x)→ 0 since x˜( j) is the optimal solution onthe support Wj.For i /∈Wj, denote ∂gi(0) = [li,ui], then Qi(x) =min{|∇i f (x)+δ | | δ ∈ [li,ui]}since xi = 0 (supp(x) = Wj). We know that ∇i f (x)’s are continuous functions,therefore it is easy to conclude that Qi(x) is continuous on Wj for i /∈Wj.To sum up, we finish the proof.Proof of Theorem 2.3.4Proof. Preliminaries:Given ∆> 0, we sort the elements of W∆ = {i1, i2, ..., im} by the time they firstenter the working set W∆ ,i.e., i1 is the first coordinate being selected and i2 is thesecond coordinate being included in W∆, etc.We denote the t-th iterate from the ∆-GCD algorithm as x(t) and the t-th iteratefrom the totally corrective greedy algorithm (TCGA) as x˜(t). Let W ]= {i˜1, i˜2, . . . , i˜k},its elements are also sorted by the time when they enter the working set.A claim: for any j ≤ k, there exist ε j > 0 such that ∀∆< ε j, the first j elementsin W∆ is the same as the first j elements in W ].Proof. We prove this claim by induction, when j = 1, ∀∆ ≤ 1, ∆-GCD and theTCGA both select the coordinate argmaxi∈[d]Qi(0) in the first iteration, thus theclaim is trivially true in the base case.Assume that the claim is true with some j > 0, then for j+1:By Proposition A.1.1, we know that Qi is continuous on Wj at x˜( j), further weassumed that argmaxi∈[d]Qi(x˜( j)) is singleton. Therefore there exist ε ′ > 0 suchthat ∀‖x− x˜( j)‖ ≤ ε ′,argmaxi∈[d]Qi(x) = i˜ j+1.By the uniqueness (recall that F is strongly convex) of x˜( j):x˜( j) := argminsupp(x)⊆W jf (x)+g(x)109and the optimality condition, we also know that ∃ δ > 0 such that ∀x ∈ Rd satisfysupp(x)⊆Wj and maxi∈Wj Qi(x)≤ δ , we have ‖x− x( j)‖ ≤ ε ′.Denote Qi(x(t)) ( recall x(t) is generated from ∆-GCD) is bounded by someconstant B ∀t > 0.Then, by setting ∆≤ (min{ε j,δ/B})2, when i j+1 first enter W∆ at some iterationt, we haveargmaxi∈WjQi(x(t))≤√∆argmaxi∈[d]Qi(x(t))≤ δBB = δ ,also by the induction assumption, we know that supp(x(t))⊆Wj. Putting these twoconditions together, we get ‖x(t)− x j‖ ≤ ε ′ and thus argmaxi∈[d]Qi(x(t)) = i˜ j+1,which implies that i j+1 = i˜ j+1. This completes the proof of this claim.Back to the proof:Following the claim, we know that ∃ εk > 0 such that for ∀∆ < εk, the first kelements in W∆ is just W ].By the nondegeneracy assumption i.e., δi > 0 ∀x∗i = 0 and continuity of Qi,∇ f ,we know that ∃ ε ′′ > 0 such that ∀‖x− x∗‖ < ε ′′ (note that x˜(k) = x∗), |∇i f (x)−∇i f (x∗)| ≤ δi ∀x∗i = 0 and this further implies Qi(x)= 0 ∀i /∈W ] (note that supp(x∗)∈W ]).Again, there exist δ ′′> 0 such that ∀x∈Rd satisfy suppW ](x) and maxi∈W ] Qi(x)≤δ ′′, we have ‖x− x∗‖ ≤ ε ′′.Thus for ∆ ≤ min{εk,δ ′′}, the first k elements in W∆ will be W ], and anycoordinate i /∈W ] can not be included in W∆. Therefore W∆ =W ].Proof of Theorem 2.3.2Proof. Given the number of iteration t, denote Zt = {i ∈ [d] | x(t′)i = 0 ∀t ′ <t}, which is the entries of x(t) that filled with 0’s. and Vt = {i ∈ [d] | x∗i =0 and |∇i f (x(t ′))−∇i f (x∗)| ≤ δi ∀t ′ ≥ t}.From Lemma 2.3.1 (in the main text), we know that any coordinates in Zt ∩Vt110will always stay at 0 and thus cannot be in W , that isW ⊂ [d]\(Zt ∩Vt) ∀t > 0⇒|W | ≤mint∈[d]{d−|Zt ∩Vt |} . (A.8)Recall the definition of the set of good steps until the t-th iteration Gt ⊂ [t]:|Vt |=d∑i=11{x∗i = 0 and |∇i f (x(t′))−∇i f (x∗)| ≤ δi ∀t ′ ≥ t}≥d∑i=11{‖∇ f (x(t ′))−∇ f (x∗)‖∞ ≤ δi ∀t ′ ≥ t}− τ(i)≥d∑i=11{L∞‖x(t ′)− x∗‖1 ≤ δi ∀t ′ ≥ t}− τ≥d∑i=11{L∞ supt ′≥t‖x(t)− x∗‖1 ≤ δi}− τ, (A.9)where (i) follows from the definition ∞-norm smoothness.By the definition of Gt in section A.1, we also have |Zt | ≥ d−|Gt |, and further|Zt ∩Vt |= |Zt |+ |Vt |− |Zt ∪Vt |≥ d−|Gt |+ |Vt |−d≥ |Vt |− |Gt |. (A.10)Plug the above result in (A.8), we get|W | ≤mint>0{d−|Vt |+ |Gt |}≤mint>0{d+ τ−d∑i=11{L∞ supt ′≥t‖x(t ′)− x∗‖1 ≤ δi}+ |Gt |}≤mint∈[d]{d+ τ−d∑i=11{L∞ supt ′≥t‖x(t ′)− x∗‖1 ≤ δi}+ t}= mint∈[d]Bt + t, (A.11)111where Bt is defined as Bt := d+τ− pδ(L∞ supi≥t{‖x(i)− x∗‖1}) in Theorem 2.3.2.Proof of Corollary 2.3.1Proof. Similar to the proof of Theorem 2.3.2, denoteZt = {i∈ [d] | x(t′)i = 0 ∀t ′< t},which is the entries of x(t) that filled with 0’s. and Vt = {i ∈ [d] | |∇i f (x(t ′))−∇i f (x∗)| ≤ δi ∀t ′ ≥ t}.From Lemma 2.3.1 (in the main text), we know that any coordinates in Zt ∩Vtwill always stay at 0 and thus cannot be in W , that isW ⊂ [d]\(Zt ∩Vt) ∀t > 0⇒|W | ≤mint∈[d]{d−|Zt ∩Vt |} . (A.12)Recall the definition of the set of good steps until the t-th iteration Gt ⊂ [t].|Vt |=d∑i=11{x∗i = 0 and |∇i f (x(t′))−∇i f (x∗)| ≤ δi ∀t ′ ≥ t}≥d∑i=11{‖∇ f (x(t ′))−∇ f (x∗)‖∞ ≤ δi ∀t ′ ≥ t}− τ(i)≥d∑i=11{L∞‖x(t ′)− x∗‖1 ≤ δi ∀t ′ ≥ t}− τ(ii)≥d∑i=11{L∞√2µ1(F(x(t))−F(x∗))≤ δi ∀t ′ ≥ t}− τ(iii)=d∑i=11{L∞√2µ1(F(x(t))−F(x∗))≤ δi}− τ(iv)= pδ(L∞√2µ1(F(x(t))−F(x∗)))− τ(v)≥ pδL∞√√√√ 2µ1|Gt |∏i=1(1− µ(τ+i−1)1L)(F(0)−F∗)− τ, (A.13)where (i) follows from the ∞-norm smoothness assumption, (ii) is from µ1 strongly112convex, (iii) is true since F(x(t)) is a decreasing sequence, (iv) is by the definitionof pδ (·), (v) directly follows from Theorem 2.3.1.By the definition of Gt , we also have |Zt | ≥ d−|Gt |, and further|Zt ∩Vt |= |Zt |+ |Vt |− |Zt ∪Vt |≥ d−|Gt |+ |Vt |−d≥ |Vt |− |Gt |. (A.14)Plug the above result in (A.12), we get|W | ≤mint>0{d−|Vt |+ |Gt |}≤mint>0d+ τ−L∞√√√√ 2µ1|Gt |∏i=1(1− µ(τ+i−1)1L)(F(0)−F∗)+ |Gt |≤mint∈[d]d+ τ−L∞√√√√ 2µ1t∏i=1(1− µ(τ+i−1)1L)(F(0)−F∗)+ t= mint∈[d]Bt + t, (A.15)where Bt is defined as Bt := d+τ− pδ(√2L2∞µ1 ∏t−1i=0(1− µ(τ+i)1L)(F(0)−F∗))inTheorem 2.3.2.A.2 Proofs for Section 2.4Proof of Proposition 2.4.1Proof. From the construction of sets A(t), it is straightforward to see thatA(1) ⊇A(2) ⊇ . . . ,113which shows that {A(t) }∞t=1 is a monotone sequence. By Rockafellar and Wets(2009, Exercise 4.3), the Painlevee´-Kuratowski set limitA(∞) = limt→∞A(t)is well-defined.First, we show that FA(M∗y∗) ⊆ A(∞). By Theorem 2.4.1, we know thatFA(M∗y∗)⊆A(t) for all t. Therefore, it follows that FA(M∗y∗)⊆A(∞).Next, we show that A(∞) ⊆FA(M∗y∗). Consider a ∈ A(∞). Since {A(t) }∞t=1 isa monotone sequence, there exist T > 0 such thata ∈ A(t), ∀t ≥ T.By the construction of A(t), we know that A(t) ⊆FA(M∗y(t),ε) for all t, and thuswe can conclude that〈a,M∗y(t)〉 ≥ σA(M∗y(t))− ε, ∀t ≥ T.Now by taking limits with respect to t to both sides of the inequality, we canconclude that〈a,M∗y∗〉 ≥ σA(M∗y∗),which implies that a ∈ FA(M∗y∗).114Appendix BAppendix for Chapter 3B.1 Standard factsB.1.1 Scalar inequalitiesFact B.1.1. For any a> 0 and b,x ∈ R, we have −ax2+bx≤ b2/4a.Fact B.1.2. e−x ≤ 1− x+ x22 for x≥ 0.Fact B.1.3. ∑ti=11√i≤ 2√t−1 for t ≥ 1.Fact B.1.4. log(x)≤ x−1 for x≥ 0.The following proposition is a variant of an inequality that is frequently used inonline learning; see, e.g., Auer et al. (2002b, Lemma 3.5), McMahan (2017, Lemma4).Proposition B.1.1. Let u> 0 and a1,a2, . . . ,aT ∈ [0,u]. ThenT∑t=1at√u+∑i<t ai≤ 2√T∑t=1at .Although it is easy to prove this inequality by induction, the following proof mayprovide more intuition. The proof is based on a generic lemma on approximatingsums by integrals.115Lemma B.1.1 (Sums with chain rule). Let S⊆ R be an interval. Let F : S→ R beconcave and differentiable on the interior of S. Let u≥ 0 and let A : {0,1, . . . ,T}→ Ssatisfy A(i)−A(i−1) ∈ [0,u] for each 1≤ i≤ T . ThenT∑i=1F ′(u+A(i−1)) · (A(i)−A(i−1)) ≤ F(A(T ))−F(A(0)).As u→ 0, the left-hand side becomes comparable to ∫ T0 F ′(A(x))A′(x)dx, anexpression that has no formal meaning since A is only defined on integers. If thisexpression existed, it would equal the right-hand side by the chain rule.Proof of Lemma B.1.1. Since F is concave, f := F ′ is non-increasing. Fix any1≤ i≤ T and observe that f (x)≥ f (A(i))≥ f (u+A(i−1)) for all x≤ A(i). Thusf (u+A(i−1)) · (A(i)−A(i−1)) ≤∫ A(i)A(i−1)f (x)dx = F(A(i))−F(A(i−1)).Summing over i, the right-hand side telescopes, which yields the result.Proof of Proposition B.1.1. Apply Lemma B.1.1 with S = R≥0, F(x) = 2√x andA(i) = ∑1≤ j≤i a j.Proposition B.1.2. Let x,y,α,β > 0.If x− y ≤ α√x+βthen x− y ≤ α√y+β +α√β +α2.Proof. The proposition’s hypothesis yieldsy+β +α24≥ x−α√x+ α24=(√x− α2)2.Taking the square root and rearranging,√x ≤√y+β +α24+α2.116Squaring both sides and rearranging,x ≤ y+α√y+β +α24+β +α22≤ y+α√y+α√β +β +α2,by subadditivity of the square root.B.1.2 Bregman divergence propertiesThe following lemma collects basic facts regarding the Bregman divergence inducedby a mirror map of the Legendre type. See Zhang (2013).Lemma B.1.2. The Bregman divergence induced by Φ satisfies the following prop-erties:• DΦ(x,y) is convex in x;• ∇Φ(∇Φ∗(z)) = z and ∇∗Φ(∇Φ(x)) = x for all x and z;• DΦ(x,y) = DΦ∗(∇Φ(y),∇Φ(x)) for all x and y.Proposition B.1.3. If Φ is ρ-strongly convex with respect to ‖ · ‖ then DΦ(x,y)≥ρ2 ‖x− y‖2.Differences of Bregman divergencesRecall that in (3.2) we defined the notationDΦ(ab ;c) := DΦ(a,c)−DΦ(b,c) = Φ(a)−Φ(b)−〈∇Φ(c), a−b〉.This has several useful properties, which we now discuss.Proposition B.1.4. DΦ(ab ; p) is linear in pˆ. In particular,DΦ(ab ;∇Φ∗(pˆ− qˆ)) = DΦ(ab ; p)+ 〈 qˆ, a−b〉 ∀qˆ ∈ Rd .Proof. Immediate from the definition.117Proposition B.1.5. For all a,b,c,d ∈ D,DΦ(ab ;d)−DΦ(ab ;c) = 〈 cˆ− dˆ, a−b〉 = DΦ(ab ;d)+DΦ(ba ;c).Proof. The first equality holds from Proposition B.1.4 with pˆ = cˆ and qˆ = cˆ− dˆ.The second equality holds since DΦ(ba ;c) =−DΦ(ab ;c).An immediate consequence is the “generalized triangle inequality for Bregmandivergence”. See Bubeck (2015, eq. (4.1)), Beck and Teboulle (2003, Lemma 4.1)or Zhang (2013, eq. (3)).Proposition B.1.6. For all a,b,d ∈ D,DΦ(a,d)−DΦ(b,d)+DΦ(b,a) = 〈 aˆ− dˆ, a−b〉Proof. Apply Proposition B.1.5 with c = a and use DΦ(a,a) = 0.Proposition B.1.7. Let a,b,c,u,v ∈ Rd satisfy γ aˆ+(1− γ)bˆ = cˆ for some γ ∈ R.ThenγDΦ(uv ;a) + (1− γ)DΦ(uv ;b) = DΦ(uv ;c).Proof. By definition of DΦ, the claimed identity is equivalent to(1− γ)(Φ(u)−Φ(v)−〈∇Φ(a), u− v〉) + γ(Φ(u)−Φ(v)−〈∇Φ(b), u− v〉)=(Φ(u)−Φ(v)−〈∇Φ(c), u− v〉).This equality holds by canceling Φ(u)−Φ(v) and by the assumption that ∇Φ(c) =(1− γ)∇Φ(a)+ γ∇Φ(b).The following proposition is the “Pythagorean theorem for Bregman diver-gence”. Recall that ΠΦX (y) = argminu∈XDΦ(u,y). Proofs may be found in Bubeck(2015, Lemma 4.1) or Zhang (2013, eq. (17)).Proposition B.1.8. Let X ⊂Rd be a convex set. Let p ∈Rd and pi =ΠΦX (p). ThenDΦ( zpi ; p) ≥ DΦ( zpi ;pi) = DΦ(z,pi) ∀z ∈ X .118A generalization of the previous proposition can be obtained by using thelinearity property.Proposition B.1.9. Let X ⊂Rd be a convex set. Let p ∈Rd and pi =ΠΦX (p). ThenDΦ( vpi ;∇Φ∗(pˆ− qˆ)) ≥ DΦ( vpi ;Φ∗(pˆi− qˆ)) ∀v ∈ X , qˆ ∈ Rd .Proof.DΦ( vpi ;∇Φ∗(pˆ− qˆ)) = DΦ( vpi ; p)+ 〈 qˆ, v−pi 〉 (by Proposition B.1.4)≥ DΦ( vpi ;pi)+ 〈 qˆ, v−pi 〉 (by Proposition B.1.8)= DΦ( vpi ;∇Φ∗(pˆi− qˆ)) (by Proposition B.1.4).B.2 Proofs for Section 3.3.1Proof of Proposition 3.3.1Proof. First we apply Proposition B.1.6 with a = x, b = x′ and d = ∇Φ∗(xˆ− qˆ) toobtainDΦ( xx′ ;w) = 〈 xˆ− dˆ, x− x′ 〉−DΦ(x′,x)= 〈 qˆ, x− x′ 〉−DΦ(x′,x)(i)≤ ‖qˆ‖∗‖x− x′‖− ρ2 ‖x− x′‖2(ii)≤ ‖qˆ‖2∗/2ρ,where (i) is from the definition of dual norm and Proposition B.1.3, (ii) is byFact B.1.1.B.3 Proofs for Section 3.3.2Proposition 3.3.3. Let a,b ∈ X and c ∈ D. Then DΦ(ab ;c)≤ Λ(a,c).119Proof of Proposition 3.3.3Proof. Since a,b ∈ X we have ‖a‖1 = ‖b‖1 = 1. ThenDΦ(ab ;c) = DKL(a,c)−DKL(b,c)=(DKL(a,c)+1−‖c‖1+ ln‖c‖1)− (DKL(b,c)+1−‖c‖1+ ln‖c‖1)= Λ(a,c)−Λ(b,c) (by definition of Λ)≤ Λ(a,c) (by Proposition 3.3.2).Proposition 3.3.4. Let a ∈ X , qˆ ∈ [0,1]d and η > 0. Then Λ(a,∇Φ∗(aˆ−η qˆ)) ≤η2〈a, qˆ〉/2.Proof. Let b = ∇Φ∗(aˆ−η qˆ). By (3.24), bi = ai exp(−η qˆi). ThenΛ(a,∇Φ∗(aˆ−η qˆ)) =d∑i=1ai ln(ai/bi)+ ln‖b‖1=d∑i=1ηaiqˆi+ ln( d∑i=1ai exp(−η qˆi))≤d∑i=1ηaiqˆi+d∑i=1ai exp(−η qˆi)−1 (by Fact B.1.4)≤d∑i=1ηaiqˆi+d∑i=1ai(1−η qˆi+ η2qˆ2i2)−1 (byFact B.1.2)≤ η2d∑i=1aiqˆi/2,using ∑di=1 ai = 1 (since a ∈ X ) and qˆ2i ≤ qˆi (since qˆ ∈ [0,1]d).B.4 Proofs for Section 3.4At many points throughout this section we will need to talk about optimality condi-tion for problems where we minimize a convex function over a convex set. Suchconditions depend on the normal cone of the set on which the optimization is takingplace.120Definition B.4.1. The normal cone to C ⊆ Rd at x ∈ Rd is the set NC(x) := {s ∈Rd | 〈s,y− x〉 ≤ 0 ∀y ∈C}.Lemma B.4.1 (Rockafellar, 1970, Theorem 27.4). Let h : C →R be a closed convexfunction such that (riC)∩ (riX ) 6= /0. Then, x ∈ argminz∈X h(z) if and only if thereis g ∈ ∂h(x) such that −g ∈ NX (x).Using the above result allows us to derive a useful characterization of pointsthat realize the Bregman projections.Lemma B.4.2. Let y ∈D and x ∈ D¯. Then x =ΠΦX (y) if and only if x ∈D∩X and∇Φ(y)−∇Φ(x) ∈ NX (x).Proof. Suppose x ∈D∩X and ∇Φ(y)−∇Φ(x) ∈ NX (x). Since ∇Φ(y)−∇Φ(x) =−∇(DΦ(·,y))(x), by Lemma B.4.1 we conclude that x ∈ argminz∈X D(z,y). Nowsuppose x = ΠΦX (y). By Lemma B.4.1 together with the definition of Bregmandivergence, this is the case if and only if there is −g ∈ ∂Φ(x) such that −(g−∇Φ(y)) ∈ NX (x). Since Φ is of Legendre type we have ∂Φ(z) = /0 for any z 6∈ D(see Rockafellar, 1970, Theorem 26.1). Thus, x ∈ D and g = ∇Φ(x) since Φ isdifferentiable. Finally, x ∈ X by the definition of Bregman projection.Before proceeding to the proof of the results from Section 3.4, we need to stateon last result about the relation of subgradients and conjugate functions.Lemma B.4.3 (Rockafellar, 1970, Theorem 23.5). Let f : X → R, let x ∈ X andlet yˆ ∈ Rd . Then yˆ ∈ ∂ f (x) if and only if x attains supx∈Rd (〈 yˆ, x〉− f (x)) = f ∗(yˆ).Proof of Proposition 3.4.1Proof. Let t ≥ 1 and let Ft : D→ R be the function being minimized on the right-hand side of (3.31). By definition we have x(t+1)=ΠΦX (y(t+1)). Using the optimalityconditions of the Bregman projection, we havex(t+1) =ΠΦX (y(t+1)) ⇐⇒ yˆ(t+1)− xˆ(t+1) ∈ NX (x(t+1)), (by Lemma B.4.2)121By further using the definitions from Algorithm 3 we getyˆ(t+1)− xˆ(t+1) = γt(xˆ(t)−ηtgt)+(1− γt)xˆ(1)− xˆ(t+1)= γt(xˆ(t)− xˆ(t+1)−ηtgt)+(1− γt)(xˆ(1)− xˆ(t+1))=−γt(∇(DΦ(·,x(t)))(x(t+1))+ηtgt)− (1− γt)∇(DΦ(·,x(1)))(x(t+1))=−∇Ft(x(t+1))Thus, we have −∇Ft(x(t+1)) ∈ NX (x(t+1)). By the optimality conditions fromLemma B.4.1 we conclude that x(t+1) ∈ argminx∈X Ft(x), as desired.Proof of Theorem 3.4.1Theorem 3.4.1 is an easy consequence of the following proposition.Proposition B.4.1. Let { ft}t≥1 with ft : X → R be a sequence of convex functionsand let η : N→ R>0 be non-increasing. Let {x(t)}t≥1 and {gt}t≥1 be as in Algo-rithm 3. Define γ [i,t] := ∏tj=i γ j for every i, t ∈ N. Then, there are {pt}t≥1 withpt ∈ NX(x(t)) for each t ≥ 1 such that ∀t ≥ 0,{x(t+1)}= argminx∈X( t∑i=1γ [i,t]〈ηigi+ pi, x〉−(γ [1,t]+t∑i=1γ [i+1,t](1−γi))〈 xˆ(1), x〉+Φ(x)).(B.1)Proof. First of all, in order to prove (B.1) we claim it suffices to prove that thereare {pt}t≥1 with pt ∈ NX(x(t)) for each t ≥ 1 such thatyˆ(t+1) =−t∑i=1γ [i,t](ηigi+ pi)+(γ [1,t]+t∑i=1γ [i+1,t](1− γi))xˆ(1), ∀t ≥ 0. (B.2)122To see the sufficiency of this claim, note thatx(t+1) =ΠΦX (y(t+1)))⇐⇒ yˆ(t+1)− xˆ(t+1) ∈ NX (x(t+1)) (Lemma B.4.2)⇐⇒ yˆ(t+1) ∈ ∂ (Φ+δ (· |X ))(x(t+1)) (∂ (δ (· |X ))(x) = NX (x))⇐⇒ x(t+1) ∈ argmaxx∈Rd(〈 yˆ(t+1), x〉−Φ(x)−δ (x |X )) (Lemma B.4.3)⇐⇒ x(t+1) ∈ argminx∈X(−〈 yˆ(t+1), x〉+Φ(x)).The above together with eq. (B.2) yields eq. (B.1). Let us now prove eq. (B.2) byinduction on t ≥ 0.For t = 0, eq. (B.2) holds trivially. Let t > 0. By definition, we have yˆ(t+1) =(1− γt)(xˆ(t)−ηtgt) + γt xˆ(1). At this point, to use the induction hypothesis, weneed to write xˆ(t) in function of yˆ(t). From the definition of Algorithm 3, we havex(t) =ΠΦX (yt). By Lemma B.4.2, the latter holds if and only if yˆ(t)− xˆ(t) ∈ NX (x(t)).That is, there is pt ∈NX (x(t)) such that xˆ(t) = yˆ(t)− pt . Plugging these facts togetherand using our induction hypothesis we haveyˆ(t+1)=γt(xˆ(t)−ηtgt)+(1− γt)xˆ(1) = γt(yˆ(t)−ηtgt − pt)+(1− γt)xˆ(1)I.H.= γt(−t−1∑i=1γ [i,t−1](ηigi+ pi)−ηtgt − pt +(γ [1,t−1]+t−1∑i=1γ [i+1,t−1](1− γi))xˆ(1))+(1− γt)xˆ(1)=−t∑i=1γ [i,t](ηigi+ pi)+(γ [1,t]+t∑i=1γ [i+1,t](1− γi))xˆ(1),and this finishes the proof of eq. (B.2).Proof (of Theorem 3.4.1). Define γ [i,t] for every i, t ∈ N as in Proposition B.4.1.If γt = 1 for all t ≥ 1, then γ [i,t] = 1 for any t, i ≥ 1. Moreover, if γt = ηt+1ηt forevery t ≥ 1, then for every t, i ∈ N with t ≥ i we have γ [i,t] = ηt+1ηi , which yields123γ [i,t](ηigi+ pi) = ηtgi+ 1ηi pi andγ [1,t]+t∑i=1γ [i+1,t](1− γi) = ηt+1η1 +t∑i=1ηt+1ηi+1(1− ηi+1ηi)=ηt+1η1+ηt+1t∑i=1( 1ηi+1− 1ηi)= 1.124Appendix CAppendix for Chapter 4C.1 Proofs for Section 4.3Proof of Proposition 4.3.2Proof. Given x1,x2 ∈ Rd , ∀i ∈ [n],‖∂ fi(x2)−∂ fi(x1)‖≤ ‖`′(hi(x2))∂hi(x2)− `′(hi(x1))∂hi(x1)‖≤ ‖`′(hi(x2))∂hi(x2)− `′(hi(x2))∂hi(x1)+ `′(hi(x2))∂hi(x1)− `′(hi(x1))∂hi(x1)‖≤ ‖`′(hi(x2))∂hi(x2)− `′(hi(x2))∂hi(x1)‖+‖`′(hi(x2))∂hi(x1)− `′(hi(x1))∂hi(x1)‖≤ ‖`′(hi(x2))(∂hi(x2)−∂hi(x1))‖+‖(`′(hi(x2))− `′(hi(x1)))∂hi(x1)‖≤ |`′(hi(x2))|‖∂hi(x2)−∂hi(x1))‖+ |`′(hi(x2))− `′(hi(x1))|‖∂hi(x1)‖≤ 2L|`′(hi(x2))|+ |hi(x2)−hi(x1)|×L(i)≤ 2L√2 fi(x2)+L2‖x2− x1‖, (C.1)where (i) follows the same argument as the proof of Proposition 4.3.1:|`′(hi(x2))|= |`′(hi(x2))− `′(0)| ≤√2(`(hi(x2))− `(0)) =√2 fi(x2).125Exchange x1 and x2, we can also get‖∂ fi(x2)−∂ fi(x1)‖ ≤ 2L√2 fi(x1)+L2‖x2− x1‖. (C.2)Combining eq. (C.2) and eq. (C.1), we finished the proof for eq.(4.5).Given x1,x2 ∈Rd , note that fi(x) is almost every differentiable by Rademacher’sTheorem, thusfi(x2) = fi(x1)+∫ 10〈∂ fi(x1+ τ(x2− x1)),x2− x1〉dτ= fi(x1)+ 〈∂ fi(x1),x2− x1〉+∫ 10〈∂ fi(x1+ τ(x2− x1))−∂ fi(x1),x2− x1〉dτ(C.3)Note that ∫ 10〈∂ fi(x1+ τ(x2− x1))−∂ fi(x1),x2− x1〉dτ≤∫ 10‖∂ fi(x1+ τ(x2− x1))−∂ fi(x1)‖‖x2− x1‖dτ≤∫ 10(2L√2 fi(x1)+L2‖τ(x2− x1)‖)‖x2− x1‖dτ≤ 2L‖x2− x1‖√2 fi(x1)+L22‖x2− x1‖2. (C.4)Plug eq. (C.4) into eq. (C.3), we finished the proof for eq. (4.6).126
Thesis/Dissertation
eng
Computer Science
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Friedlander, Michael P.
First-order methods for structured optimization
Text