Proximal Distance Algorithms: Theory and Practice.

Proximal distance algorithms combine the classical penalty method of constrained minimization with distance majorization. If f(x) is the loss function, and C is the constraint set in a constrained minimization problem, then the proximal distance principle mandates minimizing the penalized loss f ( x ) + ρ 2 dist ( x , C ) 2 and following the solution x ρ to its limit as ρ tends to ∞. At each iteration the squared Euclidean distance dist(x,C)2 is majorized by the spherical quadratic ‖x- P C (x k )‖2, where P C (x k ) denotes the projection of the current iterate x k onto C. The minimum of the surrogate function f ( x ) + ρ 2 ‖ x - P C ( x k ) ‖ 2 is given by the proximal map prox ρ -1f [P C (x k )]. The next iterate x k+1 automatically decreases the original penalized loss for fixed ρ. Since many explicit projections and proximal maps are known, it is straightforward to derive and implement novel optimization algorithms in this setting. These algorithms can take hundreds if not thousands of iterations to converge, but the simple nature of each iteration makes proximal distance algorithms competitive with traditional algorithms. For convex problems, proximal distance algorithms reduce to proximal gradient algorithms and therefore enjoy well understood convergence properties. For nonconvex problems, one can attack convergence by invoking Zangwill's theorem. Our numerical examples demonstrate the utility of proximal distance algorithms in various high-dimensional settings, including a) linear programming, b) constrained least squares, c) projection to the closest kinship matrix, d) projection onto a second-order cone constraint, e) calculation of Horn's copositive matrix index, f) linear complementarity programming, and g) sparse principal components analysis. The proximal distance algorithm in each case is competitive or superior in speed to traditional methods such as the interior point method and the alternating direction method of multipliers (ADMM). Source code for the numerical examples can be found at https://github.com/klkeys/proxdist.

Introduction

The solution of constrained optimization problems is part science and part art. As mathematical scientists explore the largely uncharted territory of high-dimensional nonconvex problems, it is imperative to consider new methods. The current paper studies a class of optimization algorithms that combine Courant’s penalty method of optimization (Beltrami, 1970; Courant, 1943) with the notion of a proximal operator (Bauschke and Combettes, 2011; Moreau, 1962; Parikh and Boyd, 2013). The classical penalty method turns constrained minimization of a function f(x) over a closed set C into unconstrained minimization. The general idea is to seek the minimum point of a penalized version f(x)+ρq(x) of f(x), where the penalty q(x) is nonnegative and vanishes precisely on C. If one follows the solution vector x as ρ tends to ∞, then in the limit one recovers the constrained solution. The penalties of choice in the current paper are squared Euclidean distances dist(x,C)2 = inf ‖x−y‖2.

The formula defines the proximal map of a function f(). Here ‖ · ‖ is again the standard Euclidean norm, and f() is typically assumed to be closed and convex. Projection onto a closed convex set C is realized by choosing f() to be the 0/∞ indicator δ() of C. It is possible to drop the convexity assumption if f() is nonnegative or coercive. In so doing, prox(y) may become multi-valued. For example, the minimum distance from a nonconvex set to an exterior point may be attained at multiple boundary points. The point in the definition (1) can be restricted to a subset S of Euclidean space by replacing f() by f() + δ(), where δ() is the indicator of S.

One of the virtues of exploiting proximal operators is that they have been thoroughly investigated. For a large number of functions f(), the map prox() for c > 0 is either given by an exact formula or calculable by an efficient algorithm. The known formulas tend to be highly accurate. This is a plus because the classical penalty method suffers from ill conditioning for large values of the penalty constant. Although the penalty method seldom delivers exquisitely accurate solutions, moderate accuracy suffices for many problems.

There are ample precedents in the optimization literature for the proximal distance principle. Proximal gradient algorithms have been employed for many years in many contexts, including projected Landweber, alternating projection onto the intersection of two or more closed convex sets, the alternating-direction method of multipliers (ADMM), and fast iterative shrinkage thresholding algorithms (FISTA) (Beck and Teboulle, 2009; Combettes and Pesquet, 2011; Landweber, 1951). Applications of distance majorization are more recent (Chi et al., 2014; Lange and Keys, 2014; Xu et al., 2017). The overall strategy consists of replacing the distance penalty dist(x,C)2 by the spherical quadratic ‖ − ‖2, where is the projection of the kth iterate onto C. To form the next iterate, one then sets The MM (majorization-minimization) principle guarantees that x decreases the penalized loss. We call the combination of Courant’s penalty method with distance majorization the proximal distance principle. Algorithms constructed according to the principle are proximal distance algorithms.

The current paper extends and deepens our previous preliminary treatments of the proximal distance principle. Details of implementation such as Nesterov acceleration matter in performance. We have found that squared distance penalties tend to work better than exact penalties. In the presence of convexity, it is now clear that every proximal distance algorithm reduces to a proximal gradient algorithm. Hence, convergence analysis can appeal to a venerable body of convex theory. This does not imply that the proximal distance algorithm is limited to convex problems. In fact, its most important applications may well be to nonconvex problems. A major focus of this paper is on practical exploration of the proximal distance algorithm.

In addition to reviewing the literature, the current paper presents some fresh ideas. Among the innovations are: a) recasting proximal distance algorithms with convex losses as concave-convex programs, b) providing new perspectives on convergence for both convex and nonconvex proximal distance algorithms, c) demonstrating the virtue of folding constraints into the domain of the loss, and d) treating in detail seven interesting examples. It is noteworthy that some our new convergence theory is pertinent to more general MM algorithms.

It is our sincere hope to enlist other mathematical scientists in expanding and clarifying this promising line of research. The reviewers of the current paper have correctly pointed out that we do not rigorously justify our choices of the penalty constant sequence ρ. The recent paper by Li et al. (2017) may be a logical place to start in filling this theoretical gap. They deal with the problem of minimizing f() subject to = through the quadratic penalized objective . For the right choices of the penalty sequence ρ, their proximal gradient algorithm achieves a O(k−1) rate of convergence for f() strongly convex. As a substitute, we explore the classical problem of determining how accurately the solution of the problem approximates the solution of the constrained problem min f(x). Polyak (1971) demonstrates that f()−f() = O(−1) for a penalty function q() that vanishes precisely on C. Polyak’s proof relies on strong differentiability assumptions. Our proof for the case q() = dist(,C) relies on convexity and is much simpler.

As a preview, let us outline the remainder of our paper. Section 2 briefly sketches the underlying MM principle. We then show how to construct proximal distance algorithms from the MM principle and distance majorization. The section concludes with the derivation of a few broad categories proximal distance algorithms. Section 3 covers convergence theory for convex problems, while Section 4 provides a more general treatment of convergence for nonconvex problems. To avoid breaking the flow of our exposition, all proofs are relegated to the Appendix. Section 5 discusses our numerical experiments on various convex and nonconvex problems. Section 6 closes by indicating some future research directions.

