Science, Faculty of
Computer Science, Department of
DSpace
UBCV
Shea, Betty
Master of Science - MSc
University of British Columbia
Many machine learning (ML) problems are solved by iteratively optimizing an objective along a single descent direction. By restricting the search for the next iterate to one direction, per iteration cost is kept low. In many important problem classes, however, it is possible to perform subspace optimization (SO) and search along additional directions nearly "for free". This can lead to finding a better next iterate and result in fewer iterations to reach a solution. In theory and in practice, however, SO is rarely used and existing theory usually assumes an exact SO. In this thesis, we show experimentally that SO could lead to more accurate solutions in less time. We also explore the effect of different SO parameter settings. The experiments are performed through a software package we developed called minFuncLinear which is available for download. On the theory side, we propose a more abstract problem class, multilinear map problems, where SO could be performed efficiently. This encompasses linear composition models, known to allow efficient SO. We give further examples of problems that fall under this framework, such as log-determinant problems. We also provide preliminary convergence analysis of gradient descent (GD) combined with SO, discuss multi-dimensional Wolfe conditions and extend initialization and interpolation methods from linesearches to subspace searches.
Subspace Optimization for Machine LearningbyBetty SheaB. Sc., Carnegie Mellon University, 1999M. Sc., Carnegie Mellon University, 2007A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFMaster of ScienceinTHE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES(Computer Science)The University of British Columbia(Vancouver)November 2021© Betty Shea, 2021The following individuals certify that they have read, and recommend to the Faculty of Graduateand Postdoctoral Studies for acceptance, the thesis entitled:Subspace Optimization for Machine Learningsubmitted by Betty Shea in partial fulfillment of the requirements for the degree of Master ofScience in Computer Science.Examining Committee:Mark Schmidt, Associate Professor, Computer Science, UBCSupervisorChen Greif, Professor, Computer Science, UBCSupervisory Committee MemberiiAbstractMany machine learning (ML) problems are solved by iteratively optimizing an objective along asingle descent direction. By restricting the search for the next iterate to one direction, per iterationcost is kept low. In many important problem classes, however, it is possible to perform subspaceoptimization (SO) and search along additional directions nearly “for free”. This can lead to findinga better next iterate and result in fewer iterations to reach a solution. In theory and in practice,however, SO is rarely used and existing theory usually assumes an exact SO. In this thesis, we showexperimentally that SO could lead to more accurate solutions in less time. We also explore the ef-fect of different SO parameter settings. The experiments are performed through a software packagewe developed called minFuncLinear which is available for download. On the theory side, wepropose a more abstract problem class, multilinear map problems, where SO could be performedefficiently. This encompasses linear composition models, known to allow efficient SO. We givefurther examples of problems that fall under this framework, such as log-determinant problems.We also provide preliminary convergence analysis of gradient descent (GD) combined with SO, dis-cuss multi-dimensional Wolfe conditions and extend initialization and interpolation methods fromlinesearches to subspace searches.iiiLay SummaryMany research problems, in computer science, economics, etc., could be distilled into finding so-lutions to mathematical formulas that represent our objectives. Usually, this is solved by an “algo-rithm” that looks for a solution. Most algorithms look for the solution in only one direction at atime. In some cases, however, searching in many directions at the same time may lead to findingbetter solutions faster. This work shows that you can search in many directions for some commonmachine learning problems. It also explains how to pick the directions and perform the search.ivPrefaceThis thesis is original, unpublished, independent work by the author, Betty Shea. The basis forthis thesis is unpublished work conceived with Mark Schmidt. The code for the experiments,minFuncLinear , is a stripped down, Julia version of Mark’s minFunc Matlab code modi-fied to exploit problem structure and to perform subspace optimizations.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiNotation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixAcronyms and Abbreviations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Thesis organization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Problem setting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.3 Related work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41.3.1 Quadratic objectives and SO . . . . . . . . . . . . . . . . . . . . . . . . . 51.3.2 Non-quadratic objectives and SO . . . . . . . . . . . . . . . . . . . . . . . 71.3.3 Efficient SO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 Problem Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112.1 Multilinear maps . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122.2 Linear maps . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132.2.1 Trace operator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132.2.2 Linear composition problems . . . . . . . . . . . . . . . . . . . . . . . . 13vi2.3 Log determinant problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172.3.1 Gaussian negative log-likelihood (NLL)s . . . . . . . . . . . . . . . . . . . 172.4 Matrix factorization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182.4.1 Matrix completion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 183 Descent directions and stepsize methods . . . . . . . . . . . . . . . . . . . . . . . . . 203.1 Descent directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 203.1.1 Primary directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 203.1.2 Secondary directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 223.2 Stepsize methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 233.2.1 Quasilinearization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 233.2.2 Solvers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 253.2.3 d-dimensional Wolfe . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 253.3 Initialization and interpolation methods . . . . . . . . . . . . . . . . . . . . . . . 253.3.1 Initialization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 263.3.2 Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 283.4 GD and exact SO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 283.4.1 Quadratic upper bound with Euclidean norm . . . . . . . . . . . . . . . . 283.4.2 Quadratic upper bound with general quadratic norm . . . . . . . . . . . . 293.4.3 Higher-order upper bound . . . . . . . . . . . . . . . . . . . . . . . . . . 314 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 324.1 minFuncLinear software . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 334.1.1 Features, usage and example . . . . . . . . . . . . . . . . . . . . . . . . . 334.1.2 Performance gain from optimizing to problem structure . . . . . . . . . . . 334.2 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 344.2.1 From linesearch to SO . . . . . . . . . . . . . . . . . . . . . . . . . . . . 354.2.2 Choice of secondary directions . . . . . . . . . . . . . . . . . . . . . . . . 374.2.3 Choice of SO descent method . . . . . . . . . . . . . . . . . . . . . . . . 384.2.4 Initialization and stopping criteria . . . . . . . . . . . . . . . . . . . . . . 384.2.5 Increasing the number of secondary directions . . . . . . . . . . . . . . . 394.2.6 Quasi-Newton methods not converging . . . . . . . . . . . . . . . . . . . 395 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45A Background and supporting material . . . . . . . . . . . . . . . . . . . . . . . . . . 49A.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49viiA.1.1 Zoutendijk’s theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49A.1.2 Linear conjugate gradient (CG) . . . . . . . . . . . . . . . . . . . . . . . . 50A.1.3 Polyak’s heavy ball (PHB) with exact SO and CG . . . . . . . . . . . . . . 52A.2 Problem structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53A.2.1 Derivatives of affine function . . . . . . . . . . . . . . . . . . . . . . . . . 53A.2.2 Gaussian models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54A.3 Descent directions and stepsize methods . . . . . . . . . . . . . . . . . . . . . . . 54A.3.1 Direction of steepest descent . . . . . . . . . . . . . . . . . . . . . . . . . 54A.3.2 Newton’s direction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54A.3.3 Quasi-Newton directions . . . . . . . . . . . . . . . . . . . . . . . . . . . 55A.3.4 Broyden Fletcher Goldfarb Shanno (BFGS) convergence . . . . . . . . . . 56A.3.5 Limited memory Broyden Fletcher Goldfarb Shanno (L-BFGS) . . . . . . . 56A.3.6 Memoryless L-BFGS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56A.3.7 Non-linear CG . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57A.3.8 Truncated Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58A.3.9 Hermite polynomial interpolation . . . . . . . . . . . . . . . . . . . . . . 58A.4 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59A.4.1 Barzilai Borwein . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59B Figures and tables from experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . 61viiiList of TablesTable 1.1 Comparison of convergence rates and whether knowledge of λ1, L or µ is re-quired. Table entries are first-order methods on a quadratic objective f (x) =12 xᵀAx+ bᵀx where A is positive definite with eigenvalues 0 < λ1 ≤ ·· · ≤ λn.k is the iteration count and εk is a sequence that goes to zero asymptotically.Assumes that arithmetic, linesearches and SOs are exact. . . . . . . . . . . . . 10Table 2.1 Number of matrix-vector multiplications per iteration of GD for objectives ofthe form f (x) = F(Ax) and second direction p2 in the case of a plane search.There is an additional matrix-vector multiplication at the first iteration for Ax0.Subsequent calculations of Axk do not require matrix-vector multiplications. . . 13Table 4.1 Binary classification of the spam dataset where * denotes the optimized version.Comparison of minimum value found, total time taken, number of iterations andstepsize search iterations and number of matrix-vector multiplications. . . . . . 34Table 4.2 Binary classification of the sido dataset with different initial solution x0. L-BFGS-0 and L-BFGS-1 correspond to x0 is a vector of zeros and ones respectively. Thelatter diverges with a linesearch but remains convergent with a SO. . . . . . . . 40Table B.1 Binary classification of the spam dataset. Comparing a linesearch to a SO. . . . 62Table B.2 Binary classification of the sido dataset. Comparing a linesearh to a SO. . . . . 62Table B.3 Binary classification of the spam dataset. Comparing momentum and truncatedNewton (TN) as a second direction. . . . . . . . . . . . . . . . . . . . . . . . . 63Table B.4 Binary classification of the sido dataset. Comparing momentum and TN as asecond direction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63Table B.5 Binary classification of the spam dataset. Comparing different descent methodsused in the SO. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64Table B.6 Binary classification of the sido dataset. Comparing different descent methodsused in the SO. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64ixTable B.7 Binary classification of the spam dataset. Comparing different initializationmethods and exit conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . 65Table B.8 Binary classification of the sido dataset. Comparing different initialization meth-ods and exit conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65Table B.9 Binary classification of the sido dataset. The effect of going from one to fourmomentum terms. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66Table B.10 Binary classification of the sido dataset. The effect of going from one to fourmomentum terms. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66xList of FiguresFigure 4.1 Binary classification of the spam dataset using various descent methods. . . . 36Figure 4.2 Binary classification of the sido dataset using various descent methods. . . . . 37Figure B.1 Binary classification of the spam dataset where * denotes the optimized ver-sion. Convergence is plotted in terms of both number of iterations and time (inseconds). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61Figure B.2 Binary classification of the spam dataset. Comparing momentum and TN as asecond direction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67Figure B.3 Binary classification of the sido dataset. Comparing momentum and TN as asecond direction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68Figure B.4 Binary classification of the spam dataset. Comparing SO methods. . . . . . . . 69Figure B.5 Binary classification of the sido dataset. Comparing SO methods. . . . . . . . 70Figure B.6 Binary classification of the spam dataset. Linesearch initialization and relativeexit conditions for the SO. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71Figure B.7 Binary classification of the sido dataset. Linesearch initialization and relativeexit conditions for the SO. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72Figure B.8 Binary classification of the spam dataset. Varying the number of momentumdirections in the SO. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73Figure B.9 Binary classification of the sido dataset. Varying the number of momentumdirections in the SO. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74xiNotationA m×n features matrixd number of descent directions to use in searchf objective functionk outer iteration countL Lipschitz constantm number of training examplesn dimension of optimization problempk single descent direction on the kth iterationPk a n×d matrix with d search directions as columnsxk kth iteratex∗ optimal solutionαk kth stepsize(s)µ strong convexity constantxiiAcronyms and AbbreviationsBB Barzilai BorweinBCD Block coordinate descentBFGS Broyden Fletcher Goldfarb ShannoCG conjugate gradientDFP Davidson Fletcher PowellFR Fletcher-ReevesGD gradient descentGED geometric descentGLM generalized linear modelHPI Hermite polynomial interpolationHS Hestenes-StiefelISM iterated-subspace minimizationL-BFGS Limited memory Broyden Fletcher Goldfarb ShannoML machine learningMLE maximum likelihood estimatorNAG Nesterov’s accelerated gradientNLL negative log-likelihoodPCA principal component analysisPHB Polyak’s heavy ballPRP Polak-Ribie`re-PolyakSESOP sequential subspace optimizationSR1 symmetric rank oneSO subspace optimizationxiiiSVM support vector machineTN truncated NewtonUGM undirected graphical modelxivAcknowledgmentsI am grateful to my supervisors Mark Schmidt and Maryam Kamgarpour for their guidance andmentorship, and to examining committee member Chen Greif for his time and suggestions. I amparticularly indebted to Mark for his many research ideas and for giving a ‘somewhat-non-standard’student an opportunity. To the many students that I am so fortunate to meet (and especially Aaron,Cathy, Fred and Wilder): Thank you for sharing your work with me and for making these two yearsmemorable. My deepest gratitude goes to Mom, Dad, Wendy and Freddy: Thank you for yoursupport and for your unwavering bullishness.xvChapter 1IntroductionMany research problems can be formulated as an optimization problem where an objective functionis minimized by an iterative method. For example, training a robot to navigate a course could bethought of as minimizing the instances it hits obstacles and the time it takes to complete the course.Learning from previous mistakes and successes, the robot hones its navigation skills. In finance,stock picking could be simplified to minimizing negative returns on a portfolio. One could startwith random allocations and repeatedly adjust the weights to find increasingly profitable portfolios.Broadly speaking, we want to progressively find better solutions where “better” is measuredby a function that captures how far we are from the goal. This abstraction hides many practicaldetails. How do we ensure that we move in the direction of our goal? How do we do so as quicklyas possible? The topic of this thesis, subspace optimization (SO), is about some aspects of thesepractical considerations.1.1 Thesis organizationSO is defined by its set of search directions and the method to determine the stepsizes for eachdirection. In the remainder of this chapter, we introduce notation and the problem setting. Wealso discuss SO in the broader context of convex minimization and list some examples of relatedwork. Chapter 2 discusses problem structures where SO can be performed efficiently. The generalclass of problems where SO can be performed efficiently are a subset of multilinear map problemsand we focus particularly on the problem class of linear map models. Assuming a large scalemachine learning (ML) setting where the main computational bottleneck are expensive matrix-vectormultiplications, we demonstrate how SO can be added “for free” to existing optimization methods.We also describe other common ML problems that fall under this broad category of multilinearmaps.Chapter 3 examines practical design decisions of performing SOs. We discuss the choice ofsearch directions, stepsize methods, and practical considerations for initializing and exiting a SO.1We conceptually divide search directions into a “primary” descent direction (for example, the gradi-ent or Newton direction) and a set of secondary directions. We give three different stepsize methodsthat include satisfying d-dimension Wolfe conditions. We also provide preliminary convergenceanalysis of adding SO to gradient descent (GD).Chapter 4 introduces software (minFuncLinear ) developed for our experiments. Our ex-periments are based on the generalized linear model discussed in chapter 2. The experiments showthat, in some cases, the common approach of using a linesearch should be replaced by using a SO.We back up our claim SOs could be performed “for free” by showing faster convergence, both innumber of iterations and in time, to more accurate solutions. The experiments also highlight prac-tical differences between different SO settings. We conclude with a discussion of future work inChapter 5.1.2 Problem settingSuppose our goal is to minimize an unconstrained function that is smooth and convex. In otherwords, the problem isminx∈Rnf (x) (1.1)where f : Rn → R. By smoothness, we mean that for scalar L > 0, Euclidean norm ‖ · ‖ and anyx,y ∈ Rn, f is differentiable and changes to the gradient ∇ f are bounded. We write‖∇ f (x)−∇ f (y)‖ ≤ L‖x− y‖, for all x and y. (1.2)(When the subscript to ‖ ·‖ is omitted, the Euclidean norm is assumed.) Sometimes we also assumethat f is µ-strongly convex or, for scalar µ > 0, we havef (y)≥ f (x)+∇ f (x)ᵀ(y− x)+ µ2‖y− x‖2, for all x and y. (1.3)A twice-differentiable f that satisfies (1.2) and (1.3) has a Hessian with eigenvalues bounded by Land µ , LI  ∇2 f (x) µI. An example of a smooth and strongly convex function is the regularizedleast squares objective f (x) = ‖Ax− y‖2+‖x‖2 where A ∈ Rm×n.The solution to (1.1) is typically through iterative methods. Starting with an initial guess x0of the optimal solution x∗, we repeat a procedure k = 1,2, . . . times where each successive stepgives us a new iterate xk+1 that is at least as accurate as our current iterate xk. In other words,f (xk+1) ≤ f (xk). Most methods fall into this framework and decrease f or an upper bound of f .The algorithm then repeats until some stopping criteria is met.The main difference between iterative methods is the way they arrive at xk+1 from xk. Onepopular strategy is a line search. We choose a “descent direction” pk and a step size αk which is theamount we move in the direction of pk. For example, GD with a fixed stepsize of 1/L is an iterative2method where the choices are pk =−∇ f (xk) and αk = 1/L. This choice guarantees convergence toa global optimal for convex and L-smooth functions (Nesterov, 2018).Throughout this thesis, we often assume that any search direction pk is a descent direction, or itsdirectional derivative is negative. The directional derivative is the derivative of the objective f alongthe search direction pk and is defined as limε→0f (xk+ε pk)ε . For f that is continuously differentiablearound xk, taking the derivative and the chain rule means that the directional derivative is∇ f (xk)ᵀpk.Cosine similarity is one way to interpret a negative directional derivative. Let θ be the angle formedby vectors ∇ f (xk) and pk. The definition of the inner product gives their cosine similaritycosθ =∇ f (xk)ᵀpk‖∇ f (xk)‖‖pk‖ . (1.4)Cosine similarity, with values between [−1,1], measures the similarity of two vectors. A value of0 indicates orthogonality, positive values imply similarity and negative values imply dissimilarity.Because the denominator is always positive, a negative directional derivative, ∇ f (xk)ᵀpk < 0, im-plies that pk is dissimilar to ∇ f (xk). If f is convex, ∇ f (x) is an ascent direction and a dissimilar pkis a descent direction. For example, if pk = −∇ f (xk) as used in GD, then cosine similarity is −1and thus pk a descent direction.1After selecting αk and pk, the next iterate is given byxk+1 = xk +αk pk (1.5)We will discuss different choices of pk in chapter 3. In terms of choosing a stepsize, we wouldideally look for the perfect step size for our chosen direction that satisfiesαk ∈ argminα>0f (xk +α pk) (1.6)In practice, (1.6) often can not be solved exactly. Instead, we choose an αk that satisfies somedesirable conditions. Popular inexact stepsize conditions include the Wolfe conditions. Given con-stants c1 ∈ (0, 12 ] and c2 ∈ (c1,1), the Wolfe conditions aref (xk +αk pk)≤ f (xk)+ c1αk∇ f (xk)ᵀpk (1.7)and∇ f (xk +αk pk)ᵀpk ≥ c2∇ f (xk)ᵀpk (1.8)Equation (1.7) is sometimes referred to as the Armijo condition or the “sufficient decrease” con-dition. It finds a stepsize that decreases the objective proportional to the directional derivative.1Zoutendijk’s theorem is a more formal way of looking at descent directions and relates it to the Wolfe conditions(1.7) and (1.8). For more information, see Appendix A.1.1.3Equation (1.8) is sometimes referred to as the curvature condition and prevents taking very smallstep sizes. The curvature condition is sometimes strengthened to|∇ f (xk +αk pk)ᵀpk| ≤ c2|∇ f (xk)ᵀpk|. (1.9)Combined with (1.7), this gives the strong Wolfe conditions. One way to view (1.9) is that themagnitude of the directional derivative at the new iterate cannot be too large. This means that thenew iterate has a slope that cannot be too large and thus is an approximate stationary point alongthe direction.Another popular line search condition is the backtracking line search with contraction factorρ ∈ (0,1) and initial positive stepsize α0. Backtracking multiplies α0 by ρ repeatedly until it findsa stepsize that satisfies the Armijo condition (1.7). For more details on the topics discussed in thissubsection, chapter 3 of Nocedal and Wright (2006) and chapter 8 of Luenberger and Ye (2008) givean overview. They also discuss alternatives to line search methods such as trust region methods.Although far less common, a linesearch can be generalized to a SO which is a search within aspace spanned by a set of directions. Analogously, we determine a set of d promising directions thatform the columns of a n×d matrix Pk. As in (1.6), ideally we would like to find stepsizes αk ∈ Rdthat minimize f along our chosen directions. In other words,αk ∈ argminαf (xk +Pkα) (1.10)Similarly, (1.10) cannot be solved exactly in practice. Thus, the design decisions of a SO are:1. which directions Pk to include in the subspace,2. the stepsize conditions to use to approximate (1.10), and3. the search method to use to find αk.To distinguish the SO subproblem (1.10) from the original optimization problem (1.1), we referto iterating through the main descent method as the outer loop of the optimization problem and tothe linesearch or SO subproblem as the inner loop of the optimization problem.1.3 Related workIn this section, we discuss how SO fits into the broader context of convex minimization. We startwith the well-known method of conjugate gradient originally introduced to solve positive definitelinear systems. Minimizing a convex quadratic objective is equivalent to solving a linear system.We use this problem class to highlight known relationships between conjugate gradient (CG), accel-eration in the sense of Nesterov and SO.41.3.1 Quadratic objectives and SOConvex quadratic objective functionsA convex quadratic objective is writtenf (x) =12xᵀAx+bᵀx+ c (1.11)where A∈Rn×n is positive semi-definite, x,b∈Rn and c∈R. By equating∇ f (x)with 0, minimizing(1.11) means setting Ax+b = 0 and thus equivalent to solving a linear system Ax =−b.Linear CG and SOThe method of CG was originally developed to solve linear systems where A is positive definite(Hestenes and Stiefel, 1952). The connections between CG and SO are many. For example, CGperforms a SO on a set of conjugate directions. This is sometimes referred to as expanding subspaceminimization. CG starts with a search direction equal to the negative gradient at the initial solutionx0, or p0 = −∇ f (x0). At every iteration k, CG chooses a new pk that is conjugate to previousdirections p0, p1, . . . , pk−1. In other words,pᵀk Ap j = 0 for all j = 0,1, . . . ,k−1 and positive definite AThe CG update step then becomesxk+1 = xk +αk pk where αk =−rᵀk pkpᵀk Apkand rk , Axk +b. (1.12)Assuming that all the eigenvalues of ∇2 f (x), 0 < λ1 ≤ ·· · ≤ λn, are known, the CG iterates satisfy2‖xk+1− x∗‖A ≤(λn−k−λ1λn−k +λ1)‖x0− x∗‖A (1.13)where the A-norm is defined as ‖x‖A , xᵀAx. This norm of the error is minimized through theiteration.Equation (1.13) is derived by showing that the eigenvalues of A produce small values of theoptimal kth degree polynomial generated at every iteration of CG. This CG convergence rate isoptimal for convex quadratics, for algorithms based on matrix-vector multiplication. For details onthe proof of (1.13), see chapters 6 and 7 of Saad (2003) and chapter 9 of Luenberger and Ye (2008).If only the largest and smallest eigenvalues, λmax and λmin, are known, a looser convergence could2The dependence on the n−kth eigenvalue reflects the way CG removes the effect of the largest remaining eigenvalueat each iteration.5be stated as‖xk− x∗‖A ≤ 2(√κ(A)−1√κ(A)+1)k‖x0− x∗‖A (1.14)where κ(A), the Euclidean condition number of A, isκ(A) = ‖A‖2‖A−1‖2 = λmax/λmin. (1.15)Equation (1.14) is derived through bounding the error of Chebyshev polynomials. For details, seeGreenbaum (1997). Note that (1.14) is faster than GD with exact linesearch, which uses κ(A) insteadof√κ(A). Thus, the SO method is much faster in theory and in practice.There are multiple ways to pick conjugate directions. In theory, one could take the eigenvectorsof A but this would be too expensive. A better choice is via a modified Gram-Schmidt algorithm thatreturns conjugate directions instead of orthogonal directions. The method of CG is a specific form ofa conjugate directions method where, in practice, one only needs the current gradient, the previousgradient and the previous direction (∇ f (xk), ∇ f (xk−1) and pk−1) to generate the next direction. Thisis a result of exactly optimizing the stepsize at each iteration and the growing subspace optimizationproperty, which makes most previous inner products equal to zero. CG also has a recursive formthat can be written as a two-dimensional subspace minimization at each step. For more details, seechapter 1 of Nesterov (2018).A key result of using conjugate directions is that, in the absence of roundoff errors, CG ter-minates in at most n steps where n is the dimension of the problem. If there are only k distincteigenvalues, then CG converges within k iterations. In modern days, the latter observation is moresignificant. If k n, then CG arrives at a solution with a small residual after relatively few itera-tions. For more details on CG as a SO, properties of CG and convergence of CG, see Appendix A.1.2and the references given there.Momentum with SOIn Miele and Cantrell (1969a); Cantrell (1969), the authors introduce “memory gradient descent”.Essentially gradient descent with memory of its previous direction, the update on the kth iteration isxk+1 = xk−αk∇ f (xk)+βk(xk− xk−1). (1.16)The term xk− xk−1 is sometimes referred to as momentum. The authors show that, for quadraticobjectives, using an exact SO to find stepsizes αk and βk in (1.16) gives the same update rule as thenon-linear CG method Fletcher-Reeves (FR). Because FR reduces to linear CG in the quadratics case(Dai and Yuan, 2000), the memory gradient method with exact SO is also equivalent to linear CG.Thus this method achieves the optimal convergence rate for quadratics.6Note that constant stepsizes in (1.16) gives the Polyak’s heavy ball (PHB) method which achievesthe convergence rate‖xk− x∗‖ ≤(√κ(A)−1√κ(A)+1+ εk)k‖x0− x∗‖where εk is a sequence that satisfies‖Ak‖ ≤ (maxi|λi(A)|+ εk)k and limk→∞εk = 0.For a detailed explanation, see Recht (2012). This shows that one can achieve CG rate without a SO.The convergence rate remains linear but is faster than the rate attained by GD. The fixed stepsize,however, requires knowing λ1(A) and λn(A) unlike CG and SO, and also requires the sequence εk.In practice, having to know the extremal eignevalues is a significant disadvantage. For convergenceproperties of PHB for different problem classes, see Ghadimi et al. (2014). For more details, seeBertsekas (1982); Nesterov (2018) and appendix A.1.3.1.3.2 Non-quadratic objectives and SOWhile the quadratic case provide a conceptual base to relate different methods, many problems arenot quadratic in practice. In this section, we discuss the role of SO in the convex objectives case. Fornon-quadratic objectives, it is typically impossible to perform SO exactly. Therefore, there is moreemphasis on search conditions such as the Wolfe conditions (1.7) and (1.9) discussed previously.Nesterov’s accelerated gradient (NAG) achieves the optimal rate for smooth convex minimiza-tion and could be considered a SO-free version of the optimal method in Nemirovsky (1982). InDrori and Taylor (2020), the authors generalize this to other problem classes. They propose a tech-nique that transform a SO method into a fixed-step method with no-worse convergence rates. Theirframework, however, assumes that methods satisfy interpolability constraints.Strongly-convex objectivesRecently, in Karimi and Vavasis (2017), the authors introduce a unifying theory for CG, NAG (Nes-terov, 1983) and geometric descent (GED) (Bubeck et al., 2015) in the context of minimizing asmooth, strongly-convex function. They give an idealized algorithm that, at every iteration, per-forms two exact SOs on a subspace that is two-dimensional on most iterations. They show that thisidealized algorithm improves the solution by a factor of (1−√1/κ) per iteration where κ = L/µin the non-quadratic setting. Their goal is to demonstrate the relationship between SO and variousalgorithms that achieve the optimal convergence rate for this problem class.7Non-linear CG and SONon-linear CG is a group of methods that attempt to generalize linear CG for the non-quadratic case.There are many nonlinear CG methods but most follow this framework. The next iterate is given by(1.5) and the directions are generated according topk+1 =−∇ f (xk+1)+βk pk (1.17)where βk follows a formula that distinguishes between different versions of nonlinear CG and αkis determined by a linesearch. The βk is generally chosen to simulate the conjugacy relationshipbetween successive search directions in linear CG. Popular non-linear CG methods, such as FR, Dai-Yuan and Polak-Ribie`re-Polyak (PRP), reduces to linear CG when the objective is quadratic. In thecase of non-quadratic functions that are thrice differentiable with a bounded Hessian, instead of n-step convergence, Cohen (1972) shows that PRP, Daniel3 and FR give r-step quadratic convergenceafter each restart where r ≥ n. For a survey of the βks corresponding to different nonlinear CGversions, see Hager and Zhang (2006).Equation (1.5) combined with (1.17) is similar in form to that of PHB and thus there are parallelsbetween momentum and nonlinear CG. For example, one could view momentum with stepsizes setvia SO as a generalization of nonlinear CG which has a prescribed stepsize βk and a second stepsizeαk found by a one-dimensional search. In Miele and Cantrell (1969a), the authors empiricallycompare their memory gradient method where stepsizes are set by a SO with FR in the convexnon-quadratic case and find that SO performs better. In Cragg and Levy (1969), the authors extend(1.16) to use the last d momentum directions in an algorithm they called the supermemory gradientmethod.In both the memory and supermemory gradient methods, the SO is performed through quasi-linearization. This technique could be viewed as a generalization of the Newton-Raphson methodwhere the unknowns are functions instead of function values. Given a general non-linear function,quasilinearization linearizes the original function and generates a sequence of functions that con-verge to the solution of the original function. For more details, see Lee (1968). We also discussquasilinearization further in Chapter 3.In Yuan and Stoer (1995), the authors use an inexact two-dimensional SO spanned by the currentgradient and previous direction to find CG-like directions for non-quadratic objectives. Their methodreduces to linear CG when the objective is quadratic and the SO is performed exactly. But theiranalysis does not address more general SO approaches where the search directions do not coincidewith those used by CG methods.Beyond convex objectives, in Conn et al. (1994), the authors introduce a general framework3The Daniel (1967) version of non-linear CG requires calculating∇2 f (x) and is different in form from most non-linearCG methods.8called iterated-subspace minimization (ISM). At the heart of ISM is an inexact SO that uses a multi-dimensional version of the Goldstein conditions to terminate. Global convergence is ensured foran objective function that is L-smooth (1.2), bounded from below and where the inexact SO returnsstepsizes that satisfy the Goldstein conditions. This essentially translates the problem of showingglobal convergence to solving an inexact SO. Because ISM is a framework, no concrete algorithm forthe search itself is given in the paper. The authors, however, provide guidance on what to include inthe subspace. For example, they advocate including at least one gradient-based direction to ensureglobal convergence and one Newton-based direction for faster convergence asymptotically. Werevisit this combination of first- and second-order directions in our experimental section.1.3.3 Efficient SOWhile above works define SO methods, a key factor is whether we can solve SO fast enough tomake it worthwhile compared to traditional linesearches. Some problem structures allow SO to beperformed accurately while keeping the computation cost low. In Narkiss and Zibulevsky (2005a),the authors propose an algorithm, sequential subspace optimization (SESOP), based on a SO at eachstep. The subspace that they consider is spanned by the current gradient ∇ f (xk) and any of thefollowing• the two Nemirovsky directions (Nemirovsky, 1982)• previous momentum directions xk− xk−1, and• previous gradients ∇ f (xk−1),∇ f (xk−2), . . . .Assuming that a SO can be performed exactly on a subspace spanned by the current gradient,Nemirovsky directions and previous steps, and if f is convex and L-smooth (but not necessarilystrongly convex), they show that SESOP iterates xk satisfyf (xk+1)− f (x∗)≤ Lk2 ‖x0− x∗‖Thus, SESOP achieves the optimal rate for L-Lipschitz and convex objectives using first-order meth-ods.4 The authors demonstrate that for some problems, SO can be done cheaply and thus highlyaccurately. Specifically, the authors consider the case of objective functions that are in the formf (Ax+b) where A is an expensive linear mapping and f is a relatively computationally cheap con-vex function. We discuss this important special case later in section 2.2.2. SESOP is further extendedto the case where the truncated Newton (TN) direction is used in place of the gradient as the primarysearch direction in (Zibulevsky, 2008), and is applied to a support vector machine (SVM) in (Narkissand Zibulevsky, 2005b).More recently, there has been renewed interest in problems where SO could be performed ef-ficiently. In Ahookhosh and Neumaler (2018), the authors introduce an algorithm that solves un-4Since strong convexity is not assumed, this is slower than the linear rate of CG, PHB and GD.9Method Rate Needs λ1? Needs µ? Needs L?CG ‖xk+1− x∗‖A ≤(λn−k−λ1λn−k+λ1)‖x0− x∗‖A No No NoPHB ‖xk+1− x∗‖ ≤(√λn−√λ1√λn+√λ1+ εk)‖xk− x∗‖ Yes No NoNon-linear CG+SO equivalent to CG No No NoGD, fixed stepsize ‖xk+1− x∗‖ ≤(λn−λ1λn+λ1)‖xk− x∗‖ No No YesTable 1.1: Comparison of convergence rates and whether knowledge of λ1, L or µ is required.Table entries are first-order methods on a quadratic objective f (x) = 12 xᵀAx+ bᵀx whereA is positive definite with eigenvalues 0 < λ1 ≤ ·· · ≤ λn. k is the iteration count and εkis a sequence that goes to zero asymptotically. Assumes that arithmetic, linesearches andSOs are exact.constrained nonsmooth convex optimization problems where the objective is also the composite ofexpensive linear operators Ai : V →Ui and non-smooth, proper and convex functions fi : Ui → Rthat are relatively cheap to computeminx∈Vf (x) :=p∑i=1fi(Aix) where V,Ui are real finite-dimensional vector spaces.As in SESOP and in this thesis, the authors rely on avoiding expensive matrix operations to performSO to a high degree of accuracy for little cost. They demonstrate empirical improvements usingSVMs. Finding directions is more difficult in non-smooth problems and the Wolfe conditions gen-erally do not work. There is, however, little theoretical analysis given within the paper on the effectof adding exact or inexact SO to the convergence rate of non-smooth methods.Similarly, in (Nesterov et al., 2019) the authors exploit problem structures that allow for compu-tationally efficient solutions to a SO, or a “small-dimensional relaxation oracle”. With this, they pro-pose an accelerated gradient method that achieves the lower bound convergence rate for both convexand non-convex objectives. Unlike existing accelerated methods, they do this without knowledgeof the Lipschitz constant L and with an exact line search. A comparison of convergence rates forquadratic objectives are given in table 1.1.10Chapter 2Problem StructureIn convex optimization, the efficiency of an iterative method is often measured by iteration com-plexity, or the number of iterations it takes to arrive at an acceptable solution. This obscures thecost of subproblems such as SO and is an even larger issue where there are further nested subprob-lems. Instead of iteration complexity, we use the oracle model first proposed in Nemirovsky andYudin (1983). Our information unit is a “computationally expensive” operation where the runtimeis dominated by the number of times this operation is performed. Depending on the optimizationproblem at hand, this could be matrix-vector multiplication, matrix inversion, etc. This allows us todiscuss problem structures where SO could be solved without adding significant cost to the overallmethod. In other words, SO can be done for a much lower cost than the computationally expensiveoperation.In this chapter, we explore models where SO could be solved efficiently. We discuss multilinearmap problems as a category that allows for efficient SO. Multilinear maps include linear compositionmodels such as SESOP, first explored in Narkiss and Zibulevsky (2005a), as well as other problemssuch as matrix factorization. The special case of linear composition models is examined in greaterdetail. An algorithm is efficient if it solves the SO subproblem without increasing the overall numberof computationally expensive operations significantly. We detail how a SO reduces to a minimizationsubproblem solved efficiently through a sequence of linesearches. As a concrete example, we showhow GD can be improved by adding a momentum term that does not require additional oracle calls.Block coordinate descent (BCD) is an iterative method that minimizes a version of the originalobjective that has lower dimensionality. Some BCD problems do allow for SO if we update a subsetof variables. For examples of these cases, see section 2.4 of Nutini et al. (2017). Although BCD isnot the focus of this thesis, this important problem class is a potential future research topic.Recall that, in the previous chapter, we introduced our overall goal of unconstrained minimiza-tion of a real-valued objective (1.1) that is L-smooth (1.2), convex and sometimes µ-strongly convex(1.3). We find a solution through an iterative method that includes a subspace search. The update at11the kth iteration generalizes (1.5) to more than one direction. This has the formxk+1 = xk +Pkαk (2.1)where the columns of Pk ∈ Rn×d contain the d search directions and αk ∈ Rd contain the step sizesfor each direction. A SO is the subproblem of setting αk by solving (1.10). Thus, it is an auxiliaryoptimization problem that increases the computational cost of each iteration of (2.1).2.1 Multilinear mapsA multilinear composition problem is a problem with an optimization objective that is a multilinearmap. This broad category includes many important ML problems. For example, special cases ofmultilinear composition problems include linear composition models (which include generalizedlinear model (GLM)s (Nelder and Wedderburn, 1972)), and matrix problems such as factorization,completion and finding the determinant. We discuss these special cases in this section and focus es-pecially on linear composition problems. For this problem class, if we consider matrix-vector mul-tiplications as an oracle call, we can essentially perform line searches and some types of subspacesearches “for free”. The cost of performing SO for linear composition problems are summarized intable 2.1.As will be discussed in detail later in this chapter, SO can be performed cheaply for linear oper-ators. Multilinear problems decompose into linear subproblems and so, in general, are cases whereSO is efficient. Given vector spaces V1, . . . ,Vk,W over a fixed field F, f : V1×·· ·×Vk→W is a multi-linear map if it is a linear map in each input vector space Vi. In other words, f (v1, . . . ,vi−1, ·,vi+1, . . . ,vk) :Vi →W is linear for all i = 1, . . . ,k. Where there are k variables, f is also referred to as k-linearmap. Special cases are 1-linear maps (or linear maps) and 2-linear maps (bilinear maps).A matrix is an example of a bilinear map. A matrix A over field F is a m×n array of elementsof F. Given standard bases e1,e2 ∈ F2, A = [ai j]m,ni, j=1 ∈ Fm×n is a bilinear mapf : Fm×Fn→ F where ai j = f (e1,e2) ∈ F.Matrix-vector multiplication Ax where vector x ∈Rn is an example of a linear map. Thus, problemsof the form F(Ax,x) are also a multilinear composition problems. Other examples include findingthe determinant of a matrix (which is an alternating k-linear map on either the columns or the rows),matrix completion, and matrix factorization.12Step size selection Mat-vec multiplications CountFixed 1/L ∇ f = Aᵀ∇F , A∇ f 2Line search, e.g. Armijo, Wolfe as above 2Plane search with momentum direction as above 2Plane search with generic direction as above and Ap2 3Table 2.1: Number of matrix-vector multiplications per iteration of GD for objectives of theform f (x) = F(Ax) and second direction p2 in the case of a plane search. There is an ad-ditional matrix-vector multiplication at the first iteration for Ax0. Subsequent calculationsof Axk do not require matrix-vector multiplications.2.2 Linear maps2.2.1 Trace operatorThe trace of a square matrix A, denoted Tr(A), is given by the sum of its diagonal entries. Tr(A) isalso equal to the sum of the eigenvalues of A. This means that Tr(A) is a linear operator, a propertythat is advantageous in a linesearch. On iteration k,Tr(A(xk +αkPk)) = Tr(Axk)+αkTr(APk).Therefore, for trial values of the stepsize αk, Tr(Axk) and Tr(APk) only need to be calculated once.2.2.2 Linear composition problemsAn important class of problems is minimizing an objective in the form f (x) = F(Ax,x). The inputsto F are a linear map and the parameter x itself, which could be used to compute a regularizationterm such as ‖x‖2.1 The function F could be non-linear. For ease of notation, we simplify byomitting the second argument to F .A linear composition problem has an objective in the formf (x) = F(Ax) (2.2)where A ∈ Rm×n is the feature matrix. In other words, the objective is the composition of a linearfunction with another function. These models form an important class of optimization problemsin machine learning. For example, classification often involves SVMs and (binary or multiclass)logistic regression.1We could use A˜ =[AI]instead of using Ax and x, but want to recognize that multiplication by A is expensive but notby I.13In large-scale problems, A could be large and dense. If F is a relatively simple calculation, thenmatrix-vector multiplications Av, for different vectors v, could be the bottleneck of the optimizationproblem. Thus, we want to minimize the number of Av operations we perform. For typical iterativemethods, we might calculate three quantities at every iteration: objective value f (x), gradient ∇ f (x)and Hessian∇2 f (x). The structure of linear composition problems allow us to minimize the numberof matrix-vector multiplications for each of the three quantities.Gradient and Hessian The gradient of linear composition problems has a convenient form∇ f (x) = Aᵀ∇F(Ax). (2.3)Calculating the gradient results in one additional matrix-vector multiplication at the end of everyiteration.2 The Hessian has the formAᵀ∇2F(Ax)A. (2.4)Depending on the form of ∇2F(Ax), the Hessian could also be efficiently computed. For moredetails on the derivatives of f (x) = F(Ax), see Appendix A.2.1.Function value Consider the case where the next iterate is determined along a single search direc-tion (1.5). At each iteration, given the search direction pk, the function value at the next iteratef (xk+1) = F(A(xk +αk pk)) = F(Axk +αkApk) = f (vk +αk p′k)where the values vk =Axk and p′k =Apk are vectors that can be stored and reused. Thus, the functionvalue of f (xk+1) could be calculated with only two matrix-vector multiplications and several callsto f for many αk values. In other words, after the first iteration, every subsequent iteration of thelinesearch requires only a scalar-vector multiplication rather than a matrix-vector multiplication.In the case where more than one search direction is used, substituting (2.1) into (2.2) gives thevalue at the next iteratef (xk+1) = F (A(xk +Pkαk)) = F(Axk +APkαk) = f (vk +P′kαk)where P′k = APk is also stored and reused. For a Pk with d generic directions, computing APk requiresd matrix-vector multiplications. For example, going from a line search to a two-dimensional planesearch would result in one additional matrix-vector multiplication. But in special cases, the extradirections may not incur any additional cost. If the additional direction is the most recent momentumdirection xk−xk−1, then calculating APk is the same as calculating Axk and A(xk−xk−1). If we store2Aᵀ may not always be available, for example, in some sparse settings arising from PDEs.14Axk−1 from the previous iteration, we can perform a two-dimensional subspace search for the samecost as a one-dimensional linesearch. The SO αk = argminα f (vk +P′kα) is now a minimizationproblem in terms of αk. There are many methods we could use, including linesearch methodswhose iterations would also be very cheap.Line search conditions In the case of inexact line searches, checking for the sufficient decreasecondition (1.7) requires no additional matrix-vector multiplications. We calculate f (xk+1) as beforeand the directional derivative ∇ f (xk)ᵀpk multiplies two known vectors. The curvature condition(1.9), however, has a directional derivative term ∇ f (xk+αk pk)ᵀ∇ f (xk) that includes the gradient atthe next iterate. Calculating ∇ f (xk +αk pk) directly would require an extra matrix-vector multipli-cation at every iteration of the line search. We avoid this by noting that the curvature condition doesnot need the gradient but the directional derivative. In other words,∇ f (xk +αk pk))ᵀ∇ f (xk) = [Aᵀ∇F(Axk +αkApk +b)]ᵀ∇ f (xk)= ∇F(Axk +αkApk +b)ᵀ [A∇ f (xk)]and we are back in the scenario where the only matrix-vector multiplications needed (Axk, Apk andA∇ f (xk)) are already computed.GD example Returning to table 2.1, we count the number of matrix-vector multiplications forwhen pk = −∇ f (x). Initially, calculating Ax0 requires one matrix-vector multiplication. Becausethe linesearch requires calculating the function value of the next iterate f (xk+1) = F(Axk+1), we canstore the vector Axk+1 and reuse in the next iteration.We first consider the case where the constant L is known and where a fixed stepsize αk = 1/Lis used. Each kth iteration requires two matrix-vector multiplications. This comes from the updatestepxk+1 = xk−αk∇ f (xk) = xk−αkAᵀ∇F(Axk)which requires one matrix-vector multiplication, Aᵀ∇F , and from calculating the function valuef (xk+1) = F(Axk+1) = F (Axk−αkA∇ f (xk))which requires one matrix-vector multiplication A∇ f .Next, we consider the case where we perform a two-dimensional SO with the first directionbeing the gradient direction and the second is a generic direction pk,2. In other words, the stepsizes15and directions are given by αk =[αk,1αk,2]and Pk = | |∇ f (xk) pk,2| |. The update becomesxk+1 = xk−Pkαk = xk−αk,1Aᵀ∇F(Axk)+αk,2 pk,2.Like the fixed stepsize case, the update step also requires one matrix-vector multiplication Aᵀ∇F .The function valuef (xk+1) = F(Axk+1) = F (Axk−αk,1A∇ f (xk)+αk,2Apk,2)requires two matrix-vector multiplications A∇ f and Apk,2. Thus, a two-dimensional SO generallytakes three instead of two operations per iteration.If pk,2 is chosen to be a momentum direction, then Apk,2 = A(xk− xk−1) = Axk−Axk−1 couldbe stored from previous iterations. Thus, in the case of GD with momentum, the number of matrix-vector multiplications per iteration is only two. This gives us PHB with an inexact but accurate SOand thus is an attempt to achieve the CG convergence rate on this problem class “for free”. Note thatwe do not rely on knowing the constant L and we return to this point in Chapter 3.To conclude the section on linear map models, we give some examples below of common mod-els that take the form of equation (2.2). We denote the ith example and label as ai ∈ Rn (the ith rowin A) and yi ∈ R respectively.Logistic regressionIn the case of binary classification, if yi ∈ {−1,1}, then the logistic loss isf (x) =m∑i=1log(1+ exp(−yixᵀai)) = log(1+ exp(−y∗Ax)) ,where f : Rn→ R, x ∈ Rn, A ∈ Rm×n and ∗ is used to denote element-wise multiplication. In themulti-class case with c classes and labels yi ∈ {1, . . . ,c}, the softmax objective isf (x) =m∑i=1−xᵀyiai+ log(c∑j=1exp(xᵀj ai))where x j and xyi denote the weights of the jth class and of true class yi respectively.16SVMsA SVM is an example of a model with a non-smooth objective f that is given byf (x) = max{0,1− y∗Ax}+ λ2‖x‖2The linear composition problem is identified as a problem where SO is efficient (Narkiss andZibulevsky, 2005a). Below, by using our multilinear composition problem perspective, we identifyother problems where SO can be done efficiently.2.3 Log determinant problemsMany log-determinant problems involve the calculation of two matrix operations: the trace and thelog-determinant. As discussed, the trace operation is linear. The log-determinant function f : Sn→R, where Sn denotes the set of symmetric n×n matrices, is given byf (A) =logdetA positive definite A 0+∞ otherwise .In some situations, these two operations could be calculated efficiently. As discussed previ-ously, the trace operation could be computed efficiently. If the goal is to iteratively improve on Xwith rank 1 updates, then the Sherman-Morrison-Woodbury formula (Sherman and Morrison, 1950;Woodbury, 1950) allows us to try different step sizes αk inexpensively. Given a square, invertiblematrix Xk, stepsize αk and column vectors uk and vk, finding Xk+1 that gives an optimal value ofdetXk+1 becomesdetXk+1 = det(Xk +αkukvᵀk ) = (1+αkvᵀk X−1k uk)detXk.In many cases, because X−1k is the gradient of logdetXk and used as a descent direction, invertingXk is a necessary calculation and its results could be stored and reused.3 If we also save the re-sults of vᵀk X−1k uk and of detXk, trying out candidate next iterate values of Xk would involve only amultiplication of two scalars.We conclude this section by giving an example of log-determinant problems below.2.3.1 Gaussian negative log-likelihood (NLL)sMultivariate Gaussians generalize Gaussian distributions to n variables by taking a linear combina-tion of independent standard univariate Gaussians. In other words, if a(i) is the ith row in A ∈Rm×n,and if a(i) is drawn independently and identically from a univariate Gaussian with mean µ ∈ Rn3For this thesis, we set aside the issue where Xk is sparse and its inverse is not easily available.17and nonsingular covariance matrix Σ ∈ Rn×n 0, then the data follows a multivariate Gaussiandistribution. The probability density is given byp(a(i)|µ,Σ) = (2pi)−n/2 detΣ−1/2 exp(−12(a(i)−µ)ᵀΣ−1(a(i)−µ)). (2.5)Analogous to the standard deviation in the univariate case, the square root of the determinant of Σabove is used as a normalization constant. The ML problem of learning about the distribution of thedata can be formulated as finding the maximum likelihood estimator (MLE) of the parameters µ andΣ. To do this, we can minimize the NLL of (2.5) with respect to µ and Σ. The NLL could be writtenas the sum of a log determinant term and a trace termf (A|µ,Σ) =− log(p(A|µ,Σ)) =−m2logdetΣ−1+12Tr(SΣ−1)where S, ∑mi=1(a(i)−µ)(a(i)−µ)ᵀ. For more details, see Murphy (2012).2.4 Matrix factorizationA goal of many machine learning problems is to represent high dimension data with low dimensionmodels. Reasons for this include lower storage requirements and better interpretability. Given datarepresented by a high dimension matrix A ∈Rm×n, one approach is to factor A into a product of twolower ranked matrices. In other words, we want to express A as the product of two rank d matrices,Z ∈ Rm×d and W ∈ Rd×n, where A = ZW and d < n. One way to view matrix factorization isthrough principal component analysis (PCA) where Z is sometimes referred to as the “part weights”or “factor loadings” and W is referred to as a “part” or a “factor”. If A is big, directly factorizingA with singular value decomposition may be infeasible. Instead, we can use (stochastic) GD on anobjective that minimizes the difference between A and its lower dimension representation,minZ,Wf (Z,W ) = ‖A−ZW‖2F where ‖ · ‖F denotes the Frobenius norm.The objective f contains the multiplication ZW which costs O(mnd) to calculate, or O(d) fora fixed A. We can also view f as a multilinear composition problem where the computationallyexpensive operation is ZW . So, if we store ZW , the SO subproblem at each iteration does not incurmuch additional cost above the initial O(d). The larger the value of d, the higher the cost savingsfor the SO subproblem. Note that although d is generally much smaller than n, d could still be largein absolute terms. For more details, see chapter 4 and the appendix of Strang (2016).182.4.1 Matrix completionSuppose A is a matrix with missing entries. The objective of matrix completion is to find reasonablevalues for the missing entries with a model that performs well on the available entries. The usualexample given for matrix completion is the Netflix competition, where A is a sparse matrix withentries that correspond to user ratings of movies. Suppose the set R contains the indices of knownentries of A. Our objective would then be to minimize the error over R, orminZ,Wf (Z,W ) = ∑(i, j)∈R(a(i, j)− z(i)w( j))2where a(i, j) is the (i, j)th entry of A, z(i) is the ith row of Z and w( j) is the jth column of W . Bysimilar reasoning, the iteration cost is O(|R|d) where |R| denotes the size of the set R, and is O(d)for a fixed R. Collaborative filtering, a variation, is essentially the above with regularization on thesize of Z and W , orminZ,Wf (Z,W ) = ∑(i, j)∈R(a(i, j)− z(i)w( j))2+λ12‖Z‖2F +λ22‖W‖2Fwhere the computation of the last two regularization terms do not change the leading order of asymp-totic complexity. For more details, see part III.5 of Strang (2019).19Chapter 3Descent directions and stepsize methodsIn this chapter, we discuss primary and secondary search directions to include in a SO and methodsto find stepsizes. Primary directions are mainly search directions from popular methods such asGD and Newton’s method. Secondary directions are those that add favorable properties to a searchmethod. For example, a primary direction equal to the negative gradient and a secondary directionequal to the most recent momentum term results in accelerated rates for convex problems.A SO is primarily concerned with finding some “good” weighting between the primary directionand secondary directions. While popular methods that involve more than one search direction withprescribed weights exist, the stepsizes are often set with the assumption that parameters to theproblem are known. An example is NAG which assumes that the Lipschitz constant (and sometimesthe strong convexity constant) are known (Nesterov, 1983). In contrast, stepsize methods discussedin this chapter do not require knowledge of these constants. We give three examples of SO methodsin this chapter, including a d-dimensional version of the Wolfe conditions.In theory and in practice, practical issues such as initialization have an impact on the qualityand cost of both linesearches and SOs (Nocedal and Wright, 2006; Schmidt, 2011). We discussextending various initialization and interpolation methods used in linesearches to SO. We concludethis chapter with preliminary convergence analysis on GD with SO.3.1 Descent directionsThe material for this section is mainly drawn from Nocedal and Wright (2006).3.1.1 Primary directionsIn the case of one descent direction pk with stepsize αk, the iterations are in the formxk+1 = xk +αk pk (3.1)20We list a few possibilities for pk below.Current gradientA natural choice for a primary direction is the negative gradient at the current iterate, or pk =−∇ f (xk). At xk, this is the direction where f decreases the most rapidly. For more information, seeAppendix A.3.1. For Barzilai Borwein (BB), see Appendix A.4.1.Newton’s directionNewton’s direction is an option for twice differentiable objectives. Here, pk =−(∇2 f (xk))−1∇ f (xk).If∇2 f (xk) 0, then Newton’s direction is a descent direction. Many variations exist to handle caseswhere ∇2 f (x) is not positive definite. For more information, see Appendix A.3.2.Quasi-NewtonBecause Newton’s direction requires computing the Hessian, which could be computationally ex-pensive, we could choose a quasi-Newton direction which uses an approximate inverse Hessian.In other words, set pk = −Hk∇ f (xk) where Hk ≈(∇2 f (xk))−1. Popular versions of quasi-Newtonmethods include symmetric rank one (SR1), Broyden Fletcher Goldfarb Shanno (BFGS) and David-son Fletcher Powell (DFP). For details, see Appendix A.3.3.Limited memory Broyden Fletcher Goldfarb Shanno (L-BFGS) is a memory efficient approxi-mation of BFGS. Instead of storing the full n× n matrix Hk, L-BFGS stores c < n vectors that areused at each iteration to approximate Hk. The tradeoff is that the rate of convergence becomes linearwhereas BFGS achieves superlinear convergence when some conditions are satisfied. For details,see A.3.4.Define sk , xk+1− xk and yk , ∇ f (xk+1)−∇ f (xk). BFGS finds the new approximation to theinverse Hessian Hk+1 with a rank two update to HkHk+1 =Vᵀk HkVk +ρksksᵀk (3.2)whereρk =1yᵀk sk, Vk = I−ρkyksᵀkL-BFGS uses the latest m pairs of {si,yi} on some choice of initial H(0)k to approximate Hk via thesame rank two update rule (3.2). The choice of H(0)k could vary between iterations. For details, seeAppendix A.3.5.A particular choice of L-BFGS is a memoryless version where H(0)k is set to the identity matrixmultiplied by a scalar and m = 1. This results in a Hk+1 that could be considered a symmetricpositive definite version of the Hestenes-Stiefel form of non-linear CG . For details, see Appendix21A.3.6.Non-linear CGNon-linear CG methods were developed to extend beyond quadratic objectives while retaining someof the desirable properties of CG. As the name suggests, the search directions of successive iterationsof CG are conjugate. This property is approximately retained in non-linear CG where the searchdirection has the formpk =−∇ f (xk)+βk pk−1 (3.3)and βk is a scalar that ensures that pk and pk−1 are conjugate. Popular forms of non-linear CG includeDai-Yuan, Fletcher-Reeves, Polak-Ribiere-Polyak and Hestenes-Stiefel. For details, see AppendixA.3.7.TNLike quasi-Newton methods, TN is an inexact Newton method where the search direction is in theform pk = −B−1k ∇ f (xk) where Bk is an approximation of ∇2 f (xk). At every iteration, TN usesCG or Lanczos to solve for Bk to some required accuracy.1 A sequence of coefficients, sometimesreferred to as a forcing sequence, gives the termination condition for this subproblem. For details,see Appendix A.3.8.3.1.2 Secondary directionsThe secondary direction could be chosen to complement the primary direction. For example, thecurrent gradient as a primary direction could be combined with a TN direction so that the searchtakes advantage of local curvature information. Conversely, Newton’s direction could be combinedwith a gradient direction to ensure global convergence or to improve convergence speed when farfrom a stationary point. The primary directions in the previous section could be used as a secondarydirection. We list some other secondary directions below.MomentumWe touched on momentum directions, x j−x j−1 where j = 1, . . . ,k, in Chapter 2. There, we showedthat we can increase the search space for objectives in the form (2.2) to include a momentum direc-tion without adding to iteration complexity.There are many theoretical reasons for using momentum as a secondary direction. Although thisthesis does not consider the stochastic setting, a momentum term can be shown to partially “cancelout” the effect of the noise term for stochastic GD (Defazio, 2020).1As CG assumes that Bk is positive definite, TN uses a modified version of CG that checks for negative curvature.22Nemirovsky directionsIn the case of quadratic objectives, CG uses a subspace spanned by the current and previous gra-dients that are mutually orthogonal. For convex L-smooth objectives, Nemirovsky (1982) relaxesthe orthogonality requirement so that the current gradient is only orthogonal to the two Nemirovksydirectionsp(1)k = xk− x0andp(2)k =k−1∑i=0wi∇ f (xi) where wi =1 i = 012 +√14 +w2i−1 i > 0Thus, at every iteration, the SO involves a two-dimensional subspace instead of the expanding sub-space used by linear CG. This gives optimal complexity for convex problems, as well as for strongly-convex problem with knowledge of µ and with restarts. Narkiss and Zibulevsky (2005a) gives anexample of using Nemirovsky directions in a subspace search.3.2 Stepsize methodsIn this section, we discuss three search methods that could be used to solve the SO subproblem. Thematerial from this section is drawn from Lee (1968); Schmidt (2005); Nocedal and Wright (2006).3.2.1 QuasilinearizationAs mentioned in Chapter 1, quasilinearization is a generalization of the Newton-Raphson methodwhere the solutions are functions instead of values. In general, it is a method to solve linear andnon-linear boundary problems, which are differential equations with a set of additional constraintsthat are referred to as boundary conditions. Having boundary conditions allows for unique solutionsto the problem. Quasilinearization is the baseline technique used in both Miele and Cantrell (1969a)and Cragg and Levy (1969) to solve their SO, to which they add steps that ensure that the stepsizesfound decrease the objective at every iteration. They accomplish this via a scaling factor ρ thatresembles the multiplier used in backtracking Armijo.Given a non-linear function, quasilinearization generates a system of linear algebraic equations.It then starts with an arbitrary initial solution and iteratively applies corrections until stopping con-ditions are met. At a very high level, the technique solves a boundary value problem in two steps:1. The problem is transformed from a boundary-value problem to a set of initial-value problems.In the general case, there are no analytical solution to these problems and so the solution iscomputed numerically.2. Then, find the integration constants of the differential equation through solving a set of linearalgebra equations determined from the boundary conditions.23Note that this method requires computing the first and second derivatives of the objective withrespect to the stepsizes. The calculations may be done analytically or numerically.As a concrete example of quasilinearization, we summarize the two-dimensional search usedat every iteration of the memory gradient method in Miele and Cantrell (1969a). Suppose our nextiterate is given by (1.16) and we want to find stepsizes that minimizes our objective at the next step.Let φ(α,β ) be in the objective in terms of stepsizes and let α0 and β 0 be initial stepsize guesses.The SO subproblem isαk,βk = argminα,βφ(α,β ) = argminα,βf (xk−α∇ f (xk)+β (xk− xk−1))Quasilinearization iteratively refines the stepsizes to some degree of accuracy. We denote the kthiterates of the subproblem with superscripts. For example, αk and β k are the k iterates in the innerloop and we omit the index to the outer loop for clarity.It can be shown that the correction to the previous stepsizes are given byαk = αk−1−ρD1D3 sign(D4D3)(3.4)andβk = βk−1−ρD2D3 sign(D4D3)(3.5)where sign(x) is +1 for positive x and -1 for negative x and ρ is the scaling factor (Miele andCantrell, 1969b). ρ starts at 1 and is multiplied by a constant between 0 and 1, such as 0.5. Similarto backtracking Armijo, we accept the first ρ that gives φ(αk+1,βk+1)< φ(αk,βk).2 Let subscriptsof functions denote partial derivatives. For example, φα and φβ are the first partial derivativesof φ . The terms Di are differences in partial derivatives calculated at αk−1 and βk−1 where D1 =φαφββ−φβφαβ , D2 = φβφαα−φαφαβ , D3 = φααφββ−φ 2αβ , and D4 = φ 2αφββ−2φαφβφαβ+φ 2βφαα .Equations (3.4) and (3.5) are repeated until improvements between iterations fall below an accuracythreshold.There is enough material on quasilinearization to cover an entire textbook. For example, it canbe proven that the technique converges monotonically and that there are relations to multiple othermethods. The choice of different stopping conditions and initial values also make a difference to theperformance of the algorithm. The details are beyond the scope of this thesis. For more informationon the quasilinearization technique, see Lee (1968). For more details of the example above, seesection 6 of Miele and Cantrell (1969a) and section 6 of Cragg and Levy (1969).2In practice, this condition is very similar to the sufficient decrease condition (1.7) because the constant c1 is generallyset to a very small number.243.2.2 SolversA SO is an optimization subproblem that needs to be solved at every iteration of a descent method.Just as we can call a black box solver on the main problem (or outer loop), we can also use a solverfor the subproblem (inner loop). This is the approach used in Narkiss and Zibulevsky (2005a,b);Zibulevsky (2008) where the SO is solved by calling the Matlab solver minFunc (Schmidt, 2005).The solver used for the SO subproblem is not required to be the same as that used for the mainproblem. In our experiments, we attempt both calling our software recursively, which allows us totake advantage of the structure of the subproblem, as well as calling a version of our software thatis structure-independent. We discuss this further in Chapter 4.3.2.3 d-dimensional WolfeThe strong Wolfe conditions could be extended from a linesearch to a SO. Suppose that, at iterationk, we have solution xk ∈ Rn, d directions stored as columns of Pk ∈ Rn×d and stepsizes αk ∈ Rd .The sufficient decrease condition (1.7) becomesf (xk +Pkαk)≤ f (xk)+ c1∇ f (xk)ᵀ(Pkαk), (3.6)and the curvature condition (1.8) becomes|∇ f (xk +Pkαk)ᵀ(Pkαk)| ≤ c2|∇ f (xk)ᵀ(Pkαk)|. (3.7)Intuitively, the conditions are now on the average direction Pkαk instead of a single direction pk.Equation (3.6) was previously used in Conn et al. (1994) but equation (3.7) is new.One challenge that needs to be addressed in order to use the above d-dimensional Wolfe con-ditions is initialization and interpolation. We attempt to address this in the next section. Anotherchallenge is whether or not the directions contained in Pk are sufficiently linearly independent. Al-though this is not addressed in this thesis, the topic is of both practical and theoretical concern andwarrants future attention.3.3 Initialization and interpolation methodsThe initial stepsize estimate and how to move from one trial stepsize to the next may appear to be anunimportant detail. Yet, suitable initialization and interpolation can have a significant impact bothin theory and in practice. In this section, we summarize some known initialization and interpolationmethods for linesearches and discuss how they might be generalized to the SO setting. Materialon initialization and interpolation methods for linesearches is drawn from Schmidt (2011); Nocedaland Wright (2006).253.3.1 InitializationDenote the initial stepsize choice for the SO subproblem at every iteration by α0 ∈ Rd . Again, weuse superscript to denote the iteration count for the inner loop. For example, αk is the stepsize foundin the kth iteration of the SO. To start, suppose d = 1 and we have a linesearch. The choice of α0depends on a few considerations, for example, the iterative method used in the outer loop, the costof calculating function values and gradients, etc. We summarize some of these considerations inthis section.For Newton and quasi-Newton methods, a stepsize of 1 is automatically accepted in the regionof quadratic convergence. Thus, for these methods, unit initialization α0 = 1 is necessary to takeadvantage of faster convergence around the solution x∗. For gradient-based methods, we may wantto include more information and arrive at a more informed α0. There are several heuristics.For example, if we assume that first-order changes are roughly equal between sequential itera-tions, then we can setαk∇ f (xk)ᵀpk = αk−1∇ f (xk−1)ᵀpk−1 (3.8)Then it makes sense to set initial estimate α0 equal to αk in (3.8) and soα0 =αk−1∇ f (xk−1)ᵀpk−1∇ f (xk)ᵀpk. (3.9)As shown in Chapter 2, for some problem classes such as linear composition problems, the direc-tional derivative∇ f (xk)ᵀpk may be inexpensive to compute even though the gradient∇ f (xk)may becostly. In this situation, it is possible to better initialize the linesearch without incurring significantcost.Another method is known as Hermite polynomial interpolation (HPI), which has been shownto potentially reduce the number of linesearch evaluations significantly (Schmidt, 2011). Unlikeother interpolation methods (such as Newton’s), HPI uses both function value and derivatives. Itis not required to have a derivative value for every point x where a function value f (x) is used inthe interpolation. To initialize the linesearch, one approach is to use information on function valuesf (xk−1) and f (xk), and on the previous directional derivative ∇ f (xk−1)ᵀpk−1. Because no newgradient calculation is necessary, this method works well for problems where gradient calculationsare more costly than function values. Let φ(α) equal the objective f (xk) = f (xk−1 +α pk−1) withthe stepsize α as input. HPI works by fitting a quadratic approximation of φ , φq, through these threepoints and then sets α0 as its minimizer. This givesα0 =2( f (xk)− f (xk−1))∇ f (xk−1)ᵀpk.For details on how α0 was derived, see Appendix A.3.9.26If d > 1, what should be the equivalent of the initialization methods above? Because the solutionto a linesearch is also the solution to a subspace search, one way to initialize a d-dimensional SOis by setting the initial weight of the primary direction to the result of a linesearch and the initialweight of all the secondary directions to zero. This ensures that the SO will produce a next iteratethat decreases the objective at least as much as that found by a traditional linesearch on the primarydirection.Are there other ways to extend the three initialization methods to higher dimensions? The equiv-alent of the unit initial stepsize is especially important for d-dim Wolfe because of its connectionto Newton and quasi-Newton methods (see Appendix A.3.4.) Setting all stepsizes to one, α0 = 1d ,would mean that the size of the initial estimated step grows as d grows. This seems incorrect. If wedo not want this, we need to somehow normalize α0 by d. A natural way to do this is to initializeall elements of α0 to 1/d but it is also unclear whether equal weights across directions make sense.For SO along d directions, the first-order change assumption in (3.8) becomes∇ f (xk)ᵀPkαk = ∇ f (xk−1)ᵀPk−1αk−1Denote Dk and Dk−1 as the d × d diagonal matrices that contain ∇ f (xk)ᵀPk and ∇ f (xk−1)ᵀPk−1along the diagonal entries respectively, and the SO version of (3.9) becomesα0 = D−1k Dk−1αk−1As in the case of a linesearch, this initialization method for SO may be computed cheaply.3In the case of HPI, the intuitive extension to α0 ∈ Rd and Pk ∈ Rn×d would be to defineφ(α) = f (xk−1 +Pk−1α) and to use interpolation conditions f (xk−1) = φ(0d), f (xk) = φ(α), andPᵀk−1∇ f (xk−1) = φ′(0d). The last condition comes from a derivative of a scalar by a vector, φ ′(α) =Pᵀk−1∇ f (xk−1+Pk−1αk−1). This produces d points instead of a single point as in the linesearch case.The challenge is that if we also want to fit a quadratic in the SO case, we have too many points. It isunclear whether we should find heuristics to combine the d points into one and proceed somewhatidentically to the linesearch case, or increase the dimensionality of the polynomial used in HPI tomatch the dimensionality of the subspace. The former may be more efficient computationally as wewould be fitting to fewer points. Also, finding a minimizer to a quadratic is much simpler than to ahigher degree polynomial. The latter approach, however, does not depend on heuristics to combined points and may be more theoretically sound.Finally, it was suggested that a least-squares and piece-wise polynomial interpolation may beuseful approaches when d gives a reasonable number of points. Although not addressed in thisthesis, it is an interesting line of inquiry for future research.3It is inexpensive to compute directional derivatives for linear composition problems and Dk is easy to invert.273.3.2 InterpolationThis section discusses how to move from α0 to α1 and subsequent trial stepsizes. We begin withSO where d = 1. At every iteration, a linesearch takes an interval that contains a satisfying stepsize,chooses a trial stepsize within that interval and either accepts the trial stepsize or reduces the size ofthe interval and repeats. For example, in the case of the Armijo condition, if αk−1 is not a satisfyingstepsize, then αk lies in the interval (0,αk−1). One popular approach for picking αk within thisinterval is backtracking, or setting αk = ραk−1 where ρ ∈ (0,1). This does not, however, use anyinformation gained from the search improve the choice of αk.As in initialization, an interpolation method such as HPI can help incorporate recent searchinformation. Let φ(α) = f (xk−1+α pk−1). The Armijo condition (1.7) could be rewritten asφ(αk)≤ φ(0)+ c1αk−1φ ′(0)Again, we can fit a quadratic approximation φq to φ and interpolate between points φ(0), φ ′(0) andφ(αk−1). Taking steps similar to Appendix A.3.9, we find the next trial stepsizeαk =− φ′(0)(αk−1)22(φ(αk−1)−φ(0)−φ ′(0)αk−1)Extending this to the d > 1 poses the same challenges as in the initialization case.For more details on initialization and interpolation, see Nocedal and Wright (2006) and Schmidt(2011).3.4 GD and exact SOThis section contains preliminary convergence analysis for applying SO to GD. As a first step,we compare GD using linesearches with GD using subspace searches. The goal is unconstrainedminimization (1.1) where the objective f is continuous, differentiable and L-smooth (1.2). Denotethe ith column (search direction) in Pk ∈ Rn×d by Pk,i ∈ Rn. We assume that one of the searchdirections in Pk is the gradient direction. In other wordsPk,i =−∇ f (xk) for some i ∈ 1, . . . ,d (3.10)where d is the number of search directions. We also assume that all columns of Pk (search directions)are independent.3.4.1 Quadratic upper bound with Euclidean normWith assumption (3.10), because any solution to a linesearch is also a solution to a subspace search,SO should do at least as well as GD. Without loss of generality, let the first column of Pk be the28negative gradient direction, or Pk,1 =−∇ f (xk). Then, SO need only set the first coordinate of αk to1/L and the other coordinates to zero to exactly match GD with an exact linesearch. In other words,f (xk+1) = minαkf (xk +Pkαk)≤ f(xk− 1L∇ f (xk)+d∑i=20 ·Pk,i)= f (xk−1/L∇ f (xk))which gives the guaranteed progress of GDf (xk+1)− f (x∗)≤ f (xk)− f (x∗)− 12L‖∇ f (xk)‖2. (3.11)For more details see Nesterov (2018). Arriving at (3.11) means that GD and exact SO does at leastas well as GD and exact linesearch. If f was also µ-strongly convex, we could combine (3.11) with‖∇ f (x)‖2 ≥ 2µ( f (xk)− f (x∗)), a consequence of (1.3), to arrive at a bound in terms of µ and L.Given (3.11), the question is whether we could do better than GD. If we follow the same ap-proach as standard convergence analysis for GD, it can be shown that SO cannot make any improve-ment. To see this, consider the quadratic upper bound based on (1.2) used to establish progressf (xk+1)≤ f (xk)+∇ f (xk)ᵀPkαk + L2‖Pkαk‖2 (3.12)Because Pk,1 =−∇ f (xk), GD would set αk with 1/L in the first coordinate and zero elsewhere into(3.12). This givesf (xk+1)≤ f (xk)− 12L‖∇ f (xk)‖2which is optimal.3.4.2 Quadratic upper bound with general quadratic normTo try to show that SO can do better than GD, one approach is to use a more general form of(1.2). We define smoothness not in terms of the Euclidean norm but a quadratic norm of the form‖x‖H =√xᵀHx where H is a positive-definite matrix‖∇ f (x)−∇ f (y)‖H−1 ≤ L‖x− y‖H , for all x and y.(The dual norm of ‖·‖H is ‖·‖H−1 .) This is the approach taken in Section 3.3 of Nutini et al. (2017).We rewrite (3.12) asf (xk+1)≤ f (xk)+αᵀk Pᵀk ∇ f (xk)+L2‖Pkαk‖2H (3.13)GD is optimal where H is the identity matrix but SO may help in other cases. Finding such an Hfor a general subspace, however, may not be straight-forward. An example of H where SO will nothelp and GD with a linesearch is optimal is H = γI where γ is a scalar. In this case, GD remains an29optimal algorithm because GD readjusts its optimal stepsize by a factor of γ . Exact SO givesαk =− 1Lγ (Pᵀk Pk)−1Pᵀk ∇ f (xk) (3.14)Substituting αk into (3.13) givesf (xk +Pkαk)≤ f (xk)− 1Lγ∇ f (xk)ᵀPk(Pᵀk Pk)−1Pᵀk ∇ f (xk)+L2(− 1Lγ∇ f (xk)ᵀ)Pk(Pᵀk Pk)−1Pᵀk γIPk(Pᵀk Pk)−1Pᵀk(− 1Lγ∇ f (xk))which reduces tof (xk+1)≤ f (xk)− 12Lγ ‖Pᵀk ∇ f (xk)‖2(Pᵀk Pk)−1 (3.15)If Pk,1 = −∇ f (xk), GD with a linesearch simplifies (3.14) to αk =[1Lγ 0 . . . 0]ᵀand (3.15)simplifies tof (xk+1)≤ f (xk)− 12Lγ ‖∇ f (xk)‖2 = f (xk)− 12L‖∇ f (xk)‖21/γwhich is optimal.Three terms, Pᵀk Pk, (Pᵀk Pk)−1 and Pk(Pᵀk Pk)−1Pᵀk , appear frequently in the analysis of SO usingquadratic upper bounds. Pᵀk Pk is a symmetric matrix where the diagonals contain ‖Pk,i‖2 and theoff-diagonals contain Pᵀk,iPk, j where i, j ∈ {1, . . . ,d}. As discussed in Chapter 1, the inner productbetween two direction vectors measures similarity. Thus, one way to think of Pᵀk Pk and (Pᵀk Pk)−1is a “covariance” and “precision” matrix respectively between search directions. Furthermore, thequantity P˜k , Pk(Pᵀk Pk)−1Pᵀk is a projection matrix because P˜2k = P˜k and P˜ᵀk = P˜k.4 This is notgenerally true if H is not in the form H = γI.It is unsurprising that the term P˜k appears because exact SO to find αk on a quadratic upperbound, or the righthand-side of (3.13), solves a least squares problem. In particular, if H = γI, theequivalent linear system to the right hand side of (3.13) isγL2(Pᵀk Pk)αk =−Pᵀk ∇ f (xk). (3.16)This could be interpreted as finding the orthogonal projection of the vector of directional derivatives,Pᵀk ∇ f (xk), onto the subspace formed by the search directions, Pᵀk Pk. The system (3.16) is satisfiablewhen Pᵀk ∇ f (xk) is contained within the subspace Pᵀk Pk. This gives an intuitive explanation of whythe current gradient should be one of the search directions, or assumption (3.10).4Assuming that d ≤ n, P˜k is a rank d matrix and if n = d then P˜k becomes the identity matrix.303.4.3 Higher-order upper boundAnother approach is to move from quadratic bounds to higher order bounds. For example, insteadof minimizing an upper first-order approximation of f as in (3.12), we could minimize an up-per second-order approximation as Nesterov and Polyak (2006) did for their analysis of Newton’smethodf (x)+∇ f (x)ᵀ(y− x)+ 12(y− x)ᵀ∇2 f (x)ᵀ(y− x)+ M6‖y− x‖3. (3.17)Minimizing (3.17) is generally more complicated. The objective f is usually assumed to be twice-differentiable with a M-smooth Hessian‖∇2 f (x)−∇2 f (y)‖ ≤M‖x− y‖ for all x and y.Furthermore, exact linesearch or SO on (3.17) means solving a system of nonlinear equations.31Chapter 4ExperimentsThis chapter describes the optimizer software minFuncLinear and provides some experimentalresults that support the analysis that a SO could find a more accurate solution in less time for someproblem classes. Linesearches are commonplace in ML while SO are used to a much less extent.Thus, much less is known about the practical aspects of performing a SO. A main goal of theexperiments is to explore the impact of different SO parameter settings. We consider the impact of• different primary and secondary directions,• increasing the dimensionality of the SO,• different initialization and stopping conditions for the SO, and• the choice of the descent method used for the SO.The experimental results suggest that SO should be initialized with a one-dimensional linesearchand terminated using relative accuracy. Furthermore, using momentum as the second direction andBB as the SO solver seem to have the best performance.As described in section 2.2.2, minFuncLinear is developed for problems that minimize anobjective function with a linear map component. We look at the performance gain from optimizingfor this structure. We also briefly explore whether including a SO reduces the risk of divergencegiven a poor choice of an initial solution.The experiments are run on two real datasets sido and spam. The sido dataset comes from theCausality Workbench website (Guyon, 2008) and contains 12,678 datapoints and 4,932 variables.The spam dataset was prepared by Carbonetto (2009) using the TREC 2005 corpus Cormack andLynam (2005) and contains 92,189 datapoints and 823,470 variables. The task is binary classifica-tion and we do not split the data into training and test set. Instead, we look exclusively at how wellwe minimize a logistic regression objective during training.324.1 minFuncLinear softwareThe minFuncLinear software developed for the thesis is heavily based on the software packageminFunc (Schmidt, 2005).1 minFuncLinear differs from minFunc because it1. has less functionality implemented (for example, fewer choices of search methods, initializa-tion and interpolation, no support for preconditioning),2. is written in Julia (Bezanson et al., 2017),3. is optimized for objectives that have a linear map component,4. allows for d-dimension SOs.minFuncLinear is available to download at https://github.com/sheaws/minFuncLinear.4.1.1 Features, usage and exampleThe descent methods implemented are GD, BB (see A.4.1), CG (Hestenes-Steifel), TN (linesearchNewton-CG) and L-BFGS. The reason that these methods are chosen is that they give at least oneeach of CG, quasi-Newton and inexact Hessian methods. Step sizes could be set according to inexactlinesearch methods (Armijo backtracking and Wolfe conditions) or inexact SO methods (calling asolver, using minFuncLinear itself or satisfying d-dimension Wolfe conditions). Most param-eters (such as search accuracy, number of L-BFGS corrections to use, etc) have default values thatcould be overridden by the user.minFuncLinear takes as input user-defined objective and gradient functions, choice of de-scent method and stepsize search and parameter settings and, if successful, returns a function valueand minimizer. example log Reg.jl gives an example of how to use minFuncLinear forbinary classification using logistic regression.4.1.2 Performance gain from optimizing to problem structureTo demonstrate the benefits from taking into account the linear structure of the objective, we ran asmall experiment using GD combined with a linesearch or a SO. The comparison is performed onbinary classification of the spam database using a midrange laptop. We look at three cases:1. GD with a linesearh using backtracking Armijo (GD,Armijo),2. GD with a linesearch satisfying the strong Wolfe conditions (GD,Wolfe), and3. GD with a two-dimensional SO that uses the most recent momentum term as the secondarydirection (GD,2d-SO).The SO is solved by calling a solver that uses GD and backtracking Armijo to return the stepsizes.The three cases are solved by minFuncLinear as well as a version of minFuncLinear thatdoes not take into account the structure of the objective function. Most notably, the non-optimizedversion does not economize on the number of matrix-vector multiplications.1https://www.cs.ubc.ca/∼schmidtm/Software/minFunc.html33Method f (x40) Time (sec) Outer Iters Inner Iters Mat-VecGD,Armijo 3,433.371 9.2 40 50 102GD,Wolfe 5,519.416 15.5 40 68 232GD,2d-SO 96.155 303.1 40 3,523 7,327GD,Armijo* 3,433.371 3.0 40 50 82GD,Wolfe* 4,307.752 4.6 40 163 82GD,2d-SO* 16.963 12.9 40 2,914 2,507Table 4.1: Binary classification of the spam dataset where * denotes the optimized version.Comparison of minimum value found, total time taken, number of iterations and stepsizesearch iterations and number of matrix-vector multiplications.The results are summarized in table 4.1 and in figure B.1. The spam dataset provides a rel-atively large features matrix A and thus the time savings are significant (over 20x faster in time).For the remaining experiments, we will only use the optimized version of minFuncLinear .We will also be more specific about SO settings. The shorthand used to describe a SO is in theform <primary direction>,< n >d-<secondary directions>-<subproblem method>. For example,GD,2d-Mm-BB indicates a two-dimensional SO performed with the current gradient and momen-tum as primary and secondary directions respectively, and where the SO subproblem is solved withBarzilai-Borwein.4.2 ExperimentsUnless otherwise specified, the following experiments use these minFuncLinear settings: Thenumber of L-BFGS corrections used is 20. The constants for the Wolfe conditions are set to c1 =10−4 and c2 = 0.9. The backtracking Armijo reduces the trial stepsize not by a fixed constant ρ butby quadratic interpolation. The descent algorithm terminates when either the maximum number ofiterations is reached (400 for spam and 1000 for sido), the value of the objective decreases by lessthan 10−12 or when the infinity norm of the gradient (the largest coordinate) has a magnitude of lessthan 10−9.There are also some items to note about the reported results. The number of linesearch iterationsfor a SO include both its outer and inner iterations. The number of matrix-vector multiplicationsinclude those in the SO although the matrices used in the subproblem are smaller (n× 2 insteadof m× n).2 GD and non-linear CG use quadratic initialization while L-BFGS and TN used Newtonstep initialization. All methods use quadratic interpolation. The y-axis of the charts display theoptimality gap, which is the difference between the objective value of the final iterate and the lowest2From the discussion in chapter 2, the number of matrix-vector multiplications for the linesearch methods Armijo andWolfe is always twice the number of iterations plus two.34objective value found across all methods. The x-axis contains either number of iterations or time inseconds. Tables display the final objective value, the time taken in seconds, the number of iterationsof the search method, the total number of linesearch iterations summed across the outer iterationsand the number of matrix-vector multiplications performed.4.2.1 From linesearch to SOThe first experiment establishes our baseline for going from a linesearch to a SO. Figures 4.1 and4.2 show the convergence paths of four descent methods: GD, non-linear CG, L-BFGS and TN.Each method sets stepsizes according to either backtracking Armijo, Wolfe conditions or a two-dimension SO. The SO uses the most recent momentum term as the second direction and usesminFuncLinear to find a combination of stepsizes that minimizes the next function value alongthe two directions. Under the hood, minFuncLinear solves the SO subproblem inexactly withdescent method BB and backtracking Armijo. Further details of this experiment are given in tablesB.1 and B.2.For methods GD and non-linear CG, going from a linesearch to a SO resulted in higher accuracywith both datasets. The improvements varied slightly depending on dataset and linesearch type butwere of the same magnitude. For example, using GD for the spam dataset, going from backtrackingArmijo and Wolfe to a SO lowered the objective value from 361.214 and 241.925 respectively to0.00150. Thus, with GD, a correct result was only obtainable via using SO. In terms of effort,SO consistently resulted in the fewest number of iterations. In terms of time, using SO is often bothfaster and more accurate. These results suggest that for binary classification with logistic regression,going from a linesearch to a SO for GD and non-linear CG could significantly improve efficiency andaccuracy.For L-BFGS and TN, the results are more mixed. With the exception of L-BFGS combined withbacktracking Armijo (which stalled at a solution far from the optimal), switching to a SO did notresult in a marked improvement with either datasets.3The combination of L-BFGS and Wolfe is already both fast and accurate, achieving the lowestobjective with only 130 iterations for spam and the second lowest objective with 170 iterations forsido. We will show later that SO can improve on this. Regardless of the impact of a SO, the choiceof the primary descent direction is of first order importance. As a comparison, GD+SO took roughly47 seconds to find a less accurate solution whereas L-BFGS+SO found the best solution in roughly2.5 seconds. This difference could also be explained by the particular combination of primary andsecondary directions.The shape of the convergence curve is a potential point of interest. When using a SO, somemethods display a stepwise decrease in the objective value. This is particularly noticeable in the GD3L-BFGS with Armijo stalled on the 25th and 48th iteration for spam and sido respectively because the function valueremained the same up to 12 decimal points.35Figure 4.1: Binary classification of the spam dataset using various descent methods.plots in the sido dataset. It could be interesting to explain this theoretically. In the quadratic case,PHB with an accurate SO results in CG iterates, which remove the effect of clusters of eigenvalues.Perhaps a similar effect could occur for non-quadratic objectives in practice.The convergence charts show some methods “flattening out” before they reach the maximumnumber of iterations. For example, in the third chart on the left of Figure 4.2, the line in bluerepresenting L-BFGS with an Armijo linesearch becomes flat after roughly 25 iterations. In theexperiments, on top of using first-order optimality conditions for terminating, we also require thata minimum progress of 10−12 in function value between iterations. Thus, a method that flattensout before reaching the maximum number of iterations is one that stopped searching due to lack ofprogress made.36Figure 4.2: Binary classification of the sido dataset using various descent methods.4.2.2 Choice of secondary directionsFor first order methods, a potential benefit of using a SO is that it allows including a secondarydirection that captures curvature information. For example, we may want to combine a GD primarydirection with Newton’s as a secondary direction. Newton’s method is not currently implementedin minFuncLinear but we could approximate this with an inexact Hessian method such as TN.In this experiment, we explore the effect of replacing momentum with TN as a secondary direction.Our experiment show that there is not much of an improvement from going to a TN secondarydirection. In some case (see GD on the spam dataset), the results are much worse. This couldbe because of the many theoretical reasons to choose a momentum direction (see chapter 3). Orperhaps an inexact Hessian is insufficient or that we need the TN approximation to be better.Further details of the experiments are given in figures B.3 and B.2 and tables B.3 and B.4.374.2.3 Choice of SO descent methodIn this experiment, we explore the effect of different descent methods used in the SO. We keep themomentum direction as the second search direction. We add descent methods GD and TN to thebenchmark BB used for the SO. Figures B.5 and B.4, tables B.5 and B.6 contain the results. Themain takeaway of this experiment is that the choice of SO parameters do matter.For these two datasets, the choice of descent method to use within the SO makes a significantimpact. Choosing a more suitable method results in significant gains both in efficiency and accuracy.For example, with the sido dataset, using GD as the outer loop method and TN as the inner loopmethod is several orders of magnitude more accurate than using GD as both the outer and inner loopmethod. The former combination also runs much faster, completing the task in 2.3 seconds insteadof 62 seconds. The same effect could be seen in the spam dataset. The difference in accuracy is ofsimilar magnitude and using TN takes 11 seconds to complete instead of the 132 seconds that it takesGD to complete. In this experiment, GD seems to be strictly worse than TN and BB. It is unclear,however, whether TN or BB is the better choice. For example, TN,2d-Mm-TN is significantly worsewith the spam dataset but better with the sido dataset.A followup topic to explore is an explanation for why, with logistic regression, using GD as theSO method performs so poorly. Conversely, we could also consider the theoretical reasons for whyTN and BB performed so much better. The difference between BB and GD is particularly interesting.Does BB work so well because, under some conditions (see A.4.1), BB achieves local superlinearconvergence? It would be interesting to come up with theoretical reasons why certain combinationsof outer and inner loop descent methods make sense and to verify theory with more classes ofoptimization problems.4.2.4 Initialization and stopping criteriaA SO is a generalization of a linesearch and thus should produce a solution that is no worse than alinesearch. One way to achieve this is through suitable initialization. For example, instead of startingeach SO with random stepsizes, we could initialize the SO with the stepsize from a linesearch on theprimary direction and zeros for the secondary direction(s). We refer to this as 1d-init.Similarly, there are potentially better exit conditions for the SO. Instead of fixed stopping condi-tions, we could start with looser conditions for satisfying the SO and progress to tighter conditions asthe estimate of the solution becomes increasingly accurate. This is simply multiplying the progressand optimality thresholds by the infinity norm of the gradient. We refer to this relative stoppingcondition as RelStop.In this experiment, we compare four settings that correspond to all combinations of 1d-init andRelStop: 2d-Mm-BB-1 uses neither 1d-init nor RelStop, 2d-Mm-BB-2 uses 1d-init but not RelStop,2d-Mm-BB-3 does not use 1d-init and uses RelStop and 2d-Mm-BB-4 uses both 1d-init and RelStop.Figures B.6 and B.7 and tables B.7 and B.8 detail the results of our experiment. For readability,38we omit all linesearches and only display results for the SOs, even though L-BFGS with a Wolfeconditions linesearch is, thus far, the best combination for the spam dataset. It is worth noting,however, that by using RelStop, both L-BFGS,2d-Mm-BB-3 and L-BFGS,2d-Mm-BB-4 outperformthe linesearch.In fact, for both datasets, turning on both 1d-init and RelStop appear to be the best combination.The resulting solutions are generally more accurate and some are found with less effort (measuredboth in time and in number of iterations). For example, in the spam dataset, the most accuratesolution comes from running L-BFGS,2d-Mm-BB-4 although this takes significantly more time at161 seconds. The second most accurate solution, however, is L-BFGS,2d-Mm-BB-3 which requiresonly 19 seconds. An observation from this experiment is that although the benefit from usingan appropriate descent method far outweigh whether or not 1d-init or RelStop are used, tweakinginitialization and stopping conditions may improve performance.In the remaining experiments, we turn on both 1d-init and RelStop.4.2.5 Increasing the number of secondary directionsIn this experiment, we look at the effect of varying the number of momentum directions used in theSO from one to four. The idea for this comes from the supermemory gradient method (Cragg andLevy, 1969). The results of this experiment are given by figures B.9 and B.8 and tables B.9 andB.10. Increasing the number of directions do not appear to have any consistently beneficial effect.4.2.6 Quasi-Newton methods not convergingNewton’s method is known to possibly diverge or have a κ2 dependence instead of κ for GD (whereκ is the condition number) when starting from a solution that is sufficiently far from a local mini-mum.In this experiment, we consider the quasi-Newton method L-BFGS as s a proxy for Newton’smethod. Using an initial solution x0 equal to a vector of all ones instead of a vector all zeros willcause the method to not converge (see table 4.2). This is the case for both backtracking Armijoand Wolfe conditions. On the other hand, using a SO meant that the method gets closer to a localsolution. This backs up the idea that we could guarantee convergence with a SO for Newton’smethod even without restarts.39Method f (x∗) Time (sec) Outer Iters Inner Iters Mat-VecL-BFGS-0,Armijo 128.496 1.017 47 983 98L-BFGS-0,Wolfe 0.000000544 3.837 169 1,057 342L-BFGS-0,2d-Mm-BB 0.000000110 4.473 39 3,209 1,939L-BFGS-1,Armijo NaN 17.199 1,000 1,000 2,002L-BFGS-1,Wolfe NaN 38.907 1,000 50,955 2,002L-BFGS-1,2d-Mm-BB 0.0000150 8.496 155 11,537 5,417Table 4.2: Binary classification of the sido dataset with different initial solution x0. L-BFGS-0 and L-BFGS-1 correspond to x0 is a vector of zeros and ones respectively. The latterdiverges with a linesearch but remains convergent with a SO.40Chapter 5ConclusionThis thesis sets out to explore subspace optimization (SO) from a practical point of view. Our goalsinclude exploring the types of problems that are amenable to efficient SO, finding guidelines for thechoice of descent directions and of search methods, and analyzing the benefits of performing SOin theory and in practice. To this end, we identified a broader class of problems where SO may beefficient, provided experimental results that show cases where SO could arrive at better solutionsfaster, and started preliminary work on extending theory around linesearches to subspace searches.We hope that the thesis convinces the reader that not all SOs are created equal and that practicalnuances such as initialization matter. It is also our hope that the reader now considers SO a viablealternative to linesearches and an interesting topic for further research.We now turn to some unanswered questions from this thesis, and from SO literature in general,that may provide interesting topics for future work. We also briefly mention some loose ends totie up from this thesis. On the practical side, it could be useful to package minFuncLinear asa standalone piece of optimization software. This could allow those who want to use SO in theirresearch to do so without having to write extensive code to solve the subproblem. This entails ex-panding minFuncLinear so that it handles more problem classes and provides more optionalityin the way the SO is solved. At the moment, minFuncLinear is optimized for binary logisticregression. This should be extended to other problem classes such as multiclass logistic regres-sion, matrix completion and log determinant problems. Furthermore, the d-dimensional subspacesearch is currently solved indirectly as a d-dimension optimization subproblem. minFuncLinearshould provide other methods for solving SO such as directly satisfying d-dimensional Wolfe con-ditions. This would also require initialization and interpolation methods that are less crude. For ex-ample, minFuncLinear should extend HPI initialization and interpolation to higher dimensionswhich could then be used within an implementation of d-dimensional Wolfe conditions. Finally,minFuncLinear is written in Julia and a Python version of it may make the software accessibleto a wider group.41On the experimental side, a significant flaw of the thesis is that only two datasets were used. Thiscalls into question the generalizability of the experimental results. Running the same experimentson more datasets, for example on classification datasets from libSVM, would lend more weightto the results. Furthermore, there are other practical aspects of SO to explore experimentally. Wesuggest three broad experimental topics: characterizing the behaviour of SO when given “problemmatrices”, exploring the characteristics of typical SO subproblems and improving the convergencerate of L-BFGS using d-dimensional Wolfe conditions.In practice, the features matrix could be sparse or ill-conditioned. Does SO still work well underthese conditions? There are known methods that deal with sparse and/or ill-conditioned matrices.Could these methods be fruitfully combined with SO? In addition, the matrix of search directionsitself could cause problems. For example, SO may not work well if some of the directions in Pk arevery close to linearly dependent. How different do the directions need to be in order for SO to beeffective?The experimental results in this thesis suggest that SO performs better when the optimizationsubproblem is solved with methods such as BB instead of methods such as GD. One guess for whythis could be the case is that the typical optimization subproblem is well-characterized by a quadraticfunction. This allows methods like BB to converge faster to a good solution to the subproblem. Thebroader implication is that the typical subproblem itself could have particular characteristics thatcould be exploited by choosing suitable SO parameters. One observation from the experiments in thethesis is that when GD is combined with momentum as a second direction, convergence proceededin a stepwise fashion. It could be interesting to verify whether n-step convergence is occurring andwhether SO is choosing to restart periodically (i.e., periodically set the weight of the momentumdirection to zero).In theory, BFGS combined with a Wolfe linesearch gives superlinear convergence under someconditions. This property, however, does not generalize to the limited memory version. The reasonis related to how well the search directions approximate the Newton direction along the directionalderivative. While it is unlikely that SO can correct this issue for L-BFGS, a potentially interestingexperiment to conduct is to see how close we can get to BFGS performance by combining L-BFGSwith a suitable secondary direction and using d-dimensional Wolfe conditions.On the theory side, this thesis shallowly explored some topics that may warrant more detailedanalysis. One topic is to formally identify the subset of multilinear map problems that allow for ef-ficient SO. Could we determine whether or not multilinear maps are the correct level of abstractionfor efficient SO methods? If so, this could give rise to other problem classes that are not currentlyassociated with SO. Tensor problems could be one such class of problems. Intuitively, given vectorspaces W and Vi where i = {1, . . . ,k}, a multilinear map f : V1× ·· ·×Vk →W could also be ex-pressed as a linear map of a tensor product of vector spaces g : V1⊗·· ·⊗Vk→W . Thus, if SO aregenerally suited to multilinear map problems, they could also be effective for tensor problems.42While the focus of the thesis uses the full gradient, real world applications generally employstochastic methods. A research question is whether or not SO could work in a stochastic setting.The problem is that SO would have to work with noisy estimates of the gradient that may not evenbe descent directions. Perhaps a remedy for the stochastic setting would be to impose additional cri-teria, such as satisfying the over-parameterized regime explored in Meng et al. (2020) and Vaswaniet al. (2019). In terms of the algorithm itself, could we come up with a stochastic version of a SOmethod inspired by, for example, techniques used in stochastic dual subspace ascent (Loizou andRichta´rik, 2017)?The convergence analysis section of this thesis leaves much to be desired. A few directionscome to mind for more rigorous analysis. Firstly, the thesis does not show that GD with additionaldirections and weights set by SO converges strictly faster than GD with a simple Armijo linesearch.Experimentally, we see that the former is several orders of magnitude faster than the latter. Thedifficulty, however, is to show theoretically that this holds for all problems within a given problemclass and using any secondary direction. In this thesis, we considered objective functions that arequadratic objectives, strongly-convex or self-concordant. Perhaps there are other classes of prob-lems that lie between these classes with characteristics favorable to SO? Furthermore, if we restrictthe analysis to specific primary and secondary directions, are there complementary characteristics ofthe search directions that could be leveraged to guarantee SO progress? As a concrete example, theintroductory chapter discussed how CG on quadratic objectives is equivalent to GD with momentumand exact SO. Could we relax some aspects of the analysis, e.g., from an exact to an inexact SO, andstill retain some meaningful result? There are aspects of non-linear CG algorithms that resemblesSO. Could non-linear CG literature inform the use of SO in this particular case? Is it possible toshow a convergence rate that deviates from CG or accelerated rates by a known quantity?Secondly, convergence analysis could be guided by the field of numerical linear algebra. Con-sider the same example of CG above. Numerical linear algebra provides a class of algorithms ofwhich CG belongs and is used in more idealized problem settings. Chapter 4 of Greenbaum (1997)looks at the effects of inexact arithmetic on the convergence of linear CG when there is error in thesearch directions, i.e., the search directions are not perfectly conjugate. Perhaps the analysis of SOcan adopt similar error analysis, e.g., when SO return stepsizes with error.Thirdly, the methods that perform SO used in this thesis is wasteful in the sense that only thefinal result of the subproblem is used. For example, for linear map problems, minFuncLineartakes advantage of the low cost of calculating directional derivatives. Therefore, in a typical SOstep, the directional derivative is calculated for a large number of trial stepsizes (solutions). Thisprocess generates curvature information about the subproblem which is unexamined and unused.A reasonable approach is to attempt to use information gathered in the process of solving the SOsubproblem to improve on the choice of search directions in the main optimization problem, andthus potentially improve the overall convergence rate.43Finally, theoretical analysis of SO could be guided by using search directions that somehowcomplement each other. It is generally the case that no single search direction performs strictlybetter in all scenarios than any other search direction. The implicit hope of this thesis is to showthat SO can adapt to an optimization problem by picking the search direction most suited to the localcharacteristics of the objective function. Essentially, this means that SO compensates any limitationsof one search direction with other available search directions, and that worst-case scenarios areless prevalent because there is no obvious Achilles heel of the combined method. As an example,consider Newton’s method which has two known weaknesses. When the Hessian is not positivedefinite, Newton’s method may generate a search direction that is an ascent direction. Furthermore,when far from a local minimum, its convergence rate could be slow compared to simpler methodssuch as GD (because of a dependency on (µ/L)2). On the other hand, while GD does not achievesuperlinear convergence close to a local solution, it does return a descent direction at every iterationand is faster than Newton’s method far from the solution (because of its dependence on µ/L instead).Could SO intelligently combine Newton’s method with GD so that we get the best of both worlds?Could we produce combined methods that work better than each method on its own? If so, then thiscould be the main appeal of SO.44BibliographyMasoud Ahookhosh and Arnold Neumaler. An optimal subgradient algorithm with subspacesearch for costly convex optimization problems. Bulletin of the Iranian Mathematical Society,2018. [Online] doi: https://doi.org/10.1007/s41980-018-0171-1.Jonathan Barzilai and Jonathan M. Borwein. Two-point step size gradient methods. IMA Journalof Numerical Analysis, 8:141–148, January 1988.Dimitri P. Bertsekas. Constrained Optimization and Lagrange Multiplier Methods. AthenaScientific, 1982.Jeff Bezanson, Alan Edelman, Stefan Karpinski, and Viral B Shah. Julia: A fresh approach tonumerical computing. SIAM Review, 59(1):65–98, 2017. doi:10.1137/141000671. URLhttps://epubs.siam.org/doi/10.1137/141000671.Stephen Boyd and Lieven Vandenberghe. Convex Optimization. Cambridge, 2004.Se´bastien Bubeck, Yin Tat Lee, and Mohit Singh. A geometric alternative to nesterov’s acceleratedgradient descent. 2015. URL https://arxiv.org/pdf/1506.08187.pdf.Joel W. Cantrell. Relation between the memory gradient method and the fletcher-reeves method.Journal of Optimization Theory and Applications, 4(1):67–71, 1969.Peter Carbonetto. New Probabilistic Inference Algorithms that Harness the Strengths ofVariational and Monte Carlo Methods. PhD thesis, University of British Columbia., 2009.Arthur I. Cohen. Rate of convergence of several conjugate gradient algorithms. SIAM J. Numer.Anal., 9(2):248–259, June 1972.A.R. Conn, Nick Gould, A. Sartenaer, and Ph. L. Toint. On iterated-subspace minimizationmethods for nonlinear optimization. AMS-IMS-SIAM, pages 50–78, 1994.Gordon V. Cormack and Thomas R. Lynam. Spam Corpus Creation for TREC, 2005. URLhttps://plg.uwaterloo.ca/∼gvcormac/treccorpus/.E. E. Cragg and A. V. Levy. Study on a supermemory gradient method for the minimization offunctions. Journal of Optimization Theory and Applications, 4(3), 1969.Y.H. Dai and Y. Yuan. Nonlinear conjugate gradient methods. Shanghai Science and TechnologyPublisher, 2000.45Aaron Defazio. Momentum via primal averaging: Theoretical insights and learning rate schedulesfor non-convex optimization. 2020. URL https://arxiv.org/pdf/2010.00406.pdf.Y. Drori and A.B. Taylor. Efficient first-order methods for convex minimization: A constructiveapproach. Math. Program., 184:183–220, 2020.Euhanna Ghadimi, Hamid Reza Feyzmahdavian, and Mikael Johansson. Global convergence ofthe heavy-ball method for convex optimization. 2014. URL https://arxiv.org/pdf/1412.7457.pdf.Anne Greenbaum. Iterative Methods for Solving Linear Systems. Frontiers in AppliedMathematics. 1997.Isabelle Guyon. Sido: A Pharmacology Dataset, 2008. URLhttp://www.causality.inf.ethz.ch/data/SIDO.html.W.W. Hager and H. Zhang. A survey of nonlinear conjugate gradient methods. Pacific Journal ofOptimization, 2:35–58, 2006.Magnus R. Hestenes and Eduard Stiefel. Methods of conjugate gradients for solving linearsystems. Journal of Research of the National Bureau of Standards, 49(6):409–436, 1952.Sahar Karimi and Stephen Vavasis. A single potential governing convergence of conjugategradient, accelerated gradient and geometric descent. 2017. URLhttps://arxiv.org/abs/1712.09498.E. Stanley Lee. Quasilinearization and Invariant Imbedding: With Applications to ChemicalEngineering and Adaptive Control. Mathematics in Science and Engineering. Elsevier, 1968.doi:https://doi.org/10.1016/S0076-5392(08)62247-X.Nicholas Loizou and Peter Richta´rik. Momentum and stochastic momentum for stochasticgradient, newton, proximal point and subspace descent methods. 2017. URLhttps://arxiv.org/pdf/1712.09677.pdf.David G. Luenberger and Yinyu Ye. Linear and Nonlinear Programming, 4th Ed. Springer, 2008.S. Meng, S. Vaswani, I. Laradji, M. Schmidt, and S. Lacoste-Julien. Fast and furious convergence:Stochastic second order methods under interpolation. In The 23rd International Conference onArtificial Intelligence and Statistics, AISTATS, 2020.A. Miele and J.W. Cantrell. Study on a memory gradient method for the minimization of functions.Journal of Optimization Theory and Applications, 3(6), 1969a.A. Miele and J.W. Cantrell. Gradient methods in mathematical programming, part 2, memorygradient method. Technical Report 56, Rice University, Aero-Astonautics Report, 1969b.Kevin P. Murphy. Machine Learning: a Probabilistic Perspective. MIT Press, 2012.Guy Narkiss and Michael Zibulevsky. Sequential subspace optimization method for large-scaleunconstrained problems. Technical report, Technion - Israel Institute of Technology, 2005a.46Guy Narkiss and Michael Zibulevsky. Support vector machine via sequential subspaceoptimization. Technical report, Technion - Israel Institute of Technology, 2005b.J.A. Nelder and R.W.M. Wedderburn. Generalized linear models. Journal of the Royal StatisticalSociety, Series A, 135:370–384, 1972.A.S. Nemirovsky. The unit-vector method of smooth convex minimization. EngineeringCybernetics, Soviet journal of computer and systems science, 2, 1982. Orth-method for smoothconvex optimization (in russian).A.S. Nemirovsky and D.B. Yudin. Problem Complexity and Method Efficiency in Optimization.Wiley, 1983.Yurii Nesterov. A method for unconstrained convex minimization problem with the rate ofconvergence o(1/k2). In Doklady AN USSR, pages 543–547, 1983.Yurii Nesterov. Lectures on Convex Optimization, 2nd Ed. Springer Optimization and itsApplications. Springer, 2018.Yurii Nesterov and Boris Polyak. Cubic regularization of newton method and its globalperformance. Math. Program., 108:177–205, 2006.Yurii Nesterov, Alexander Gasnikov, Sergey Guminov, and Pavel Dvurechensky. Primal-dualaccelerated gradient methods with small-dimensional relaxation oracle. 2019. URLhttps://arxiv.org/pdf/1809.05895.pdf.Jorge Nocedal and Stephen J. Wright. Numerical Optimization, 2nd Ed. Springer, 2006.Julie Nutini, Issam Laradji, and Mark Schmidt. Let’s make block coordinate descent go fast:Faster greedy rules, message-passing, active-set complexity, and superlinear convergence. 2017.URL https://arxiv.org/abs/1712.08859.Barnaba´s Po´czos and Ryan Tibshirani. Convex optimization. Lecture notes for 10-725/ 36-725Fall 2013, 2013. URL https://www.stat.cmu.edu/∼ryantibs/convexopt-F13/.Marcos Raydan. The barzilai and borwein gradient method for the large scale unconstrainedminimization problem. SIAM J. Optim., 7(1):26–33, February 1997.Benjamin Recht. CS726 - Lyapunov analysis and the Heavy Ball Method, 2012. URLhttp://pages.cs.wisc.edu/∼brecht/cs726docs/HeavyBallLinear.pdf.Yousef Saad. Iterative Methods for Sparse Linear System, 2nd Ed. SIAM Publications, 2003.Mark Schmidt. minfunc: unconstrained differentiable multivariate optimization in matlab, 2005.URL https://www.cs.ubc.ca/∼schmidtm/Software/minFunc.html.Mark Schmidt. Practical Session on Convex Optimization, MLSS, 2011. URLhttps://www.cs.ubc.ca/∼schmidtm/MLSS.Jack Sherman and Winfred J. Morrison. Adjustment of an inverse matrix corresponding to achange in one element of a given matrix. Ann. Math. Statist., 21:124–127, 1950.47Gilbert Strang. Introduction to Linear Algebra, Fifth Edition. 2016.Gilbert Strang. Linear Algebra and Learning from Data. 2019.S. Vaswani, A. Mishkin, I. Laradji, M. Schmidt, G. Gidel, and S. Lacoste-Julien. Painlessstochastic gradient: Interpolation, line-search, and convergence rates. NeurIPS, 2019.M.A. Woodbury. Inverting modified matrices. Technical Report 42, Statistical Research Group,Princeton University, 1950.Y. Yuan and J. Stoer. A subspace study on conjugate gradient algorithms. Zamm-zeitschrift FurAngewandte Mathematik Und Mechanik, 75:69–77, 1995.Michael Zibulevsky. Sesop-tn: Combining sequential subspace optimization with truncatednewton method. 2008. URL https://ie.technion.ac.il/∼mcib/sesoptn paper1.pdf.48Appendix ABackground and supporting materialA.1 IntroductionA.1.1 Zoutendijk’s theoremThis section describes Zoutendijk’s theorem and the theorem’s conclusion, which is often referredto as the Zoutendijk conditions. The material comes from Nocedal and Wright (2006); Hager andZhang (2006).Suppose updates are in the form of (1.5) where pk is a descent direction and αk satisfies theWolfe conditions (or Goldstein conditions or strong Wolfe conditions), the objective f is continu-ously differentiable, bounded below and L-smooth, then∑k≥0cos2 θk‖∇ f (xk)‖2 < ∞ (A.1)where θk is the angle between pk and the negative gradient (1.4). (A.1) is often referred to as theZoutendijk condition. Intuitively, the theorem gives how far a descent direction pk can deviate fromthe negative gradient direction and still, with a well chosen steplength αk, result in an iteration thatis globally convergent. (A.1) is used to show global convergence for linesearch methods and also innon-linear CG convergence proofs. Because the sum is convergent, (A.1) also implies thatcos2 θk‖∇ f (xk)‖2→ 0.Note that global convergence means that the solution is a stationary point, or limk→0 ‖∇ f (xk)‖= 0,and not that the method finds a global minimum.49A.1.2 Linear CGThis section gives more details on the properties, convergence rates and many interpretations ofCG. We assume exact arithmetic, linesearches and SOs. The material comes from the textbooks ofLuenberger and Ye (2008); Nocedal and Wright (2006); Greenbaum (1997); Saad (2003); Nesterov(2018); Bertsekas (1982).CG as a SOOne interpretation of CG on a quadratic objective, sometimes referred to as the expanding manifoldproperty, states that the method generates iterates that are minimizers over an expanding linearKrylov subspace. Define residual rk , Axk +b. CG generates a sequence {xk} according to the rulexk = argmin{ f (x) : x ∈ x0+Lk},k ≥ 1whereLk = Lin{r0,Ar0, . . . ,Ak−1r0}(A.2)is the linear Krylov subspace of the initial residual.1CG is an instance of a conjugate direction method, meaning that the algorithm minimizes over asubspace of search directions Lk = Lin{p0, p1, . . . , pk−1}. In other words, search directions pk andresiduals rk are contained in the kth linear Krylov subspace for the initial residual. Furthermore,it can be shown that Lk is equivalently a linear subspace spanned by previous gradients, or Lk =Lin{∇ f (x0), . . . ,∇ f (xk−1)}.Properties of CGEigenvalue clusters For any x0 ∈ Rn, the sequence {xk} generated by the CG update step (1.12)converges to the solution x∗ of the linear system Ax = b in at most n steps. If A has only c distincteigenvalues then the CG iteration will terminate in at most c iterations. In general, if the eigenvaluesof A occur in c clusters, then CG will find the solution in roughly c steps.Conjugate directions CG iterates over a set of search directions that are mutually conjugate. Inother words, pᵀk Ap j = 0 for all j = 0, . . . ,k−1. A key reason for using conjugate directions is thatminimizing over such a subspace is equivalent to minimizing the objective along each directionsequentially. One way to view this is to first consider an A in the quadratic objective (1.11) that isdiagonal. This means that ∇2 f (x) = A is diagonal and minimizing each coordinate consecutively1Note that Ax∗ =−b where x∗ is the optimal solution. Thus, r0 = Ax0−b = A(x0− x∗) and Lk could also be writtenas the Krylov space spanned by the initial error Lk = Lin{A(x0− x∗), . . . ,Ak(x0− x∗)} for k ≥ 1.50will lead to the solution x∗. For a non-diagonal, positive definite A that is diagonalizable, we canredefine the quadratic objective asfˆ (xˆ), f (Sxˆ) = 12xˆᵀ(SᵀAS)xˆ− (Sᵀb)ᵀxˆ where xˆ = S−1xwhere S is a matrix that contains all the conjugate directions p0, . . . , pn−1 as columns. The conjugacyrequirement for the directions pk, k = 0, . . . ,n−1 ensure that SᵀAS is a diagonal matrix.CG generates a new direction pk that is a linear combination of the negative residual, or negativegradient, and the previous direction pk−1. Therefore, new direction has the formpk =−rk +βk pk−1 (A.3)where βk =rᵀk Apk−1pᵀk−1Apk−1.2 This means that the only storage requirement for CG is pk−1.Orthogonal residuals At each iteration, the residual rk is orthogonal to the residuals of previousiterations as well as to previous directions.3 In other words,rᵀk r j = 0 and rᵀk p j = 0 for all j = 0,1, . . . ,k−1Optimal polynomials Denote Pk(A) as a polynomial that takes as an argument the square matrixA. It can be shown that CG solves the quadratic objective (1.11) by generating the optimal degree kpolynomial P∗k at every iteration k restricted to the subspace (A.2). Polynomial P∗k is optimal in thatit minimizes the distance between the next iterate and the solution when distance is measured by A.In other words,xk+1 = x0+P∗k (A)r0whereP∗k = minPk‖x0+Pk(A)r0− x∗‖A.Informally, CG chooses the Pk that has roots at the k largest eigenvalues λn, . . . ,λn−k and at themidpoint between the smallest eigenvalue λ1 and λn−k.CG convergenceAssuming a quadratic objective given by (1.11), a comparison of the convergence rates of CG withGD shows that the former could be considered a generalization of the latter. Let A be a positive2We get this βk from left multiplying (A.3) by pᵀk−1A and imposing conjugacy on directions, i.e. set pᵀk−1Apk = 0.3For a quadratic objective, the residual is equal to the gradient, or rk = ∇ f (xk). This means that despite the name, CGis an orthogonal gradients method with conjugate search directions.51definite matrix and κ(A) be its Euclidean condition number given by (1.15).The update step of GD with an exact line search isxk+1 = xk−(∇ f (xk)ᵀ∇ f (xk)∇ f (xk)ᵀA∇ f (xk))∇ f (xk)and this results in the convergence rate‖xk+1− x∗‖A ≤(λn−λ1λn+λ1)‖xk− x∗‖A (A.4)where 0< λ1 ≤ ·· · ≤ λn are the eigenvalues of A. Recall the convergence rate of CG given by (1.13),which is rewritten below for convenience.‖xk+1− x∗‖A ≤(λn−k−λ1λn−k +λ1)‖xk− x∗‖A.The form of the two bounds are similar with a main difference. At every iteration, GD reduces thedistance to the solution by a fixed multiple that depends on the fixed ratio of the largest to the small-est eigenvalues. The convergence of GD therefore depends heavily on κ(A) (1.15). On the otherhand, the convergence rate of CG depends less on κ(A) as CG sequentially removes from consider-ation the largest eigenvalues. Conceptually, the problem becomes increasingly well-conditioned.Assumptions for global convergence A well-known property of linear CG on quadratic objectivesis that the method finds the optimal solution in at most n steps where n is the dimensionality. Thisproperty, however, only holds if we assume exact arithmetic. If exact conjugate directions are usedbut the search along that direction for the stepsize is inexact, by the conjugacy property, the errorintroduced will not be corrected in future iterations. Alternatively, if the directions are not perfectlyconjugate, the method may terminate in more than n steps. Finally, CG may not converge to theglobal minimum if the first step is not in the direction of the negative gradient.A.1.3 PHB with exact SO and CGFor quadratic objectives (1.11), PHB with exact SO give identical iterates to linear CG. A slightlydifferent proof is in Cantrell (1969) where the author shows that the memory gradient method isequivalent to FR for quadratic objectives.PHB with exact SO means thatxk+1 = xk−αk∇ f (xk)+βk pk−1 where pk−1 = xk− xk−1 and αk,βk = argminα,βf (xk+1)52Rewrite (1.11) in terms of the stepsizes and minimizeminφ(α,β ) = minα,β12‖xk−α∇ f (xk)+β pk−1‖2A+bᵀ(xk−α∇ f (xk)+β pk−1)+ cSet the partial derivatives to zeroφα = xᵀk A∇ f (xk)+β pᵀk−1A∇ f (xk)+∇ f (xk)ᵀb−α∇ f (xk)ᵀA∇ f (xk) = 0φβ = xᵀk Apk−1+β pᵀk−1Apk−1+ pᵀk−1b−α∇ f (xk)ᵀApk−1 = 0Use fact that ∇ f (x) = Ax+b and that A is symmetricα∇ f (xk)ᵀA∇ f (xk) = f (xk)ᵀ∇ f (xk)+β pᵀk−1A∇ f (xk)α∇ f (xk)ᵀApk−1 = pᵀk−1∇ f (xk)+β pᵀk−1Apk−1Because of the quadratic structure, f (Axk+1) = f (Axk−αA∇ f (xk)+βApk) the partial derivativesare just the directional derivatives −∇ f (xk)ᵀA∇ f (xk) and ∇ f (xk)ᵀApk−1. In an exact SO, these areequal to zero.4 This gives the stepsizesα =f (xk)ᵀ∇ f (xk)∇ f (xk)ᵀA∇ f (xk)and β =− pᵀk−1∇ f (xk)pᵀk−1Apk−1which gives the exact same update steps as CG. A.2 Problem structureA.2.1 Derivatives of affine functionGiven A ∈ Rm×n with columns [a11 . . .a1m]ᵀ , . . . , [an1 . . .anm]ᵀ and x,b ∈ Rn, then Ax+b is a m×1vector [∑nk=1 a1kxk +b1 . . . ∑nk=1 amkxk +bk]ᵀ. The ith component of Ax+b has a derivative that isa n× 1 vector ∂ ∑nk=1 aikxk+bk∂x = [ai1 . . .ain]ᵀ. This means that ∂Ax+b∂x = [ai j]m,ni, j=1 = A. Therefore, thegradient of a function f (x) = F(Ax+b) has the form ∇ f (x) = Aᵀ∇F(Ax+b).If ∇ f (x) is differentiable, the Hessian ∇2 f (x) could be viewed as its derivative. Applying thechain rule again gives ∇2 f (x) = Aᵀ∇2F(Ax+ b)A. In the case of logistic regression, this couldbe useful because ∇2F(Ax+ b) could be represented by a diagonal matrix D and ∇2 f (x) = AᵀDAwould be computationally efficient. For details, see Boyd and Vandenberghe (2004).4Intuitively, ∇ f (xk)ᵀApk−1 = 0 corresponds to one of the conjugacy properties of CG.53A.2.2 Gaussian modelsIn the case of undirected graphical model (UGM)s, A is an unknown model parameter that we es-timate. If, however, A 0 is known and suppose we are able to find its Cholesky decompositionA = LLᵀ where L is a lower triangular matrix with positive diagonal entries. Then finding the logdeterminant of A reduces to a trace operationlogdetA = logdet(LLᵀ) = log(detLdetLᵀ) = log((detL)2)= 2logdetL= 2log(∏iLii)= 2∑ilogLii = 2Tr (logL)where Tr(·) denotes the trace of a matrix. For details, see Murphy (2012).A.3 Descent directions and stepsize methodsThe material for this section comes from Nocedal and Wright (2006), Nesterov (2018) and Po´czosand Tibshirani (2013).A.3.1 Direction of steepest descentGiven search direction p, stepsize α and twice differentiable f , Taylor’s theorem means that forsome t ∈ (0,α),f (x+α p) = f (x)+α∇ f (x)ᵀp+12α2 p∇2 f (x+ t p)pᵀpThe rate of change of f along direction p at x is the directional derivative ∇ f (x)ᵀp. Solve for theunit direction of the most rapid decrease ormin‖p‖=1∇ f (x)ᵀp≡ min‖p‖=1‖p‖‖∇ f (x)‖cosθ ≡min‖∇ f (x)‖cosθwhere θ is the angle between p and ∇ f (x). This is minimal when cosθ =−1 and so p =− ∇ f (x)‖∇ f (x)‖and the negative gradient at point x is the direction of the steepest descent.A.3.2 Newton’s directionNewton’s direction minimizes the second order Taylor series approximation to f (xk +αk pk) wherethe step size αk is set to one. In other words,pk = argminpf (xk)+∇ f (xk)ᵀp+12∇2 f (xk)pᵀp54Assume that ∇2 f (xk) 0, we solve by setting the derivative with respect to p to zero and this givesthe Newton directionpk =−(∇2 f (xk))−1∇ f (xk)If ∇2 f (xk) 0, then the directional derivative is negative∇ f (xk)ᵀpk =−∇ f (xk)ᵀ∇2 f (xk)pk∇ f (xk) =−‖∇ f (xk)‖∇2 f (xk) < 0and Newton’s direction is a descent direction.A.3.3 Quasi-Newton directionsIntuitively, quasi-Newton methods, also known as variable metric methods, approximate ∇2 f (xk)with Bk ∈ Rn×n using the secant equationsBk(xk− xk−1)≈ ∇ f (xk)−∇ f (xk−1) (A.5)The intuition comes from quadratic functions 1.11 where the above holds exactly∇ f (y)−∇ f (x) = A(y− x) for any x,y ∈ RnA quasi-Newton method could be viewed as a method that lies between a gradient method whereBk is assumed to be I and Newton’s method where Bk is the true Hessian. In practice, quasi-Newtonalgorithms directly approximate the inverse Hessian Hk ≈(∇2 f (x))−1. These algorithms choose aHk that is symmetric, like the Hessian itself, and update Hk+1 with a low rank adjustment to Hk.Define ∆Hk = Hk+1−Hk, γk = ∇ f (xk+1)−∇ f (xk), δk = xk+1− xk. Then these update rulessatisfy the quasi-Newton relation1. Rank-one correction scheme∆Hk =(δk−Hkγk)(δk−Hkγk)ᵀ(δk−Hkγk)ᵀγk2. DFP:∆Hk =δkδ ᵀkγᵀk δk− Hkγkγᵀk Hk(Hkγk)ᵀγk3. BFGS:∆Hk = βkδkδ ᵀkγᵀk δk− Hkγkδᵀk +δkγᵀk Hkγᵀk δkwhere βk = 1+(Hkγk)ᵀγk/γᵀk δk.DFP is a conjugate direction method and gives CG directions while constructing Hk. BFGS is likeDFP except it updates Bk, and not Hk, and then uses Sherman-Morrison twice to get Hk.55A.3.4 BFGS convergenceSuperlinear convergence is guaranteed when some conditions are met.1. The method uses a quasi-Newton update;2. The stepsize αk satisfies either Wolfe or strong Wolfe;3. The approximation to the Hessian, Bk, in (A.5) is positive definite on every iteration;4. The linesearch at every iteration is initialized with a unit stepsize αk = 1;5. The constant used in the sufficient decrease condition c1 ≤ 12 .It can be shown that BFGS updates produce positive definite approximations Bk if the initial approx-imation B0 is positive definite and if (xk+1− xk)ᵀ(∇ f (xk−1)−∇ f (xk))> 0. In general, condition 3is guaranteed by most rank two type quasi-Newton methods. The intuition behind these conditionsis that the sequence of matrices that approximate the Hessian are ‘close enough’ to the Hessianalong the descent direction. This property is not satisfied for the limited memory version L-BFGS.In other words, L-BFGS iterates do not converge to Newton iterates as we get closer to the solution.For further detail, see Nocedal and Wright (2006).A.3.5 L-BFGSTo avoid storing a large (inverse) Hessian, L-BFGS stores instead a small number c n of vectorpairs δk = xk+1−xk and γk =∇ f (xk+1)−∇ f (xk). These pairs are referred to as corrections and ap-proximate the inverse Hessian Hk used in the BFGS update (3.2). In practice, the product Hk∇ f (xk)is calculated directly by starting with an initial approximation of Hk, H0k , and repeatedly applyingthe same BFGS update equation (3.2).The choice of H0k could vary between outer iterations and there are different ways to choose H0k .For example, we could define the initial approximation of Hessian B0k and solve for the inner loopestimate of Hk, r, in B0kr = ∇ f (xk) to directly get H0k∇ f (xk). This gives an implicit choice of H0k .Another choice of H0k is to set H0k =δᵀk−1γk−1γᵀk−1γk−1I. This choice works well in practice because it uses ascaling factor that tries to estimate the size of the true Hessian matrix along the most recent searchdirection. For more details, see sections 6.1 and 7.2 of Nocedal and Wright (2006).A.3.6 Memoryless L-BFGSA particular choice of L-BFGS is a memoryless version where the initial approximation to the Hes-sian H(0)k is set to the identity matrix and the number of corrections c = 1. This results in a Hk+1that could be considered (nearly) a symmetric positive definite version of the Hestenes-Stiefel (HS)form of non-linear CG (A.7). Combined with an exact linesearch, however, the memoryless L-BFGSupdate gives exactly linear CG. This reveals a relationship between quasi-Newton methods and CGmethods.As in Appendix A.3.3, define γk = ∇ f (xk+1)−∇ f (xk) and δk = xk+1− xk. The stepsize of the56momentum term used in HS non-linear CG gives us the descent directionpk+1 =−∇ f (xk+1)+ ∇ f (xk+1)ᵀγkpᵀk γkpkwhich we can rewrite aspk+1 =−Hˆk+1∇ f (xk+1) where Hˆk+1 ,(I− δkγᵀkγᵀk δk).We cannot use Hˆk+1 directly as an approximate Hessian because it is neither positive definite norsymmetric. Instead we use something similar to symmetric positive definite Hˆᵀk+1Hˆk+1 that alsosatisfies the secant equationHk+1 =(I− δkγᵀkγᵀk δk)(I− γkδᵀkγᵀk δk)+δkδ ᵀkγᵀk δk(A.6)which is exactly a BFGS update on I, or a memoryless L-BFGS. It can be shown that (A.6) with anexact linesearch gives linear CG updates. For more details, see Nocedal and Wright (2006).A.3.7 Non-linear CGThe search direction of non-linear CG has the form (3.3) which ensures that successive search di-rections are conjugate. Most non-linear CG methods differ only in their choice of the stepsize of themomentum term βk. Popular choices of βk are1. Dai-Yuan (1999)βk =‖∇ f (xk+1)‖2(∇ f (xk+1)−∇ f (xk))ᵀpk2. Fletcher-Reeves (1964)βk =−‖∇ f (xk+1)‖2‖∇ f (xk)‖23. Polak-Ribiere-Polyak (1969)βk =−∇ f (xk+1)ᵀ(∇ f (xk+1)−∇ f (xk))‖∇ f (xk)‖2Most variations of non-linear CG, including the three listed above, reduce to linear CG for quadraticobjectives. It is also possible to use the original HS momentum term stepsize introduced in theoriginal linear CG paper.5βk =∇ f (xk+1)ᵀ(∇ f (xk+1)−∇ f (xk))pᵀk (∇ f (xk+1)−∇ f (xk))(A.7)5In fact, the non-linear CG option in minFuncLinear implements this HS version.57Unlike linear CG on quadratics, non-linear CG on convex objectives generally does not terminatewithin n steps. Therefore, in most practical implementations, there is a restarting strategy that setsβk = 0. This is usually employed every n steps, ensures global convergence by taking a gradient stepand is sometimes also referred to as a spacer step. Although non-linear CG does not attain the fasterconvergence rate of Newton or quasi-Newton methods, the search directions could be computedeasily and also has very low storage requirements.For a comparison of non-linear CG methods, see Hager and Zhang (2006).A.3.8 Truncated NewtonTN, also referred to as linesearch Newton-CG, is an inexact Newton method that uses CG to obtaina search direction pk that is an inexact Newton step. TN avoids direct calculation of the Hessian byapproximately solving ∇2 f (xk)pk = −∇ f (xk). A common stopping criteria for the inner iterativesolver is based on the residual rk. The CG iterations terminate when‖rk‖ ≤ ηk‖∇ f (xk)‖ (A.8)where rk = ∇2 f (xk)pk +∇ f (xk). ηk with 0 < ηk < 1 for all k is referred to as the forcing sequence{ηk}. One possible choice of the forcing sequence is ηk = min(0.5,√‖∇ f (xk)‖). The choice ofthe forcing sequence is important in the convergence analysis of inexact Newton methods.In practice, TN uses a version of CG that accounts for systems that are not positive definite.When xk is not close to a solution, ∇2 f (xk) may have negative eigenvalues. TN terminates theCG iterations as soon as a direction of negative curvature is generated. When pk is a direction ofnegative curvature, it satisfies pᵀk∇2 f (xk)pk < 0 and is a descent direction.For details on local convergence of inexact Newton methods, see theorems 7.1 and 7.2 in No-cedal and Wright (2006).A.3.9 Hermite polynomial interpolationIn the linesearch case, let φ(α) = f (xk−1+α pk−1). We want to fit a quadratic approximation of φ ,φq, that interpolates three points, f (xk−1), f (xk) and ∇ f (xk−1)ᵀpk−1. Rewriting these three pointsin terms of φ givesf (xk−1) = φ(0) , f (xk) = φ(α) and ∇ f (xk−1)ᵀpk−1 = φ ′(0).58The interpolation conditions, φq(0) = φ(0), φ ′q(0) = φ ′(0) and φq(α) = φ(α), giveφq(αq) =(φ(α)−φ(0)−αφ ′(0)α2)α2q +φ′(0)αq+φ(0)=f (xk)− f (xk−1)−α∇ f (xk−1)ᵀpk−1α2α2q +∇ f (xk−1)ᵀpk−1αq+ f (xk−1)Setting the derivative with respect to αq to zero gives∇αqφq = 0 =⇒ 2(f (xk)− f (xk−1)−α∇ f (xk−1)ᵀpk−1α2)αq+∇ f (xk−1)ᵀpk−1 = 0=⇒ αq =− ∇ f (xk−1)ᵀpk−1α22( f (xk)− f (xk−1)−α∇ f (xk−1)ᵀpk−1)Set α0 = α = αq gives−α0∇ f (xk−1)ᵀpk−1 = 2(f (xk)− f (xk−1)−α0∇ f (xk−1)ᵀpk−1)and soα0 =2( f (xk)− f (xk−1))∇ f (xk−1)ᵀpk−1Note that both the numerator and the denominator are negative and so the stepsize is positive, andthat, in the special case where pk =−∇ f (xk), we get back the result in Schmidt (2011).A.4 ExperimentsA.4.1 Barzilai BorweinIn the paper Barzilai and Borwein (1988), the authors introduce a method that, like GD, uses thenegative gradient at the current iterate as the search direction and a specific stepsize. A search inthe BB direction is equivalent to one in the direction of the current gradient direction because theydiffer only by a scalar. Therefore, BB was not included in the list of descent directions in chapter 3.The BB update step isxk+1 = xk− 1αk∇ f (xk) where αk =(xk− xk−1)ᵀ(∇ f (xk)−∇ f (xk−1))‖xk− xk−1‖2 .Like many other situations, a quadratic objective provide an intuitive way to look at BB. The ideais similar to that of quasi-Newton methods: that the secant equations (A.5) should be satisfied atevery iteration. Recall that the quasi-Newton rule is to choose a symmetric positive definite matrix,59Hk+1 = Hᵀk+1 0, that approximates the inverse Hessian (∇2 f (xk+1))−1 such thatHk+1(∇ f (xk+1)−∇ f (xk)) = xk+1− xk.Satisfying the secant equations in the quadratic case is equivalent to finding the analytical solutionto the stepsize that minimizes the objective at the next step. Thus, in the quadratic case, BB achievessuperlinear convergence.For non-quadratics, there needs to be a stepsize strategy to achieve global convergence. InRaydan (1997), the authors demonstrated that the addition of a non-monotonic linesearch achievesthis. The non-monotonicity is needed because BB may increase the objective on some iterations.Instead of satisfying Armijo conditions (1.7)6, at every iteration, the non-monotonic linesearch hasa relaxed checkf (xk+1)≤ max0≤ j≤Mf (xk− j)+ρ∇ f (xk)ᵀ(xk+1− xk) where ρ ∈ (0,1).6Using Armijo would reduce BB to GD which is known to be slow. Using Armijo directly also does not make sensebecause the stepsize choice should not be related to the function value at the next step but to the eigenvalues of the Hessianat the minimizer.60Appendix BFigures and tables from experimentsFigure B.1: Binary classification of the spam dataset where * denotes the optimized version.Convergence is plotted in terms of both number of iterations and time (in seconds).61Method f (x∗) Time (sec) Outer Iters Inner Iters Mat-VecGD,Armijo 361.214 25.865 400 410 802GD,Wolfe 241.925 42.478 400 1,964 802GD,2d-Mm-BB 0.00150 35.838 145 4,776 5,071CG(nl),Armijo 313.809 37.622 400 591 802CG(nl),Wolfe 0.0603 66.098 400 2,832 802CG(nl),2d-Mm-BB 0.000869 39.225 168 2,172 4,251L-BFGS,Armijo 112.283 7.874 25 82 52L-BFGS,Wolfe 0.00000318 56.574 130 829 262L-BFGS,2d-Mm-BB 0.000738 42.766 88 1,261 1,995TN,Armijo 0.000144 36.154 399 84 800TN,Wolfe 0.000129 34.783 362 540 726TN,2d-Mm-BB 0.0000674 42.790 400 1,914 2,445Table B.1: Binary classification of the spam dataset. Comparing a linesearch to a SO.Method f (x∗) Time (sec) Outer Iters Inner Iters Mat-VecGD,Armijo 259.818 20.290 1000 1,025 2,002GD,Wolfe 269.327 24.286 1000 2,331 2,002GD,2d-Mm-BB 0.000141 59.293 933 45,439 22,877CG(nl),Armijo 32.722 20.815 1000 1,445 2,002CG(nl),Wolfe 0.000119 21.935 804 4,369 1,610CG(nl),2d-Mm-BB 0.0000150 14.044 161 18,406 6,669L-BFGS,Armijo 128.496 1.651 48 983 98L-BFGS,Wolfe 0.000000544 5.205 170 1,057 342L-BFGS,2d-Mm-BB 0.000000110 3.147 40 3,209 1,939TN,Armijo 0.00000258 9.831 476 61 954TN,Wolfe 0.00000293 8.345 393 266 788TN,2d-Mm-BB 0.0000507 22.292 1000 3,306 3,477Table B.2: Binary classification of the sido dataset. Comparing a linesearh to a SO.62Method f (x∗) Time (sec) Outer Iters Inner Iters Mat-VecGD,Wolfe 241.925017 45.851 400 1,964 802GD,2d-Mm-BB 0.0015 38.411 144 4,776 5,071GD,2d-TN-BB 847.895 38.343 400 2,717 1,446CG(nl),Wolfe 0.0603 55.750 400 2,832 802CG(nl),2d-Mm-BB 0.000869 31.979 167 2,172 4,251CG(nl),2d-TN-BB 689.397 38.252 400 3,058 1,475L-BFGS,Wolfe 0.00000318 49.424 129 829 262L-BFGS,2d-Mm-BB 0.000738 37.858 87 1,261 1,995L-BFGS,2d-TN-BB 0.000000137 19.069 99 3,336 1,368TN,Wolfe 0.000129 32.729 361 540 726TN,2d-Mm-BB 0.0000674 44.198 400 1,914 2,441TN,2d-TN-BB 35.200 71.519 400 3,526 2,162Table B.3: Binary classification of the spam dataset. Comparing momentum and TN as a sec-ond direction.Method f (x∗) Time (sec) Outer Iters Inner Iters Mat-VecGD,Wolfe 269.326556 19.441 1,000 2,331 2,002GD,2d-Mm-BB 0.000141 48.378 932 45,439 22,877GD,2d-TN-BB 0.136457 28.876 1,000 19,099 9,276CG(nl),Wolfe 0.000119 19.466 803 4,369 1,610CG(nl),2d-Mm-BB 0.0000146 12.058 160 18,406 6,669CG(nl),2d-TN-BB 0.00000769 19.006 1,000 6,195 4,562L-BFGS,Wolfe 0.000000544 4.688 169 1,057 342L-BFGS,2d-Mm-BB 0.00000011 2.699 39 3,209 1,939L-BFGS,2d-TN-BB 0.0000165 17.520 1,000 2,373 3,828TN,Wolfe 0.00000293 6.457 392 266 788TN,2d-Mm-BB 0.0000507 17.574 1,000 3,306 3,473TN,2d-TN-BB 0.000000001 10.753 534 4,235 2,484Table B.4: Binary classification of the sido dataset. Comparing momentum and TN as a seconddirection.63Method f (x∗) Time (sec) Outer Iters Inner Iters Mat-VecGD,Wolfe 241.925 43.021 400 1,964 802GD,2d-Mm-BB 0.00150 35.283 144 4,776 5,071GD,2d-Mm-GD 0.106 131.652 400 21,212 43,127GD,2d-Mm-TN 0.000557 10.675 51 315 2,021CG(nl),Wolfe 0.0603 55.485 400 2,832 802CG(nl),2d-Mm-BB 0.000869 32.639 167 2,172 4,251CG(nl),2d-Mm-GD 0.139 127.540 400 19,776 41,873CG(nl),2d-Mm-TN 0.000436 28.564 172 693 4,737L-BFGS,Wolfe 0.00000318 46.428 129 829 262L-BFGS,2d-Mm-BB 0.000738 35.813 87 1,261 1,995L-BFGS,2d-Mm-GD 0.148 220.248 400 20,037 42,851L-BFGS,2d-Mm-TN 0.000410 33.485 87 257 1,961TN,Wolfe 0.000129 31.909 361 540 726TN,2d-Mm-BB 0.0000674 38.888 400 1,914 2,441TN,2d-Mm-GD 0.543 53.628 400 4,444 9,327TN,2d-Mm-TN 84.833 36.424 400 375 1,365Table B.5: Binary classification of the spam dataset. Comparing different descent methodsused in the SO.Method f (x∗) Time (sec) Outer Iters Inner Iters Mat-VecGD,Wolfe 269.327 19.864 1,000 2,331 2,002GD,2d-Mm-BB 0.000141 46.584 932 45,439 22,877GD,2d-Mm-GD 0.000256 61.981 1,000 60,676 105,307GD,2d-Mm-TN 0.00000403 2.285 49 385 2,361CG(nl),Wolfe 0.000119 18.260 803 4,369 1,610CG(nl),2d-Mm-BB 0.0000146 11.577 160 18,406 6,669CG(nl),2d-Mm-GD 0.000212 52.311 801 38,876 82,369CG(nl),2d-Mm-TN 0.0000194 3.921 104 527 3,387L-BFGS,Wolfe 0.000000544 4.570 169 1,057 342L-BFGS,2d-Mm-BB 0.000000110 2.334 39 3,209 1,939L-BFGS,2d-Mm-GD 0.0000543 35.939 643 30,043 63,791L-BFGS,2d-Mm-TN 0.0000235 1.508 41 226 1,437TN,Wolfe 0.00000293 6.558 392 266 788TN,2d-Mm-BB 0.0000507 17.127 1,000 3,306 3,473TN,2d-Mm-GD 0.00261 19.174 1,000 4,828 9,597TN,2d-Mm-TN 0.0000109 16.904 1,000 276 3,715Table B.6: Binary classification of the sido dataset. Comparing different descent methods usedin the SO.64Method f (x∗) Time (sec) Outer Iters Inner Iters Mat-VecGD,2d-Mm-BB-1 0.00150 42.745 144 4,776 5,071GD,2d-Mm-BB-2 0.00371 87.006 400 5,390 8,277GD,2d-Mm-BB-3 0.00202 39.272 138 4,796 5,105GD,2d-Mm-BB-4 0.00183 95.026 400 7,598 10,173CG(nl),2d-Mm-BB-1 0.000869 36.573 167 2,172 4,251CG(nl),2d-Mm-BB-2 0.000265 80.561 400 2,152 6,051CG(nl),2d-Mm-BB-3 0.00148 40.321 179 2,472 4,655CG(nl),2d-Mm-BB-4 0.00285 82.579 400 2,898 7,047L-BFGS,2d-Mm-BB-1 0.000738 42.019 87 1,261 1,995L-BFGS,2d-Mm-BB-2 0.0000288 191.088 400 1,504 5,031L-BFGS,2d-Mm-BB-3 0.000000322 18.996 38 1,122 1,465L-BFGS,2d-Mm-BB-4 0.000000076 160.950 326 2,108 4,363TN,2d-Mm-BB-1 0.0000674 47.782 400 1,914 2,441TN,2d-Mm-BB-2 0.0000527 51.688 399 1,764 2,653TN,2d-Mm-BB-3 0.0000115 21.492 144 1,749 1,675TN,2d-Mm-BB-4 0.0000109 53.644 399 1,763 3,071Table B.7: Binary classification of the spam dataset. Comparing different initialization meth-ods and exit conditions.Method f (x∗) Time (sec) Outer Iters Inner Iters Mat-VecGD,2d-Mm-BB-1 0.000141 51.642 932 45,439 22,877GD,2d-Mm-BB-2 0.274 54.695 1000 50,325 23,545GD,2d-Mm-BB-3 0.000159 63.644 1000 88,347 36,297GD,2d-Mm-BB-4 0.000386 58.382 1000 69,624 33,911CG(nl),2d-Mm-BB-1 0.0000146 12.557 160 18,406 6,669CG(nl),2d-Mm-BB-2 0.0000652 40.816 1000 5,396 12,517CG(nl),2d-Mm-BB-3 0.00000860 15.872 213 19,860 8,371CG(nl),2d-Mm-BB-4 0.00000375 41.410 1000 4,378 12,417L-BFGS,2d-Mm-BB-1 0.000000110 3.021 39 3,209 1,939L-BFGS,2d-Mm-BB-2 0.000000033 2.442 32 2,432 1,515L-BFGS,2d-Mm-BB-3 0.000000013 3.182 51 2,512 2,001L-BFGS,2d-Mm-BB-4 0.00000002 2.222 34 1,775 1,307TN,2d-Mm-BB-1 0.0000507 21.623 1,000 3,306 3,473TN,2d-Mm-BB-2 0.00000184 25.403 1,000 3,106 4,655TN,2d-Mm-BB-3 0.000000505 8.182 259 5,810 2,587TN,2d-Mm-BB-4 0.000000032 7.701 273 3,493 2,177Table B.8: Binary classification of the sido dataset. Comparing different initialization methodsand exit conditions.65Method f (x∗) Time (sec) Outer Iters Inner Iters Mat-VecGD,2d-Mm-BB-1 0.000141 49.862 932 45,439 22,877GD,2d-Mm-BB-2 0.305 55.012 1,000 58,599 25,085GD,2d-Mm-BB-3 0.0307 57.769 1,000 58,962 27,312GD,2d-Mm-BB-4 0.654 62.793 1,000 74,167 30,842CG(nl),2d-Mm-BB-1 0.0000146 12.514 160 18,406 6,669CG(nl),2d-Mm-BB-2 0.0000817 8.989 173 8,330 5,457CG(nl),2d-Mm-BB-3 0.000134 14.991 279 14,453 8,266CG(nl),2d-Mm-BB-4 0.0000718 19.179 332 22,321 11,255L-BFGS,2d-Mm-BB-1 0.000000110 2.728 39 3,209 1,939L-BFGS,2d-Mm-BB-2 0.00000142 3.263 54 3,071 2,138L-BFGS,2d-Mm-BB-3 0.0000135 4.213 74 3,391 2,584L-BFGS,2d-Mm-BB-4 0.000000014 2.425 32 2,688 2,019TN,2d-Mm-BB-1 0.0000507 20.121 1,000 3,306 3,473TN,2d-Mm-BB-2 0.00001832 18.565 1,000 4,212 4,141TN,2d-Mm-BB-3 0.000000003 8.038 294 7,456 3,470TN,2d-Mm-BB-4 0.475 22.941 1,000 12,545 7,082Table B.9: Binary classification of the sido dataset. The effect of going from one to fourmomentum terms.Method f (x∗) Time (sec) Outer Iters Inner Iters Mat-VecGD,2d-Mm-BB-1 0.00150 35.018 144 4,776 5,071GD,2d-Mm-BB-2 0.00233 26.916 114 3,823 4,254GD,2d-Mm-BB-3 0.00150 21.973 86 3,383 3,648GD,2d-Mm-BB-4 0.00107 24.458 107 3,204 3,876CG(nl),2d-Mm-BB-1 0.000869 30.129 167 2,172 4,251CG(nl),2d-Mm-BB-2 0.00108 33.319 185 2,215 4,449CG(nl),2d-Mm-BB-3 0.00142 31.290 169 2,273 4,286CG(nl),2d-Mm-BB-4 0.00258 41.432 219 3,127 5,764L-BFGS,2d-Mm-BB-1 0.000000110 2.728 39 3,209 1,939L-BFGS,2d-Mm-BB-2 0.00000142 3.263 54 3,071 2,138L-BFGS,2d-Mm-BB-3 0.0000135 4.213 74 3,391 2,584L-BFGS,2d-Mm-BB-4 0.000000014 2.425 32 2,688 2,019L-BFGS,2d-Mm-BB-1 0.000738 33.822 87 1,261 1,995L-BFGS,2d-Mm-BB-2 0.00000907 17.221 40 1,485 1,662L-BFGS,2d-Mm-BB-3 0.000191 23.003 56 1,383 1,764L-BFGS,2d-Mm-BB-4 0.0000298 16.884 36 1,863 1,937Table B.10: Binary classification of the sido dataset. The effect of going from one to fourmomentum terms.66Figure B.2: Binary classification of the spam dataset. Comparing momentum and TN as asecond direction.67Figure B.3: Binary classification of the sido dataset. Comparing momentum and TN as a sec-ond direction.68Figure B.4: Binary classification of the spam dataset. Comparing SO methods.69Figure B.5: Binary classification of the sido dataset. Comparing SO methods.70Figure B.6: Binary classification of the spam dataset. Linesearch initialization and relativeexit conditions for the SO.71Figure B.7: Binary classification of the sido dataset. Linesearch initialization and relative exitconditions for the SO.72Figure B.8: Binary classification of the spam dataset. Varying the number of momentum di-rections in the SO.73Figure B.9: Binary classification of the sido dataset. Varying the number of momentum direc-tions in the SO.74
Thesis/Dissertation
eng
Computer Science
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Schmidt, Mark
Subspace optimization for machine learning
Text