An Inexact Augmented Lagrangian Framework for Nonconvex Optimization with Nonlinear Constraints
Mehmet Fatih Sahin Armin Eftekhari Ahmet Alacaoglu Fabian Latorre Volkan Cevher


Unifying Thermodynamic Uncertainty Relations
    Jean-Charles Delvenne
    Received: date / Accepted: date
============================================









We propose a practical inexact augmented Lagrangian method (iALM) for nonconvex problems with nonlinear constraints. 

We characterize the total computational complexity of our method subject to a verifiable geometric condition, which is closely related to the Polyak-Lojasiewicz and Mangasarian-Fromowitz conditions. 



In particular, when a first-order solver is used for the inner iterates, we prove that iALM  finds a first-order stationary point with 𝒪̃(1/ϵ^4) calls to the first-order oracle. If, in addition, the problem is smooth and a second-order solver is used for the inner iterates, iALM  finds a second-order stationary point with 𝒪̃(1/ϵ^5) calls to the second-order oracle, which matches the known theoretical complexity result in the literature. 



We also provide strong numerical evidence on large-scale machine learning problems, including the Burer-Monteiro factorization of  semidefinite programs, and a novel nonconvex relaxation of the standard basis pursuit template. For these examples, we also show how to verify our geometric condition.


 This paper was updated on April 21, 2022, due to an error in the proof of Corollary (<ref>). The gradient complexity for obtaining a first order stationary point is now 𝒪̃( 1/ϵ^4) instead of 𝒪̃( 1/ϵ^3). The error was due to an order mistake while using the complexity result of APGM from <cit.>.









§ INTRODUCTION



We study the  nonconvex optimization problem


    x∈ℝ^dmin   f(x)+g(x)    A(x) = 0,


where f:ℝ^d→ℝ is  a continuously-differentiable nonconvex function and A:ℝ^d→ℝ^m is a  nonlinear operator. 


We assume that g:ℝ^d→ℝ∪{∞} is a proximal-friendly convex function <cit.>. 





A host of problems in computer science <cit.>, machine learning <cit.>, and signal processing <cit.> naturally fall under the template (<ref>),  including max-cut, clustering, generalized eigenvalue decomposition, as well as the quadratic assignment problem (QAP) <cit.>.

To solve (<ref>), we 

propose an intuitive and easy-to-implement augmented Lagrangian algorithm, and provide its total iteration complexity under an interpretable geometric condition. Before we elaborate on the results, let us first motivate (<ref>) with an application to semidefinite programming (SDP): 




  
Vignette: Burer-Monteiro splitting.
A powerful convex relaxation for max-cut, clustering, and many others is provided by the trace-constrained SDP 

    X∈𝕊^d × dmin⟨ C, X ⟩  B(X) = b,   tr(X)≤α,    X ≽ 0,


where C∈^d× d, X is a positive semidefinite d× d matrix, 

and B: 𝕊^d× d→ℝ^m is a linear operator. If the unique-games conjecture is true, the SDP (<ref>) obtains the best possible approximation for the underlying discrete problem <cit.>.

Since d is often large, many first- and second-order methods for solving such SDP's are immediately ruled out, not only due to their high computational complexity, but also due to their storage requirements, which are 𝒪(d^2). 

A contemporary challenge in optimization is therefore to solve SDPs using little space and in a scalable fashion. The recent homotopy conditional gradient method, which is based on linear minimization oracles (LMOs), can solve (<ref>) in a small space via sketching <cit.>. However, such LMO-based methods are extremely slow in obtaining accurate solutions. 




A different approach for solving (<ref>), dating back to <cit.>, is the so-called Burer-Monteiro (BM) factorization X=UU^⊤, where U∈ℝ^d× r and r is selected according to the guidelines in <cit.>, which is tight <cit.>. 


The BM factorization leads to the following  nonconvex problem in the template (<ref>):

    U∈ℝ^d × rmin⟨ C, UU^⊤⟩ B(UU^⊤) = b,    U_F^2 ≤α,


The BM factorization does not introduce any extraneous local minima <cit.>. Moreover, <cit.> establishes the connection between the local minimizers of the factorized problem (<ref>) and the global minimizers for (<ref>). 










To solve (<ref>), the inexact Augmented Lagrangian method (iALM) is widely used <cit.>, due to its cheap per iteration cost and its empirical success. 

Every (outer) iteration of iALM calls a solver to solve an intermediate augmented Lagrangian subproblem to near stationarity. The choices include first-order methods, such as the proximal gradient descent <cit.>, or second-order methods, such as the trust region method and BFGS <cit.>.[BFGS is in fact a quasi-Newton method that emulates second-order information.]

Unlike its convex counterpart <cit.>, the convergence rate and the complexity of iALM for  (<ref>) are not well-understood, see Section <ref> for a review of the related literature. Indeed, addressing this important theoretical gap is one of the contributions of our work. In addition,



▹ 
We derive the convergence rate of iALM 
to first-order optimality for solving (<ref>)  or second-order optimality for solving (<ref>) with g=0, and find the total iteration complexity of iALM using different solvers for the augmented Lagrangian subproblems. 
We provide an extensive comparison with the existing complexity results in optimization, see Section <ref>.  


▹ Our iALM framework is future-proof in the sense that different subsolvers can be substituted. 

▹ We propose a geometric condition that simplifies the algorithmic analysis for iALM, and clarify its connection to well-known Polyak-Lojasiewicz <cit.> and Mangasarian-Fromovitz <cit.> conditions. 

We also verify this condition for key problems in Appendices <ref> and <ref>. 

















§ PRELIMINARIES 





  
Notation.
We use the notation ⟨· ,·⟩ and · for the standard inner product and the norm on ^d. For matrices, · and ·_F denote the spectral and the Frobenius norms, respectively.  

For the convex function g:^d→, the subdifferential set at x∈^d is denoted by ∂ g(x) and  we will occasionally use the notation ∂ g(x)/=̱{ z/:̱ z∈∂ g(x)}. 
When presenting iteration complexity results, we often use O(·) which suppresses the logarithmic dependencies.







We denote δ_𝒳:^d→ as the indicator function  of a set 𝒳⊂^d. 