Derivation

The derivation of our proximal distance algorithms exploits the majorization-minimization (MM) principle (Hunter and Lange, 2004; Lange, 2010). In minimizing a function f(), the MM principle exploits a surrogate function g( | ) that majorizes f() around the current iterate . Majorization mandates both domination g( | ) ≥ f() for all feasible and tangency g( | ) = f() at the anchor . If minimizes g( | ), then the descent property f() ≤ f() follows from the string of inequalities and equalities Clever selection of the surrogate g( | ) can lead to a simple algorithm with an explicit update that requires little computation per iterate. The number of iterations until convergence of an MM algorithm depends on how tightly g( | ) hugs f(). Constraint satisfaction is built into any MM algorithm. If maximization of f() is desired, then the objective f() should dominate the surrogate g( | ) subject to the tangency condition. The next iterate is then chosen to maximize g( | ). The minorization-maximization version of the MM principle guarantees the ascent property.

The constraint set C over which the loss f() is minimized can usually be expressed as an intersection of closed sets. It is natural to define the penalty using a convex combination of the squared distances. The neutral choice is one we prefer in practice. Distance majorization gives the surrogate function for an irrelevant constant c. If we put , then by definition the minimum of the surrogate g(x | x) occurs at the proximal point We call this MM algorithm the proximal distance algorithm. The penalty q() is generally smooth because at any point where the projection P() is single valued (Borwein and Lewis, 2006; Lange, 2016). This is always true for convex sets and almost always true for nonconvex sets. For the moment, we will ignore the possibility that P() is multi-valued.

For the special case of projection of an external point z onto the intersection C of the closed sets C, one should take . The proximal distance iterates then obey the explicit formula Linear programming with arbitrary convex constraints is another example. Here the loss is f() = , and the update reduces to If the proximal map is impossible to calculate, but f() is L-smooth (∇f() is Lipschitz with constant L), then one can substitute the standard majorization for f(). Minimizing the sum of the loss majorization plus the penalty majorization leads to the MM update This is a gradient descent algorithm without an intervening proximal map.

In moderate-dimensional problems, local quadratic approximation of f() can lead to a viable algorithm. For instance, in generalized linear statistical models, Xu et al. (2017) suggest replacing the observed information matrix by the expected information matrix. The latter matrix has the advantage of being positive semidefinite. In our notation, if ≈ d2f(), then an approximate quadratic surrogate is The natural impulse is to update x by the Newton step This choice does not necessarily decrease f(). Step halving or another form of backtracking restores the descent property.

A more valid concern is the effort expended in matrix inversion. If is dense and constant, then extracting the spectral decomposition of reduces formula (4) to which can be implemented as a sequence of matrix-vector multiplications. Alternatively, one can take just a few terms of the series when ρ is sufficiently large. For a generalized linear model, parameter updating involves solving the linear system for a diagonal matrix with positive diagonal entries. This task is equivalent to minimizing the least squares criterion In the unweighted case, extracting the singular value decomposition = facilitates solving the system of equations (5). The svd decomposition is especially cheap if there is a substantial mismatch between the number rows and columns of . For sparse , the conjugate gradient algorithm adapted to least squares (Paige and Saunders, 1982b) is subject to much less ill conditioning than the standard conjugate gradient algorithm. Indeed, the algorithm LSQR and its sparse version LSMR (Fong and Saunders, 2011) perform well even when the matrix is ill conditioned.

The proximal distance principle also applies to unconstrained problems. For example, consider the problem of minimizing a penalized loss ℓ()+p(). The presence of the linear transformation in the penalty complicates optimization. The strategy of parameter splitting introduces a new variable and minimizes ℓ() + p() subject to the constraint = . If P(z) denotes projection onto the manifold then the constrained problem can be solved approximately by minimizing the function for large ρ. If P(z) consists of two subvectors and corresponding to and , then the proximal distance updates are

Given the matrix is n × p, one can attack the projection by minimizing the function This leads to the solution If n < p, then the Woodbury formula reduces the expense of matrix inversion.

Traditionally, convex constraints have been posed as inequalities C = { : a() ≤ t}. Parikh and Boyd (2013) point out how to project onto such sets. The relevant Lagrangian for projecting an external point amounts to with λ ≥ 0. The corresponding stationarity condition can be interpreted as a[prox()] = t. One can solve this one-dimensional equation for λ by bisection. Once λ is available, = prox() is available as well. Parikh and Boyd (2013) note that the value a[prox()] is decreasing in λ. One can verify their claim by implicit differentiation of equation (7). This gives and consequently the chain rule inequality

Convergence: Convex Case

In the presence of convexity, the proximal distance algorithm reduces to a proximal gradient algorithm. This follows from the representation involving the penalty q(). Thus, the proximal distance algorithm can be expressed as In this regard, there is the implicit assumption that q() is 1-smooth. This is indeed the case. According to the Moreau decomposition (Bauschke and Combettes, 2011), for a single closed convex set C where is the Fenchel conjugate of the indicator function Because proximal operators of closed convex functions are nonexpansive (Bauschke and Combettes, 2011), the result follows for a single set. For the general penalty q() with m sets, the Lipschitz constants are scaled by the convex coefficients α and added to produce an overall Lipschitz constant of 1.

It is enlightening to view the proximal distance algorithm through the lens of concaveconvex programming. Recall that the function is closed and convex for any nonempty closed set C. Danskin’s theorem (Lange, 2016) justifies the directional derivative expression This equality allows us to identify the subdifferential ∂s() as the convex hull convP(x). For any ∈ ∂s(), the supporting hyperplane inequality entails where d is a constant not depending on . The same majorization can be generated by rearranging the majorization when is the convex combination of vectors p from P(). These facts demonstrate that the proximal distance algorithm minimizing is a special case of concave-convex programming when f() is convex. It is worth emphasizing that is often strongly convex regardless of whether f() itself is convex. If we replace the penalty dist(,C)2 by the penalty dist(,C)2 for a matrix , then the function s() is still closed and convex, and minimization of can also be viewed as an exercise in concave-convex programming.

In the presence of convexity, the proximal distance algorithm is guaranteed to converge. Our exposition relies on well-known operator results (Bauschke and Combettes, 2011). Proximal operators in general and projection operators in particular are nonexpansive and averaged. By definition an averaged operator is a convex combination of a nonexpansive operator N() and the identity operator . The averaged operators on with α ∈ (0,1) form a convex set closed under functional composition. Furthermore, M() and the base operator N() share their fixed points. The celebrated theorem of Krasnosel’skii (1955) and Mann (1953) says that if an averaged operator M() = α + (1 − α)N() possesses one or more fixed points, then the iteration scheme = M() converges to a fixed point.

