A simple Newton method for local nonsmooth optimization
Adrian Lewis Calvin Wylie
math.OC cs.NA math.NA 90C25, 65K05, 49M15







 
















thmTheorem[section]
cor[thm]Corollary
prop[thm]Proposition
lem[thm]Lemma
defn[thm]Definition
rem[thm]Remark
exa[thm]Example
ass[thm]Assumption
alg[thm]Algorithm
[2]
Proof  


 
R










R^n

x̅
y̅

 
 
 
 




 1pt →-12pt
           2.8pt → 
myequation
   
    thm1myeqnarray
        ||‖‖{}
        
        A simple Newton method for local nonsmooth optimization
        
            
        A.S. Lewis
        ORIE, Cornell University, Ithaca, NY 14853, USA.
        
        Research supported in part by National Science Foundation Grant DMS-1613996.
        C.J.S. Wylie
        ORIE, Cornell University, Ithaca, NY 14853, USA.
            December 30, 2023
        ============================================================================================================================================================================================================
        
        
        Superlinear convergence has been an elusive goal for black-box nonsmooth optimization.  Even in the convex case, the subgradient method is very slow, and while some cutting plane algorithms, including traditional bundle methods, are popular in practice, local convergence is still sluggish.  Faster variants depend either on problem structure or on analyses that elide sequences of “null” steps.  Motivated by a semi-structured approach to optimization and the sequential quadratic programming philosophy, we describe a new bundle Newton method that incorporates second-order objective information with the usual linear approximation oracle.  One representative problem class consists of maxima of several smooth functions, individually inaccessible to the oracle.  Given as additional input just the cardinality of the optimal active set, we prove local quadratic convergence.  A simple implementation shows promise on more general functions, both convex and nonconvex, and suggests first-order analogues.
         Key words: nonsmooth optimization, bundle method, Newton method, 
        
        sequential quadratic programming, nonconvex
         AMS 2010 Subject Classification: 90C25, 65K05, 49M15
        
        § INTRODUCTION
        
        Accurate deterministic minimization of a smooth objective function f → constitutes the foundation of the classical optimization literature <cit.>.  Typically we assume a black-box oracle, returning, at any point x ∈, the objective value f(x), gradient ∇ f(x), and perhaps the Hessian ∇^2 f(x).  Fast local convergence — superlinear or quadratic — relies ultimately on Newton's method.
        
        By contrast, nonsmooth optimization, even for convex objectives, is more challenging, at least in the analogous black-box model <cit.>.  Armed with the right tools for the circumstances — some manner of global structural knowledge — reasonably accurate nonsmooth convex minimization may be tractable.  The contemporary optimization literature highlights various examples: the classical proximal point method and its various splitting-type extensions, when we can suitably decompose the objective into smooth and prox-friendly ingredients <cit.>; basic techniques for smooth minimization (along the lines surveyed in <cit.>) extended to nonsmooth objectives amenable to suitable smoothing techniques <cit.>; polynomial-time interior-point techniques for semidefinite-representable minimization <cit.>.
        
        Outside such structured convex optimization realms, however, we fall back on black-box techniques, using function values and subgradients but blind to any global structure beyond perhaps convexity <cit.>.  The simplicity and all-purpose appeal of such algorithms is tempered by their local convergence in practice.  At large scale, the classical subgradient algorithm is the notoriously slow method of last resort <cit.>.  In the convex case, modern cutting plane methods based on analytic centers <cit.> or the volumetric barrier <cit.>, can be highly effective at moderate scale, but convergence is linear.  
        
        Bundle methods, which originate with <cit.>, have proved quite successful in large-scale practice, are supported by convergence theory <cit.>, and  extend to the nonconvex case (although implementation is delicate — see <cit.>).  However, classical bundle methods converge sublinearly, and while more recent black-box versions such as <cit.> are faster, they too are complex to implement, and the convergence analysis considers only “serious” steps rather than oracle calls.  Sagastizábal's 2018 ICM lecture <cit.> includes an elegant and comprehensive survey.  Impractically large subproblems disadvantage some otherwise promising analogous developments, like level bundle methods <cit.>. 
        
        The quest for superlinear acceleration described in <cit.> has for two decades also driven important “𝒱𝒰” and “partly smooth” theory, pioneered by Lemaréchal, Mifflin, Oustry, and Sagastizábal <cit.> and paralleled in <cit.>.  However, still missing in black-box nonsmooth optimization, nonconvex or convex, is the fast local convergence of Newton's method.  Our aim here is just such a black-box Newton-type local optimization method, supported by a superlinear convergence theory on a range of interesting nonsmooth functions, and promising in practice.  
        
        As a step in that direction, we take an unconventional  semi-structured approach.  As motivation, we first consider how to minimize a common type of structured nonsmooth objective function, specifically a pointwise maximum of finitely-many smooth component functions, but using a black-box oracle for the objective that  cannot access the component functions individually.  The resulting algorithm becomes a model-based black-box “bundle Newton” method (terminology coming from <cit.>) that converges locally quadratically on nonsmooth functions of this max-type, but that also shows promise in practice on more general nonsmooth functions, both convex and nonconvex.  
        
        Like any local method, including Newton's method itself, globalizing the algorithm is an important challenge, and a topic for future work, as is the first-order analogue we describe briefly at the conclusion of this work.  Nonetheless, this Newton-type local approach seems a promising new ingredient for nonsmooth optimization.
        
        
        
        § AN ALGORITHM FOR NONSMOOTH MINIMIZATION
        
        We seek a local method for finding a minimizer ∈ of a continuous nonsmooth  objective function f.
        We assume, around every point in some set, 𝒟⊂, that the objective f is twice continuously differentiable.   Given any point s ∈𝒟, an oracle returns the value f(s), gradient ∇ f(s), and Hessian ∇^2 f(s), allowing us to build the corresponding linear and quadratic approximations, l_s and q_s, to f.  Beyond that, we call on no further information about the function f.
        
        Since the objective f is nonsmooth, we cannot typically hope to find a point s ∈𝒟 where the gradient ∇ f(s) is small.  Instead we seek a finite set (or  bundle) S ⊂𝒟 with small diameter
        
            S  = max{|s-s'| :s,s' ∈ S }
        
        (where |·| denotes the Euclidean norm),
        and small  optimality measure
            Θ(S)  = min|(∇ f(S))|.
        
        We call the elements of the bundle  reference points.  In terms of the simplex
        
            Δ_S  = {λ∈^S_+ : ∑_s ∈ Sλ_s = 1 },
        
        the optimality measure is
        Θ(S)  = min_λ∈Δ_S| ∑_s ∈ Sλ_s ∇ f(s) |.
        
        If the set of gradients ∇ f(S) is affinely independent, then the minimization problem (<ref>) has a unique optimal solution λ.  This vector, which plays an important role in the algorithm, we call the  Lagrange multiplier estimate.  
        
        We describe an algorithm that iteratively updates a bundle S one reference point at a time, stopping if it encounters a point outside the set 𝒟.  Thus S has a fixed size (cardinality) throughout:  a judicious choice of this integer parameter, denoted k, is a crucial feature of the algorithm that distinguishes it from more standard cutting plane and bundle methods.  Too small a value of k causes the method to fail because the optimality measure Θ(S) remains bounded away from zero;  too large a value may lead to an ill-defined Lagrange multiplier estimate.
        
        
        
          §.§.§ The role of convexity
        
        The core algorithm below is aimed at strongly convex objectives f, even though its statement makes sense without convexity and our interest in it is purely local.  Indeed, much of the analysis that follows is independent of convexity, and a straightforward modification of the core method results in our culminating nonconvex algorithm.
        
        Nonetheless, simple examples show that the unmodified method can fail for nonconvex objectives.  We present this algorithm first because it is easier to motivate, our aim in Section <ref>.  Most importantly, the final step in our key quadratic convergence proof (Theorem <ref>) depends crucially on a convexity argument.  Without more ado, here is an informal description of the core algorithm.
        
        [k-bundle Newton method for strongly convex f]
        
        
        Notice that the case k=1 is just the classical Newton method.
        
        
        
        § MOTIVATING THE METHOD
        
        
         §.§ The optimality measure
        
        Basic aspects of the algorithm coincide with traditional cutting plane methods.  The method aims to construct a sequence of bundles S converging to the minimizer  (and hence in particular with diameter converging to zero) in such a way that the optimality measure δ = Θ(S)
        also converges to zero:  in that case we say that the algorithm  succeeds.  In practice we deduce approximate optimality when both these measures are small.
        
        The algorithm computes an optimal vector λ (the Lagrange multiplier estimate) for the optimality measure expression (<ref>).  If we define the weighted average of the reference points
        s̅ = ∑_s ∈ Sλ_s s,
         
        then f(s̅) is certainly an upper bound on the optimal value min f.  To obtain an approximate lower bound, we use convexity to note that each linear approximation
        
            l_s(x) = f(s) + ∇ f(s)^T(x-s)      (x ∈),
        
        minorizes f, and hence so does their weighted average ∑_s λ_s l_s.  Denoting by L a Lipschitz constant for f on some convex set containing all the reference points, we have 
        |∇ f(s)| ≤ L for all s ∈ S.  Then for all points x ∈ we have
        
            f(s̅) - L (S)  ≤  f(s̅) - ∑_s ∈ Sλ_s L|s-s̅| 
             ≤  f(s̅) + ∑_s ∈ Sλ_s ∇ f(s)^T(s̅ - s)
            
                   = ∑_s ∈ Sλ_s l_s(s̅)  ≤ ∑_s ∈ Sλ_s l_s(x) + δ|x-s̅|
             ≤  f(x) + δ|x-s̅|,
        
        so in conclusion we have 
        min f  ≤  f(s̅)  ≤ min{f+δ |· - s̅|} +  L (S).
        
        Thus if both the diameter of the bundle S and the optimality measure δ are small, then the objective value at the point s̅ lies in the small interval between the minimum values of the objective function and a slightly perturbed function.  In this sense, the current bundle constitutes an approximate certificate of optimality.  
        
        As a consequence of this argument, success of the algorithm requires a certain lower bound on the bundle size k, as we discuss next.
        
        
        
         §.§ A lower bound on bundle size:  Carathéodory number
        
        Given any set Γ⊂ containing zero in its convex hull, define the  Carathéodory numberΓ, to be the minimum size of a subset whose convex hull contains zero.  By Carathéodory's theorem, we see
        
            1  ≤  Γ ≤  1 + ( Γ).
        
        Suppose that the set 𝒟 has full measure.  When 𝒟 is simply the set of all points where f is twice continuously differentiable, this assumption typically holds in practice, and in particular if the objective f is semi-algebraic <cit.>.   Define the limiting gradient set
        Γ = {lim_r ∇ f(x_r) : lim_r x_r = ,  x_r ∈𝒟  r=1,2,…}.
        
        Since f is locally Lipschitz, being continuous and convex, we have (see <cit.>)
        
        0 ∈∂ f() =  Γ.
        
        For the algorithm to succeed, the optimality measure Θ(S) must converge to zero for some sequence of bundles S converging to the minimizer , so zero is a convex combination of k elements of the set Γ, and hence the lower bound
        
            k ≥ Γ
        
        must hold.
        
        
        
         §.§ The active subspace
        
        Having computed the optimality measure δ = Θ(S), the method next seeks a new reference point. If the current bundle S is close to the minimizer , then the cutting plane model 
        f̃→ defined by
        
            f̃(x)  = max_s ∈ S l_s(x),
        
        minorizes the objective f, and approximates it around .  Furthermore, at every point x on the  active subspace where all the linear approximations are equal,
        
            M  = {x ∈ : l_s(x)   s ∈ S},
        
        the cutting plane model has subdifferential
        
            ∂f̃(x) = ( ∇ f(S) )
        
        and hence nonsmooth slope (the fastest rate of decrease) equal to δ.  Thus when the optimality measure is small, the cutting plane model is approximately minimized throughout the active subset, so there is where the method seeks  a new reference point.
        
        
        
         §.§ An upper bound on bundle size: affine independence
        
        For algorithmic stability, the minimizer  should be close to the active subspace M.  Since the bundle S is close to , the values of the linear approximations l_s() are all close to 
        f().  Equivalently, therefore, we ask that the point 
        (,f()) ∈× should be close to the affine subspace
        
            { (x,t) ∈× : l_s(x) = t    s ∈ S }.
        
        As we have observed, the residual in the linear system defining this subspace is small at the point (,f()). Standard linear algebra shows that this point is therefore close the solution set (which in particular is nonempty) providing that the smallest singular value of the matrix for the system is bounded below by some fixed tolerance 
        σ > 0.  That (n+1)-by-k matrix has columns ∇ f(s)1 (for s ∈ S), so this lower bound amounts to uniform affine independence of the gradients. 
        
        Assuming that this condition holds, success of the algorithm then requires an upper bound on the bundle size k, since by taking a convergent subsequence of the matrices above we arrive at a limiting matrix with k linearly independent columns of the form g1, where each vector g lies in the limiting gradient set Γ in (<ref>), and hence in the subdifferential ∂ f().  Thus the function f has at least k affine-independent subgradients at , from which we deduce the upper bound
        
            k ≤ 1+(∂ f()).
        
        
        As we shall see, the uniform affine independence  property holds automatically in our convergence proof for functions of max-type.  For more general versions of the algorithm, however, we verify the property as follows.  For any vectors 
        g_i ∈ indexed by i in a list I of length k, we denote by
        
            σ_I{g_i : i ∈ I}
        
        the kth largest singular value of a matrix with k columns 
        g_i1 (for i ∈ I).  This nonnegative number is zero exactly when the list is affine dependent.  Using this notation, we fix a parameter σ > 0 at the outset, and add the following check at the beginning of each iteration.
        
        
         §.§ Choosing the bundle size
        
        To summarize, if the algorithm succeeds, the bundle size k must satisfy both upper and lower bounds involving the Carathéodory number of the limiting gradient set and the dimension of the subdifferential:
        
            Γ ≤  k  ≤  1+(∂ f()).
        
        In general, these lower and upper bounds may be far apart.  For example, for the Euclidean norm 
        f = |·| at the point = 0, with the set 𝒟 = ∖{0}, the bounds become
        
            2 ≤ k ≤ n+1.
        
        However, in the following case of particular interest to us in this work, the two bounds are equal.
        
        [Convex max functions] 
        Consider a nonsmooth function of the form
        
            f(x)  = max_i=1,…,k f_i(x)    (x ∈),
        
        for smooth convex functions f_i ^n →, for i=1,2,…,k.  
        At the point ∈, suppose that the function values f_i() are all equal, with gradients 
        g_i = ∇ f_i(), so we have
        
            ∂ f()  = { g_i : i=1,2,…,k}.
        
        Assuming that the list of gradients {g_i} is affinely independent, we have 
        
            k  =  1+(∂ f()).
        
        Furthermore, assuming that  is a minimizer, so 0 ∈∂ f(), the system
        
            ∑_i λ_i g_i = 0,   ∑_i λ_i = 1,   λ∈^k_+
        
        must then have a unique solution λ̂∈^k.  Choose 𝒟 to be the set of points x ∈ for which the maximizing index set _i f_i(x) is a singleton, so the limiting gradient set (<ref>) is
        
            Γ = {g_i : i=1,2,…,k}.
        
        Hence the Carathéodory number Γ is the number of nonzero components of the vector 
        λ̂, which is k exactly when  is in fact a  nondegenerate minimizer, meaning that zero lies in the relative interior of the subdifferential ∂ f().  
        
        
        In general, estimating the lower bound on the bundle size, the Carathéodory number Γ, seems challenging in practice.  On the other hand, with respect to the upper bound, global nonsmooth optimization methods — the various methods we discussed in the introduction, including the subgradient method, the proximal point and proximal gradient methods, other splitting methods, bundle and level bundle methods, along with nonsmooth BFGS <cit.>, and gradient sampling <cit.>, for example —  typically suggest subdifferential dimension information as they progress.  Given any finite set of points Ω⊂𝒟 near the minimizer x̅, we can use equation (<ref>) to estimate
        
            ∂ f (x̅)  ≈ (∇ f(Ω)).
        
        The dimension of the set on the right-hand side is the rank of a matrix with columns
        ∇ f(x)1 (for x ∈Ω).
        This suggests that a reasonable estimate of the dimension of ∂ f() is the approximate rank — the number of singular values larger than some tolerance — of this same matrix.  This approximate rank then might serve as the integer k in our Newton method.
        
        
        
         §.§ The quadratic subproblem
        
        At the end of each iteration we update the bundle S by substituting a new reference point x̂∈ for that point in S whose deletion minimizes the resulting optimality measure.  The Newtonian flavor of the algorithm arises from the choice of x̂∈, which solves a simple, linearly-constrained quadratic program that we discuss next.
        
        To succeed, the algorithm must generate bundles S that cluster tightly.  By assumption, the corresponding list of gradients ∇ f(S) is robustly affinely independent, so this gradient information is inconsistent with any smooth model for the objective function f:   we instead must seek a simple, well-behaved, nonsmooth model.  We use a max function model, motivated by Example <ref>, and consider twice continuously differentiable functions f_s → satisfying 
        
            f_s(s) = f(s),   ∇ f_s(s) = ∇ f(s),   ∇^2 f_s(s) = ∇^2 f(s),    f_s(s') < f(s')
        
        for distinct points s,s' ∈ S.  The precise form of these functions is immaterial to the algorithm.
        It is nonetheless reassuring to note that such functions always exist.  We could for example define
        
            f_s(x)  =  q_s(x) - α|x-s|^4    (x ∈),
        
        for a sufficiently large constant α > 0.  On the other hand, if in fact f is a max-function, we could simply consider each f_s as one of the functions comprising the pointwise maximum.
        
        Now we consider the (unknown) function f̃→ defined by
        
            f̃(x) = max_s ∈ S f_s(x)     (x ∈),
        
        as our working model of the function f:  it agrees with f up to second order at each of the reference points s ∈ S.  We can minimize this model via the classical nonlinear program
        
            {[                     t;            f_s(x) - t          ≤          0    (s ∈ S);                  x ∈,                  t ∈. ].
        
        Since the functions f_s are unknown, we cannot solve this problem exactly.
        Instead, we consider a feasible solution (ŝ,f(ŝ)), for any point ŝ∈ S, and follow (loosely) a classical sequential quadratic programming approach to improve it.  We remark that SQP techniques have some history in the nonsmooth optimization literature <cit.>.  
        
        A standard SQP approach <cit.> would proceed in two steps, the first of which estimates the Lagrange multipliers.  Taking, as a first approximation, each inequality constraint to be active, we seek to solve 
        
            min_λ∈Δ_S| ∑_s ∈ Sλ_s ∇ f_s(ŝ) |
         
        Approximating the point ŝ by the point s in each summand leads to the optimality measure in the algorithm:
        
            Θ(S) = min_λ∈Δ_S| ∑_s ∈ Sλ_s ∇ f(s) |.
         
        
        Fixing the resulting Lagrange multiplier estimate λ, the Lagrangian for the nonlinear program is
        
            (x,t) ↦∑_s ∈ Sλ_s f_s(x).
        
        The second SQP step then aims to reduce its quadratic model at the feasible solution 
        (ŝ,f(ŝ)) over a feasible region defined by linearized constraints.
        We approximate the quadratic model by the function
        
            (x,t) ↦∑_s ∈ Sλ_s q_s(x),
        
        and we approximate the linearized feasible region, using the active subspace from Section <ref>, as
        
            {( x,t) : l_s(x)=t   s ∈ S }.
        
        We hence arrive at exactly the quadratic subproblem in the algorithm.  The subproblem is feasible, as we saw in Section <ref>, and bounded below by our assumption of strong convexity.
        
        This loose explanation can be tightened.  In fact, the proof of Theorem <ref> will show that, if the point x̃∈ minimizes the model f̃, and the quantity ν = max |S - x̃| is small, then, under reasonable conditions the solution x̂ of the quadratic subproblem satisfies |x̂ - x̃| = O(ν^2).  In other words, the algorithm computes a good approximation of the minimizer of the model function as the next reference point.
        
        
        
        § A SEQUENTIAL QUADRATIC PROGRAMMING TOOL
        
        For the philosophy underlying the bundle Newton method, the tool we describe in this section is central.  It is a slight variant of a standard sequential quadratic programming technique <cit.>— for completeness, we prove it directly.  Convexity plays no role, throughout this section.
        
        Given functions h_i ^n →, for i=1,2,…,k, we consider an equality-constrained optimization problem of the form
        
            (P)     {[                             c^T y;                            h_i(y)                = 0    (i=1,2,…,k);                                 y                ∈              ^n. ].
        
        For our purposes, a linear objective function c^T y (for some vector c ∈) suffices.  
        
        We say that a point y̅∈^n satisfies the  strong second-order sufficient conditions if h_i(y̅) = 0 for each i (feasibility), each h_i is twice continuously differentiable around 
        y̅, the list of constraint gradients
        
            T  = {∇ h_i(y̅) : i=1,2,…,k }
        
        is linearly independent, and there exists a Lagrange multiplier vector λ̅∈^k (necessarily unique) satisfying
        
            ∑_i λ̅_i ∇ h_i(y̅) = -c         ∑_i λ̅_i ∇^2 h_i(y̅)    T^⊥.
        
        The point y̅ is then a strict local minimizer for the problem (P).
        
        We next consider a Lagrange multiplier estimate λ∈^k close to λ̅.  The traditional approach of sequential quadratic programming linearizes the constraints around a trial point close to the minimizer y̅, and replaces the objective by the corresponding quadratic approximation to the Lagrangian c^T y + ∑_i λ_i h_i(y).  Instead, we use a  different reference point 
        y_i ∈^n (near y̅) for each constraint, forming the corresponding linear approximations 
        p_i ^n → defined by
        
            p_i(y)  =  h_i(y_i) + ∇ h_i(y_i)^T (y-y_i)    (i=1,2,…,k).
         
        We denote by Y the reference points (y_1,y_2,…,y_k) in the product space (^n)^k, which is close to 
        Y̅ = (y̅,y̅,…,y̅), and consider the following quadratic program, parametrized by Y and λ:
        
            (QP)     {[                                                                   c^T y + 
            ∑_i λ_i ( p_i(y) + 1/2 (y-y_i)^T ∇^2 h_i(y_i)(y - y_i) );                                                                                                           p_i(y) = 0    (i=1,2,…,k);                                                                                                                              y ∈^n. ].
        
        The new list of constraint gradients {∇ h_i(y_i) } is also linearly independent, so any minimizer y ∈^n for the quadratic program (QP) must satisfy the conditions for a  stationary point:
        
            p_i(y)    =    0    (i=1,2,…,k) 
            
            c + 
            ∑_i λ_i ( ∇ h_i(y_i) +∇^2 h_i(y_i)(y - y_i) )
                =    ∑_i μ_i ∇ h_i(y_i)
        
        for some multiplier vector μ∈^k.  
        
        
        Consider a point y̅∈^n satisfying the strong second-order sufficient conditions for the problem 
        (P).  Then for all Y ∈ (^n)^k near Y̅, and any multiplier vector 
        λ = λ̅+ O(Y-Y̅) in ^k, the quadratic program (QP) has a unique stationary point ŷ = y̅ + O(Y-Y̅^2), which furthermore is the unique minimizer.
        
        We can write the stationary point conditions as a linear system:
        
            (M(Y,λ))(y,μ)  =  b(Y,λ),
        
        for a linear operator M(Y,λ) on ^n ×^k and a vector 
        b(Y,λ) ∈^n ×^k, both depending continuously on the parameter (Y,λ).
        When (Y,λ) = (Y̅,λ̅), the corresponding homogeneous system is 
        
            ∇ h_i(y̅)^T y    =    0    (i=1,2,…,k) 
            ∑_i λ̅_i  ∇^2 h_i(y̅)y 
                =    ∑_i μ_i ∇ h_i(y̅).
        
        By the second-order sufficient conditions, this system has only the trivial solution.  Hence the operator M(Y̅,λ̅) is invertible.
        
        As γ = Y-Y̅→ 0 with |λ - λ̅|  = O(γ), we have
        
            p_i(y̅) = O(γ^2)    (i=1,2,…,k)
        
        and
        
            c + ∑_i λ_i ( ∇ h_i(y_i) +∇^2 h_i(y_i)(y̅ - y_i) ) 
               =   ∑_i (λ_i - λ̅_i) ∇ h_i(y̅) + O(γ^2)  
               =   ∑_i (λ_i - λ̅_i) ∇ h_i(y_i) + O(γ^2).
        
        We deduce 
        
            (M(Y,λ))(y̅,λ - λ̅)  -  b(Y,λ)
             = 
            O(γ^2).
        
        The norm of the inverse of M(Y,λ) is uniformly bounded for (Y,λ) near 
        (Y̅,λ̅), so
        
            (y̅,λ - λ̅)  - (M(Y,λ))^-1( b(Y,λ) )
             = 
            O(γ^2).
        
        So there exists a unique stationary point ŷ = y̅ + O(Y-Y̅^2).
        
        Providing γ is sufficiently small, the point ŷ is in fact a strict minimizer for the quadratic program, since it also satisfies the sufficient condition
        that ∑_i λ_i ∇^2 h_i(y_i) is positive definite on the subspace
        
            {
            z ∈ : 
            ∇ h_i(y_i)^T z = 0   i
            }.
        
        Otherwise there would exists sequences of points
        
            { y_i^r : i=1,2,…,k },
        
        and corresponding sequences of multipliers λ_i^r
        satisfying  y_i^r →y̅ for each i 
        and λ_i^r →λ̅ as r →∞, and unit vectors z^r ∈ satisfying
        
            (z^r)^T ( ∑_i λ_i^r ∇^2 h_i(y_i^r) ) z^r ≤ 0     ∇ h_i(y_i^r)^T z^r = 0   i.
        
        In that case, after taking a subsequence, we can suppose that z^r converges to a unit vector z ∈ satisfying
        
            z^T ( ∑_i λ̅_i ∇^2 h_i(y̅) ) z ≤ 0     ∇ h_i(y̅)^T z = 0   i,
        
        in contradiction to the second-order sufficient conditions.
        
        
        § THE MAX FUNCTION CASE
        
        In this section we analyze carefully how the bundle Newton method, Algorithm <ref>, behaves when minimizing a max function f →, as in Example <ref>.  Thus the function f has the form
        
        f(x) = max_i=1,…,k f_i(x),    (x ∈),
        
        for some  unknown functions f_i ^n → that are now assumed to be twice continuously differentiable for i=1,2,…,k.  
        
        Like the last section, the development in this section does not rely on convexity, with one key exception.  The proof of the final quadratic convergence result, Theorem <ref> depends critically on strong convexity.
        
         Note.
        We emphasize a crucial point.  Our interest in max functions is as local models for more general objectives, and as a test bed for the general-purpose Algorithm <ref>.  We could easily minimize an explicitly described objective of the form (<ref>) by applying a standard algorithm to the equivalent inequality-constrained optimization problem
        
            (IP)    {[                                 t;                        f_i(x) - t                ≤ 0    (i=1,2,…,k);                               x ∈                              t ∈. ].
        
        However, we seek to minimize the objective function f using an oracle  with no access to the individual functionsf_i.  
        
        
        We consider Algorithm <ref> for the objective f, on a neighborhood of a point 
        .  Corresponding to the fixed (but implicit) representation (<ref>) of f, we assume  the strong second-order conditions defined below.
        
        Given a max function representation of the form (<ref>), we say that a point 
        ∈ satisfies the  strong second-order conditions when the following properties hold.
        
          *  Full activity:  the values f_i(x̅) are equal for all i.
        
          *  Independence:  the gradients {∇ f_i(x̅) : i=1,2,…,k }
        are affinely independent.
        
          *  Stationarity:  
        There exists a Lagrange multiplier vector λ∈^k_+ (necessarily unique) satisfying 
        ∑_i λ_i = 1 and 
        ∑_i λ_i ∇ f_i(x̅) = 0.
        
          *  Strict complementarity:
        λ_i > 0 for all i.
        
          *   Second-order sufficiency:
        The Lagrangian ∑_i λ_i ∇^2 f_i(x̅) is positive definite on the subspace
        {
        z ∈ : 
        ∇ f_i(x̅)^T z   i}.
        
        
        These assumptions are closely related both to problem (IP) and to Example <ref>.  If we consider the feasible point (,f()) for (IP), full activity amounts to all the constraints being active, independence amounts to the usual linear independence constraint qualification, and the remaining three conditions correspond exactly to the analogous conditions for (IP).  Classically, these conditions are sufficient for the point (,f()) to be a strict local minimizer for (IP), and hence for  to be a strict local minimizer for the max function f.  We imposed the first four conditions in Example <ref>:  there, we referred to the strict complementarity assumption as nondegeneracy, since, assuming the first three conditions, it amounts to 
        0 ∈(∂ f()) (even in the nonconvex case).
        
        We refer to the disjoint open sets
        
            𝒟_i  = { x ∈ : f_i(x) > f_j(x)  (j  i) },
        
        as  activity regions:  the values of the functions f and f_i coincide on 𝒟_i, as do their gradients ∇ f and ∇ f_i, and their Hessians ∇^2 f and ∇^2 f_i.  At any point in the open set
        
            𝒟 = ⋃_i=1^k 𝒟_i,
        
        we suppose that an oracle returns the value of f along with its gradient and Hessian.  Our challenge in minimizing f is that we only have access to precise information about the value of each component function f_i, its gradient and Hessian, on the region 𝒟_i:  elsewhere in  we only know the bound f_i ≤ f.
        
        The algorithm stops if it encounters a point outside 𝒟. This is a reasonable assumption in practice, since the complement 𝒟^c is typically a small set.  In particular, around the local minimizer , since the gradients ∇ f_i() are all distinct (being affinely independent), 𝒟^c is contained in the union of the manifolds 
        (f_i - f_j)^-1(0) for i  j.  Thus 𝒟 is a dense open set around .
        
        We consider a closed ball B_γ() of small radius γ > 0 around .  At the outset of the algorithm, we consider a  full bundle S ⊂ B_γ(),
        meaning that it contains exactly one reference point in each of the k activity regions 𝒟_i. 
        We can therefore write
        
        S  = { x_1,x_2,…,x_k },      x_i ∈𝒟_i    i=1,2,…,k.
        
        We will prove that this property is maintained as the algorithm proceeds.  
        
        We denote by σ̅ a certain kth largest singular value: 
        σ̅ = σ_k
        ( [
        [ ∇ f_1() ∇ f_2()       ⋯ ∇ f_k();       1       1       ⋯       1 ]] ).
        
        The affine independence assumption implies σ̅> 0.
        
        
        For any tolerance σ∈ (0,σ̅), full bundles S close to  always satisfy the robust affine independence condition σ_S(∇ f(S)) > σ.
        
        If the result fails, then for each index i=1,2,…,k there exists a sequence of points (x_i^r) in 
        the activity region 𝒟_i converging to the minimizer , such that 
        
            σ_k ( [
            [ ∇ f_1(x_1^r) ∇ f_2(x_2^r)            ⋯ ∇ f_k(x_k^r);            1            1            ⋯            1 ]] )
            ≤σ < σ̅
        
        for infinitely many r.
        But that contradicts the continuity of the kth largest singular value.
        
        
        We now consider Algorithm <ref> with the choice of tolerances ϵ̅= 0,
        δ̅= 0, so that the optimality checks never cause the method to stop.  (The affine independence check discussed in Section <ref> will not stop the algorithm either, providing we fix the tolerance σ∈ (0,σ̅).)  We begin each iteration by forming the corresponding linear and quadratic approximations:
        
            l_i(·)    =    f_i(x_i) + ∇ f_i(x_i)^T(· - x_i),  
            
            q_i(·)    =    l_i(·) + 1/2(· - x_i)^T ∇^2f_i(x_i)(· - x_i),
        
        and estimate the Lagrange multiplier vector.  We have the following result.  
        
        
        For any small radius γ > 0, there exists a unique minimizer for the problem
        
            min_λ∈^k{| ∑_i λ_i ∇ f_i(x_i) | : ∑_i λ_i = 1 },
        
        and it satisfies λ = λ̅+ O(γ).
        
        A vector λ∈^k solves the problem if and only if there exists a number α∈ such that
        
            ∑_i=1^k λ_i = 1     α + ∇ f_j(x_j)^T ( ∑_i=1^k λ_i ∇ f_i(x_i) ) = 0     (j=1,2,…,k).
        
        This square linear system is defined by an operator depending smoothly on its parameters, the points x_i (for i=1,2,…,k).
        In the limit, when x_i = x̅ for each i, the system becomes
        
            ∑_i=1^k λ_i = 1     α + ∇ f_j()^T ( ∑_i=1^k λ_i ∇ f_i(x̅) ) = 0     (j=1,2,…,k).
        
        The corresponding homogeneous system has only the trivial solution, by affine independence of the set 
        {∇ f_i(x̅) }, so the defining operator is invertible.  The unique solution of this limiting system is clearly (λ̅, 0), and the result now follows.
        
        
        As a consequence of this result and the strict complementarity assumption, the Lagrange multiplier estimates are all positive for small radius γ.
        We next turn to the computation of the new reference point x̂.
        
        Being a fully active strict local minimizer of the inequality-constrained problem (IP), the point (,f()) is also a local minimizer for the more restrictive problem
        
            (P')     {[                                 t;                        f_i(x) - t                = 0    (i=1,2,…,k);                               x ∈                              t ∈. ].
        
        Furthermore, as is easy to verify, it satisfies the strong second-order sufficient conditions with the same Lagrange multiplier vector λ̅.  Hence we can apply Theorem <ref>.  The corresponding quadratic program is
        
            {[                                 t + 
            ∑_i λ_i ( q_i(x) - t );                               l_i(x) - t = 0    (i=1,2,…,k);                                                 x ∈,   t ∈, ].
        
        or equivalently, exactly the quadratic program in Algorithm <ref>:
        
            {[                        ∑_i λ_i q_i(x);                    l_i(x)   i=1,2,…,k;                                  x ∈. ].
        
        By Theorem <ref>, this quadratic program has a unique minimizer x̂∈, which furthermore satisfies x̂ = x̅ + O(γ^2).  
        
        The final step in the iteration substitutes the new reference point x̂ for one of the existing reference points.  Since x̂ is also close to the point , this substitution produces another full bundle.  To see this, the following idea is the key tool.  
        
        
        Near , the optimality measures of full bundles are always strictly less than those of bundles that are not full.
        
        If the result fails, then there exists a sequence of full bundles S_r, and a sequence of not full bundles S'_r, both shrinking to , with 
        
            Θ(S_r) ≥Θ(S'_r)      (r=1,2,3,…).
         
        The left-hand side converges to zero, because
        
            0  ∈ {∇ f_i () : i=1,2,…,k},
        
        and hence so does the right-hand side.  After taking a subsequence, we can suppose that there is an index j such that S'_r ∩𝒟_j = ∅, and hence
        
            lim inf_r Ω(S'_r)  ≥ min | {∇ f_i () : i  j}|.
        
        But the right-hand side is strictly positive, because
        
            0  ∉ {∇ f_i () : i  j}.
        
        This contradiction completes the proof.
        
        Since the reference points stay close to the point , we next deduce that we maintain full bundles as the algorithm progresses, as follows.
         
        
        For any full bundle S near , and any new reference point x̂∈𝒟 near , there is a unique reference point s ∈ S minimizing the optimality measure of the new bundle 
        S' = (S ∖{s}) ∪{x̂}, and S' is then also a full bundle.
        
        The previous result shows that when x̂ lies in the activity region 𝒟_i, the unique optimal choice of s is the unique reference point in 𝒟_i.  The result  then follows.
        
        
        In summary, we have proved the following result.  
        
        
        Given a max function representation of the objective
        
            f(x) = max_i=1,…,k f_i(x)      (x ∈),
        
        suppose the point  satisfies the strong second-order conditions in Definition <ref>.  Then there exists a constant M>0 such that the k-bundle Newton method (Algorithm <ref>), with the tolerances ϵ̅= 0, δ̅= 0, has the following property.  
        Any iteration starting with a full bundle S sufficiently close to 
         generates a new reference point x̂ satisfying
        
            |x̂ - |  ≤  M max_s ∈ S |s-|^2,
        
        and assuming x̂∈𝒟, then
        generates a new full bundle by substituting x̂ for the unique reference point in S from the same activity region.
        
        
        While this is a suggestive result, it does not yet guarantee convergence.  
        To ensure that the sequence of bundles shrinks to x̅, we finally call on our assumption that the objective f is strongly convex.  (In the next section we discuss how to modify the algorithm to handle nonconvex objectives.)
        
        We first develop a simple tool.  As usual, for vectors z ∈^k we define z_max = max_j |z_j|.
        
        
        Given constants ϵ,M > 0, consider any sequence of vectors in the orthant _+ with the property that, for each successive pair z,z' in the sequence, there exists an index i such that z'_j = z_j for all j  i, and furthermore
        
            z_i ≥ϵz_max        
            z'_i ≤ M z_max^2.
        
        Then providing the initial vector is sufficiently small, the sequence converges to zero at a k-step quadratic rate.
        
        By induction we see that z_max is nondecreasing as the vector z evolves along the sequence, providing that the initial vector is sufficiently small. Suppose z = z_ at the outset of some iteration, and set 
        θ = z__max.  At this and every subsequent iteration, the updated component z_i is always set to a new value in the interval (0,Mθ^2].  Each updated component therefore cannot be updated again until we have 
        
            max_j=1,2,…,k z_j  ≤ Mθ^2/ϵ.
        
        This inequality must therefore hold after at most k iterations, at which point, 
        if z = z_, we have
        
            z__max ≤ M/ϵz_^2_max,
        
        which completes the proof.
        
        
        We can now prove our main result.
        
        
        
        Given a max function representation of the objective
        
            f(x) = max_i=1,…,k f_i(x)      (x ∈),
        
        suppose that the point  satisfies the strong second-order conditions in Definition <ref>, with each Hessian ∇^2 f_i() positive definite.  Then, given the tolerances 
        ϵ̅= 0, δ̅= 0, 
        
        the k-bundle Newton method (Algorithm <ref>) starting from any full bundle in a neighborhood of , either stops at a point outside the set 𝒟, or generates a sequence of full bundles that converge k-step quadratically to .
        
        Given any small radius γ > 0, each function f_i is ρ-strongly convex on B_γ(), for some constant ρ > 0.  In fact, the proof that follows takes place entirely in the ball B_γ(), so we lose no generality in assuming that each function f_i is convex.  
        
        According to Theorem <ref>, there exists a constant M > 0 such that, for small enough 
        γ > 0, starting from any full bundle in the ball B_γ(), Algorithm <ref> produces a sequence of full bundles
        
            S  = { x_1,x_2,…,x_k },      
            x_i ∈𝒟_i ∩ B_γ()    i=1,2,…,k,
        
        and at each iteration replaces a reference point, x_i for some index i with a new reference point
        x̂∈  D_i ∩ B_γ() which furthermore satisfies
        
        |x̂ - |  ≤  M max_j=1,…,k |x_j-|^2.
        
        By the construction, the new reference point x̂ satisfies
        
            l_j(x̂) = l_i(x̂)     j=1,2,…,k,
        
        and since the functions f_i are twice continuously differentiable, there exists a constant R > 0 such that
        
            f_i(x)  ≤  l_i(x) + R/2|x-x_i|^2     x ∈ B_γ().
        
        On the other hand, by strong convexity we have
        
            f_j(x)  ≥  l_j(x) + ρ/2|x-x_j|^2     x ∈ B_γ(),  j=1,2,…,k.
        
        Therefore we deduce
        
            l_i(x̂) + R/2|x̂-x_i|^2 
             ≥  
            f_i(x̂)
             =  
            f_j(x̂)
             ≥ 
            l_j(x̂) + ρ/2|x̂-x_j|^2
             = 
            l_i(x̂) + ρ/2|x̂-x_j|^2,
        
        and setting α = √(ρ/R) > 0 gives
        
            |x̂-x_i|  ≥ α |x̂-x_j|    j=1,2,…,k.
        
        
        Let β = max_j |x_j - |, so by inequality (<ref>) we have
        |x̂ - | ≤ Mβ^2.  The preceding inequality implies, for each j=1,2,…,k, the inequality
        
            α |x_j - | 
               ≤   α |x_j - x̂| + α |x̂ - | 
             ≤  
            |x_i - x̂| + Mαβ^2 
               ≤    
            |x_i - | + | - x̂| + Mαβ^2
            ≤ |x_i - | +  M(α+1)β^2
        
        Maximizing over j implies
        
            |x_i - | ≥αβ - M(α+1)β^2 ≥α/2β
        
        providing the radius γ is sufficiently small.  So we conclude that the old reference point x_i and the new reference point x̂ satisfy the two key inequalities
        
            |x_i - | ≥α/2max_j |x_j - |
                    
            |x̂ - | ≤ Mmax_j |x_j - |^2.
        
        We now define
        
            z_j = |x_j - | ∈ (0,γ)       (j=1,2,…,k),
         
        and apply Lemma <ref> to complete the proof.
          
        
        Notice that the final assumption in the strong second-order conditions — that the Hessian of the Lagrangian is positive definite on the tangent subspace — is in fact superfluous here, since we are assuming that each Hessian ∇^2 f_i() is positive definite.
        
        
        
        § MINIMIZING SMOOTH-NONSMOOTH SUMS
        
        With easy modifications, we can use the same Newton bundle idea as in Algorithm <ref> to minimize a function of the form F = f + r, where the first component f → is strongly convex but nonsmooth as before, and the second component r → is nonconvex but smooth (twice continuously differentiable).  
        
        Using the natural choice of optimality measure for a bundle S ⊂𝒟, namely
        Θ(S)  = min| ( ∇ F(S) ) |.
        
        the new algorithm proceeds as follows.
        
        [k-bundle Newton method for F=f+r]
        
        
        
        We make some immediate comments in comparison with Algorithm <ref> for minimizing f alone.  Important to notice is that the affine functions l_s (for s ∈ S), whose equality defines the subproblem constraints, are linear approximations for the component f, rather than for the objective F:  this subtlety is important for the analysis, which relies on strong convexity of f.  In contrast, the functions q_s appearing in the subproblem objective must be quadratic approximations for F to ensure quadratic convergence.  
        
        In this new setting we should consider the possibility of an unbounded quadratic subproblem, since the quadratic approximations are no longer necessarily convex.  For the time being, to emphasize the parallel with Algorithm <ref>, we omit this check, but include it in the final version (still to come).
        
        As in  Algorithm <ref>, we might also check affine independence of the gradients:
        whether we use the gradients of the component f or the objective F is immaterial, since for bundles of small diameter the two differ by roughly a constant due to the continuity of the gradient of the other component r.  When f is a max function, the analogue of Proposition <ref> still holds:  providing we choose a sufficiently small affine dependence parameter, full bundles close to 
         always satisfy the check.  Like the unboundedness check, we omit this check for now, but include it in the final version.
        
        The optimality measure (<ref>) is the natural extension of the earlier version, and leads to the precise analogue of the approximate optimality property (<ref>), as well as to our strategy for substituting the new reference point at the end of each iteration.
        
        All of our discussion of the choice of bundle size applies in this new case too, simply replacing the old objective f by the new objective F, and noting the simple relationship between the subdifferentials
        
            ∂ F(x) = ∂ f(x) + ∇ r(x)      (x ∈),
        
        where we understand the left-hand side in the Clarke sense <cit.>.  To motivate the quadratic subproblem, we follow exactly the same argument as in Section <ref>, but modifying the underlying nonlinear program to
        
            {[              r(x) + t;            f_s(x) - t          ≤          0    (s ∈ S);                  x ∈,                  t ∈. ].
        
        We note that the constraints are not changed:  they do not involve the function r.  Hence, as we noted, neither does the active subspace in the algorithm.
        
        We can then mimic the analysis of the original bundle-Newton method for max functions.  For notational simplicity, we proved the sequential quadratic programming tool, Theorem <ref>, for a linear objective, but the case we now need, involving a nonlinear objective, is almost identical.  
        
        We thus arrive at the main analysis, Section <ref>.  Given a representation of the nonsmooth component function f in the objective as
        
            f(x)  = max_i=1,…,k f_i(x)      (x ∈),
        
        we assume that the point ∈ satisfies the  strong second-order conditions for the problem of minimizing the objective f+r.  By this, we  mean that the feasible point (,f()) for the appropriate modification of the problem (IP), namely,
        
            {[                          r(x) + t;                        f_i(x) - t                ≤ 0    (i=1,2,…,k);                               x ∈                              t ∈, ].
        
        has all constraints active, and satisfies the linear independence constraint qualification, stationarity, strict complementarity, and the second-order sufficient conditions. Once again we note that the constraints are unchanged:  they do not involve the function r.  With a virtually identical proof , we arrive at the following generalization of Theorem <ref>.
        
        
         
        Given a max function representation of the nonsmooth function
        
            f(x) = max_i=1,…,k f_i(x)      (x ∈),
        
        and a twice continuously differentiable function r →, suppose that the point 
        ∈ satisfies the strong second-order conditions for minimizing the objective f+r, with each Hessian ∇^2 f_i() positive definite.  Then, given the tolerances 
        ϵ̅= 0 and δ̅= 0, 
        
        the k-bundle Newton method (Algorithm <ref>) starting from any full bundle in a neighborhood of , either stops at a point outside the set 𝒟, or generates a sequence of full bundles that converge k-step quadratically to .
        
        
        § MINIMIZING WEAKLY CONVEX FUNCTIONS
        
        Consider a  weakly convex function F →, by which we mean that, for some  weak convexity parameter value
        η, the function f = F + η/2|·|^2 is convex.  By increasing η if necessary, f becomes strongly convex.  Assuming that F is also twice continuously differentiable around every point in the set 𝒟, we can then define the smooth function r = -η/2|·|^2 and arrive at the following version of Algorithm <ref> for minimizing F, using the optimality measure:
        
            Θ(S)  = min| ( ∇ F(S) ) |.
        [k-bundle Newton minimization for weakly convex F]
        
        
        
        We note that the algorithm is almost identical to the original convex version, Algorithm <ref>, the only change (other than termination checks) being to the definition of the linear approximations defining the active subspace.  As the weak convexity parameter η grows large, the active subspace converges to the subspace of all points equidistant from all the reference points in the bundle.
        
        Theorem <ref> then specializes to our culminating result.  Recall that the constant σ̅> 0 is defined by equation (<ref>).
        
        
         
        Given a max function representation of the objective
        
            F(x) = max_i=1,…,k f_i(x)      (x ∈),
        
        suppose that the point ∈ satisfies the strong second-order conditions given in Definition <ref>.  Then, given the tolerances 
        ϵ̅= 0, δ̅= 0, small σ > 0, and sufficiently large weak convexity parameter η,
        the k-bundle Newton method (Algorithm <ref>) starting from any full bundle in a neighborhood of , either stops at a point outside the set 𝒟, or generates a sequence of full bundles that converge k-step quadratically to .
         Note:
        Any choice of parameters σ∈ (0,σ̅) and η strictly larger than the largest eigenvalues of each negative Hessian -∇^2 f_i() in fact suffices.
        
        This follows directly from Theorem <ref>, once we observe that, since each function f_i is twice continuously differentiable, with the given choice of weak convexity parameter 
        η the functions f_i + η/2|·|^2 are all strongly convex on a neighborhood of , and this local property suffices for the proof.
        
        
        § NUMERICAL EXPERIMENTS
        
        
        We illustrate the local bundle Newton method on several nonsmooth objective functions.  This handful of simple experiments is meant as a proof of concept rather than comprising any algorithmic recommendations.  Nonetheless, the results appear clearly promising enough to invite future research.
        
        
        
         §.§ Practical considerations
         
        We implemented none of the stopping criteria, simply terminating the algorithm manually when rounding error prevented any further progress.
        
        
        
          §.§.§ Choosing an initial bundle
        
        In each experiment, we ran a standard global nonsmooth optimization method to
        generate a finite set of points Ω⊂𝒟 near a minimizer x̅, and used the corresponding gradients to estimate the dimension of the subdifferential ∂ f(x̅) and 
        hence choose the bundle size k, as we discussed in Section <ref>.
        We then ran a heuristic subset selection procedure <cit.> to choose a set of
        k points in Ω with robustly affinely independent gradients to form the initial bundle.
        
        For convex problems, we implemented the simple “Bundle Method with Multiple Cuts”<cit.>, which we reproduce below in our notation.
        In the nonconvex case, we implemented the nonsmooth BFGS method <cit.>.
        
        [Multiple cut bundle method to minimize convex f]
        
        
        For the bundle method, we chose Ω to be the set of points whose cutting planes were strongly active in the final iteration.  That is,
        
            Ω =  s ∈ S : α_s > 0,
        
        where α_s is the dual variable associated with cutting plane l_s(·).
        For BFGS, we chose the set Ω to be the final 2n iterates.
        
        
        
          §.§.§ Solving the quadratic subproblems
        
        The algorithm involves two quadratic programming subproblems.  
        The first involves computing the optimality measure Θ(S) which amounts to projecting
        0 onto the convex hull of vectors ∇ f(s) : s ∈ S.  
        We implemented this as a quadratic program, solved in Gurobi.
        For the equality-constrained quadratic programs, 
        
            {[                ∑_s λ_s q_s(x)                              ;                    l_s(x) - t              =   0    (s ∈ S);                          t ∈,                          x ∈. ].
        
        we simply solve the (linear) optimality conditions directly:
        
            ∑_s ∈ Sλ_s ∇^2 f(s)(x - s) + ∑_s ∈ Sμ_s ∇ f(s)
                 = 0, 
            ∑_s ∈ Sμ_s    = 1, 
            
               l_s(x) - t    = 0    (s ∈ S).
        
        The x variable of the solution is then our Newton iterate x̂.
        
        
        
         §.§ Illustrative Examples
        
        
          §.§.§ A Strongly Convex Problem
        
        
        Our first experiment is to minimize max functions of the form
        
            f(x) = max_i = 1, … k*g_i^T x + 1/2 x^T H_i x + c_i/24 ||x||^4
        
        for 1 ≤ k ≤ n+1.
        We randomly generate positive constants c_i, symmetric positive definite matrices H_i, and 
        affinely independent random vectors g_i satisfying ∑_i λ_i g_i = 0
        for some λ randomly sampled in 
        λ > 0 : ∑_i λ_i = 1.
        Then 
        
            0 ∈∂ f(0) = g_i : 1 ≤ i ≤ k
        
        so f is nonsmooth at the minimizer of 0.
        The structure of f is unknown to the algorithms, whose access is limited to a 
        black box that returns function values, gradients, and Hessians.
        
        In random trials for dimension n=50, we applied the bundle method, Algorithm <ref>, in a first phase, with parameters ρ = 1 and β = 10^-5 and starting point z=(1, …, 1). The stopping tolerance was set to 10^-6, at which point we initialized and switched to 
        the bundle Newton method, Algorithm <ref>.  
        Results for a number of random trials are shown in Figures <ref> and <ref>.
        For the bundle method phase, we observed a roughly linear rate of convergence of function values to zero, proportional to the degree of nonsmoothness k.
        Switching to the bundle Newton method results in much more rapid convergence in accordance with the rates predicted by Theorem <ref>.
        
        
        
          §.§.§ A Nonconvex Problem
        
        To test the nonconvex version of the algorithm, we used Euclidean sum functions of the form
        
            f(x) = ∑_i=1^k | g_i^T x + 1/2 x^T H_i x + c_i/24 ||x||^4 |
        
        for 1 ≤ k ≤ n+1.  The constants c_i, vectors g_i and matrices H_i were randomly generated as in the previous experiment.  As usual, access to f was limited to a black box that returns function values, gradients, and Hessians.
        
        
        
        In random trials for dimension n=50, we applied nonsmooth BFGS in a first phase
        until a breakdown occurred due to numerical instability (as usual with this method <cit.>).   At this point we switched to the Algorithm <ref> with weak convexity parameter dynamically chosen as
        
            η = max_s ∈ Sλ_max(-∇^2 f(s) )
        
        at each iteration.
        
        
        
         §.§ Partly Smooth Functions
        
        
        While the theoretical results of this paper were limited to objective functions with finite max structure, experiments suggest that variants of the bundle Newton method may be effective much more broadly.  However, one particular implementation hurdle arises even for simple nonsmooth functions like the Euclidean norm: solving the system 
        (<ref>) directly will be numerically unstable,
        due to the ill-conditioning of the Hessians ∇^2 f(s).
        
        However, for the broad class of partly smooth functions in the sense of <cit.>, this ill-conditioning is often highly structured: the nonsmoothness is associated with a certain subspace 𝒱 spanned by the eigenvectors corresponding to large eigenvalues of the Hessian ∇^2 f(s).  Along an orthogonal manifold, the function behaves smoothly, with well conditioned Hessian.
        Motivated by this idea, we follow a simple strategy for solving the system (<ref>), similar to reduced system approaches for nonlinear programming described in standard texts <cit.>, and avoiding full Hessian computations.
        
        
        
          §.§.§ A reduced system approach
        
        
        Let G and b be a matrix and vector satisfying
        
            x : Gx = b = x : l_s(x)  equal for all s ∈ S,
         
        so that we can write the optimality conditions (<ref>) as
        
            ∑_s ∈ Sλ_s ∇^2 f(s)(x - s) + G^T ν   = -∑_s ∈ Sλ_s ∇ f(s), 
            
              Gx    = b.
        
        for ν∈^k-1.
        Suppose that we have found matrices U and V such that the matrix
        [ U V ]∈^n × n is full rank 
        and GU = 0 (via a QR factorization of G^T, for example).
        The columns of U are then a basis for the space (G), and 
        we can write any solution of (<ref>) as x = U x_u + V x_v.
        The constraint Gx = b then implies GV x_v = b,
        which can be solved for x_v, since we assume G (and hence GV) 
        is full rank.  We deduce
        
            x : Gx = b = (U) + p,
        
        where p is the particular solution V(GV)^-1 b.
        Substituting this into the stationarity condition and multiplying through by U^T
        yields the reduced system
        
            ∑_s ∈ Sλ_s U^T ∇^2 f(s)(Ux_u + p  - s) = -∑_s ∈ Sλ_s U^T∇ f(s).
        
        In a slight modification to the algorithm, if we project each reference point onto the active subspace we arrive at the linear system
        
            ∑_s ∈ Sλ_s U^T ∇^2 f(s)U x_u 
              = ∑_s ∈ Sλ_s [ (U^T ∇^2 f(s)U) U^T(s - p) - U^T∇ f(s) ].
        
        This system only involves the projected HessiansU^T ∇^2 f(s)U, which remain well conditioned if the span of V is close to the subspace 𝒱, a property that we have experimentally observed to hold in practice.
        
        
        
          §.§.§ An Eigenvalue Problem
        
        
        Our final experiment to illustrate the reduced systems approach is an eigenvalue problem.  Specifically, given
        symmetric matrices A_0, …, A_n ∈^m × m we seek to minimize
        
            f(x) = λ_max(A_0 + ∑_i=1^n x_i A_i ),
        
        where λ_max(·) is the largest eigenvalue function.
        Typically minimizers occur at points where λ_max has multiplicity t > 1, necessitating nonsmooth minimization 
        techniques. Under reasonable conditions, the set of points x ∈ for which λ_max has fixed multiplicity t is a manifold of codimension t(t+1)/2, relative to which f is partly smooth <cit.>.
        
        For illustration, in Figure <ref> we show convergence of the bundle method, BFGS, and bundle Newton method on a typical trial for this problem using random data.  All algorithms were run without termination conditions until numerical issues prevented any further progress.
        In this example for n=50 matrices in ^25 × 25,
        the optimal eigenvalue multiplicity was 6, and we again observe fast 
        convergence of the bundle Newton method once the subdifferential dimension 
        t(t+1)/2 - 1 = 20 can be identified.
        
        
        
        (Note that since the optimal objective value is unknown, we instead used the best value found after running the algorithms with a large number of random starting points. This introduces a slight bias in the accuracy reported for the bundle Newton method.)
        In Figure <ref>, we observe that the bundle Newton method
        achieves an eigenvalue clustering several orders of magnitude better than is possible with a bundle method or BFGS.
        
        
        
        Using the active manifold to accelerate eigenvalue optimization is not new <cit.>.  What is remarkable is that the bundle Newton method, combined with a first phase
        such as a traditional bundle method, rapidly convergences to the minimizer
        without any structural knowledge of the function.
        
        
        
          §.§.§ First-order analogues
        
        The Newton philosophy that we explore in this work is suggestive even in the more usual setting where Hessians are unavailable.  One straightforward first-order analogue of Algorithm 2,2 replaces the Hessians by suitably tuned multiples of the identity matrix.  Simple implementations seem effective on max functions: a broader investigation is the topic of ongoing work.
        
        
        
        
        ''10atkinson-vaidya
        D.S. Atkinson and P.M. Vaidya.
         A cutting plane algorithm for convex programming that uses analytic
          centers.
         Math. Programming, 69:1–43, 1995.
        
        bubeck
        S. Bubeck.
         Convex optimization: algorithms and complexity.
         Foundations and Trends in Machine Learning, 8:231–357, 2015.
        
        rob_grad_samp
        J.V. Burke, A.S. Lewis, and M.L. Overton.
         A robust gradient sampling algorithm for nonsmooth, nonconvex
          optimization.
         SIAM J. Optim., 15:751–779, 2005.
        
        clarke
        F.H. Clarke.
         Optimization and Nonsmooth Analysis.
         Wiley Interscience, New York, 1983.
        
        Coste-semi
        M. Coste.
         An Introduction to Semialgebraic Geometry.
         RAAG Notes, 78 pages, Institut de Recherche Mathématiques de
          Rennes, October 2002.
        
        du-ruszcynski
        Y. Du and A. Ruszczyński.
         Rate of convergence of the bundle method.
         J. Optim. Theory Appl., 173:908–922, 2017.
        
        goffin-vial
        J.-L. Goffin and J.-P. Vial.
         On the computation of weighted analytic centers and dual ellipsoids
          with the projective algorithm.
         Math. Programming, 60:81–92, 1993.
        
        golub-vanloan
        G.H. Golub and C.F. Van Loan.
         Matrix Computations.
         Johns Hopkins University Press, Baltimore, MD, fourth edition, 2013.
        
        kiwiel
        K. C. Kiwiel.
         Efficiency of proximal bundle methods.
         J. Optim. Theory Appl., 104:589–603, 2000.
        
        kiwiel-book
        K.C. Kiwiel.
         Methods of Descent for Nondifferentiable Optimization, volume
          1133 of  Lecture Notes in Mathematics.
         Springer-Verlag, Berlin, 1985.
        
        lee-sidford-wong
        Yin Tat Lee, A. Sidford, and Sam Chiu-wai Wong.
         A faster cutting plane method and its implications for combinatorial
          and convex optimization.
         In  FOCS 2015, pages 1049–1065. IEEE Computer Soc., Los
          Alamitos, CA, 2015.
        
        lemarechal-bundle
        C. Lemaréchal.
         An extension of Davidon methods to non differentiable problems.
         Math. Programming Stud., 3:95–109, 1975.
        
        preamb
        C. Lemaréchal, A. Nemirovskii, and Y. Nesterov.
         New variants of bundle methods.
         Math. Programming, 69:111–147, 1995.
        
        LOS
        C. Lemaréchal, F. Oustry, and C. Sagastizábal.
         The U-lagrangian of a convex function.
         Transactions of the American Mathematical Society,
          352:711–729, 2000.
        
        LS97
        C. Lemaréchal and C. Sagastizábal.
         Practical aspects of the Moreau-Yosida regularization:
          theoretical preliminaries.
         SIAM J. Optim., 7:367—385, 1997.
        
        Lewis-active
        A.S. Lewis.
         Active sets, nonsmoothness, and sensitivity.
         SIAM J. Optim., 13:702–725, 2002.
        
        BFGS
        A.S. Lewis and M.L. Overton.
         Nonsmooth optimization via quasi-Newton methods.
         Math. Program., 141:135–163, 2013.
        
        luksan-vlcek
        L. Lukšan and J. Vlček.
         A bundle-Newton method for nonsmooth unconstrained minimization.
         Math. Programming, 83:373–391, 1998.
        
        mifflin-semismooth
        R. Mifflin.
         An algorithm for constrained optimization with semismooth functions.
         Math. Oper. Res., 2:191–207, 1977.
        
        MC02
        R. Mifflin and C. Sagastizábal.
         Proximal points are on the fast track.
         Journal of Convex Analysis, 9:563—579, 2002.
        
        MC05
        R. Mifflin and C. Sagastizábal.
         A VU-algorithm for convex minimization.
         Math. Program., 104:583–608, 2005.
        
        miller-malick
        S.A. Miller and J. Malick.
         Newton methods for nonsmooth convex minimization: connections among
          𝒰-Lagrangian, Riemannian Newton and SQP methods.
         Math. Program., 104:609–633, 2005.
        
        nem-yud
        A. S. Nemirovsky and D. B. Yudin.
         Problem Complexity and Method Efficiency in Optimization.
         John Wiley & Sons, Inc., New York, 1983.
         Translated from the Russian and with a preface by E. R. Dawson.
        
        nesterov-cutting
        Y. Nesterov.
         Complexity estimates of some cutting plane methods based on the
          analytic barrier.
         Math. Programming, 69:149–176, 1995.
        
        Nes04b
        Y. Nesterov.
         Introductory Lectures on Convex Optimization.
         Kluwer Academic, Dordrecht, The Netherlands, 2004.
        
        nest_nonsm
        Y. Nesterov.
         Smooth minimization of non-smooth functions.
         Math. Program., 103:127–152, 2005.
        
        nesterov-nemirovskii
        Y. Nesterov and A. Nemirovskii.
         Interior-Point Polynomial Algorithms in Convex Programming.
         SIAM, Philadelphia, PA, 1994.
        
        nocedal_wright
        J. Nocedal and S.J. Wright.
         Numerical Optimization.
         Springer, New York, second edition, 2006.
        
        oustry99
        F. Oustry.
         The 𝒰-Lagrangian of the maximum eigenvalue function.
         SIAM J. Optim., 9:526–549, 1999.
        
        overton88
        M.L. Overton.
         On minimizing the maximum eigenvalue of a symmetric matrix.
         In  Linear algebra in signals, systems, and control, pages
          150–169. SIAM, Philadelphia, PA, 1988.
        
        parikh-boyd
        N. Parikh and S. Boyd.
         Proximal algorithms.
         Foundations and Trends in Optimization, 1:123–231, 2013.
        
        robinson-perturbed
        S.M. Robinson.
         Perturbed Kuhn-Tucker points and rates of convergence for a class
          of nonlinear-programming algorithms.
         Math. Programming, 7:1–16, 1974.
        
        ryu-boyd
        E.K. Ryu and S. Boyd.
         A primer on monotone operator methods.
         Appl. Comput. Math., 15:3–43, 2016.
        
        claudia-icm
        C. Sagastizábal.
         A VU-point of view of nonsmooth optimization.
         In  Proceedings of the International Congress of Mathematicians,
          Rio de Janeiro, volume 3, pages 3785–3806, 2018.
        
        shapiro-fan
        A. Shapiro and M.K.H. Fan.
         On eigenvalue optimization.
         SIAM J. Optim., 5:552–569, 1995.
        
        vaidya
        P. M. Vaidya.
         A new algorithm for minimizing convex functions over convex sets.
         Math. Programming, 73:291–341, 1996.