The distance function from a point x to 𝒳 is denoted by (x,𝒳) = min_z∈𝒳x-z. 
For integers k_0 ≤ k_1, we use the  notation [k_0:k_1]={k_0,…,k_1}. For an operator A:^d→^m with components {A_i}_i=1^m,  DA(x) ∈ℝ^m× d denotes the Jacobian of A, where the ith row of DA(x) is the  vector ∇ A_i(x) ∈^d. 






  
Smoothness.
We assume smooth f:^d→ and A:^d→^m; i.e., there exist λ_f,λ_A≥ 0 s.t.



    ∇ f(x) - ∇ f(x')≤λ_f x-x',    DA(x) - DA(x') ≤λ_A x-x',
      ∀ x, x' ∈^d .






  
Augmented Lagrangian method (ALM).
ALM is a classical algorithm, which first appeared in <cit.> and extensively studied afterwards in <cit.>.
For solving (<ref>), ALM suggests solving the problem



    min_xmax_y   ℒ_β(x,y) + g(x),



where, for penalty weight >̱0,  ℒ_$̱ is the corresponding augmented Lagrangian, defined as 

    ℒ_β(x,y) := f(x) + ⟨ A(x), y ⟩ +  β/2A(x)^2.



The minimax formulation in (<ref>) naturally suggests the following algorithm for solving (<ref>):

    x_k+1∈x  ℒ_β(x,y_k)+g(x),

    y_k+1   = y_k+_k A(x_k+1),

where the  dual step sizes  are denoted as{_k}_k.
However, computingx_k+1above requires solving the nonconvex problem (<ref>) to optimality, which is typically intractable. Instead, it is  often easier to find an approximate first- or second-order stationary point of (<ref>).


Hence, we argue that by gradually improving the stationarity precision and  increasing the penalty weight$̱ above, we can reach a stationary point of the main problem in (<ref>), as detailed  in Section <ref>. 



  
Optimality conditions.
First-order necessary optimality conditions for (<ref>) are well-studied. Indeed, x∈^d is a first-order stationary point of (<ref>) if there exists y∈^m such that 











    -∇_x ℒ_β(x,y) ∈∂ g(x),      A(x) = 0,

which is in turn the necessary optimality condition for (<ref>). 





















Inspired by this, we say that x is an (ϵ_f,)̱ first-order stationary point of (<ref>) if there exists a y ∈^m such that 

    (-∇_x ℒ_β(x,y), ∂ g(x)) ≤ϵ_f,      A(x) ≤ϵ_f,

for ϵ_f≥ 0. 
In light of (<ref>), a metric for evaluating the stationarity of a pair (x,y)∈^d×^m is 

    (-∇_x ℒ_β(x,y), ∂ g(x) ) + A(x) ,

which we use as the first-order stopping criterion. 

As an example, for a convex set 𝒳⊂^d, suppose that g = δ_𝒳 is the indicator function on 𝒳.
Let also T_𝒳(x) ⊆^d denote the tangent cone to 𝒳 at x, and with P_T_𝒳(x):^d→^d we denote the orthogonal projection onto this tangent cone. Then, for u∈^d, it is not difficult to verify that 

    (u, ∂ g(x) ) =  P_T_𝒳(x)(u) .
 

When g = 0, a first-order stationary point x∈^d of (<ref>)   is also second-order stationary if 



    λ_min(∇ _xxℒ_β(x,y))≥ 0,



where ∇_xxℒ_$̱ is the Hessian ofℒ_$̱ with respect to x, and λ_min(·) returns the smallest eigenvalue of its argument.
Analogously, x is an (ϵ_f, ϵ_s,)̱ second-order stationary point if, in addition to (<ref>), it holds that  



    λ_min(∇ _xxℒ_β(x,y)) ≥ -ϵ_s,

for  ϵ_s≥ 0. 


Naturally, for second-order stationarity, we  use λ_min(∇ _xxℒ_β(x,y)) as the stopping criterion.































  
Smoothness lemma. This next result controls the smoothness of Ł_(̱·,y) for a fixed y. The proof is standard but nevertheless is included in Appendix <ref> for completeness. 

 For fixed y∈^m and ρ,ρ'≥ 0, it holds that 
 
    ∇_x ℒ_β(x, y)-  ∇_x ℒ_β(x', y) ≤λ_x-x',

 for every x,x' ∈{ x”: x”≤ρ, A(x”) ≤ρ'}, where
 
    λ_β≤λ_f + √(m)λ_A y +  (√(m)λ_Aρ'  +  d λ'^2_A )=̱: λ_f + √(m)λ_A y +   λ”(A,ρ,ρ') .̱

Above, λ_f,λ_A were defined in (<ref>) and 

    λ'_A := max_x≤ρDA(x).











§ ALGORITHM 


To solve the equivalent formulation of (<ref>)  presented in (<ref>), we  propose the inexact ALM (iALM), detailed in Algorithm <ref>. 
At the k^th iteration, Step 2 of Algorithm <ref> calls a solver that finds an approximate stationary point of the augmented Lagrangian Ł__̱k(·,y_k) with the accuracy of ϵ_k+1, and this accuracy gradually increases in a controlled fashion. 
The increasing sequence of penalty weights {_̱k}_k and the dual update (Steps 4 and 5) are responsible for continuously enforcing the constraints in (<ref>). The appropriate choice for {_̱k}_k will be specified in Corrollary  Sections  <ref> and <ref>. 


The particular choice of the dual step sizes {_k}_k in Algorithm <ref> ensures that the dual variable y_k remains bounded.













§ CONVERGENCE RATE 

This section presents the total iteration complexity of Algorithm <ref> for finding first and second-order stationary points of problem (<ref>). 
All the proofs are deferred to Appendix <ref>.
Theorem <ref> characterizes the convergence rate of Algorithm <ref> for finding stationary points in the number of outer iterations.
(convergence rate)
  For integers 2 ≤ k_0≤ k_1, consider the interval K=[k_0:k_1], and let {x_k}_k∈ K be the output sequence of Algorithm <ref> on the interval K.[The choice of k_1 = ∞ is valid here too.]  Let also ρ:=sup_k∈ [K]x_k.[If necessary, to ensure that ρ<∞, one can add a small factor of x^2 to ℒ_ in (<ref>). Then it is easy to verify that the iterates of Algorithm <ref> remain bounded, provided that the initial penalty weight β_0 is large enough, sup_x ∇ f(x)/x< ∞,  sup_x A(x) < ∞, and sup_x D A(x) <∞.  ] Suppose that f and A satisfy (<ref>) and  let 

    λ'_f = max_x≤ρ∇ f(x),    λ'_A = max_x≤ρDA(x),

be the (restricted) Lipschitz constants of f and A, respectively. 

With ν>0, assume that 
 
    νA(x_k)   ≤( -DA(x_k)^⊤ A(x_k) , ∂ g(x_k)/_̱k-1),

for every k∈ K. We consider two cases:



  * If a first-order solver is used in Step 2, then x_k is an (ϵ_k,f,_̱k) first-order stationary point of (<ref>) with 

    ϵ_k,f    =  1/β_k-1(2(λ'_f+λ'_A y_max) (1+λ_A' σ_k)/ν+1) =: Q(f,g,A,_1)/β_k-1,

for every k∈ K, where y_max(x_1,y_0,_1):=y_0+ c  A(x_1).

  * If a second-order solver is used in Step 2, then x_k is an (ϵ_k,f, ϵ_k,s,_̱k) second-order stationary point of (<ref>) with ϵ_k,s specified above and with 

    ϵ_k,s   = ϵ_k-1 + σ_k √(m)λ_A  2λ'_f +2 λ'_A y_max/ν_̱k-1 = ν + σ_k √(m)λ_A 2λ'_f +2 λ'_A y_max/ν_̱k-1  =: Q'(f,g,A,_1)/β_k-1.










Theorem <ref> states that Algorithm <ref> converges to a (first- or second-) order stationary point of (<ref>) at the rate of 1/_̱k, further specified in 
Corollary <ref> and Corollary <ref>.

A few remarks are in order about Theorem <ref>.




  
Regularity.  The key geometric condition in Theorem <ref> is (<ref>) which, broadly speaking, ensures that the primal updates of Algorithm <ref> reduce the feasibility gap as the penalty weight _̱k grows. We will verify this condition for several examples in Appendices <ref> and <ref>.


This condition in (<ref>) is closely related to those in the existing literature. In the special case where g=0 in (<ref>), (<ref>) reduces to;

    DA(x)^⊤ A(x)≥νA(x).

Polyak-Lojasiewicz (PL) condition <cit.>. Consider the  problem with λ_f̃-smooth objective,

    min_x∈^df̃(x).

f̃(x) satisfies the PL inequality if the following holds for some μ > 0,

    1/2∇f̃(x) ^2 ≥μ (f̃(x) - f̃^*),  ∀ x
    PL inequality

This inequality implies that gradient is growing faster than a quadratic as we move away from the optimal. Assuming that the feasible set { x: A(x) = 0} is non-empty, it is easy to verify that <ref> is equivalent to 
 the PL condition for minimizing f̃(x) = 1/2A(x)^2 with ν = √(2 μ)  <cit.>. 

PL condition itself is a special case of Kurdyka-Lojasiewicz with θ = 1/2, see <cit.>. When g=0, it is also easy to see that (<ref>) is  weaker than the Mangasarian-Fromovitz (MF) condition in nonlinear optimization <cit.>. 
Moreover, when g is the indicator on a convex set, (<ref>) is a consequence of the basic constraint qualification in <cit.>, which itself generalizes the MF condition to the case when g is an indicator function of a convex set. 


We may think of  (<ref>) as a local condition, which should hold within a neighborhood of the constraint set {x:A(x)=0} rather than everywhere in ℝ^d. 
Indeed, the iteration count k appears in (<ref>) to reflect this local nature of the condition. 
Similar kind of arguments on the regularity condition also appear in <cit.>. 
There is also a constant complexity algorithm in <cit.> to reach so-called “information zone”, which supplements Theorem <ref>. 





 































  
Penalty method.
A classical algorithm to solve (<ref>) is the penalty method, which is characterized by the absence of the dual variable (y=0) in (<ref>). Indeed, ALM can be interpreted as an adaptive penalty or smoothing method with a variable center  determined by the dual variable. It is worth noting that, with the same proof technique, one can establish the same convergence rate of Theorem <ref> for the penalty method. However, while both methods have the same convergence rate in theory, we ignore the uncompetitive penalty method since it is significantly outperformed by iALM in practice. 




  
Computational complexity. Theorem <ref> specifies the number of (outer) iterations that Algorithm <ref> requires to reach a near-stationary point of problem (<ref>) with a  prescribed precision and, in particular, specifies the number of calls made to the solver in Step 2.  In this sense, Theorem <ref> does not fully capture the computational complexity of Algorithm <ref>, as it does not take into account the computational cost of the solver in Step 2. 

To better understand the total iteration complexity of Algorithm <ref>, we consider two scenarios in the following. In the first scenario, we take the solver in Step 2 to be the Accelerated Proximal Gradient Method (APGM), a well-known first-order algorithm <cit.>. In the second scenario, we will use the second-order trust region method developed in <cit.>.
We have the following two corollaries showing the total complexity of our algorithm to reach first and second-order stationary points. 
Appendix <ref> contains the proofs and more detailed discussion for the complexity results.




For b>1, let β_k =b^k for every k. If we use APGM from <cit.> for Step 2 of Algorithm <ref>, the algorithm finds an (ϵ_f,_̱k) first-order stationary point of (<ref>), 
 after T calls to the first-order oracle, where


    T = 𝒪( Q^3 ρ^2/ϵ^4log_b( Q/ϵ)) =   𝒪̃( Q^3ρ^2/ϵ^4).





For Algorithm <ref> to reach a near-stationary point with an accuracy of ϵ_f in the sense of (<ref>) and with the lowest computational cost, we therefore need to perform only one iteration of Algorithm <ref>, with _̱1 specified as a function of ϵ_f by (<ref>) in Theorem <ref>. In general, however, the constants in (<ref>) are unknown and this approach is thus not feasible. Instead, the homotopy approach taken by Algorithm <ref> ensures achieving the desired accuracy by gradually increasing the penalty weight. This homotopy approach increases the computational cost of Algorithm <ref> only by a  factor logarithmic in the ϵ_f, as detailed in the proof of Corollary <ref>.


For b>1, let β_k =b^k for every k.
We assume that 

    ℒ_β(x_1, y) - min_xℒ_β(x, y) ≤ L_u,    ∀β.

If we use the trust region method from <cit.> for Step 2 of Algorithm <ref>, the algorithm finds an ϵ-second-order stationary point of (<ref>) in T calls to the second-order oracle where


    T = 𝒪( L_u Q'^5/ϵ^5log_b( Q'/ϵ)) = 𝒪( L_u Q'^5/ϵ^5).





  
Remark. 
These complexity results for first and second-order are stationarity with respect to (<ref>).
We note that second order complexity result matches <cit.> and <cit.>.
However, the stationarity criteria and the definition of dual variable in these papers differ from ours.
We include more discussion on this in the Appendix.



  
Effect of  β_k in <ref>.
 We consider two cases, when g is the indicator of a convex set (or 0), the subdifferential set will be a cone (or 0), thus β_k will not have an effect. On the other hand, when g is a convex and Lipschitz contiunous function defined on the whole space, subdifferential set will be bounded <cit.>. This will introduce an error term in  <ref> that is of the order (1/β_k). One can see that b^k choice for β_k causes a linear decrease in this error term. In fact, all the examples in this paper fall into the first case. 














§ RELATED WORK 
 


ALM has a long history in the optimization literature, dating back to <cit.>.
In the special case of (<ref>) with a convex function f and a linear operator A, standard,  inexact, and linearized versions of ALM have been extensively studied <cit.>.

Classical works on ALM focused on the general template of (<ref>) with nonconvex f and nonlinear A, with arguably stronger assumptions and required exact solutions to  the subproblems of the form (<ref>), which appear in Step 2 of Algorithm <ref>, see for instance  <cit.>.

A similar analysis was conducted in <cit.> for the general template of (<ref>).
The authors considered inexact ALM and proved convergence rates for the outer iterates, under specific assumptions on the initialization of the dual variable. 
However, in contrast, the authors did not analyze how to solve the subproblems inexactly and  did not provide total complexity results with verifiable conditions.


Problem (<ref>) with similar assumptions to us is also studied in <cit.> and <cit.> for first-order and second-order stationarity, respectively, with explicit iteration complexity analysis.
As we have mentioned in Section <ref>, our second order iteration complexity result matches these theoretical algorithms with a simpler algorithm and a simpler analysis.
In addition, these algorithms require setting final accuracies since they utilize this information in the algorithm while our Algorithm <ref> does not set accuracies a priori.


<cit.> also considers the same template (<ref>) for first-order stationarity with a penalty-type method instead of ALM.
Even though the authors show 𝒪(1/ϵ^2) complexity, this result is obtained by assuming that the penalty parameter remains bounded. We note that such an assumption can also be used to improve our complexity results to match theirs.



<cit.> studies the general template (<ref>) with specific assumptions involving local error bound conditions for the (<ref>).
These conditions are studied in detail in <cit.>, but their validity for general SDPs (<ref>) has never been established. This work also lacks the total iteration complexity analysis presented here. 

Another work <cit.> focused on solving (<ref>) by adapting the primal-dual method of Chambolle and Pock <cit.>.
The authors proved the convergence of the method and provided convergence rate by imposing error bound conditions on the objective function that do not hold for standard SDPs.

<cit.> is the first work that proposes the splitting X=UU^⊤ for solving SDPs of the form (<ref>).
Following these works, the literature on Burer-Monteiro (BM) splitting for the large part focused on using ALM for solving the reformulated problem (<ref>).
However, this proposal has a few drawbacks: First, it requires exact solutions in Step 2 of Algorithm <ref> in theory, which in practice is replaced with inexact solutions. Second, their results only establish convergence without providing the rates.  In this sense, our work provides a theoretical understanding of the BM splitting with inexact solutions to Step 2 of Algorithm <ref> and complete iteration complexities.

<cit.> are among the earliest efforts to show convergence rates for BM splitting, focusing on the special case of SDPs without any linear constraints.
For these specific problems, they prove the convergence of gradient descent to global optima with convergence rates, assuming favorable initialization.
These results, however, do not apply to general SDPs of the form (<ref>) where the difficulty arises due to the linear constraints.


Another popular method for solving SDPs are due to <cit.>, focusing on the case where the constraints in (<ref>) can be written as a  Riemannian manifold after BM splitting.
In this case, the authors apply the Riemannian gradient descent and Riemannian trust region methods for obtaining first- and second-order stationary points, respectively. 
They obtain 𝒪(1/ϵ^2) complexity for finding first-order stationary points  and 𝒪(1/ϵ^3) complexity for finding second-order stationary points. 

While these complexities appear better than ours, the smooth manifold requirement in these works is indeed restrictive. In particular, this requirement holds for max-cut and generalized eigenvalue problems, but it is not satisfied for other important SDPs such as quadratic programming (QAP), optimal power flow and clustering with general affine constraints. 
In addition, as noted in <cit.>, per iteration cost of their method for max-cut problem is an astronomical 𝒪(d^6).

Lastly, there also exists a line of work for solving SDPs in their original convex formulation, in a storage efficient way <cit.>. These works have global optimality guarantees by their virtue of directly solving the convex formulation. On the downside, these works require the use of eigenvalue routines and exhibit significantly slower convergence as compared to nonconvex approaches <cit.>.






§ NUMERICAL EVIDENCE 


We first begin with a caveat:  It is known that quasi-Newton methods, such as BFGS and lBFGS, might not converge for nonconvex problems <cit.>. For this reason, we have used the trust region method as the second-order solver in our  analysis in Section <ref>, which is well-studied for nonconvex problems <cit.>.  Empirically, however, BFGS and lBGFS are extremely successful and we have therefore opted for those solvers in this section since the subroutine does not affect Theorem <ref> as long as the subsolver performs well in practice. 



 §.§ Clustering

Given data points {z_i}_i=1^n, the entries of the corresponding Euclidean distance matrix D ∈^n× n are D_i, j =   z_i - z_j^2. 
Clustering is then the problem of finding a co-association matrix Y∈^n× n such that Y_ij = 1 if points z_i and z_j are within the same cluster and Y_ij = 0 otherwise. In <cit.>, the authors provide  a SDP relaxation of the clustering problem, specified as 

    Y ∈^nxnmintr(DY)  Y1 = 1,
     tr(Y) = s, 
     Y≽ 0, Y ≥ 0,

where s is the number of clusters and Y is both positive semidefinite and has nonnegative entries. 
Standard  SDP solvers do not scale well with the number of data points n, since they often require projection onto the semidefinite cone with the complexity of 𝒪(n^3). We instead use the BM factorization to solve (<ref>), sacrificing convexity to reduce the computational complexity. More specifically, we solve the program

    V ∈^n× rmintr(DVV^⊤) 
    VV^⊤1 = 1,  V_F^2 ≤  s, 
      V  ≥ 0,

where 1∈^n is the vector of all ones.
Note that Y ≥ 0 in (<ref>) is replaced above by the much stronger but easier-to-enforce constraint V ≥ 0 in (<ref>), see <cit.> for the reasoning behind this relaxation.
Now, we can cast (<ref>) as an instance of (<ref>). Indeed, for every i≤ n, let x_i ∈^r denote the ith row of V. We next form x ∈^d with d = nr by expanding the factorized variable V, namely, 
x := [x_1^⊤, ⋯, x_n^⊤]^⊤∈^d,
and then set 

    f(x) =∑_i,j=1^n D_i, j⟨ x_i, x_j ⟩,
         g = δ_C,     
    A(x) = [x_1^⊤∑_j=1^n x_j -1, ⋯, x_n^⊤∑_j=1^n x_j-1]^⊤,

where C is the intersection of the positive orthant in ^d with the Euclidean  ball of radius √(s). In Appendix <ref>, we verify that Theorem <ref> applies to (<ref>) with f,g,A specified above. 




In our simulations, we use two different solvers for Step 2 of  Algorithm <ref>, namely, APGM and lBFGS. APGM is a  solver for nonconvex problems of the form (<ref>) with convergence guarantees to first-order stationarity, as discussed in Section <ref>. lBFGS is a limited-memory version of BFGS algorithm in <cit.> that approximately leverages the second-order information of the problem. 
We compare our approach against the following convex methods:


  * HCGM: Homotopy-based Conditional Gradient Method in <cit.> which directly solves (<ref>). 

  * SDPNAL+: A second-order augmented Lagrangian method for solving SDP's with nonnegativity constraints <cit.>. 




As for the dataset, our experimental setup is similar to that described by <cit.>. We use the publicly-available fashion-MNIST data in <cit.>, which is released as a possible replacement for the MNIST handwritten digits. Each data point is a 28× 28 gray-scale image, associated with a label from ten classes, labeled from 0 to 9. 
First, we extract the meaningful features from this dataset using a simple two-layer neural network with a sigmoid activation function. Then, we apply this neural network to 1000 test samples from the same dataset, which gives us a vector of length 10 for each data point, where each entry represents the posterior probability for each class.  Then, we form the ℓ_2 distance matrix D from these probability vectors. 
The solution rank for the template (<ref>) is known and  it is equal to number of clusters k <cit.>. As discussed in <cit.>, setting rank r>k leads more accurate reconstruction in expense of speed. 
Therefore, we set the rank to 20. 
For iAL lBFGS, we used β_1 = 1 and σ_1 = 10 as the initial penalty weight and dual step size, respectively. For HCGM, we used β_0 = 1 as the initial smoothness parameter. We have run SDPNAL+ solver with 10^-12 tolerance. 
The results are depicted in Figure <ref>.
We implemented 3 algorithms on MATLAB and used the software package for SDPNAL+ which contains mex files. 
It is predictable that the performance of our nonconvex approach would even improve by using mex files.





 §.§ Additional demonstrations

We provide several additional experiments in Appendix <ref>. 
Section  <ref> discusses a novel nonconvex relaxation of the standard basis pursuit template which performs comparable to the state of the art convex solvers. 
In Section <ref>, we provide fast numerical solutions to the generalized eigenvalue problem. 
In Section <ref>, we give a contemporary application example that our template applies, namely,  denoising with generative adversarial networks. 
Finally, we provide improved bounds for sparse quadratic assignment problem instances in Section <ref>. 









§ CONCLUSIONS


In this work, we have proposed and analyzed an inexact augmented Lagrangian method for solving nonconvex optimization problems with nonlinear constraints.
We prove convergence to the first and second order stationary points of the augmented Lagrangian function, with explicit complexity estimates.
Even though the relation of stationary points and global optima is not well-understood in the literature, we find out that the algorithm has fast convergence behavior to either global minima or local minima in a wide variety of numerical experiments.







§ ACKNOWLEDGEMENTS



The authors would like to thank Nicolas Boumal and Nadav Hallak for the helpful suggestions. 
 
 

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement n^∘ 725594 - time-data) and was supported by the Swiss National Science Foundation (SNSF) under  grant number 200021_178865 / 1. This project was also sponsored by the Department of the Navy, Office of Naval Research (ONR)  under a grant number N62909-17-1-2111 and was supported by Hasler Foundation Program: Cyber Human Systems (project number 16066). This research was supported by the PhD fellowship program of  the Swiss Data Science Center (SDSC) under grant lD number P18-07. 


abbrv      
   



























































































§ COMPLEXITY RESULTS




 §.§ First-Order Optimality 



Let us first consider the case where the solver in Step 2 is is the first-order algorithm APGM, described in detail in  <cit.>. At a high level, 
APGM makes use of ∇_x ℒ_β(x,y) in (<ref>), the proximal operator prox_g, and the classical Nesterov acceleration <cit.> to reach first-order stationarity for the subproblem in (<ref>).
 


Suppose that g=δ_𝒳 is the indicator function on a bounded convex set 𝒳⊂^d and let 

    ρ= max_x∈𝒳x,

be the radius of a ball centered at the origin that includes 𝒳. 
 Then, adapting the results in <cit.> to our setup, APGM reaches x_k in Step 2 of Algorithm <ref> after 

    𝒪 ( λ_β_k^2 ρ^2/ϵ_k+1^2)

(inner) iterations, where λ_β_k denotes the Lipschitz constant of ∇_xℒ_β_k(x, y), bounded in (<ref>). For the clarity of the presentation, we have used a looser bound in (<ref>) compared to <cit.>. 
Using (<ref>), we derive the following corollary, describing the total iteration complexity of Algorithm <ref> in terms of the number calls made to the first-order oracle in APGM. 


For b>1, let β_k =b^k for every k. If we use APGM from <cit.> for Step 2 of Algorithm <ref>, the algorithm finds an (ϵ_f,_̱k) first-order stationary point, 
 after T calls to the first-order oracle, where




    T = 𝒪( Q^3 ρ^2/ϵ^4log_b( Q/ϵ)) =   𝒪̃( Q^3ρ^2/ϵ^4).




Let K denote the number of (outer) iterations of Algorithm <ref> and let ϵ_f denote the desired accuracy of Algorithm <ref>, see (<ref>). Recalling Theorem <ref>, we can then write that 

    ϵ_f = Q/_̱K,

or, equivalently, _̱K = Q/ϵ_f. 

We now count the number of total (inner) iterations T of Algorithm <ref> to reach the accuracy ϵ_f. From (<ref>) and for sufficiently large k, recall that λ__̱k≤λ”_̱k is the smoothness parameter of the augmented Lagrangian. Then, from (<ref>) ad by summing over the outer iterations, we bound the total number of (inner) iterations of Algorithm <ref>  as


    T    = ∑_k=1^K𝒪 ( λ_β_k-1^2 ρ^2 /ϵ_k^2) 
        = ∑_k=1^K𝒪 (β_k-1^4 ρ^2  ) 
        (Step 1 of Algorithm <ref>)
       ≤𝒪(Kβ_K-1^4 ρ^2 )
        ( {_̱k}_k  is increasing)
     
       ≤𝒪( K Q^3ρ^2/ϵ_f^4).
         (see (<ref>))

In addition, if we specify β_k=b^k for all k, we can further refine T. Indeed, 

    β_K = b^K  ⟹   K = log_b ( Q/ϵ_f),

which, after substituting into (<ref>) gives the final bound in Corollary <ref>.









 §.§ Second-Order Optimality 


Let us now consider the second-order optimality case where the solver in Step 2 is the the trust region method developed in <cit.>.
Trust region method minimizes a quadratic approximation of the function within a dynamically updated trust-region radius.
Second-order trust region method that we consider in this section makes use of Hessian (or an approximation of Hessian) of the augmented Lagrangian in addition to first order oracles.

As shown in <cit.>, finding approximate second-order stationary points of convex-constrained problems is in general NP-hard. For this reason, we focus in this section on the special case of (<ref>) with g=0. 


















Let us compute the total computational complexity of Algorithm <ref> with the trust region method in Step 2, in terms of the number of calls made to the second-order oracle.  By adapting the result in <cit.> to our setup, we find that the number of (inner) iterations required in Step 2 of Algorithm <ref> to produce x_k+1 is


    𝒪 ( λ_β_k, H^2 (ℒ_β_k(x_1, y) - min_xℒ_β_k(x, y))/ϵ_k^3),



where λ_β, H is the Lipschitz constant of the Hessian of the augmented Lagrangian, which is of the order of β, as can be proven similar to Lemma <ref> and x_1 is the initial iterate of the given outer loop.
In <cit.>, the term ℒ_β(x_1, y) - min_xℒ_β(x, y) is bounded by a constant independent of ϵ.
We assume a uniform bound for this quantity for every β_k, instead of for one value of β_k as in <cit.>. Using (<ref>) and Theorem <ref>, we arrive at the following: 







For b>1, let β_k =b^k for every k.
We assume that 

    ℒ_β(x_1, y) - min_xℒ_β(x, y) ≤ L_u,    ∀β.

If we use the trust region method from <cit.> for Step 2 of Algorithm <ref>, the algorithm finds an ϵ-second-order stationary point of (<ref>) in T calls to the second-order oracle where


    T = 𝒪( L_u Q'^5/ϵ^5log_b( Q'/ϵ)) = 𝒪( L_u Q'^5/ϵ^5).





Before closing this section, we note that the remark after Corollary <ref> applies here as well. 



 §.§ Approximate optimality of (<ref>).


Corollary <ref> establishes the iteration complexity of Algorithm <ref> to reach approximate first-order stationarity for the equivalent formulation of (<ref>) presented in (<ref>). Unlike the exact case, approximate first-order stationarity in  (<ref>) does not immediately lend itself to approximate stationarity in (<ref>), and the study of approximate stationarity for the penalized problem (special case of our setting with dual variable set to 0) has also precedent in <cit.>. 


For a precedent in convex optimization for relating the convergence in augmented Lagrangian to the constrained problem using duality, see <cit.>.
For the second-order case, it is in general not possible to establish approximate second-order optimality for (<ref>) from Corollary <ref>, with the exception of linear constraints. <cit.> provides an hardness result by showing that checking an approximate second-order stationarity is NP-hard.





§ PROOF OF THEOREM <REF> 


For every k≥2, recall from (<ref>) and Step 2 of Algorithm <ref> that x_k satisfies 

    (-∇ f(x_k) - DA(x_k)^⊤ y_k-1
            - _̱k-1 DA(x_k)^⊤ A(x_k) ,∂ g(x_k) )  
        = (-∇_x Ł__̱k-1 (x_k ,y_k-1) ,∂ g(x_k) ) ≤ϵ_k.

With an application of the triangle inequality, it follows that 

    ( -_̱k-1 DA(x_k)^⊤ A(x_k) , ∂ g(x_k) ) 
           ≤∇ f(x_k ) +  DA(x_k)^⊤ y_k-1 + 
    ϵ_k,

which in turn implies that 

    ( -DA(x_k)^⊤ A(x_k) , ∂ g(x_k)/ _̱k-1  ) 
       ≤∇ f(x_k )/_̱k-1 +  DA(x_k)^⊤ y_k-1/_̱k-1 + 
    ϵ_k/_̱k-1
       ≤λ'_f+λ'_A y_k-1+ϵ_k/_̱k-1 ,

where λ'_f,λ'_A were defined in (<ref>). 






We next translate (<ref>) into a bound on the feasibility gap  A(x_k). Using the regularity condition (<ref>),  the left-hand side of (<ref>) can be bounded below as 

    ( -DA(x_k)^⊤ A(x_k) , ∂ g(x_k)/ _̱k-1  ) ≥νA(x_k) .
        (see (<ref>))






By substituting (<ref>) back into (<ref>), we find that 

    A(x_k)≤λ'_f + λ'_A y_k-1 + ϵ_k/ν_̱k-1.

In words, the feasibility gap is directly controlled by the dual sequence {y_k}_k. We next establish that the dual sequence is bounded. Indeed, for every k∈ K, note that 

    y_k    =  y_0 + ∑_i=1^k_i A(x_i)   (Step 5 of Algorithm <ref>)
       ≤y_0+ ∑_i=1^k _i A(x_i)    (triangle inequality)
       ≤y_0+ ∑_i=1^k A(x_1)log^2 2 / k log^2(k+1)  (Step 4)
       ≤y_0+ c  A(x_1) log^2 2 =: y_max,
  
where 

    c ≥∑_i=1^∞1/k log^2 (k+1).

Substituting (<ref>) back into (<ref>), we reach

    A(x_k)   ≤λ'_f + λ'_A y_max + ϵ_k/ν_̱k-1
       ≤ 2λ'_f +2 λ'_A y_max/ν_̱k-1 ,

where the second line above holds if k_0 is large enough, which would in turn guarantees that ϵ_k=1/_̱k-1 is sufficiently small since {_̱k}_k is increasing and unbounded. 




















It remains to control the first term in (<ref>). To that end, after recalling Step 2 of Algorithm <ref> and applying the triangle inequality, we can write that 

    ( -∇_x Ł__̱k-1 (x_k,y_k),    ∂ g(x_k) )  
       ≤( -∇_x Ł__̱k-1 (x_k,y_k-1) , ∂ g(x_k) ) 
         + ∇_x Ł__̱k-1 (x_k,y_k)-∇_x Ł__̱k-1 (x_k,y_k-1)   .

The first term on the right-hand side above is bounded by ϵ_k, by Step 5 of  Algorithm <ref>. 
For the second term on the right-hand side of (<ref>), we write that 

    ∇_x Ł__̱k-1 (x_k,y_k)-∇_x Ł__̱k-1 (x_k,y_k-1)   
        =  DA(x_k)^⊤ (y_k - y_k-1)     (see (<ref>))
       ≤λ'_A y_k- y_k-1    (see (<ref>))
        = λ'_A _k A (x_k)     (see Step 5 of Algorithm <ref>)
       ≤2λ'_A _k /ν_̱k-1( λ'_f+ λ'_Ay_max)  .
        (see (<ref>))

By combining (<ref>,<ref>), we find that 

    ( ∇_x Ł__̱k-1 (x_k,y_k),    ∂ g(x_k) )   
       ≤2λ'_A _k /ν_̱k-1( λ'_f+ λ'_Ay_max)  +  ϵ_k.

By combining (<ref>,<ref>), we find that 

    ( -∇_x Ł__̱k-1(x_k,y_k),∂ g(x_k)) +  A(x_k)
       ≤( 2λ'_A _k /ν_̱k-1( λ'_f+ λ'_Ay_max)  +  ϵ_k ) 
             + 2( λ'_f + λ'_A y_max/ν_̱k-1).

Applying  _k≤_1, we find that 

    ( -∇_x Ł__̱k-1(x_k,y_k),∂ g(x_k)) +  A(x_k)
       ≤ 2λ'_A_1 + 2/ν_̱k-1 ( λ'_f+λ'_A y_max) + ϵ_k.

For the second part of the theorem, we use the Weyl's inequality and Step 5 of Algorithm <ref> to write

    λ_min   (∇_xxℒ_β_k-1(x_k, y_k-1)) ≥λ_min (∇_xxℒ_β_k-1(x_k, y_k)) 
       -  σ_k ∑_i=1^m A_i(x_k)  ∇^2 A_i(x_k) .

The first term on the right-hand side is lower bounded by -ϵ_k-1 by Step 2 of Algorithm <ref>. We next bound the second term on the right-hand side above as 

    σ_k ∑_i=1^m A_i(x_k) ∇^2 A_i(x_k)   
       ≤σ_k √(m)max_i A_i(x_k)∇^2 A_i(x_k)
       ≤σ_k √(m)λ_A  2λ'_f +2 λ'_A y_max/ν_̱k-1,

where the last inequality is due to (<ref>,<ref>). 
Plugging into (<ref>) gives

    λ_min(∇_xxℒ_β_k-1(x_k, y_k-1))
       ≥ -ϵ_k-1 - σ_k √(m)λ_A  2λ'_f +2 λ'_A y_max/ν_̱k-1,

which completes the proof of Theorem <ref>. 

































































































§ PROOF OF LEMMA <REF>


Note that

    ℒ_β(x,y) = f(x) + ∑_i=1^m y_i A_i (x)  + /2∑_i=1^m (A_i(x))^2,

which implies that 

    ∇_x ℒ_β(x,y) 
        = ∇ f(x) + ∑_i=1^m y_i  ∇ A_i(x) + /2∑_i=1^m A_i(x) ∇ A_i(x) 
        = ∇ f(x) + DA(x)^⊤ y  + ḎA(x)^⊤ A(x),

where DA(x) is the Jacobian of A at x. By taking another derivative with respect to x, we reach 

    ∇^2_x ℒ_β(x,y)     = ∇^2 f(x) + ∑_i=1^m ( y_i + A̱_i(x) ) ∇^2 A_i(x) 
            +∑̱_i=1^m ∇ A_i(x) ∇ A_i(x)^⊤.

It follows that 

    ∇_x^2 ℒ_β(x,y)
       ≤∇^2 f(x)  +  max_i ∇^2 A_i(x) (y_1+A(x)_1 ) 
            +β∑_i=1^m ∇ A_i(x)^2 
       ≤λ_h+  √(m)λ_A  (y+A(x)) + DA(x)^2_F.

For every x such that x≤ρ and A(x)≤ρ, we conclude that 

    ∇_x^2 ℒ_β(x,y)   ≤λ_f + √(m)λ_A (y + ρ̱' )+ m̱a̱x̱_x≤ρDA(x)_F^2,

which completes the proof of Lemma <ref>. 



































































§ CLUSTERING 



We only verify the condition in (<ref>) here. 
Note that

    A(x) = VV^⊤1-  1,


    DA(x)     = 
    [
    [  w_1,1 x_1^⊤            ⋯  w_1,n x_1^⊤;            ⋮;   w_n,1x_n^⊤            ⋯ w_n,n1 x_n^⊤ ]] 
        = [ 
    [ V ⋯ V ]]
    + 
    [
    [ x_1^⊤      ;           ⋱      ;             x_n^⊤ ]],


where w_i.i=2 and w_i,j=1 for i j. In the last line above, n copies of V appear and the last matrix above is block-diagonal. For x_k, define V_k accordingly and  let x_k,i be the ith row of V_k. 
Consequently,

    DA(x_k)^⊤ A(x_k)     = 
    [
    [ (V_k^⊤ V_k - I_n) V_k^⊤1;                        ⋮; (V_k^⊤ V_k - I_n) V_k^⊤1 ]] 
            + 
    [
    [ x_k,1 (V_k V_k^⊤1-   1)_1;                         ⋮; x_k,n (V_k V_k^⊤1-   1)_n ]],

where I_n∈^n× n is the identity matrix. 
Let us make a number of simplifying assumptions. First, we assume that x_k< √(s) (which can be  enforced in the iterates by replacing C with (1-ϵ)C for a small positive ϵ in the subproblems).  Under this assumption, it follows that 

    (∂ g(x_k))_i =
    
    0     (x_k)_i > 0
    {a: a≤ 0}    (x_k)_i = 0,
         i≤ d.

Second, we assume that V_k has nearly orthonormal columns, namely, V_k^⊤ V_k ≈ I_n. This can also be  enforced in each iterate of Algorithm <ref> and naturally corresponds to well-separated clusters. While a more fine-tuned argument can remove these assumptions, they will help us simplify the presentation here. Under these assumptions, the (squared) right-hand side of (<ref>) becomes

    ( -DA(x_k)^⊤ A(x_k) , ∂ g(x_k)/_̱k-1)^2 
        = ( -DA(x_k)^⊤ A(x_k) )_+^2
         (a_+ = max(a,0))
     
        = 
    [
    [ x_k,1  (V_k V_k^⊤1-   1)_1;                          ⋮; x_k,n  (V_k V_k^⊤1-   1)_n ]]
    ^2 
         (x_k∈ C ⇒ x_k≥ 0)
    
        = ∑_i=1^n  x_k,i^2 (V_kV_k^⊤1-1)_i^2 
       ≥min_i  x_k,i^2
    ·∑_i=1^n  (V_kV_k^⊤1-1)_i^2 
        = min_i  x_k,i^2
    · V_kV_k^⊤1-1^2.

Therefore, given a prescribed ν, ensuring min_i x_k,i≥ν guarantees (<ref>). When the algorithm is initialized close enough to the constraint set, there is indeed no need to separately enforce (<ref>). In practice, often n exceeds the number of true clusters and a more intricate analysis is required to establish (<ref>) by restricting the argument to a particular subspace of ^n. 


























































§ ADDITIONAL EXPERIMENTS




 §.§ Basis Pursuit


Basis Pursuit (BP) finds sparsest solutions of an under-determined system of linear equations by solving

    min_zz_1 
    Bz = b,

where B ∈^n × d and b ∈^n. 

Various primal-dual convex optimization algorithms are available in the literature to solve BP, including <cit.>.  

We compare our algorithm against state-of-the-art primal-dual convex methods for solving (<ref>), namely, Chambole-Pock <cit.>, ASGARD <cit.> and ASGARD-DL <cit.>. 





Here, we take a different approach and cast (<ref>) as an instance of (<ref>). Note that any z ∈^d can be decomposed as z = z^+ - z^-, where z^+,z^-∈^d are the positive and negative parts of z, respectively. Then consider the change of variables z^+ = u_1^∘ 2 and  z^-= u_2^∘ 2∈^d, where ∘ denotes element-wise power. Next, we concatenate u_1 and u_2 as x := [ u_1^⊤, u_2^⊤ ]^⊤∈^2d and define B := [B, -B] ∈^n × 2d. Then, (<ref>) is equivalent to (<ref>)  with 

    f(x) =   x^2,    g(x) = 0,
    A(x) = Bx^∘ 2- b.


We draw the entries of B independently from a zero-mean and unit-variance Gaussian distribution.  For a fixed sparsity level k, the support of z_*∈^d and its nonzero amplitudes are also drawn from the standard Gaussian distribution. 
Then the measurement vector  is created as b = Bz + ϵ, where ϵ is the noise vector with entries drawn independently from the zero-mean Gaussian distribution with variance σ^2=10^-6. 

The results are compiled in Figure <ref>. Clearly,  the performance of Algorithm <ref> with a second-order solver for BP is comparable to the rest. 

It is, indeed, interesting to see that these type of nonconvex relaxations gives the solution of convex one and first order methods succeed.







  
Discussion: 
The true potential of our reformulation is in dealing with more structured norms rather than ℓ_1, where computing the proximal operator is often intractable. One such case is the latent group lasso norm <cit.>, defined as 

    z_Ω = ∑_i=1^I  z_Ω_i,

where {Ω_i}_i=1^I are (not necessarily disjoint) index sets of {1,⋯,d}. Although not studied here, we believe that the nonconvex framework presented in this paper can serve to solve more complicated problems, such as the latent group lasso. We leave this research direction for future work.




  
Condition verification:
In the sequel, we verify that Theorem <ref> indeed applies to (<ref>) with the above f,A,g. 
Note that 


    DA(x) = 2 Bdiag(x),

where diag(x)∈^2d× 2d is the diagonal matrix formed by x. The left-hand side of (<ref>) then reads as 

    dist(  -DA(x_k)^⊤ A(x_k) , ∂ g(x_k)/_̱k-1) 
        = dist(  -DA(x_k)^⊤ A(x_k) , {0})      (g≡ 0)
        = DA(x_k)^⊤ A(x_k)   
        =2  diag(x_k) B^⊤ ( Bx_k^∘ 2 -b) .
        (see (<ref>))

To bound the last line above, let x_* be a solution of (<ref>) and note that B x_*^∘ 2 = b by definition.  Let also z_k,z_*∈^d denote the vectors corresponding to x_k,x_*. Corresponding to x_k, also define u_k,1,u_k,2 naturally and let |z_k| = u_k,1^∘ 2 + u_k,2^∘ 2∈^d be the vector of amplitudes of z_k. To simplify matters, let us assume also that B is full-rank. 
We then rewrite the norm in the last line of (<ref>) as 

    diag(x_k) B^⊤ ( Bx_k^∘ 2 -b) ^2 
        =   diag(x_k) B^⊤B (x_k^∘ 2 -x_*^∘ 2) ^2 
         (B x_*^∘ 2 = b) 
    
        =  diag(x_k)B^⊤ B (x_k - x_*) ^2 
        =  diag(u_k,1)B^⊤ B (z_k - z_*) ^2 
            + diag(u_k,2)B^⊤ B (z_k - z_*) ^2 
        =  diag(u_k,1^∘ 2+ u_k,2^∘ 2) B^⊤ B (z_k - z_*) ^2 
        =  diag(|z_k|) B^⊤ B (z_k - z_*) ^2 
       ≥η_n (   B diag(|z_k|) )^2
     B(z_k - z_*) ^2 
        =
    η_n (   B diag(|z_k|) )^2
     B z_k -b ^2
         ( Bz_* = B x^∘2_* = b) 
       ≥min_|T|=nη_n(B_T)·  |z_k,(n)|^2 Bz_k - b^2,

where η_n(·) returns the nth largest singular value of its argument. In the last line above, B_T is the restriction of B to the columns indexed by T of size n. Moreover, z_k,(n) is the nth largest entry of z in magnitude. 
Given  a prescribed ν,  (<ref>) therefore holds if 

    |z_k,(n) | ≥ν/2 √(min_|T|=nη_n(B_T)) ,

for every iteration k. If Algorithm <ref> is initialized close enough to the solution z^* and the entries of z^* are sufficiently large in magnitude, there will be no need to directly enforce (<ref>). 




 §.§ Generalized Eigenvalue Problem







Generalized eigenvalue problem has extensive applications in machine learning, statistics and data analysis <cit.>.
The well-known nonconvex formulation of the problem is <cit.> given by

    x∈ℝ^nmin x^⊤ C x 
     
    x^⊤ B x = 1,

where B, C ∈^n × n are symmetric matrices and B is positive definite, namely, B ≻ 0.

The generalized eigenvector computation is equivalent to performing principal component analysis (PCA) of C in the norm B. It is also equivalent to computing the top eigenvector of symmetric matrix S = B^-1/2CB^1/2 and multiplying the resulting vector by B^-1/2. However, for  large values of n, computing B^-1/2 is extremely expensive. 
The natural convex SDP relaxation for (<ref>) involves lifting Y = xx^⊤ and removing the nonconvex rank(Y) = 1 constraint, namely, 

    Y ∈^n × nmintr(CY)
    tr(BY) = 1,    X ≽ 0.

Here, however, we opt to directly solve (<ref>) because it  fits into our template with

    f(x) =    x^⊤ C x,    g(x) = 0,
    
    A(x) =    x^⊤ B x - 1.

We compare our approach against three different methods: manifold based Riemannian gradient descent and Riemannian trust region methods in <cit.> and  the linear system solver in <cit.>, abbrevated as GenELin.  We have used Manopt software package in <cit.> for the manifold based methods. For GenELin, we have utilized Matlab's backslash operator as the linear solver. The results are compiled in Figure <ref>.



  
Condition verification:
Here, we verify the regularity condition in (<ref>) for problem (<ref>). Note that


    DA(x) = (2Bx)^⊤.

Therefore,

    ( -DA(x_k)^⊤ A(x_k) , ∂ g(x_k)/_̱k-1)^2     = ( -DA(x_k)^⊤ A(x_k) , {0})^2 
         (g≡ 0) 
     
        =  DA(x_k)^⊤ A(x_k) ^2 
    
        = 2Bx_k (x_k^⊤ Bx_k - 1)^2 
        (see (<ref>))
        = 4 (x_k^⊤ Bx_k - 1)^2Bx_k^2 
    
        = 4Bx_k^2 A(x_k)^2 
        (see (<ref>))
       ≥η_min(B)^2x_k^2 A(x_k)^2,

where η_min(B) is the smallest eigenvalue of the positive definite matrix B. 
Therefore, for a prescribed ν, the regularity condition in (<ref>) holds with x_k≥ν /η_min for every k. If the algorithm is initialized close enough to the constraint set, there will be again no need to directly enforce this latter condition.




 §.§ ℓ_∞ Denoising with a Generative Prior

The authors of <cit.> have proposed to
project onto the range of a Generative Adversarial network (GAN)
<cit.>, as a way to defend against adversarial examples. For a
given noisy observation x^* + η, they consider a projection in the
ℓ_2 norm. We instead propose to use our augmented Lagrangian method to
denoise in the ℓ_∞ norm, a much harder task:

    [       x, zmin               x^* + η - x_∞;          s.t.       x=G(z). ]




We use a pretrained generator for the MNIST dataset, given by a standard
deconvolutional neural network architecture <cit.>. We compare the
succesful optimizer Adam <cit.> and gradient Descent against our
method. Our algorithm involves two forward and one backward pass through the
network, as oposed to Adam that requires only one forward/backward pass. For
this reason we let our algorithm run for 2000 iterations, and Adam and GD for
3000 iterations. Both Adam and gradient descent generate a sequence of feasible
iterates x_t=G(z_t). For this reason we plot the objective evaluated at the
point G(z_t) vs iteration count in figure <ref>. Our
method successfully minimizes the objective value, while Adam and GD do not.



 §.§ Quadratic assginment problem


Let K, L be n × n symmetric metrices. QAP in its simplest form can be written as


    maxtr(KPLP),        subject to   P   be a permutation matrix


A direct approach for solving (<ref>) involves a combinatorial search. 
To get the SDP relaxation of (<ref>), we will first lift the QAP to a problem involving a larger
matrix. Observe that the objective function takes the form


    tr((K⊗ L) (vec(P)vec(P^⊤))),

 where ⊗ denotes the Kronecker product. Therefore, we can recast (<ref>) as 
 
 
    tr((K⊗ L) Y)     subject to    Y = vec(P)vec(P^⊤),

where P is a permutation matrix. We can relax the equality constraint in (<ref>) to a semidefinite constraint and write it in an equivalent form as

 
    X = 
    [        1 vec(P)^⊤;   vec(P)        Y ]≽ 0    for a symmetric X ∈𝕊^(n^2+1) × (n^2+1)


We now introduce the following constraints such that 

    B_k(X) = b_k,       b_k∈^m_k

to make sure X has a proper structure. Here, B_k is a linear operator on X and the total number of constraints is m = ∑_k m_k. Hence, SDP relaxation of the quadratic assignment problem takes the form,  








    max      ⟨ C, X ⟩
    
     subject to       P1 = 1,    1^⊤ P = 1,   P≥ 0 
       trace_1(Y) = I    trace_2(Y) = I 
       vec(P) = diag(Y) 
       trace(Y) = n   [        1 vec(P)^⊤;   vec(P)        Y ]≽ 0,

where trace_1(.) and trace_2(.) are partial traces satisfying,

    trace_1(K ⊗ L) = trace(K)L       and     trace_2(K⊗ L)=  K trace(L)


    trace_1^*(T) = I ⊗ T        and     trace_2^*(T)=  T ⊗ I

1st set of equalities are due to the fact that permutation matrices are doubly stochastic. 2nd set of equalities are to ensure permutation matrices are orthogonal, i.e., PP^⊤ = P^⊤ P=I. 3rd set of equalities are to enforce every individual entry of the permutation matrix takes either 0 or 1,  i.e., X_1, i = X_i,i   ∀ i ∈ [1, n^2+1]. . Trace constraint in the last line is to bound the problem domain. 
By concatenating the B_k's in (<ref>), we can rewrite (<ref>) in standard SDP form as 


    max      ⟨ C, X ⟩
    
     subject to       B(X) = b,    b∈^m  
       trace(X) = n+1 
        X_ij≥ 0,      i,j   𝒢
        X≽0,

where 𝒢 represents the index set for which we introduce the nonnegativities. When 𝒢 covers the wholes set of indices, we get the best approximation to the original problem. However, it becomes computationally undesirable as the problem dimension increases. Hence, we remove the redundant nonnegativity constraints and enforce it for the indices where Kronecker product between K and L is nonzero. 

We penalize the non-negativity constraints and add it to the  augmented Lagrangian objective since a projection to the positive orthant approach in the low rank space as we did for the clustering does not work here. 

We take <cit.> as the baseline. This is an SDP based approach for solving QAP problems containing a sparse graph. 
We compare against the best feasible upper bounds reported in <cit.> for the given instances. Here, optimality gap is defined as 

    %Gap = |bound - optimal|/optimal× 100

We used a (relatively) sparse graph data set from the QAP library.
We run our low rank algorithm for different rank values. r_m in each instance corresponds to the smallest integer satisfying the Pataki bound <cit.>. 
Results are shown in Table <ref>. 
Primal feasibility values except for the last instance esc128 is less than 10^-5 and we obtained bounds at least as good as the ones reported in <cit.> for these problems.

For esc128, the primal feasibility is ≈ 10^-1, hence, we could not manage to obtain a good optimality gap.













