These results immediately apply to minimization of the penalized loss Given the choice , the algorithm map = prox−1() is an averaged operator, being the composition of two averaged operators. Hence, the Krasnosel’skiiMann theorem guarantees convergence to a fixed point if one or more exist. Now z is a fixed point if and only if for all . In the presence of convexity, this is equivalent to the directional derivative inequality for all , which is in turn equivalent to z minimizing h(). Hence, if h() attains its minimum value, then the proximal distance iterates converge to a minimum point.

Convergence of the overall proximal distance algorithm is tied to the convergence of the classical penalty method (Beltrami, 1970). In our setting, the loss is f(), and the penalty is . Assuming the objective f() + ρq() is coercive for all ρ ≥ 0, the theory mandates that the solution path is bounded and any limit point of the path attains the minimum value of f() subject to the constraints. Furthermore, if f() is coercive and possesses a unique minimum point in the constraint set C, then the path converges to that point.

Proximal distance algorithms often converge at a painfully slow rate. Following Mairal (2013), one can readily exhibit a precise bound.

Proposition 1 Suppose C is closed and convex and f() is convex. If the point z minimizes , then the proximal distance iterates satisfy

The O(ρk−1) convergence rate of the proximal distance algorithm suggests that one should slowly send ρ to ∞ and refuse to wait until convergence occurs for any given ρ. It also suggests that Nesterov acceleration may vastly improve the chances for convergence. Nesterov acceleration for the general proximal gradient algorithm with loss ℓ() and penalty p() takes the form where L is the Lipschitz constant for ∇p(x) and d is typically chosen to be 3. Nesterov acceleration achieves an O(k−2) convergence rate (Su et al., 2014), which is vastly superior to the O(k−1) rate achieved by proximal gradient descent. The Nesterov update possesses the further desirable property of preserving affine constraints. In other words, if = and = , then = as well. In subsequent examples, we will accelerate our proximal distance algorithms by applying the algorithm map M() given by equation (2) to the shifted point z of equation (10), yielding the accelerated update = M(z). Algorithm 1 provides a schematic of a proximal distance algorithm with Nesterov acceleration. The recent paper of Ghadimi and Lan (2015) extends Nestorov acceleration to nonconvex settings.

In ideal circumstances, one can prove linear convergence of function values in the framework of Karimi et al. (2016).

Proposition 2 Suppose C is closed and convex and f() is L-smooth and μ-strongly convex. Then possesses a unique minimum point , and the proximal distance iterates xk satisfy

We now turn to convergence of the penalty function iterates as the penalty constants ρ tends to ∞. To simplify notation, we restrict attention to a single closed constraint set S. Let us start with a proposition requiring no convexity assumptions.

Proposition 3 If f() is continuous and coercive and S is compact, then the proximal distance iterates are bounded and the distance to the constraint set satisfies for some constant c. If in addition f() is continuously differentiable, then for some further constant d. Similar claims hold for the solutions yk of the penalty problem except that the assumption that S is compact can be dropped.

As a corollary, if the penalty sequence ρ tends to ∞, then all limit points of must obey the constraint. Proposition 3 puts us into position to prove the next important result.

Proposition 4 If f() is continuously differentiable and coercive and S is convex, then the penalty function iterates defined by satisfy where y attains the constrained minimum and d is the constant identified in Proposition 3.

Convergence: General Case

Our strategy for addressing convergence in nonconvex problems fixes ρ and relies on Zangwill’s global convergence theorem (Luenberger and Ye, 1984). This result depends in turn on the notion of a closed multi-valued map (). If converges to ∞ and ∈ N() converges to ∞, then for () to be closed, we must have ∞ ∈ N(∞). The next proposition furnishes a prominent example.

Proposition 5 If S is a closed nonempty set in , then the projection operator P() is closed. Furthermore, if the sequence is bounded, then the set ∪P() is bounded as well.

Zangwill’s global convergence theorem is phrased in terms of an algorithm map M() and a real-valued objective h(). The theorem requires a critical set Γ outside which M() is closed. Furthermore, all iterates ∈ M() must fall within a compact set. Finally, the descent condition h() ≤ h() should hold for all ∈ M(), with strict inequality when ∉ Γ. If these conditions are valid, then every convergent subsequence of tends to a point in Γ. In the proximal distance context, we define the complement of Γ to consist of the points with for all ∈ M(). This definition plus the monotonic nature of the proximal distance algorithm force the satisfaction of Zangwill’s final requirement. Note that if f() is differentiable, then a point belongs to Γ whenever 0 ∈ ∇f(x) + ρx − ρP(x).

In general, the algorithm map M(x) is multi-valued in two senses. First, for a given z ∈ P(), the minimum may be achieved at multiple points. This contingency is ruled out if the proximal map of f() is unique. Second, because S may be nonconvex, the projection may be multi-valued. This sounds distressing, but the points where this occurs are exceptionally rare. Accordingly, it makes no practical difference that we restrict the anchor points z to lie in P() rather than in convP().

Proposition 6 If S is a closed nonempty set in , then the projection operator P() is single valued except on a set of Lebesgue measure 0.

In view of the preceding results, one can easily verify the next proposition.

Proposition 7 The algorithm map M(x) is everywhere closed.

To apply Zangwill’s global convergence theory, we must in addition prove that the iterates x = M(x) remain within a compact set. This is true whenever the objective is coercive since the algorithm is a descent algorithm. As noted earlier, the coercivity of f() is a sufficient condition. One can readily concoct other sufficient conditions. For example, if f() is bounded below, say nonnegative, and S is compact, then the objective is also coercive. Indeed, if S is contained in the ball of radius r about the origin, then which proves that dist(,S) is coercive. The next proposition summarizes these findings.

Proposition 8 If S is closed and nonempty, the objective is coercive, and the proximal operator prox−1f(x) is everywhere nonempty, then all limit points of the iterates k+1 ∈ M() of the proximal distance algorithm occur in the critical set Γ.

This result is slightly disappointing. A limit point could potentially exist with improvement in the objective for some but not all ∈ convP(). This fault is mitigated by the fact that P() is almost always single valued. In common with other algorithms in nonconvex optimization, we also cannot rule out convergence to a local minimum or a saddlepoint. One can improve on Proposition 8 by assuming that the surrogates g( | ) are all μ-strongly convex. This is a small concession to make because ρ is typically large. If f() is convex, then g( | ) is ρ-strongly convex by definition. It is also worth noting that any convex MM surrogate g( | ) can be made μ-strongly convex by adding the viscosity penalty majorizing 0. The addition of a viscosity penalty seldom complicates finding the next iterate x and has little impact on the rate of convergence when μ > 0 is small.

Proposition 9 Under the μ-strongly convexity assumption on the surrogates g( | ), the proximal distance iterates satisfy limk→∞‖xk+1 −x‖ = 0. As a consequence, the set of limit points is connected as well as closed. Furthermore, if each limit point is isolated, then the iterates converge to a critical point.

Further progress requires even more structure. Fortunately, what we now pursue applies to generic MM algorithms. We start with the concept of a Fréchet subdifferential (Kruger, 2003). If h() is a function mapping into , then its Fréchet subdifferential at ∈ dom f is the set The set ∂h() is closed, convex, and possibly empty. If h() is convex, then ∂h() reduces to its convex subdifferential. If h() is differentiable, then ∂h() reduces to its ordinary differential. At a local minimum , Fermat’s rule 0 ∈ ∂h(x) holds.

Proposition 10 In an MM algorithm, suppose that h() is coercive, that the surrogates g( | ) are differentiable, and that the algorithm map M() is closed. Then every limit point z of the MM sequence is critical in the sense that 0 ∈ ∂ (−h)(z).

We will also need to invoke Łojasiewicz’s inequality. This deep result depends on some rather arcane algebraic geometry (Bierstone and Milman, 1988; Bochnak et al., 2013). It applies to semialgebraic functions and their more inclusive cousins semianalytic functions and subanalytic functions. For simplicity we focus on semialgebraic functions. The class of semialgebraic subsets of is the smallest class that: A function a : is said to be semialgebraic if its graph is a semialgebraic set of . The class of real-valued semialgebraic functions contains all polynomials p() and all 0/1 indicators of algebraic sets. It is closed under the formation of sums, products, absolute values, reciprocals when a() 6≠ 0, nth roots when a(x) ≥ 0, and maxima max{a(), b()} and minima min{a(), b()}. For our purposes, it is important to note that dist(,S) is a semialgebraic function whenever S is a semialgebraic set.

contains all sets of the form { : q() > 0} for a polynomial q() in p variables,

is closed under the formation of finite unions, finite intersections, and set complementation.

Łojasiewicz’s inequality in its modern form (Bolte et al., 2007) requires a function h() to be closed (lower semicontinuous) and subanalytic with a closed domain. If z is a critical point of h(), then for all x ∈ B(z)∩dom∂h satisfying h() > h() and all in ∂h(). Here the exponent θ ∈ [0,1), the radius r, and the constant c depend on z. This inequality is valid for semialgebraic functions since they are automatically subanalytic. We will apply Łojasiewicz’s inequality to the limit points of an MM algorithm. The next proposition is an elaboration and expansion of known results (Attouch et al., 2010; Bolte et al., 2007; Cui et al., 2018; Kang et al., 2015; Le Thi et al., 2009).

Proposition 11 In an MM algorithm suppose the objective h() is coercive, continuous, and subanalytic and all surrogates g( | ) are continuous, μ-strongly convex, and satisfy the L-smoothness condition on the compact set {x : h() ≤ h(0)}. Then the MM iterates k+1 = argming( | ) converge to a critical point.

The last proposition applies to proximal distance algorithms. The loss f() must be subanalytic and differentiable with a locally Lipschitz gradient. Furthermore, all surrogates should be coercive and μ-strongly convex. Finally, the constraints sets S should be subanalytic. Semialgebraic sets and functions will do. Under these conditions and regardless of how the projected points P() are chosen, the MM iterates are guaranteed to converge to a critical point.

Examples

The following examples highlight the versatility of proximal distance algorithms in a variety of convex and nonconvex settings. Programming details matter in solving these problems. Individual programs are not necessarily long, but care must be exercised in projecting onto constraints, choosing tuning schedules, folding constraints into the domain of the loss, implementing acceleration, and declaring convergence. All of our examples are coded in the Julia programming language. Whenever possible, competing software was run in the Julia environment via the Julia module MathProgBase (Dunning et al., 2017; Lubin and Dunning, 2015). The sparse PCA problem relies on the software of Witten et al. (Witten et al., 2009), which is coded in R. Convergence is tested at iteration k by the two criteria where ϵ1=10−6 and ϵ2=10−4 are typical values. The number of iterations until convergence is about 1000 in most examples. This handicap is offset by the simplicity of each stereotyped update. Our code is available as supplementary material to this paper. Readers are encouraged to try the code and adapt it to their own examples.

Linear Programming

Two different tactics suggest themselves for constructing a proximal distance algorithm. The first tactic rolls the standard affine constraints = into the domain of the loss function vx. The standard nonnegativity requirement ≥ 0 is achieved by penalization. Let x be the current iterate and = ()+ be its projection onto . Derivation of the proximal distance algorithm relies on the Lagrangian One can multiply the corresponding stationarity equation by and solve for the Lagrange multiplier in the form assuming has full row rank. Inserting this value into the stationarity equation gives the MM update where − = ()−1 is the pseudo-inverse of .

The second tactic folds the nonnegativity constraints into the domain of the loss. Let denote the projection of onto the affine constraint set = . Fortunately, the surrogate function splits the parameters. Minimizing one component at a time gives the update x with components The projection can be computed via where − is again the pseudo-inverse of .

Table 1 compares the accelerated versions of these two proximal distance algorithms to two efficient solvers. The first is the open-source Splitting Cone Solver (SCS) (O’Donoghue et al., 2016), which relies on a fast implementation of ADMM. The second is the commercial Gurobi solver, which ships with implementations of both the simplex method and a barrier (interior point) method; in this example, we use its barrier algorithm. The first seven rows of the table summarize linear programs with dense data , , and v. The bottom six rows rely on random sparse matrices with sparsity level 0.01. For dense problems, the proximal distance algorithms start the penalty constant ρ at 1 and double it every 100 iterations. Because we precompute and cache the pseudoinverse − of , the updates (12) and (13) reduce to vector additions and matrix-vector multiplications.

Table 1:

CPU times and optima for linear programming. Here m is the number of constraints, n is the number of variables, PD1 is the proximal distance algorithm over an affine domain, PD2 is the proximal distance algorithm over a nonnegative domain, SCS is the Splitting Cone Solver, and Gurobi is the Gurobi solver. After m = 512 the constraint matrix is initialized to be sparse with sparsity level s = 0.01.

Dimensions		Optima				CPU Times (secs)
m	n	PD1	PD2	SCS	Gurobi	PD1	PD2	SCS	Gurobi
2	4	0.2629	0.2629	0.2629	0.2629	0.0142	0.0010	0.0034	0.0038
4	8	1.0455	1.0457	1.0456	1.0455	0.0212	0.0021	0.0009	0.0011
8	16	2.4513	2.4515	2.4514	2.4513	0.0361	0.0048	0.0018	0.0029
16	32	3.4226	3.4231	3.4225	3.4223	0.0847	0.0104	0.0090	0.0036
32	64	6.2398	6.2407	6.2397	6.2398	0.1428	0.0151	0.0140	0.0055
64	128	14.671	14.674	14.671	14.671	0.2117	0.0282	0.0587	0.0088
128	256	27.116	27.125	27.116	27.116	0.3993	0.0728	0.8436	0.0335
256	512	58.501	58.512	58.494	58.494	0.7426	0.1538	2.5409	0.1954
512	1024	135.35	135.37	135.34	135.34	1.6413	0.5799	5.0648	1.7179
1024	2048	254.50	254.55	254.47	254.48	2.9541	3.2127	3.9433	0.6787
2048	4096	533.29	533.35	533.23	533.23	7.3669	17.318	25.614	5.2475
4096	8192	991.78	991.88	991.67	991.67	30.799	95.974	98.347	46.957
8192	16384	2058.8	2059.1	2058.5	2058.5	316.44	623.42	454.23	400.59

For sparse problems the proximal distance algorithms update ρ by a factor of 1.5 every 50 iterations. To avoid computing large pseudoinverses, we appeal to the LSQR variant of the conjugate gradient method (Paige and Saunders, 1982b,a) to solve the linear systems (11) and (14). The optima of all four methods agree to about 4 digits of accuracy. It is hard to declare an absolute winner in these comparisons. Gurobi and SCS clearly perform better on low-dimensional problems, but the proximal distance algorithms are competitive as dimensions increase. PD1, the proximal distance algorithm over an affine domain, tends to be more accurate than PD2. If high accuracy is not a concern, then the proximal distance algorithms are easily accelerated with a more aggressive update schedule for ρ.

Constrained Least Squares

Constrained least squares programming subsumes constrained quadratic programming. A typical quadratic program involves minimizing the quadratic subject to x ∈ C for a positive definite matrix Q. Quadratic programming can be reformulated as least squares by taking the Cholesky decomposition Q = LL of Q and noting that The constraint x ∈ C applies in both settings. It is particularly advantageous to reframe a quadratic program as a least squares problem when Q is already presented in factored form or when it is nearly singular (Bemporad, 2018). To simplify subsequent notation, we replace by the rectangular matrix and −1 by . The key to solving constrained least squares is to express the proximal distance surrogate as as in equation (6). As noted earlier, in sparse problems the update can be found by a fast stable conjugate gradient solver.

Table 2 compares the performance of the proximal distance algorithm for least squares estimation with probability-simplex constraints to the open source nonlinear interior point solver Ipopt (Wächter and Biegler, 2005, 2006) and the interior point method of Gurobi. Simplex constrained problems arise in hyperspectral imaging (Heylen et al., 2011; Keshava, 2003), portfolio optimization (Markowitz, 1952), and density estimation (Bunea et al., 2010). Test problems were generated by filling an n×p matrix and an n-vector with standard normal deviates. For sparse problems we set the sparsity level of to be 10/p. Our setup ensures that has full rank and that the quadratic program has a solution. For the proximal distance algorithm, we start ρ at 1 and multiply it by 1.5 every 200 iterations. Table 2 suggests that the proximal distance algorithm and the interior point solvers perform equally well on small dense problems. However, in high-dimensional and low-accuracy environments, the proximal distance algorithm provides much better scalability.

Table 2:

CPU times and optima for simplex-constrained least squares. Here , PD is the proximal distance algorithm, IPOPT is the Ipopt solver, and Gurobi is the Gurobi solver. After n = 1024, the predictor matrix A is sparse.

Dimensions		Optima			CPU Times
n	p	PD	IPOPT	Gurobi	PD	IPOPT	Gurobi
16	8	4.1515	4.1515	4.1515	0.0038	0.0044	0.0010
32	16	10.8225	10.8225	10.8225	0.0036	0.0039	0.0010
64	32	29.6218	29.6218	29.6218	0.0079	0.0079	0.0019
128	64	43.2626	43.2626	43.2626	0.0101	0.0078	0.0033
256	128	111.7642	111.7642	111.7642	0.0872	0.0151	0.0136
512	256	231.6455	231.6454	231.6454	0.1119	0.0710	0.0619
1024	512	502.1276	502.1276	502.1276	0.2278	0.4013	0.2415
2048	1024	994.2447	994.2447	994.2447	1.2575	2.3346	1.1682
4096	2048	2056.8381	2056.8381	2056.8381	1.3253	15.2214	7.4971
8192	4096	4103.4611	4103.4611	4103.4611	3.0289	146.1604	49.7411
16384	8192	8295.2136	8295.2136	8295.2136	6.8739	732.1039	412.3612

Closest Kinship Matrix

In genetics studies, kinship is measured by the fraction of genes two individuals share identical by descent. For a given pedigree, the kinship coefficients for all pairs of individuals appear as entries in a symmetric kinship matrix . This matrix possesses three crucial properties: a) it is positive semidefinite, b) its entries are nonnegative, and c) its diagonal entries are unless some pedigree members are inbred. Inbreeding is the exception rather than the rule. Kinship matrices can be estimated empirically from single nucleotide polymorphism (SNP) data, but there is no guarantee that the three highlighted properties are satisfied. Hence, it helpful to project to the nearest qualifying matrix.

This projection problem is best solved by folding the positive semidefinite constraint into the domain of the Frobenius loss function . As we shall see, the alternative of imposing two penalties rather than one is slower and less accurate. Projection onto the constraints implied by conditions b) and c) is trivial. All diagonal entries x of are reset to , and all off-diagonal entries x are reset to max{x,0}. If P() denotes the current projection, then the proximal distance algorithm minimizes the surrogate where c is an irrelevant constant. The minimum is found by extracting the spectral decomposition of and truncating the negative eigenvalues. This gives the update = + in obvious notation. This proximal distance algorithm and its Nesterov acceleration are simple to implement in a numerically oriented language such as Julia. The most onerous part of the calculation is clearly the repeated eigen-decompositions.

Table 3 compares three versions of the proximal distance algorithm to Dykstra’s algorithm (Boyle and Dykstra, 1986). Higham proposed Dykstra’s algorithm for the related problem of finding the closest correlation matrix Higham (2002). In Table 3 algorithm PD1 is the unadorned proximal distance algorithm, PD2 is the accelerated proximal distance, and PD3 is the accelerated proximal distance algorithm with the positive semidefinite constraints folded into the domain of the loss. On this demanding problem, these algorithms are comparable to Dykstra’s algorithm in speed but slightly less accurate. Acceleration of the proximal distance algorithm is effective in reducing both execution time and error. Folding the positive semidefinite constraint into the domain of the loss function leads to further improvements. The data matrices in these trials were populated by standard normal deviates and then symmetrized by averaging opposing off-diagonal entries. In algorithm PD1 we set ρ = max{1.2,222}. In the accelerated versions PD2 and PD3 we started ρ at 1 and multiplied it by 5 every 100 iterations. At the expense of longer compute times, better accuracy can be achieved by all three proximal distance algorithms with a less aggressive update schedule.

Table 3:

CPU times and optima for the closest kinship matrix problem. Here the kinship matrix is n × n, PD1 is the proximal distance algorithm, PD2 is the accelerated proximal distance, PD3 is the accelerated proximal distance algorithm with the positive semidefinite constraints folded into the domain of the loss, and Dykstra is Dykstra’s adaptation of alternating projections. All times are in seconds.

Size	PD1		PD2		PD3		Dykstra
n	Loss	Time	Loss	Time	Loss	Time	Loss	Time
2	1.64	0.36	1.64	0.01	1.64	0.01	1.64	0.00
4	2.86	0.10	2.86	0.01	2.86	0.01	2.86	0.00
8	18.77	0.21	18.78	0.03	18.78	0.03	18.78	0.00
16	45.10	0.84	45.12	0.18	45.12	0.12	45.12	0.02
32	169.58	4.36	169.70	0.61	169.70	0.52	169.70	0.37
64	837.85	16.77	838.44	2.90	838.43	2.63	838.42	4.32
128	3276.41	91.94	3279.44	18.00	3279.25	14.83	3279.23	19.73
256	14029.07	403.59	14045.30	89.58	14043.59	64.89	14043.46	72.79

Projection onto a Second-Order Cone Constraint

Second-order cone programming is one of the unifying themes of convex analysis (Alizadeh and Goldfarb, 2003; Lobo et al., 1998). It revolves around conic constraints of the form { : ‖Au + b‖ ≤ c + d}. Projection of a vector onto such a constraint is facilitated by parameter splitting. In this setting parameter splitting introduces a vector , a scalar r, and the two affine constraints = + and r = + d. The conic constraint then reduces to the Lorentz cone constraint ‖w‖ ≤ r, for which projection is straightforward (Boyd and Vandenberghe, 2009). If we concatenate the parameters into the single vector and define L = { : ‖w‖ ≤ r} and M = { : = + and r = + d}, then we can rephrase the problem as minimizing subject to ∈ L ∩ M. This is a fairly typical set projection problem except that the and r components of y are missing in the loss function.

Taking a cue from Example 5.1, we incorporate the affine constraints in the domain of the objective function. If we represent projection onto L by then the Lagrangian generated by the proximal distance algorithm amounts to This gives rise to a system of three stationarity equations Solving for the multipliers and θ in equations (16) and (17) and substituting their values in equation (15) yield This leads to the MM update The updates = + and r = + d follow from the constraints.

Table 4 compares the proximal distance algorithm to SCS and Gurobi. Echoing previous examples, we tailor the update schedule for ρ differently for dense and sparse problems. Dense problems converge quickly and accurately when we set ρ0 = 1 and double ρ every 100 iterations. Sparse problems require a greater range and faster updates of ρ, so we set ρ0 = 0.01 and then multiply ρ by 2.5 every 10 iterations. For dense problems, it is clearly advantageous to cache the spectral decomposition of + as suggested in Example 5.2. In this regime, the proximal distance algorithm is as accurate as Gurobi and nearly as fast. SCS is comparable to Gurobi in speed but notably less accurate.

Table 4:

CPU times and optima for the second-order cone projection. Here m is the number of constraints, n is the number of variables, PD is the accelerated proximal distance algorithm, SCS is the Splitting Cone Solver, and Gurobi is the Gurobi solver. After m = 512 the constraint matrix is initialized with sparsity level 0.01.

Dimensions		Optima			CPU Seconds
m	n	PD	SCS	Gurobi	PD	SCS	Gurobi
2	4	0.10598	0.10607	0.10598	0.0043	0.0103	0.0026
4	8	0.00000	0.00000	0.00000	0.0003	0.0009	0.0022
8	16	0.88988	0.88991	0.88988	0.0557	0.0011	0.0027
16	32	2.16514	2.16520	2.16514	0.0725	0.0012	0.0040
32	64	3.03855	3.03864	3.03853	0.0952	0.0019	0.0094
64	128	4.86894	4.86962	4.86895	0.1225	0.0065	0.0403
128	256	10.5863	10.5843	10.5863	0.1975	0.0810	0.0868
256	512	31.1039	31.0965	31.1039	0.5463	0.3995	0.3405
512	1024	27.0483	27.0475	27.0483	3.7667	1.6692	2.0189
1024	2048	1.45578	1.45569	1.45569	0.5352	0.3691	1.5489
2048	4096	2.22936	2.22930	2.22921	1.0845	2.4531	5.5521
4096	8192	1.72306	1.72202	1.72209	3.1404	17.272	15.204
8192	16384	5.36191	5.36116	5.36144	13.979	133.25	88.024

With a large sparse constraint matrix , extraction of its spectral decomposition becomes prohibitive. If we let = (ρ−1 c), then we must solve a linear system of equations defined by the Gramian matrix = . There are three reasonable options for solving this system. The first relies on computing and caching a sparse Cholesky decomposition of . The second computes the QR decomposition of the sparse matrix . The R part of the QR decomposition coincides with the Cholesky factor. Unfortunately, every time ρ changes, the Cholesky or QR decomposition must be redone. The third option is the conjugate gradient algorithm. In our experience the QR decomposition offers superior stability and accuracy. When is very sparse, the QR decomposition is often much faster than the Cholesky decomposition because it avoids forming the dense matrix . Even when only 5% of the entries of are nonzero, 90% of the entries of can be nonzero. If exquisite accuracy is not a concern, then the conjugate gradient method provides the fastest update. Table 4 reflects this choice.

Copositive Matrices

A symmetric matrix is copositive if its associated quadratic form is nonnegative for all x ≥ 0. Copositive matrices find applications in numerous branches of the mathematical sciences (Berman and Plemmons, 1994). All positive semidefinite matrices and all matrices with nonnegative entries are copositive. The variational index is one key to understanding copositive matrices (Hiriart-Urruty and Seeger, 2010). The constraint set S is the intersection of the unit sphere and the nonnegative cone . Projection of an external point onto S splits into three cases. When all components of are negative, then P() = e, where y is the least negative component of , and is the standard unit vector along coordinate direction i. The origin 0 is equidistant from all points of S. If any component of is positive, then the projection is constructed by setting the negative components of equal to 0, and standardizing the truncated version of to have Euclidean norm 1.

As a test case for the proximal distance algorithm, consider the Horn matrix (Hall and Newman, 1963) The value μ(M) = 0 is attained for the vectors , , and equivalent vectors with their entries permuted. Matrices in higher dimensions with the same Horn pattern of 1’s and −1’s are copositive as well (Johnson and Reams, 2008). A Horn matrix of odd dimension cannot be written as a positive semidefinite matrix, a nonnegative matrix, or a sum of two such matrices.

The proximal distance algorithm minimizes the criterion and generates the updates It takes a gentle tuning schedule to get decent results. The choice ρ = 1.2 converges in 600 to 700 iterations from random starting points and reliably yields objective values below 10−5 for Horn matrices. The computational burden per iteration is significantly eased by exploiting the cached spectral decomposition of . Table 5 compares the performance of the proximal distance algorithm to the Mosek solver on a range of Horn matrices. Mosek uses semidefinite programming to decide whether can be decomposed into a sum of a positive semidefinite matrix and a nonnegative matrix. If not, Mosek declares the problem infeasible. Nesterov acceleration improves the final loss for the proximal distance algorithm, but it does not decrease overall computing time.

Table 5:

CPU times (seconds) and optima for approximating the Horn variational index of a Horn matrix. Here n is the size of Horn matrix, PD is the proximal distance algorithm, aPD is the accelerated proximal distance algorithm, and Mosek is the Mosek solver.

Dimension	Optima			CPU Seconds
n	PD	aPD	Mosek	PD	aPD	Mosek
4	0.000000	0.000000	feasible	0.5555	0.0124	2.7744
5	0.000000	0.000000	infeasible	0.0039	0.0086	0.0276
8	0.000021	0.000000	feasible	0.0059	0.0083	0.0050
9	0.000045	0.000000	infeasible	0.0055	0.0072	0.0082
16	0.000377	0.000001	feasible	0.0204	0.0237	0.0185
17	0.000441	0.000001	infeasible	0.0204	0.0378	0.0175
32	0.001610	0.000007	feasible	0.0288	0.0288	0.1211
33	0.002357	0.000009	infeasible	0.0242	0.0346	0.1294
64	0.054195	0.000026	feasible	0.0415	0.0494	3.6284
65	0.006985	0.000026	infeasible	0.0431	0.0551	2.7862

Testing for copositivity is challenging because neither the loss function nor the constraint set is convex. The proximal distance algorithm offers a fast screening device for checking whether a matrix is copositive. On random 1000×1000 symmetric matrices , the method invariably returns a negative index in less than two seconds of computing time. Because the vast majority of symmetric matrices are not copositive, accurate estimation of the minimum is not required. Table 6 summarizes a few random trials with lower-dimensional symmetric matrices. In higher dimensions, Mosek becomes non-competitive, and Nesterov acceleration is of dubious value.

Table 6:

CPU times and optima for testing the copositivity of random symmetric matrices. Here n is the size of matrix, PD is the proximal distance algorithm, aPD is the accelerated proximal distance algorithm, and Mosek is the Mosek solver.

Dimension	Optima			CPU Seconds
n	PD	aPD	Mosek	PD	aPD	Mosek
4	−0.391552	−0.391561	infeasible	0.0029	0.0031	0.0024
8	−0.911140	−2.050316	infeasible	0.0037	0.0044	0.0045
16	−1.680697	−1.680930	infeasible	0.0199	0.0272	0.0062
32	−2.334520	−2.510781	infeasible	0.0261	0.0242	0.0441
64	−3.821927	−3.628060	infeasible	0.0393	0.0437	0.6559
128	−5.473609	−5.475879	infeasible	0.0792	0.0798	38.3919
256	−7.956365	−7.551814	infeasible	0.1632	0.1797	456.1500

Linear Complementarity Problem

The linear complementarity problem (Murty and Yu, 1988) consists of finding vectors x and y with nonnegative components such that = 0 and = x + for a given square matrix A and vector b. The natural loss function is . To project a vector pair (,) onto the nonconvex constraint set, one considers each component pair (u, v) in turn. If u ≥ max{v,0}, then the nearest pair (,) has components (x, y) = (u,0). If v ≥ max{u,0}, then the nearest pair has components (x, y) = (0, v). Otherwise, (x, y) = (0,0). At each iteration the proximal distance algorithm minimizes the criterion where (, ) is the projection of (, ) onto the constraint set. The stationarity equations become Substituting the value of from the second equation into the first equation leads to the updates The linear system (19) can be solved in low to moderate dimensions by computing and caching the spectral decomposition of and in high dimensions by the conjugate gradient method. Table 7 compares the performance of the proximal gradient algorithm to the Gurobi solver on some randomly generated problems.

Table 7:

CPU times (seconds) and optima for the linear complementarity problem with randomly generated data. Here n is the size of matrix, PD is the accelerated proximal distance algorithm, and Gurobi is the Gurobi solver.

Dimension	Optima		CPU Seconds
n	PD	Mosek	PD	Mosek
4	0.000000	0.000000	0.0230	0.0266
8	0.000000	0.000000	0.0062	0.0079
16	0.000000	0.000000	0.0269	0.0052
32	0.000000	0.000000	0.0996	0.4303
64	0.000074	0.000000	2.6846	360.5183

Sparse Principal Components Analysis

Let be an n × p data matrix gathered on n cases and p predictors. Assume the columns of are centered to have mean 0. Principal component analysis (PCA) (Hotelling, 1933; Pearson, 1901) operates on the sample covariance matrix . Here we formulate a proximal distance algorithm for sparse PCA (SPCA), which has attracted substantial interest in the machine learning community (Berthet and Rigollet, 2013b,a; D’Aspremont et al., 2007; Johnstone and Lu, 2009; Journée et al., 2010; Witten et al., 2009; Zou et al., 2006). According to a result of Ky Fan (Fan, 1949), the first q principal components (PCs) 1,…, can be extracted by maximizing the function tr() subject to the matrix constraint = , where is the ith column of the p×q matrix . This constraint set is called a Stiefel manifold. One can impose sparsity by insisting that any given column have at most r nonzero entries. Alternatively, one can require the entire matrix U to have at most r nonzero entries. The latter choice permits sparsity to be distributed non-uniformly across columns.

Extraction of sparse PCs is difficult for three reasons. First, the Stiefel manifold and both sparsity sets are nonconvex. Second, the objective function is concave rather than convex. Third, there is no simple formula or algorithm for projecting onto the intersection of the two constraint sets. Fortunately, it is straightforward to project onto each separately. Let denote the projection of onto the Stiefel manifold. It is well known that can be calculated by extracting a partial singular value decomposition = of and setting P() = (Golub and Van Loan, 2012). Here and are orthogonal matrices of dimension p×q and q×q, respectively, and Σ is a diagonal matrix of dimension q × q. Let denote the projection of onto the sparsity set Because operates column by column, it suffices to project each column vector to sparsity. This entails nothing more than sorting the entries of by magnitude, saving the r largest, and sending the remaining p−r entries to 0. If the entire matrix must have at most r nonzero entries, then can be treated as a concatenated vector during projection.

The key to a good algorithm is to incorporate the Stiefel constraints into the domain of the objective function (Kiers, 1990; Kiers and ten Berge, 1992) and the sparsity constraints into the distance penalty. Thus, we propose decreasing the criterion at each iteration subject to the Stiefel constraints. The loss can be majorized via because is positive semidefinite. The penalty is majorized by up to an irrelevant constant c since the squared Frobenius norm satisfies the relation on the Stiefel manifold. It now follows that f() is majorized by up to an irrelevant constant. Accordingly, the Stiefel projection provides the next MM iterate.

Figures 1 and 2 compare the proximal distance algorithm to the SPC function from the R package PMA (Witten et al., 2009). The breast cancer data from PMA provide the data matrix . The data consist of p = 19,672 RNA measurements on n = 89 patients. The two figures show computation times and the proportion of variance explained (PVE) by the p × q loading matrix . For sparse PCA, PVE is defined as , where = ()−1 (Shen and Huang, 2008). When the loading vectors of are orthogonal, this criterion reduces to the familiar definition tr()/tr() of PVE for ordinary PCA. The proximal distance algorithm enforces either matrix-wise or column-wise sparsity. In contrast, SPC enforces only column-wise sparsity via the constraint ‖u‖1 ≤ c for each column of . We take c = 8. The number of nonzeroes per loading vector output by SPC dictates the sparsity level for the column-wise version of the proximal distance algorithm. Summing these counts across all columns dictates the sparsity level for the matrix version of the proximal distance algorithm.

Figure 1:

Proportion of variance explained by q PCs for each algorithm. Here PD1 is the accelerated proximal distance algorithm enforcing matrix sparsity, PD2 is the accelerated proximal distance algorithm enforcing column-wise sparsity, and SPC is the orthogonal sparse PCA method from PMA.

Figure 2:

Computation times for q PCs for each algorithm. Here PD1 is the accelerated proximal distance algorithm enforcing matrix sparsity, PD2 is the accelerated proximal distance algorithm enforcing column-wise sparsity, and SPC is the orthogonal sparse PCA method from PMA.

Figures 1 and 2 demonstrate the superior PVE and computational speed of both proximal distance algorithms versus SPC. The type of projection does not appear to affect the computational performance of the proximal distance algorithm, as both versions scale equally well with q. However, the matrix projection, which permits the algorithm to more freely assign nonzeroes to the loadings, attains better PVE than the more restrictive column-wise projection. For both variants of the proximal distance algorithm, Nesterov acceleration improves both fitting accuracy and computational speed, especially as the number of PCs q increases.

Discussion

The proximal distance algorithm applies to a host of problems. In addition to the linear and quadratic programming examples considered here, our previous paper (Lange and Keys, 2014) derives and tests algorithms for binary piecewise-linear programming, ℓ0 regression, matrix completion (Cai et al., 2010; Candès and Tao, 2010; Chen et al., 2012; Mazumder et al., 2010), and sparse precision matrix estimation (Friedman et al., 2008). Other potential applications immediately come to mind. An integer linear program in standard form can be expressed as minimizing cx subject to + = , ≥ 0, and . The latter two constraints can be combined in a single constraint for which projection is trivial. The affine constraints should be folded into the domain of the objective. Integer programming is NP hard, so that the proximal distance algorithm just sketched is merely heuristic. Integer linear programming includes traditional NP hard problems such as the traveling salesman problem, the vertex cover problem, set packing, and Boolean satisfiability. It will be interesting to see if the proximal distance principle is competitive in meeting these challenges. Our experience with the closest lattice point problem (Agrell et al., 2002) and the eight queens problem suggests that the proximal distance algorithm can be too greedy for hard combinatorial problems. The nonconvex problems solved in this paper are in some vague sense easy combinatorial problems.

The behavior of a proximal distance algorithm depends critically on a sensible tuning schedule for increasing ρ. Starting ρ too high puts too much stress on satisfying the constraints. Incrementing ρ too quickly causes the algorithm to veer off the solution path guaranteed by the penalty method. Given the chance of roundoff error even with double precision arithmetic, it is unwise to take ρ all the way to ∞. Trial and error can help in deciding whether a given class of problems will benefit from an aggressive update schedule and strict or loose convergence criteria. In problems with little curvature such as linear programming, more conservative updates are probably prudent. The linear programming, closest kinship matrix, and SPCA problems document the value of folding constraints into the domain of the loss. In the same spirit it is wise to minimize the number of constraints. A single penalty for projecting onto the intersection of two constraint sets is almost always preferable to the sum of two penalties for their separate projections. Exceptions to this rule obviously occur when projection onto the intersection is intractable. The integer linear programming problem mentioned previously illustrates these ideas.

Our earlier proximal distance algorithms ignored acceleration. In many cases the solutions produced had very low accuracy. The realization that convex proximal distance algorithms can be phrased as proximal gradient algorithms convinced us to try Nesterov acceleration. We now do this routinely on the subproblems with ρ fixed. This typically forces tighter path following and a reduction in overall computing times. Our examples generally bear out the contention that Nesterov acceleration is useful in nonconvex problems (Ghadimi and Lan, 2015). It is noteworthy that the value of acceleration often lies in improving the quality of a solution as much as it does in increasing the rate of convergence. Of course, acceleration cannot prevent convergence to an inferior local minimum.

On both convex and nonconvex problems, proximal distance algorithms enjoy global convergence guarantees. On nonconvex problems, one must confine attention to subanalytic sets and subanalytic functions. This minor restriction is not a handicap in practice. Determining local convergence rates is a more vexing issue. For convex problems, we review existing theory for a fixed penalty constant ρ. The classical results buttress an O(ρk−1) sublinear rate for general convex problems. Better results require restrictive smoothness assumptions on both the objective function and the constraint sets. For instance, when f() is L-smooth and strongly convex, linear convergence can be demonstrated. When f() equals a difference of convex functions, proximal distance algorithms reduce to concave-convex programming. Le Thi et al. (2009) attack convergence in this setting.

We hope readers will sense the potential of the proximal distance principle. This simple idea offers insight into many existing algorithms and a straightforward path in devising new ones. Effective proximal and projection operators usually spell the difference between success and failure. The number and variety of such operators is expanding quickly as the field of optimization relinquishes it fixation on convexity. The current paper research leaves many open questions about tuning schedules, rates of convergence, and acceleration in the face of nonconvexity. We welcome the contributions of other mathematical scientists in unraveling these mysteries and in inventing new proximal distance algorithms.